Tutorial Computational Methods For Aeroacoustics
Tutorial Computational Methods For Aeroacoustics
Tutorial Computational Methods For Aeroacoustics
Sheryl M. Grace Dept. of Aerospace and Mechanical Engineering Boston University [email protected]
Opening comments
I refer to ANY computational method focussing on the computation of the sound associated with a fluid flow as computational aeroacoustics (CAA). The CAA methods are strongly linked to CFD CAA methods use specific techniques to resolve wave behavior well which makes this different than general computational fluid dynamics (CFD).
Kinds of applications
Rotors Helicopter noise, wind turbine noise Propulsion systems Stator/rotor, jet noise, combustion noise, propellers (underwater) Airframes Cavity noise, high lift wings Automobiles Cavities, mirrors HVAC and piping systems Fans, duct acoustics
Aeroacoustic domains
Solid body
Computational boundary
Direct numerical, large eddy, and detached eddy simulationsDNS/LES/DES Useful for problems where the sound is from turbulence Euler and Linearized Euler Equations (LEE) LEE very popular when viscous effects can be considered 2nd order as a source of sound Splitting methods Based on LEE, applied to specific unsteady fluid-stucture interaction problems. Integral approaches -- need input from something (all above + CFD) Near field computation using some method above, acoustic field computed using an appropriate form of Ffowcs-Williams and Hawkings or the Kirchhoff method Other acoustic propagation methods Solve a wave equation or associated Euler equation numerically
Governing equations for the different approaches Application references (at end of talk) Implementation of methods for CAA applied to LEE and other CFD like methods Dispersion relation preserving method, Pad methods Boundary conditions Integral approaches Ffowcs-Williams and Hawkings, Kirrchoff method Other propagators
DNS
Two - dimensional governing equations in conservative form Cartesian co-ordinate system for a perfect gas
- density p - pressure u - velocity vector eT - total energy - viscous stress tensor T - temperature R - gas constant - Prandtl number - ratio of specific heats q - thermal conduction
Governing equations
DNS LES DES CFD Euler Linearized Euler Splitting method (at end)
Euler Equation
Two - dimensional governing equations in conservative form Cartesian co-ordinate system for a perfect gas
Moving to the conservative form of the energy equation, and noticing that the terms in the energy equation that involve can be replaced by u p, one derives a form of the equation that is useful when defining the linearized Euler equations:
LEE
Two - dimensional governing equations in conservative form Cartesian co-ordinate system
S is comprised of mass, momentum, and heat sources.
The matrix A has 3 distinct eigenvalues. The repeated eigenvalue is associated with an eigenvector describing the entropy wave field, and an eigenvector describing the vortical wave field. The other 2 eigenvalues are associated with the acoustic wave field.
DRP scheme
If one wants to match the dispersion relation in the numerical simulation, then one must match the Fourier transform. The finite difference representation of a derivative takes the form:
Defining the Fourier transform as The transform of the finite difference expression becomes
So we see the approximation to the wave number is DRP schemes formed by minimizing the wave number error, i.e. minimizing
j=0 0 0 0 0
1.45
Waves with wavelengths longer than 4.33 grid spacings will be adequately approximated
DRP 6th order Central 4th order Central 2nd order Central
Pad/compact methods
So the approximate wave number is given by Tridiagonal schemes, 3 pt. stencil l.h.s., 5 pt. stencil r.h.s
Scheme Pad Pad a1 0.25 0.33333333 a2 0 0 b1 0.6666666 4.6666666 b2 0 0.11111111 b3 0 0 order 4th 6th
DRP (4th order) 6th order Central 4th order Central 2nd order Central
Time discretization
One can use a time discretization scheme constructed using the DRP methodology. (Tam and Web) The method requires the addition of some artificial damping because of the existence of spurious short waves. Tam and Dong: J. of Comp. Acoustics, 1:1-30, 1993; Add artificial damping such that the damping is confined to the high wavenumber range. Then the short waves are damped leaving the long waves basically unaffected. Zhuang and Chen: AIAA J. 40(3):443-449, 2002; Use high-order optimized upwind schemes that damp out the spurious short waves automatically. Hu et al, J, of Comp. Phys, 124:177-191, 1996; Use low-dissipation low-dispersion, low-storage Runge-Kutta schemes.
LDDRK (Hu)
Consider the time discretization
applied to the one-dimensional wave equation where from the spatial transform discussed earlier so that
LDDRK cont.
Dissipation error
Dispersion error
Classical Rung-Kutta matches the expansion for e-i So that 1 = 1, 2 = 1/2!, 3 = 1/3!, 4 = 1/4! would give a fourth order approximation
Specifies range of
Here m chosen to minimize and to satisfy the stability limit |r| <= 1
LDDRK (cont.)
3
0.1629978 0.166558 1/3!
4
0.0407574 0.0395041 1/4!
5
0.00781071 0.00781005
6
0.00132141
Order
4 5 6
accuracy limit
stability limit
Hu et al. discusses the equivalence of this method to minimizing the error in the dispersion error. Also, the implementation of boundary conditions is discussed. Bogey, Bailly, J. of Comp. Phys, 194:194-219, 2004 is another good example of applying this method to define a LDDRK.
1-|r|
0.6
0.75
2nd order 5 stage LDDRK 2nd order 4 stage LDDRK 2nd order 4stage classical RK Exact (dashed)
Boundary conditions
One must set computational boundaries that draw a line between a flow region of interest and other regions that are to be neglected. These boundaries must not produce unrealistic reflections or spurious solutions BC type Characteristicbased Asymptotic ** Pro Straightforward and robust Con Inaccurate for wave angles not perpendicular to the boundary BC surface must be in far-field, not always applicable Large zone, may produce some reflections at interface Stability has been an issue, seems to be better now. Set up for linearized eqs.
Buffer zone
Quite familiar to CFD community Absorbs well with smaller zone, no reflections
BCs cont.
Consider the asymptotic boundary conditions (and BCs that will apply be applied in the linearized region of the flow field) For inviscid, nonheat-conducting, calorically perfect gas, one can decouple the equations into equations that govern the vorticity fluctuation, the entropy fluctuation, and the pressure fluctuations.
Radiation boundary condition radiation boundary condition outflow boundary condition
acoustic wave source
inflow
Good reviews: Colonius, Annu. Rev. Fluid Mech. 36:315-45, 2004 Givoli, J. Comp. Phys, 94(1):1-29, 1991 Hagstrom , Acta. Numerica, 8:47-106, 1999
Asymptotic BCs
Tam and Web : LEE, uniform flow in x-direction only, simple starting point.
Transformed system
Already mentioned that A has 3 distinct eigenvalues giving rise to 4 eigenvectors, X1, X2, X3, X4 where X1 can be associated with the entropy wave, X2 can be associateed with the vorticy wave, and X3 & X4 can be associated with the two modes of the acoustic waves.
Radiation BC
(Tam/Web)
can be solved in terms of the eigenvectors: Focus on the two acoustic modes, in particular, the outgoing mode at a boundary. Transform back to space and time using
Radiation (asymptotic) BC
Other conditions that are similar Bayliss, Turkel, J. Comp. Phys. 48:182-199, 1982 Hagstrom, Hariharan, SIAM J. Sci. Comput., 25(3):1088:1101, 2003 (high order!)
Outflow (asymptotic) BC
At the outflow, entropy, vorticity, and acoustic waves must traverse the boundary. Use the same process of evaluating the behavior of the entropic and vortical parts of the solution. The density perturbation is associated with the entropic mode and the vortical mode is associated with the perturbation velocity vector.
BCs cont.
Tam and Dong: J. of Comp. Acoustics. 4(2):175-201, 1996; Extended the work of Tam and Web to multidimensional nonuniform mean flow. The boundary condition equations have the same form excep that
PML
Computational domain
Add absorbing layer at end of computational domain where the absorption quantities are based on the plane wave solutions of the linearized Euler equations. Create the absorption quantities such that all three wave types are absorbed in the appropriate area. PML differs from the regular buffer zone technique in that the equations used in the added region will not cause ANY reflection when entering the region at any frequency and angle of incidence. Hu uses the nondimensional form of LEE, but follows the same process of forming the plane wave solutions (based on the dispersion relations) discussed previously. One forms but in nondimensional quantities, is 1.
and then find the eigenvalues of A which when set to zero give the dispersion relations.
PML cont.
Combined acoustic eigenvalues give y
wave angle
entropy
vortical
A single Fourier/Laplace component of U is still formed from where X_j are the eigenvectors
PML cont.
Absorption coefficients are then introduced using a splitting method
absorption coefficients
Introduction of scaled spatial parameters, allows one to write the acoustic mode for a given Fourier/Laplace component in such a way that the damping is clear.
There are similar expressions for the entropy and vortical modes
PML cont.
It was shown that there is an instability arising due to convective acoustic waves that have a positive group velocity but a negative phase velocity in the x-direction.
Another spatial transformation is used to overcome this instability, and the final equation that one uses in the PML becomes
0 in vertical layer
0 in horizontal layer
vertical layer x=0
PML cont.
The absorption coefficients can be varied gradually, for example (given in Hus paper): m x = 2 (x is the grid size) =2
location where PML domain starts
Outer edge of the PML can use characteristic, asymptotic, or even very simple reflective type boundary conditions.
yl
xl
Flow field quantities are known in a region near the source, use the integral approaches to find the propagation of the acoustics to the farfield.
Lighthills Eq.
creation of sound generation of vorticity refraction, convection, attenuation, known a priori Lighthill stress tensor mean speed of sound mean density
Colonius, Freund, AIAA J. 38(2):368-370. Applied to jet flow, coupled to a DNS. Oberai et al. (AIAA J., 40(11):2206-2216, 2002) - airfoil selfnoise, coupled to FEM LES
FWH
Differential form of the FWH Eq.
Turbulence is moving Two distinct regions of fluid flow Solid boundaries in the flow
FWH cont.
Integral form of the FWH Eq. Quadrupole
If S shrinks to the body dipole = fluctuating surface forces monopole = aspiration through the surface
There is a formulation for permeable (or porous) surfaces developed for use with CFD where the surface has to be placed
quite close to the body, but not on the body Paper to appear in AIAA J (in near future) - A. Morgans et al. CFD + permeable surface FWH for transonic helicopter noise
Brentner, Farrassat, Progress in Aerospace Sciences 39:83-120, 2003 Great review/overview of use of FWH in rotor noise studies Gloerfelt et al., JSV 266:119-146, 2003, FWH and porous FWH for 2D
cavity problem coupled to DNS Acoustic part computed in the frequency domain (no need for retarded time variable this way.)
All sources of sound and nonuniform flow regions must be inside the
surface of integration. Integration surface must be placed in the linear region of the flow.
FWH is same if the surface is chosen as it is for the Kirchhoff method FWH superior Based on the governing equation of motion (not wave equation) Valid in the nonlinear region Both methods when used in the linear regime, may capture a lower
maximum frequency (CFD method may use grid-stretching).
Brentner, Farrassat AIAA J. 36(8):1379-1386, 1998 Compare FWH and Kirchhoff FWH does separate contributions to the noise (if surface is placed
closed to or on body)
Gloerfelt et al. also compares the two methods (for the cavity problem) Lyrintzis gives a great review of coupling CFD to FWH and Kirchhoff
in Int. J. of Aeroacoustics 2(2):95-128, 2003
Uzun et al., AIAA Paper No. 2004-0517, shows LES coupled to open
FWH and Kirchhoff surfaces for jet (meaning jet outflow not enclosed by the surface - see figure above)
Rahier et al. Aero. Sci. & Tech. 8:453-467,2004, open vs. closed
surfaces for jet analysis. Open surface makes more sense.
Alternative couplings
Solve the linearized Euler equations (with an additional term) using near-field DNS or LES as boundary information. Additional term drives the density towards the Navier-Stokes value. Applied to jet M = 0.9, JFM 438:277-305,2001)
DNS examples
Gloerfelt, Bailly, Juve (JSV 266:119-146, 2003) - subsonic cavity Use DRP to discretize equations Use non-reflecting boundary conditions + absorbing layer Couple to an integral approach
Colonius, Freund, Lele (AIAA J, 38:2023, 2000) - supersonic jet Use Pade methods for discretization Use non-reflecting boundary conditions
LES examples
Bogey, Bailly, Juve (Theor. Comp. Fluid Dyn, 16:273-297,2003) - jet Use DRP to discretize equations Use non-reflecting boundary conditions Uzun, Lynrinztis, Blaisdell (AIAA Paper No. 2004-0517) - jet Use DRP to discretize equations Use non-reflecting boundary conditions Couple to an integral approach Sheen, Meecham (ASME Fluids Div, Sum. Mtg, 2:651-657, 2003) - jet Coupled to an integral approach Oberai et.al. (AIAA J., 40(11):2206-2216, 2002) - incompressible airfoil Use finite element incompressible LES Use non-reflecting boundary conditions Coupled to an integral approach
Euler, examples
Main applications found in literature -- rotor type simulations where flow disturbance is periodic and dominant in the creation of sound.
LEE, examples
Florea, Hall (AIAA J, 39(6):1047-1056,2001) - cascade/gust Low-order discretization, finite volume Bailly, Juve (AIAA J, 38(1):22-29,2000) - apps. DRP scheme Longatte, et al (AIAA J, 38(3):389-394,2000) - sheared ducted flow Lim, et al (JSV, 268(2):385-401,2003) - diffraction from impedance barriers High order discretization Ozyorok, et al (JSV, 270(4-5):933-950,2004) - turbofan noise Chen, et al (JSV, 270(3):573-586,2004) - sound from unflanged duct Mankbadi et al (AIAA J, 36(2):140-147, 1998) - jet
Splitting (LEE-based)
The vortical part of the velocity still satisfies exactly as it did in the original splitting, but the boundary conditions are now defined so that there is no singularity
the boundary condition along the surface is the jump of the potential velocity in a wake must be 0 far upstream:
Splitting, examples
9 The vortical part is first solved analytically or numerically, and then the potential part is found numerically. 9 Most often these problems are computed in the frequency domain. 9 As shown on last slide, valid for subsonic, nonswirling flows
Scott, Patrick/Grace, Atassi, (JCP 119(1):75-93, 1995, AIAA J. 31(1):12-19, 1993) - airfoil/gust Low order-finite difference Coupled to an integral approach Fang, Atassi, (JFE 115:573-579, 1993) - cascade/gust Low order-finite difference Novel non-reflecting boundary conditions Verdon, Hall (AIAA J. 29(9):1463-1471, 1991)- cascade/gust Peake et al. - (JFM 463(25):25-52, 2002) cascade/gust Golubev, Atassi (AIAA J 38(7):1142-1158,2000) - cascades/swirling flow Vortical velocity is no longer the solution to a homogeneous equation
Summary
Weve consider the main features of computational aeroacoustic methods Governing equations - hierarchy of approximations Discretization schemes
Many choices are problem dependent -- makes it difficult to incorporate good acoustic calculations in general CFD type codes If one is implementing these methods, it is good to use the CAA benchmark problems as preliminary method checks.
The End
More questions?
LES
Two - dimensional governing equations in conservative form Cartesian co-ordinate system for a perfect gas Spatially filtered (overbar), Favre (or density weighted) average (tilde) Smagorinksy turbulence model
- density p - pressure u - velocity vector eT - total energy - viscous stress tensor T - temperature R - gas constant - ratio of specific heats q - thermal conduction t - turbulent viscosity T - subgrid scale stress tensor Cs - Smagorinsky constant
DES
Bisseseur et. al. (Aerospace Science Mtg. Proc. pg. 1673-1685, 2004) Use high-order compact difference scheme Couple to an integral approach
CFD, examples
Turbulence is modelled, usually high level of dispersion in discretization
Kim, et al (AIAA Pap. 2003-3202) - general apps FLUENT for near-field Coupled to an integral approach Hendriana, et al (AIAA Pap. 2003-01-1361) - sideview mirror FLUENT for near-field Coupled to an integral approach Grace, Curtis (ASME, IMECE, NCAD 26:103-108, 1999) - cavity FLUENT (URANS) for near-field Coupled to solution of wave equation
Splitting Technique
Ideally (uniform mean flow) the unsteady velocity can be split into solenoidal (vortical/entropic) and irrotational (potential) parts with separate governing equations. The components are linked through boundary conditions. For realistic flows (no shocks or swirl), the velocity components cannot be split as such. Potential governed by a single inhomogeneous, non-constant coefficients, convective wave equation forced by the solenoidal component Vortical part governed by homogeneous, non-constant coeffeicient, convective wave equation Entropic part governed by energy equation.
Used when the disturbance is vortical or entropic not acoustic. Upstream where flow is uniform:
So entropy and velocity upstream are the boundary conditions imposed on the flow.
Splitting cont.
Equations for the split variables are derived from the nonconservative form of the governing equations with the energy equation expressed in terms of entropy.
One can show that the solution to this set of equations can be written as: where
The components of the argument are Lagrangian coordinates of the steady flow fluid particles. The components of X are defined as independent integrals of
Splitting cont.
Solenoidal part of flow is defined and then the single equation:
must be solved, subject to the boundary condition on the surface: and far upstream:
This splitting above leads to singular behavior of the solenoidal part along the solid boundary.
Time discretization
Consider the time discretization:
3 bs chosen so that the Taylors series are satisfied to 3rd order. Leaving one free parameter b0 to minimize dispersion error.
b0 chosen to minimize E1. allows one to adjust the degree of emphasis placed on wave propagation (real part) or damping characteristics (imaginary part). Tam uses a value of 0.36
Time cont.
~ is not unique based on . Spurious numerical solutions exist because
Optimization range has been selected based on the behavior of the ~t is well behaved for values approximated frequency. In particular, less than 0.6, and the optimization ranges from -0.5 to 0.5. For stability, the entire computation must be restricted to the range of t from -0.6 to 0.6. Time step must be chosen to ensure numerical stability:
M is the mean flow Mach number, c0 is the mean flow speed of sound, 0.4 is the value under which all of the roots of are damped. Numerical damping due to the small imaginary part of the approximate frequency. This gives a more stringent requirement on the numerator above. (Details in Tam)
Acoustic field pressure fourth power of velocity Acoustic Power eighth power of velocity
FWH Eq.
Need a more general solution when:
Turbulence is moving Two distinct regions of fluid flow Solid boundaries in the flow Define a surface S by f = 0 that encloses sources and boundaries (or separates regions of interest) Surface moves with velocity V Heavy side function of f : H(f) Rule:
Curle Eq.
When the surface is stationary the equation reduces to Curles Equation
Comments on analogy
Analogy is based on the fact that one never knows the fluctuating fluid flow very accurately Get the equivalent sources that give the same effect Insensitivity of the ear as a detector of sound obviates the need for highly accurate predictions Just use good flow estimates Alternative wave operators that include some of the refraction etc. effects that can occur due to flow nonuniformity near the source have been derived: Phillips eq. , Lilleys eq. Using these is getting close to the direct calculation of sound.
FWH cont.
Multiply continuity and Navier Stokes equation by H Rearrange terms, add and subtract appropriate quantities Take the time derivative of the continuity and combine it with the divergence of the NS equation
Differential form of the FWH Eq. Dipole type term
i Quadrupole
Curles Eq.
Ex. cont.
far-field expansion integration of pressure is lift interchange space and time derivatives
Ex. cont.
our CL vs. t Curle acoustic calc using our CL
Curle acoustic calc using our CL Purely analytic acoustic calc (based on analytic CL
because
When control surface moves with the coordinate system becomes ratio of the area elements of the surface S in the two spaces
Square brackets indicate evaluation at the retarded time e > 1 for approaching subsonic source < 1 for receding subsonic source
Doppler shift
Far-field form
From before. Pressure goes as fourth power of velocity and power as eighth power of velocity
Moving turbulence
Pressure goes as scaled fourth power of velocity Power goes as eighth power scaled by
Example
Method based on : Unsteady CFD -> forced wave equation solved numerically Goal : Make use of existing CFD through a hybrid method for computational aeroacoustics (i.e. no integral formulation for the acoustics) viscous acoustic
t=0.5168 s
t=0.5247 s
t=0.5326 s
t=0.5405 s