OpenFOAM

How to compute acoustic power generated by isotropic turbulence

• 2018/05/14

Commercial CFD codes like Fluent use a very simple equation that calculates the "acoustic power" due to vorticity: 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:   