Acoustics PDE¶

Governing equations¶

The governing equation for acoustic simulations is the linear homogeneous wave equation for the acoustic pressure $p_{\mathrm a}$ $$\left( \frac{1}{c_0^2}\frac{\partial^2 }{\partial t^2} - \Delta \right) p_{\mathrm a} = 0 \,, \label{eq:acoupressure}$$ Thereby, $c_0$ is the speed of sound, e.g. approximately $343\,$m/s in air at $20^ \circ$C. For establishing the weak formulation a test funtion $\varphi$ is multiplied to the linear wave equation \eqref{eq:acoupressure} and it is integrated over the computation domain $\Omega$. In a next step Green's first integration theorem is applied to reduce the order of spatial derivatives. This gives an additional integral term over the boundary surfaces $\partial \Omega = \Gamma$ $$\label{eq:weak_formulation} \int_{\Omega} \frac{1}{c_0^2}\, \varphi \frac{\partial^2 p_{\mathrm a}}{\partial t^2} {\mathrm {d}} \Omega + \int_{\Omega} \nabla \varphi \cdot \nabla p_{\mathrm a} {\mathrm {d}} \Omega - \int_{\Gamma} \varphi \nabla p_{\mathrm a} \cdot {\mathbf n} {\mathrm {d}} \Gamma = 0$$ with ${\mathbf n}$ as the surface normal vector. Hence, the integral over the domain surface $\Gamma$ is used to apply boundary conditions onto simulations.

For some cases, an acoustic potential formulation is appropriate, which also solves the acoustic wave equation $$\left( \frac{1}{c_0^2}\frac{\partial^2 }{\partial t^2} - \Delta \right) \psi_{\mathrm a} = 0 \, , \label{eq:1}$$

with $\psi_{\mathrm a}$ as the acoustic scalar potential

$$p_{\mathrm a} = \rho_0 \frac{\partial \psi_{\mathrm a}}{\partial t} \,.$$

Boundary conditions¶

In acoustics, two main boundary conditions exist. (1) Soundhard, a homogeneous Neumann boundary condition

$$\nabla p_{\textrm{a}} \cdot n = 0 \, ,$$

and (2) soundsoft, a homogeneous Dirichlet boundary condition

$$p_{\textrm{a}} = 0 \, .$$

A sound soft boundary condition can be variably defined in openCFS. Please note that all surfaces that are not manually applied to any boundary condition are automatically defined as sound hard (natural boundary condition).

They are defined in the simulation-xml file within <bcsAndLoads> section.

        <bcsAndLoads>
<pressure name="" value=""/>
<potential name="" value="" />
<impedance name="" volumeRegion="" impedanceId=""/>
<normalVelocity name="" value="" volumeRegion=""/>
<normalAcceleration name="" value="" volumeRegion=""/>
<soundSoft name=""/>

• pressure / potential: Can only be defined on surface elements.
• impedance: needs an <impedanceList> defining the acoustic impedance for the corresponding surface, e,g.,
        <impedanceList>
<interpolImpedance id="imID">
<dataName_real>Datei_A.txt</dataName_real>
<dataName_imag>Datei_B.txt</dataName_imag>
</interpolImpedance>
<fctImpedance id="fctImId">
<fctReal>10*x^2</fctReal>
<fctImag>10*f^2</fctImag>
</fctImpedance>
</impedanceList>

• normalVelocity: Normal velocity of a surface $v_n = \mathbf v \cdot \mathbf n$. Has to be defined on elements of a surface region (can not be defined on a node!). Can only be prescired for potential formulation and for the harmonic pressure formulation of PDE (Does not work for SoSatLaplace PDE).
• normalAcceleration: Normal acceleration of a surface $a_n = \mathbf a \cdot \mathbf n$. Has to be defined on elements of a surface. Can only be prescribed for pressure formulation.
• soundsoft: corresponds to an homogeneous Dirichlet boundary condition, e.g. $p_\mathrm{a}=0$.

In case of free radiation, an absorbing boundary condition (ABC) or a perfectly matched layer (PML) can be applied.

Sources¶

For computing acoustic simulations a right hand side is applied to the homogeneous wave equation. Thereby, it is possible to apply a synthetic source term, to certain nodes/elements or surfaces, or to use aeroacoustic source terms (1 2), computed using an aeroacoustic analogy . Consequently, the acousticPDE can be used for computing pure acoustic, as well as aeroacoustic simulations that relay on the hybrid aeroacoustic workflow 3.

For applying synthethic acoustic sources 4 different RHS loads can be applied, both in the time and frequency domain. (Note: If a load is definied in the freuqnecy domain, a phase can be defined as well)

    <bcsAndLoads>
<rhsValues name="" value=""/>
<rhsDensity name="" value=""/>
<rhsDensityVector name="" value=""/>

• rhsValues: Has to be applied to nodes and is not avaiable for SoSatLaplace. This source term replaces the whole RHS of the weak formulation of the inhomogeneous wave equation. If aeroacoustic source terms are used and they are interpolated conservatively, they have to be prescribed using rhsValues
• rhsDensity: Has to be defined on nodes.
• rhsDensityVector: Has to be defined on nodes. Using Lighthill's analogy, a vectoral source term $\nabla \cdot \mathbf [ \mathbf T \mathbf ]$ is derived. For using this source term rhsDensityVector has to be used.

Material¶

A standard material file for an acoustic field computation with openCFS can be seen below. It is required to define the density $\rho$ and the compression modulus $K$. During the simulation process openCFS uses those material properties for computing the speed of sound $c_0$ $$c_0 = \sqrt{\frac{K}{\rho}} \, .$$ Additionally, it is possible to define a Rayleigh damping in the material file. Thereby, a numerical damping matrix $\mathbf [ \mathbf C \mathbf ]$, based on a combination of the mass matrix $\mathbf [ \mathbf M \mathbf ]$ and the linear stiffness matrix $\mathbf [ \mathbf K \mathbf ]$ is defined

$$\mathbf [ \mathbf C \mathbf ]= \alpha_{\mathrm M} \mathbf [ \mathbf M \mathbf ] + \alpha_{\mathrm K} \mathbf [ \mathbf K \mathbf ] \, .$$

Thereby, the lossTangensDelta and the measuredFreq contribute to $\alpha_{\mathrm M}$ and $\alpha_{\mathrm K}$

$$\tan(\delta_i) = \frac{\alpha_{\mathrm M}+\alpha_{\mathrm k} \omega_i ^2}{\omega_i} \, .$$

For more details look into 5

  <material name="air_approx">
<acoustic>
<density>
<linear>
<real> 1.184 </real>
</linear>
</density>
<compressionModulus>
<linear>
<real> 141610 </real>
</linear>
</compressionModulus>
<damping>
<rayleigh>
<lossTangensDelta>100</lossTangensDelta>
<measuredFreq>100</measuredFreq>
</rayleigh>
</damping>
</acoustic>
</material>



To model e.g MPPs (Micro-Perforated-Plates) the wave equation in the pressure formulation is transformed into the frequency domain (Helmholtz Equation) 6

$$\frac{\omega ^2}{K_0(\omega)} \hat{p}_{\mathrm a} + \nabla \cdot \frac{1}{\rho_0(\omega)} \nabla \hat{p}_{\mathrm a} = 0 \, ,$$

delivering the frequency depending density $\rho_0(\omega)$ and bulk modulus $K_0(\omega$). For defining such parameters e.g. the Delany-Bazley model can be used. The frequency depending parameters can be stored in text files and loaded into the material file as shown below.


<material name="Porous1">
<acoustic>
<density>
<linear>
<real> 1.2 </real>
</linear>
</density>
<densityComplex>
<linear>
<real> sample1D('porous1RhoR.txt',f,1) </real>
<imag> sample1D('porous1RhoI.txt',f,1) </imag>
</linear>
</densityComplex>
<compressionModulus>
<linear>
<real>  138720.0 </real>
</linear>
</compressionModulus>
<compressionModulusComplex>
<linear>
<real> sample1D('porous1CompR.txt',f,1) </real>
<imag> sample1D('porous1CompI.txt',f,1) </imag>
</linear>
</compressionModulusComplex>
</acoustic>
</material>



Analysis Types¶

Since in general we are dealing with a time-dependent PDE, we can distinguish three analysis types:

• Transient: $\frac{\partial }{\partial t}\neq 0$
    <analysis>
<transient>
<numSteps>10</numSteps>
<deltaT>0.1</deltaT>
</transient>
</analysis>

• Harmonic: $\partial / \partial t (\cdot) = j \omega (\cdot)$
<analysis>
<harmonic>
<numFreq>10</numFreq>
<startFreq>10</startFreq>
<stopFreq>3000</stopFreq>
<sampling>linear</sampling>
</harmonic>
</analysis>

• Eigenfrequency: $\partial / \partial t (\cdot) = j \omega (\cdot)$ and $RHS = 0$
    <analysis>
<eigenFrequency>
<numModes>5</numModes>
<freqShift>0</freqShift>
<writeModes>yes</writeModes>
</eigenFrequency>
</analysis>


Postprocessing results¶

The following results can be computed, using openCFS and the acoustic PDE

Defined on nodes

• Acoustic pressure: Can only be computed as nodal result using the acoustic pressure formulation of the PDE, which is the standard formulation. So <acoustic formulation="acouPressure"> and <acoustic> in the <pdeList> is the same!
<nodeResult type="acouPressure">

• Acoustic potential: Can only be defined for the acoustic potential formulation of the PDE defined in the xml-file by <acoustic formulation="acouPotential">.
<nodeResult type="acouPotential">

• 1st order time derivatives of acoustic potential: First time derivative of acoustic potential formulation, which corresponds with the acoustic pressure. Can only be defined when using <acoustic formulation="acouPotential">.
<nodeResult type="acouPotentialD1">

• 2nd order time derivatives of acoustic potential: First time derivative of acoustic potential formulation, which corresponds with the acoustic pressure. Can only be defined when using <acoustic formulation="acouPotential">.
<nodeResult type="acouPotentialD2">

• Nodal loads: right hand side source term in FE formulation as nodal loads
<nodeResult type="acouRhsLoad">


Defined on elements

• Acoustic pressure: in case of acoustic potential formulation when defining <acoustic formulation="acouPotential">
<elemResult type="acouPressure"/>

• Acoustic particle velocity
<elemResult type="acouVelocity">

• Acoustic intensity $$\mathbf I_\mathrm{a} = p_\mathrm{a} \mathbf v_\mathrm{a}$$
<elemResult type="acouIntensity">

<elemResult type="acouRhsLoadDensity"/>

• Acoustic speed of sound: e.g., when the it depends on temperature
<elemResult type="acouSpeedOfSound"/>

• Acoustic density: acoustic density computedt via stete equation with acoustic pressure
<elemResult type="density"/>

• PML damping factor: Can only be defined when a PML is used
<elemResult type="pmlDampFactor"/>


Defined on surface elements

• Acoustic particle velocity normal component
<elemResult type="acouNormalVelocity">

• Acoustic surface intensity
<elemResult type="acouSurfIntensity">


Defined on surface region

• Acoustic power on surface region (surface region result) $$P_\mathrm{a} = \int\limits_\Gamma \mathbf I_\mathrm{a} \cdot \mathbf n\, \mathrm{d}\Gamma'$$
 <surfRegionResult type="acouPower">


Defined on regions

• Acoustic energy:
 <regionResult type="acouEnergy">

• Acoustic kinetic energy:
<regionResult type="acouKinEnergy"/>

• Acoustic potential energy:
<regionResult type="acouPotEnergy"/>


Aeroacoustic Simulations¶

Regarding aeroacoustic propagation simulations 1, two different wave propagation models will be explained here:

1) Lighthill's aeroacoustic analogy:

Lighthills aeroacoustic analogy is an exact reformulation of the conservation equations of mass and momentum and considers all aerodynamic and aeroacoustic pressure production and propagation mechanisms. For propagation simulation Lighthill's inhomogenous wave equation is considered

\left( \frac{1}{c_0^2} \frac{\partial^2}{\partial t ^2}-\Delta \right)(c_0^2(\rho-\rho_0))= \nabla \cdot \nabla \cdot \mathbf [ \mathbf T \mathbf ] \, .

The right hand side of Lighthill's equation is called Lighthill's source term and has to be derived from a CFD simulation. For establishing the weak formulation of Lighthill's equation, the some methodology as presented above is used.

$$\begin{split} \frac{1}{c_0^2} \int_\Omega \varphi \frac{\partial ^2}{\partial t^2} & (p - p_0) {\mathrm {d}} \Omega + \int_ \Omega \nabla \varphi \cdot \nabla (p - p_0) {\mathrm {d}} \Omega = \\ & \int_\Gamma ( \nabla(p - p_0) \cdot \mathbf n + \nabla \cdot \mathbf [ \mathbf T \mathbf ] \cdot \mathbf n )\varphi {\mathrm {d}} \Gamma - \int_\Omega \nabla \cdot \mathbf [ \mathbf T \mathbf ] \cdot \nabla \varphi {\mathrm {d}} \Omega \, , \end{split}$$

For computing a simulation with Lighthill's analogy the pressure formulation of the wave equation has to be used.

     <pdeList>
<acoustic formulation="acouPressure" >
.
.
.
</acoustic>
</pdeList>


2) Pertubed Convetcive Wave Equation (PCWE) 1 2:

The PCWE rests upon a generic splitting of physical quantities of incompressible Navier-Stokes equations. Because of this splitting, it is possible to derive pure acoustic variables. Consequently, perturbation equations can be used to evaluate acoustic results even in flow regions, where aerodynamic and acoustic quantities are superposed in the CFD simulation. The PCWE is also in inhomogeneous wave equation

\frac{1}{c_0^2} \frac{D^2 \psi^{\mathrm a}}{D t^2} - \nabla \cdot \nabla \psi^{\mathrm a} = -\frac{1}{\rho_0 c_0^2} \frac{D p^{\mathrm{ ic}}}{D t} \, .

The source term of this inhomogeneous wave equation is the substantial derivative of the incompressible pressure considering the mean flow and is also derived from CFD simulations, or from CFSdat. Deriving the weak formulation is done the same way as presented above.

For computing a simulation with the PCWE formulation the acoustic potential formulation of the wave equation has to be used.

     <pdeList>
<acoustic formulation="acouPotential" >
.
.
.
</acoustic>
</pdeList>


Using Aeroacoustic Source Terms

As stated before, aeroacoustic source terms have to be applied as rhsDensity, rhsDensityVectors oder rhsValues. Therefore, the source terms have to be stored at the loaded input mesh.

  <rhsDensity name="">
<grid>
<defaultGrid quantity="acouRhsDensity" dependtype="GENERAL">
<globalFactor>1</globalFactor>
</defaultGrid>
</grid>
</rhsDensity>


As name, the name of the region on which the source terms are located has to be stated. If there are source terms on multiple regions, the rhsDensity as above has to be added multiple times to the xml file. As global factor, spatial and temporal blending functions can be defined. A temporal blending function is especially important for transient simulations, to smoothly couple in source terms. The spatial blending function is especially important when non-conforming interfaces are used. Directly at the non-conforming interface no source terms should be applied 1 7.

In all blue regions the blending function set's the source terms to zero, whereat yellow denotes the regions where the full source term is applied. Blending functions can be superposed by multiplying them. A good working temporal blending function would be

 <rhsValues name="">
<grid>
<globalFactor>
((t lt 1e-3)? (1-cos(0.5*pi/5e-4*(t-5e-4) )^2) : 1) * ((t lt 5e-4)? 0 : 1)
</globalFactor>
</defaultGrid>
</grid>
</rhsValues>


For the very first 0.5 ms all source terms are suppressed. This is to remove any numerical artifacts from the beginning of the simulations, like time derivatives with not enough time steps for the complete derivation stencil. In the second 0.5 ms, the sources are increased smoothly to their actual value 7.

Background flow¶

openCFS allows to consider background flows 8. Therefore, a background flow has to be defined

<flowList>
<flow name="Flow">
<comp dof="x" value=""/>
<comp dof="y" value=""/>
<comp dof="z" value=""/>
</flow>
</flowList>


and added to the region List of the defined PDE

<acoustic formulation="">
<regionList>
<region name="" flowId="Flow"/>
</regionList>
</acoustic>


This added flow velocity $\bar{\mathbf v}$ is considered by the substiantial time derivative

$$\frac{D}{Dt} = \frac{\partial}{\partial t} + \bar{\mathbf v} \cdot \nabla \, .$$

Such background flows can also be read from an external grid, or the grid on which the source terms are defined.

Temperature fields¶

Just like background flows it is also possible to add temperature fields that changes the speed of sound that is indirectly defined in the material file. The temperature field can be defined for each PDE formulation.

<temperatureList>
<temperature name="Temp" value="150">
</temperature>
</temperatureList>


The defined name id of the temperature field has to be added to the region List of the PDE again

<acoustic formulation="acouPressure">
<regionList>
<region name="" temperatureId="Temp"/>
</regionList>
</acoustic>


Again the temperature field can be read from an external grid, or the grid on which the source terms are defined. However, the temperature field can only be defined for the pressure formulation.

Tutorials¶

Acoustics Cylindrical Wave

References¶

1. Stefan Schoder, Clemens Junger, and Manfred Kaltenbacher. Computational aeroacoustics of the eaa benchmark case of an axial fan. Acta Acustica, 4(5):22, 2020.

2. Stefan Schoder, Michael Weitz, Paul Maurerlehner, Alexander Hauser, Sebastian Falk, Stefan Kniesburges, Michael Döllinger, and Manfred Kaltenbacher. Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation. The Journal of the Acoustical Society of America, 147(2):1179–1194, 2020.

3. Stefan Schoder and Manfred Kaltenbacher. Hybrid aeroacoustic computations: state of art and new achievements. Journal of Theoretical and Computational Acoustics, 27(4):1950020, 2019.

4. Michael Weitz, Stefan Schoder, and Manfred Kaltenbacher. Numerical investigation of the resonance behavior of flow-excited helmholtz resonators. PAMM, 19(1):e201900033, 2019.

5. M. Kaltenbacher. Numerical Simulation of Mechatronic Sensors and Actuators: Finite Elements for Computational Multiphysics. Springer Berlin Heidelberg, 2015. ISBN 978-3-642-40169-5. URL: https://www.springer.com/de/book/9783642401695

6. Manfred Kaltenbacher and Sebastian Floss. Nonconforming finite elements based on nitsche-type mortaring for inhomogeneous wave equation. Journal of Theoretical and Computational Acoustics, 26:1850028, 2018. doi:10.1142/S2591728518500287

7. Clemens Junger. Computational aeroacoustics for the characterization of noise sources in rotating systems. PhD thesis, Technische Universität Wien, Wien Austria, August 2019.

8. Manfred Kaltenbacher and Stefan Schoder. Computational aeroacoustics based on a helmholtz-hodge decomposition. Technical Report, SAE Technical Paper, 2018.