機能
How to compute acoustic power generated by isotropic turbulence
- 2018/05/14
- 日本ESI
Commercial CFD codes like Fluent use a very simple equation that calculates the "acoustic power" due to vorticity:
[source http://www.afs.enea.it/project/neptunius/docs/fluent/html/th/node237.htm#eq-acoustic-power ]
Suppose you have launched a steady analysis with simpleFoam using the k-epsilon turbulence model. How to compute the acoustic power P [W/m³] directly in OpenFOAM? The formula is fairly easy to implement:
P = alpha * rho * epsilon * ( sqrt(2 * k) / c ) ^ 5
where:
alpha: a rescaled constant = 0.1
rho: air density
k / epsilon: turbulence fields given from simpleFoam
c: speed of sound
We can perform its computation by inserting the C++ custom code we need directly into a function object in controlDict, as the following listing shows:
application simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 2000;
deltaT 1;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
acousticPower
{
// Load the library containing the 'coded' functionObject
libs ("libutilityFunctionObjects.so");
type coded;
// Name of on-the-fly generated functionObject
name acousticPower;
codeEnd
#{ // from here we can insert the C++ source code we need to compute the acoustic power
// alpha constant [dimensionless]
scalar alpha = 0.1;
// rho air density [kg/m³]
dimensionedScalar rho
(
"rho",
dimensionSet(1, -3, 0, 0, 0, 0 ,0),
1.225
);
// c speed of sound [m/s]
dimensionedScalar c
(
"c",
dimensionSet(0, 1, -1, 0, 0, 0 ,0),
340.0
);
// Pref reference acoustic power [W/m^3]
dimensionedScalar Pref
(
"Pref",
dimensionSet(1, -1, -3, 0, 0, 0 ,0),
10e-12
);
// Lookup k
Info<< "Looking up field k\n" << endl;
const volScalarField& k = mesh().lookupObject<volScalarField>("k");
// Lookup epsilon
Info<< "Looking up field epsilon\n" << endl;
const volScalarField& epsilon = mesh().lookupObject<volScalarField>("epsilon");
volScalarField acousticPower("acousticPower", alpha * rho * epsilon * pow( (sqrt(2*k)/c) , 5) );
Info<<"Writing acousticPower to " << acousticPower.objectPath()
<< endl;
acousticPower.write();
// Pa = alpha * rho * epsilon * pow( (sqrt(2*k)/c) , 5)
// LP = 10*log(Pa/Pref) = 10.0*log( alpha * rho * epsilon * pow( (sqrt(2*k)/c) , 5)/Pref)
// we do a transformation of base log() = ln(10)*log10() = 2.302585093 * log10() to avoid "log" keyword conflict
// LP = 10*log(Pa/Pref) = 10.0*ln(10)*log10( alpha * rho * epsilon * pow( (sqrt(2*k)/c) , 5)/Pref) // where ln(10) = 2.302585093
volScalarField acousticPowerDB("acousticPowerDB", 10.0*2.302585093*log10( alpha * rho * epsilon * pow( (sqrt(2*k)/c) , 5)/Pref) );
Info<<"Writing acousticPower in dB to " << acousticPowerDB.objectPath()
<< endl;
acousticPowerDB.write();
#};
}
}
// ************************************************************************* //
In this way, once the steady analysis is completed, OpenFOAM will write two additional fields, acousticPower [W/m³] and acousticPowerDB [in decibel], that you can postprocess in Paraview.
Here are for example the results obtained by applying the code to the pitzdaily tutorial:
tutorials/incompressible/simpleFoam/pitzDaily
[NOTE: before running the code check that in constant/turbulenceProperties the RASModel is set to "kEpsilon;" ]
Let's apply the code to the motorbike tutorial and compare 20m/s vs 40m/s:
Another example comparing round vs square pillars. Vorticity and acoustic power are higher with the square shape: