How to compute acoustic power generated by isotropic turbulence

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("k");

             // Lookup epsilon
             Info<< "Looking up field epsilon\n" << endl;
             const volScalarField& epsilon = mesh().lookupObject("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:

お問い合わせ