A Frequency-Domain Linearized Euler Model for Noise Radiation

like in the case of realistic engine geometries, no exact solution can be found and the problem must be solved numerically.


Introduction
Aeroacoustics is the branch of Fluid Mechanics studying the mechanism of generation of noise by fluid motions and its propagation. The noise generation is associated with turbulent and unsteady vortical flows, including the effects of any solid boundary in the flow. Experimental studies in this field are very difficult, requiring anechoic wind tunnels and very sensitive instruments able to capture high frequency, low amplitude, pressure fluctuations. Computational Aeroacoustics (CAA) can be a powerful tool to simulate the aerodynamic noise associated to complex turbulent flow fields. As sound production represents only a very minute fraction of the energy associated to the flow motion, CAA methods for acoustic propagation have to be more accurate compared to the solution schemes normally used in Computational Fluid Dynamics (CFD). A direct approach to aeroacoustic problems would imply to solve numerically the full Navier-Stokes equations for three-dimensional, unsteady, compressible flows. Sound would then be that part of the flow field which dominates at large distances from the region characterized by intense hydrodynamic fluctuations, propagating at the local sound speed. The direct noise simulation is hardly achievable in practice, except for very simple configurations of academic interest, and in a limited region of space. Even if the small amplitude of the fluctuations allow a linearization, the equations of acoustic disturbances on an arbitrary base flow are very complicated and their solution is not straightforward.
To solve aeroacoustic problems of practical interest some simplifying approximations are necessary. One way to obtain realistic solutions is provided by the hybrid approach. It decouples the computation of the flow field from the computation of the acoustic field. When using a hybrid method the aeroacoustic problem is solved in two steps: in the first step, the hydrodynamic flow field is solved using a CFD method, then the noise sources are identified and the acoustic field is obtained. Extraction of noise sources from the fluid dynamic field can be done using an aeroacoustic theory such as the Lighthill's analogy [Crighton (1975); Goldstein (1976)]. A hybrid approach is based on the fundamental assumption that there is a one-way coupling of mean flow and sound, i.e., the unsteady mean flow generates sound and modifies its propagation, but sound waves do not affect the mean flow in any significant way. This assumption is not so restrictive, because acoustic feedback is possible only when the mechanical energy in the unsteady mean flow is weak enough to be influenced by acoustic disturbances. This occurs principally in the vicinity of a starting point for flow instability (for instance, upstream edges of cavities or initial areas of shear layers). Since the fluid-dynamic 6 www.intechopen.com field and the acoustic field are computed separately, numerical accuracy for the mean flow simulations used as an input of hybrid methods is less critical than in direct computation. Simpler, more flexible and lower-resolution schemes are applicable provided that numerical dissipation is carefully controlled to prevent the artificial damping of high-frequency source components. Incompressible flow solutions can be adequate for evaluating acoustic source terms based on the low Mach numbers approximation. Time-accurate turbulence simulation approaches such as DNS, LES, DES and unsteady RANS methods can be used to compute the space-time history of the flow field, from which acoustic sources are extracted. Because of the high computational cost of the time-accurate simulations, there have been efforts to use steady RANS calculations in conjunction with a statistical model to generate the turbulent acoustic terms. Once the acoustics sources have been evaluated, the generated noise has to be propagated in the surrounding region with linearized propagation models. The main focus of the present chapter is the description of a computational method for noise propagation in turbomachinery applications. In the next Section a linearized model is presented. Section 3 describes the numerical algorithm based on a Discontinuous Galerkin approximation on unstructured grids, and in Section 4 several applications are presented.

Governing equations
In principle the propagation of acoustic waves could be directly studied using the equations of the fluid motion, i.e. the Navier Stokes equations. However, it is possible to introduce some approximations in the Navier Stokes equations in order to obtain equations more suitable for aeroacoustics. At frequencies of most practical interest, viscous effects are negligible in the acoustic field because the pressure represents a far greater stress field than the viscous stresses. Moreover, these disturbances are always small, also for very loudly acoustic waves. The threshold of pain, i.e. the maximum Sound Pressure Level (SPL) which a human can endure for a very short period of time without the risk of permanent ear damage, is equal to 140 dB, which corresponds to pressure fluctuations of amplitude equal to where p ref is the reference pressure corresponding to the threshold of hearing at 1 kHz for a typical human hear. For sound propagating in gases it is equal to p ref = 2 × 10 −5 Pa. The atmospheric pressure of the standard air is equal to p 0 = 101325 Pa, which is 10 3 greater than the pressure variation associated with an acoustic wave at the threshold of pain, i.e., p ′ /p 0 = O 10 −3 , where the superscript (.) ′ denotes acoustic quantities and the subscript (.) 0 denotes mean flow quantities. The corresponding density fluctuations of a progressive plane wave are also of the order of 10 −3 , because in air ρ 0 c 2 0 /p 0 = γ = c p /c v = 1.4 with c 0 being the speed of sound. These estimates demonstrate that the flow perturbations involved in acoustic waves are very small compared to the mean-flow quantities: the acoustic field can be considered as a small perturbation of the mean flow field. Therefore it is possible to linearize the equations of motion. Considering acoustic waves as a perturbation of the mean flow field, defining and assuming small perturbations, it is possible to obtain the equations for the propagation of the sound waves, i.e., the Linearized Euler Equations (LEE). For a two-dimensional problem and a steady mean flow field, LEE are formulated as T is the acoustic perturbation vector, F x and F y are the fluxes along x and y directions respectively, H contains the mean flow derivatives and S represents the acoustic sources. The fluxes, F x and F y , and the term H have the following expressions It is evident from Eqs.(5) that, in order to solve the LEE, the mean flow field must be known in advance. For turbomachinery tonal noise propagation it is better to express the LEE in a cylindrical coordinate system. Given the Cartesian coordinate system (x, y, z), the cylindrical system (r, z, θ) is defined as With respect to this reference frame, the LEE for a three-dimensional problem read where T is the acoustic perturbation vector expressed in the cylindrical coordinate system, i.e. u ′ , v ′ , and w ′ are the velocity components in (z, r, θ) directions respectively. F AX z , F AX r , F AX θ are the fluxes along z, r, and θ directions respectively, H AX contains the terms due to the cylindrical reference frame and to the mean flow derivatives and S AX represents the acoustic sources. In Section 4.3 it will be shown that a generic turbomachinery tonal wave can be expanded in a sum of complex duct modes, having the formf (z, r, t) · exp (Imθ) where m is an integer number which identifies the azimuthal mode and I the imaginary unit. Therefore, using the dependence of the acoustic field on θ, the problem can be reduced, from a three-dimensional problem, to a two-dimensional one in (r, z). For a single duct-mode the LEE become where the superscript( .) reminds that the variable comes from a mode expansion. Assuming that the mean flow is axial-symmetric, i.e. the azimuthal component of the mean flow velocity is zero, w 0 = 0, the fluxes,F AX z ,F AX r , andF AX θ have the following expressionŝ whereas the termĤ AX is given bŷ

Frequency domain approach
The linearized Euler equations, beside acoustic waves, support also instability waves that, for a mean flow with shear-layers, are the well-known Kelvin-Helmholtz instabilities. In the complete physical problem this instabilities are limited and modified by non-linear and viscous effects. Indeed, in the linearized Euler equations, these two effects are not present. Therefore when solving LEE in presence of a shear-layer type mean flow, Kelvin-Helmholtz instabilities can grow indefinitely as they propagate down-stream from the point of introduction and the acoustic solution may be obscured by the non-physical instabilities [Agarwal et al. (2004); Özyörük (2009)]. By using a Fourier decomposition of the acoustics sources and solving the linearized Euler equations in the frequency domain one can, in principle, avoid the unbounded growth of the shear-layer type instability, since the acoustic and instability modes correspond to different values of complex frequency [Rao & Morris (2006)]. However, this could be accomplished in practice only if the discretized form of the equations is solved using a direct solver. The use of iterative techniques to solve the resulting global matrix has been discussed by Agarwal et al. (2004). It is proved that the use of any iterative technique to solve the global matrix is equivalent to a pseudo-time marching method, and hence, produces an instability wave solution. Therefore, the solution of the global matrix needs to be sought by using direct methods such as Gaussian elimination or LU decomposition techniques.

GTS-like approximation
In order to reduce computational time and memory requirements, the pressure gradients in the momentum equations are neglected. A similar approximation, termed Gradient Terms Suppression (GTS), is often used to overcome instability problems that prevent convergence of time domain algorithms for the LEE [Tester et al. (2008); Zhang et al. (2003)]. While the original GTS approximation suppresses all mean-flow gradients, which are likely to be small in the considered subsonic flows, in the present case, being interested in reducing the computational time, only the density mean-flow gradients in momentum equations are neglected. This allows to decouple the continuity equation and to solve only momentum and energy equations. For an axial-symmetric problem the number of total unknowns is thus reduced by a factor T = 5/4, whereas the non-zero terms in the coefficient matrix of the linear system associated with the discretized form of Eqs. (8) is reduced by a factor T 2 . Indeed, a smaller linear system can be solved faster, and, more important, its resolution requires less memory.

Boundary conditions
When a problem is solved numerically, the governing equations must be solved only within the domain, whereas on its borders appropriate conditions, called boundary conditions, must be imposed. In many situations the boundary conditions associated with the continuous problem do not completely supply the discrete problem, and numerical boundary conditions must be added.

Rigid Walls
If walls are assumed impermeable and acoustically rigid, no flow passes through the boundary and acoustic waves are totally reflected. Assuming that the mean flow satisfies the slip flow boundary condition, an analogous slip flow condition must be imposed on the velocity fluctuations u ′ · n = 0 , where u ′ is the acoustic velocity and n is the normal vector to the wall. To apply this condition, Eq. (11) is used to express one of the velocity components in terms of the others.

Axial symmetry
When dealing with axial-symmetric problems, the equations could be solved only for r ≥ 0 if an appropriate boundary condition is applied on the symmetry axis. Along that boundary, the acoustic velocity should be aligned with the r = 0 axis, this can be achieved applying a wall type boundary condition.

Far-field boundary
One of the major issues in CAA is to truncate the far-field domain preserving a physically meaning solution. This leads to the necessity to have accurate and robust non-reflecting far-field boundary conditions. A large number of families of non-reflecting boundary conditions has been derived in literature. The most widely used for the Euler equations are the characteristics-based boundary conditions [Giles (1990); Thompson (1990)]. These methods are derived applying the one-dimensional characteristic-variable splitting in the boundary-normal direction. This technique is usually efficient and robust. The main drawback is that reflections are prevented only for waves that are traveling in the boundary-normal direction. Not negligible reflections can be seen for waves that hit the boundary with other angles. Another class of non-reflecting boundary conditions is based on the asymptotic solutions of the wave equation [Bayliss & Turkel (1980); Tam & Webb (1993)]. In this case, the governing equations are replaced in the far field by an analytic solution obtained imposing an asymptotic behavior to the system. These conditions can be very accurate. Unfortunately, the asymptotic solution can be achieved only in a limited number of cases, reducing the applicability of this model to test cases. Another family of non-reflective boundary conditions is composed by the buffer zone technique [Bodony (2006); Hu (2004)]. In this case, an extra zone is added to damp the reflected waves. The damping can be introduced as a low-pass filter, grid stretching or accelerating the mean flow to supersonic speed. The main drawback of these techniques is represented by the increase of the computational cost, as the thickness of the buffer zone could be important to achieve a good level of accuracy. More recently, the Perfectly Matched Layer (PML) technique has been developed as a new class of non-reflective boundary conditions. The basic idea of the PML approach is to modify the governing equations in order to absorb the out-going waves in the buffer region. The advantage of this technique is that the absorbing layer is theoretically capable to damp waves of any direction and frequency, resulting in thinner layers with respect to other buffer zone approaches, with benefits on the efficiency and the accuracy of the solution. Originally proposed by Berenger (1994) for the solution of the Maxwell equations, the PML technique was extended to CAA applying the split physical variable formulation to the linearized Euler equations with uniform mean flow [Hu (1996)]. It was shown that the PML absorbing zone is theoretically reflectionless to the acoustic, vorticity and entropy waves. Nonetheless, numerical instability arises in this formulation, and in Tam et al. (1998) the presence of instability waves is demonstrated. In Hu (2001), it was shown that the instability of the split formulation is due to an inconsistency of the phase and group velocities of the acoustic waves in presence of a mean flow, and a stable PML formulation for the linearized Euler equation was proposed, based on an unsplit physical variable formulation. The PML technique can be seen as a change of variable in the frequency domain, for example, considering the vertical layer, this change of variable can be written as where σ x > 0 is the absorption coefficient and x 0 is the location of the PML/LEE interface.
To avoid instabilities, a proper space-time transformation must be used before applying the PML change of variable, so that in the transformed coordinates all linear waves supported by the LEE have consistent phase and group velocities. Assuming that the mean flow in the absorbing layer is uniform and parallel to the x axis, the proper space-time transform involves a transformation in time of the form [Hu (2001)] where M 0 = u 0 /c 0 . A stable PML formulation for the two-dimensional LEE can be obtained applying the space-time transformation Eq. (13) to Eqs.(8) and then using the PML change of variable Eq. (12) in the transformed coordinates. Expressing the formulation in the original (x, y) coordinates the PML formulation becomes The terms F PML x , F PML y , and H PML of Eq. (14) are defined as follow where α x = 1 + σ x iω , α y = 1 + σ y iω andF x ,F y , andH are the terms of Eq. (5) under the assumption that the mean flow is uniform and parallel to the x axis. The damping constants σ x and σ y have the following expressions where D x and D y are the widths of the absorbing layers in the x and y directions respectively and x l and y l are the positions of the interfaces between the PML region and the physical domain. The maximum value of the damping σ max is usually taken as 2c 0 /Δx and the coefficient β is set to 2 [Hu (2001)]. At the end of the PML domain, no special boundary conditions are needed except those that are necessary to maintain the numerical stability of the scheme. For this reason at the external boundary of the absorbing layer wall boundary conditions are applied.
Applying the same PML formulation to the LEE written for the turbomachinery duct modes, the following system is obtained The , and H PML of Eq. (17) are defined as follow where α z and α r are defined as in the two-dimensional case andF z ,F r ,F θ , andH are the terms of Eq. (9) and Eq. (10) under the assumption that the mean flow is uniform and parallel to the z axis.

Acoustic inlet
The PML formulation is also used to impose incoming waves at acoustic inlet boundaries. On those boundaries incoming waves should be specified, but at the same time outgoing waves should leave the computational domain without reflections. This can be achieved applying the PML equations to the reflected wave, u re [Özyörük (2009)], which can be expressed as the total acoustic field, u, minus the incoming prescribed acoustic wave, u in Considering the two-dimensional problem and substituting Eq. (20) into Eq. (14) the equation for the inlet PML domain reads Since the linearized Euler flux functions are linear, Eq. (21) becomes The same procedure can be used for the acoustic inlet boundaries of the axial-symmetric LEE.

Numerical methods
The numerical solution of the LEE requires highly accurate and efficient algorithms able to mimic the non-dispersive and non-diffusive nature of the acoustic waves propagating over long distances. One of the most popular numerical scheme in CAA is the Dispersion Relation Preserving (DRP) algorithm originally proposed by Tam & Webb (1993). The DRP scheme is designed for Cartesian or highly regular curvilinear coordinates. However, in many practical applications, complex geometries must be considered and unstructured grids may be necessary. One of the most promising numerical scheme able to fulfill all the above requirements is the Discontinuous Galerkin method (DGM or DG method).
The DGM was firstly proposed in the early seventies by Reed and Hill in the frame of the neutron transport [Reed & Hill (1973)]. Since then, the method has found its use in many different computational models. In the last years, in the context of CFD, DGM has gained an increasing popularity because of its superior properties with respect to more traditional schemes in terms of accuracy and intrinsic stability [Cockburn et al. (2000)].
The DG method displays many interesting properties. It is compact: regardless of the order of the element, data are only exchanged between neighboring elements. It is well suited for complex geometries because the expected dispersion and dissipation properties are retained also on unstructured grids. Furthermore in the framework of DGM it is straightforward to implement the boundary conditions, since only the flux needs to be specified at the boundary. The main disadvantage of the DGM is its computational cost. Because of the discontinuous character, there are extra degrees of freedom at cell boundaries in comparison to the continuous finite elements, demanding more computational resources. This drawback can be partially reduced with a static condensation technique and with a parallel implementation of the algorithm, operations which are made easier by the compactness of the scheme [Bernacki et al. (2006)].

Discontinuous Galerkin formulation
The DGM will be initially presented for the scalar problem of finding the solution u of the hyperbolic conservation equation where F (u) is the flux vector, H is the source term and S is the forcing term. Defining a test function vector space, W, the weak form of the problem (23) over the domain Ω consists in finding u ∈ W such that The discontinuous Galerkin formulation is based on the idea of discretizing the domain Ω into a set of E non-overlapping elements Ω e . Introducing the notations the weak form can be rewritten as To obtain an expression which explicitly contains the flux at the element interfaces, the divergence term in Eq. (26) is integrated by parts where n is the outward-pointing normal versor referred to each element edge. For interfaces on the domain borders, the normal flux vector is evaluated using appropriate boundary conditions. In the general case a boundary condition defines the normal flux as F (u) · n = F BC (u) + G BC . On internal interfaces, F (u) · n is evaluated from the values of u. In order for the formulation to be consistent, the normal flux vector evaluated on right side of an internal interface must be equal to minus the normal flux vector evaluated on the left side of the same interface. Since one of the key feature of the DGM is the discontinuity of the solution among the elements, the consistency is not automatically guaranteed by the formulation. Therefore the normal flux F (u) · n is replaced by a numerical flux F R (u) which is uniquely defined no matter of the side on which it is evaluated (see section 3.2). For ease of notation it is convenient to introduce the following definition where F ∂ + G ∂ is equal to F BC + G BC for interfaces on the domain borders and is equal to to F R for internal ones. Furthermore, assuming that the flux vector is a linear function of the unknown, yields where A and A ∂ are two matrices representing the Jacobian of the physical flux and the Jacobian of the numerical flux respectively. Using Eq. (28) and Eqs. (29), the weak formulation reads Given Eq. (30), the discontinuous Galerkin approximation is obtained considering a finite element space, W h , to approximate W. On each element, a set of points called nodes or degrees of freedom is identified. The number and the position of the nodes depend on the type of approximation used. The set of nodes is chosen to be the same on each element, in this way, on element's borders, there is a direct correspondence among the nodes defined on neighboring elements. The nodes are numbered globally using the index j glob = 1, 2, . . . , n glob with n glob being the global number of degrees of freedom. Beside the global numbering, there is a local numbering. On each element the nodes are identified using the index j e loc = 1, 2, . . . , n e loc where n e loc is the number of degrees of freedom of the e-th element. The correspondence between local node numbers and global node numbers can be expressed through a matrix called connectivity matrix j glob = C e, j e loc .
The nodes of the discretization are used to define the finite element space W h : the vector space W h is generated by the Lagrangian polynomials defined on the nodes of the discretization. The variable u ∈ W is therefore approximated in the W h space with an interpolation of its nodal values where u j (t) is the value of u in the j-th global node x j , y j at the time t and Φ j is the Lagrangian polynomial defined on the j-th global node with the property Φ i x j , y j = δ ij i, j = 1, 2, . . . , n glob .
Although in this work Lagrangian interpolation functions are used, other types of interpolation are possible. Considering the vector space W h , the discrete weak form of problem (23) consists in finding u h ∈ W h such that Substituting Eq. (32) into the discrete weak form (34) leads to This equation must hold for every admissible choice of weight functions w h , therefore it is sufficient to test it for the n glob linearly independent functions of a base of W h . In this way it is possible to obtain n glob independent algebraic equations to solve for the n glob unknowns u j . The vector space W h is defined as the space formed by the Lagrangian polynomials Φ i , therefore the functions Φ i form a base for W h . The i-th algebraic equation is obtained Taking the Fourier transform of Eq. (36), the weak formulation associated with the l-th mode can be written as where( .) (l) j is the l-th component of the Fourier transform of (.), ω (l) is the angular frequency of the l-th Fourier mode and I is the imaginary unit. Equation (37) represents the weak-form discontinuous Galerkin model for a scalar hyperbolic problem in the frequency domain. It can also be written in matrix notation as where Solving this linear system it is possible to obtain the nodal values of the l-th Fourier mode, The same formulation can be applied to the vectorial problem of finding the solution u of the system of equations ∂u ∂t where u is a vector of n vars unknowns, F (u) is the flux tensor, H is the source term, and S is the forcing term. For the vector problem the weak formulation consists in finding u ∈ W such that Assuming a linear flux function, i.e., F (u) = Au and F (u) · n = A ∂ u + G ∂ , the Discontinuous Galerkin approximation in the frequency domain leads to the following linear system where u (l) contains the vectorial nodal values of the l-th Fourier mode and the system is defined as with I being the identity matrix.

Interface flux
The flux through an interface has to be uniquely computed, but, due to the discontinuous function approximation, flux terms are not uniquely defined at element interfaces. Therefore, to evaluate the flux at element interfaces, a technique traditionally used in finite volume schemes is borrowed by the discontinuous Galerkin formulation: the flux function F (u) · n of the vector weak form, Eq. (42), is replaced by a numerical flux function, called Riemann flux, F R (u). Arbitrarily designating one element of the interface to be on the left , l, and the other to be on the right, r, the numerical flux depends only on the internal interface state, u l , on the neighboring element interface state, u r , and on the direction n normal to the interface, i.e. F R (u) = F R (u r , u l , n). In order to guarantee the formal consistency of the scheme, F R is required to satisfy the relations which are the consistency and the conservative conditions respectively. In the present work, the Riemann flux F R is approximated by the Osher flux. This approach is based on the diagonalization of the Jacobian matrix [Toro (1999)]. Assuming a linear dependence of the flux function on the unknown u, the flux along the interface normal direction can be written as where A n = A · n. The numerical method is applied to a hyperbolic system, i.e., LEE, which has a diagonalizable Jacobian matrix A n , that is where K is the non-singular matrix whose columns are the right eigenvectors of A n K = K 1 ; K 1 ;...;K n vars ; , and Λ is the diagonal matrix formed by the eigenvalues λ i Given the diagonalization of A n it is convenient to introduce the diagonal matrix formed by the absolute eigenvalues, |Λ|, and the corresponding absolute flux matrix Using the Osher approach the numerical flux F R can be written as Assuming that the Jacobian matrix does not depend upon the unknown and using the hypothesis of linear fluxes, Eq. (47), the numerical flux becomes

Numerical integration
Integrals of Eq. (39) and Eq. (40) can be evaluated numerically for every element of the mesh using Gauss quadrature formulae. However, it is not convenient to evaluate the integrals directly on the generic element: it is easier to transform (or map) every element of the finite element mesh, Ω e , into a reference element,Ω, called master element and perform the numerical integration on this master element. The transformation between Ω e andΩ is accomplished by a coordinate transformation from the physical coordinates (x, y) to the reference coordinates (ξ, η) where m is the number of parameters used to identify the transformation, (x i , y i ) are the global coordinates of the points of the element used in the transformation and Ψ e i denote the interpolation functions used in the transformation. It is important to point out that the functions Ψ e i used for the approximation of the geometry differ from the functions Φ e i used for the interpolation of the dependent variables. In this work linear interpolation functions Ψ e i are used: on each element the number of parameters used to identify the transformation is equal to the number of vertexes, n = n vertex , and the (x i , y i ) points used in the transformation are the vertexes of the element. Once integrands are expressed on the master elementΩ, numerical integration is performed using Gauss quadrature formulae, in the form where M denotes the number of quadrature points, (ξ i , η i ) are the Gauss points and W i denotes the corresponding Gauss weights.

Interpolation functions
The basis {Φ i } are also evaluated over the master elements: they are the Lagrangian polynomials defined on the node set T p = {x i ; i = 1, . . . , N}, where N is the number of nodes in the node set. For rectangular elements the basis are obtained as the tensor product of the corresponding one-dimensional Lagrangian polynomials defined on the Gauss-Lobatto nodes. Given the one-dimensional polynomials φ l (ξ) with l = 1, . . . , N ξ and φ r (η) with r = 1,...,N η , the two-dimensional ones are defined as For triangular elements the Lagrangian polynomials are constructed on a set of nodes which is defined in such a way that the internal-node positions are the solutions of a steady state, minimum energy electrostatics problem, whereas the nodes along the edges are specified as one-dimensional Gauss-Lobatto quadrature points [Hesthaven (1998)].

Static condensation
One of the main disadvantages of using the DGM for solving LEE in frequency domain is the requirement of a huge amount of memory. The method leads to a linear system of equations which, as explained above, has to be solved with a direct solver, thus requiring a great amount of memory. To partially overcome this problem, a static condensation method can be applied. Static condensation allows to assemble and solve a system matrix which contains only the degrees of freedom associated with the element boundary nodes [Karniadakis & Sherwin (2005)]. Distinguishing between the boundary and interior components of the vectors u e and f e using u e b , u e i and f e b , f e i respectively, that is the DGM linear system (38) can be written as In this decomposition the block K b corresponds to the global assembly of the elemental boundary-boundary basis interaction, K c1 and K c2 correspond to the global assembly of the elemental boundary-interior coupling and K i corresponds to the interior-interior coupling. The static condensation of internal degrees of freedom consists in performing a block elimination by a pre-multiplication of the system by the matrix leading to The elemental boundary unknowns can therefore be evaluated solving the linear system From equation (61) it is evident that, using the static condensation, it is possible to assembly and solve a system that contains only the degrees of freedom associated to the boundary nodes. The therm K b − K c1 K −1 i K c2 is the Schur complement of the full system matrix and can be globally assembled starting from the Schur complements of the elemental matrix where the superscript (.) e denotes elemental matrices and M e is a block diagonal matrix which has been formed by the local matrices M e with e = 1, . . . , E. A b is the matrix that performs the scattering from the global boundary degrees of freedom to the boundary degrees of freedom, that is where u e b contains the components of u e b in element e. Similarly, A T b is the matrix which performs the assembly process from local to global degrees of freedom. Once the linear system of Eq. (61) is solved and the elemental boundary solution is known the solution for the interior elemental nodes is given by the second row of Eq. (60), i.e., Since Eq. (60) involves matrix-vector product of known quantities, it can be evaluated locally within every element, leading to

Linear system solver
The discrete problem leads to a complex matrix system where the complex Fourier coefficients of the acoustic fluctuations are the unknowns. As stated in Section 2.1, this system must be solved with a direct method in order to avoid the Kelvin-Helmholtz instabilities. For this purpose the MUMPS (MUltifrontal Massively Parallel Solver) package Amestoy et al. (2006) will be adopted. MUMPS uses a direct method based on a multifrontal approach which performs a direct factorization K = LU or K = LDL t depending on the symmetry of the matrix. In the multifrontal method the factorization of a sparse matrix is achieved through the partial factorization of many, smaller dense matrices (called frontal matrices).

Multi-geometry scattering problem
A typical test-case to assess the ability of an aeroacustical code to resolve complex geometries is the two-dimensional scattering of sound generated by a spatially distributed monopole source from two rigid circular cylinders, as defined in the Fourth Computational Aeroacoustics (CAA) Workshop on Benchmark Problems [Scott & Sherer (2004)]. The scattering problem is presented here in terms of non-dimensional quantities. Assuming a mean flow at rest, variables can be non-dimensionalized using the mean flow pressure p 0 , density ρ 0 , and speed of sound c 0 . To generate a time-harmonic monopole, only the source term in the energy equation has to be different from zero, i.e. S = [0, 0, 0, S e ]. The forcing term S e is a Gaussian function and can be written in a source-centered coordinate system as where ω = 8π, b = 0.2, ǫ = 0.4. The cylinders have unequal diameters (D 1 = 1.0, D 2 = 0.5), with the source located on the x−axis and equidistant from the center of each cylinder. In the (x S , y S )-coordinate system centered on the source, the locations of the cylinders are given by L 1 = (−4, 0), and L 2 = (4, 0). Considering the symmetry of the problem, only the y ≥ 0 half-domain can be considered if an appropriate symmetry boundary condition is applied on the x−axis. To obtain such a symmetry boundary condition it is sufficient to consider the x−axis as an acoustically rigid wall. The physical domain extends for x ∈ [−10, 10], y ∈ [0, 10] and is surrounded by a PML region with a thickness equal to 0.75. The domain is discretized with an unstructured grid, figure (1), of about 27, 000 elements (both triangles and quadrangles) and on each element Lagrangian basis of degree p = 4 are used.  (3) and (4), the RMS of the fluctuating pressure is plotted along the center line and on the surface of the cylinders. Figure (4) shows a very good agreement between the numerical and the analytical solution [Scott & Sherer (2004)]. Computational time is about 5 minutes using 4 cores of a dual Intel Xeon quad-core computer and the calculation requires 5Gb of RAM.

Sound propagation around a high-lift airfoil
To assess the ability of the present method to simulate realistic geometries, the scattering of a monopole source from an high-lift airfoil is considered here . The airfoil is a three element airfoil based on the RA16SC1   Computational time using 6 th order polynomials is less then 10 minutes using 4 cores of a dual Intel Xeon quad-core computer and the calculation requires about 6.5Gb of RAM. To compare the solutions obtained using different polynomial orders, RMS pressure values are extracted on a circle with radius of 0.7m centered in (0, 0), see figure (7(a)). In figure (7(b)) it can be shown that as the interpolation polynomial order increases the directivity converges to the reference solution.

Turbomachinery noise
A turbomachine produces two type of noise: a broad-band noise, associated with vortex shedding, wake turbulence and blade vibrations, and a tonal noise, related to steady aerodynamic blade loading and to blade thickness effects. Indeed, this last kind of noise is generated because the rotor blade rotation and the aerodynamic interaction of rotor blades with stationary vane wakes generates a periodic pressure field. This field is the source of a narrow banded spectrum type acoustic field, known as spinning mode tones. For a duct having a constant and geometrically simple cross section, spinning mode tones can be studied analytically with the duct mode theory, see appendix (6). However, when the duct has a non-constant or non-simple cross section, like in the case of realistic engine geometries, no exact solution can be found and the problem must be solved numerically. While recent research programs brought to significant progress in reducing both the turbomachinery noise generation and the radiation of noise from the intake, there is still a lack of knowledge about the exhaust noise radiation problem and the need to develop accurate models for its prediction. This problem represents a challenge for CAA, due to the fact that radiated sound propagates through the shear layers separating core, bypass and free-stream fields.

Circular duct propagation and radiation
In order to validate the propagation model, an idealized case for the propagation of sound waves inside the exhaust nozzle is studied. This idealized problem studies the propagation of a sound wave inside a semi-infinite circular cylinder and the subsequent radiation of the wave outside the pipe. This idealized problem is the acoustic diffraction by a sound wave propagating out of a rigid semi-infinite cylindrical duct ( figure 8(a)). The radius of the cylinder is equal to r 1 = 1.212 m and the duct wall has no thickness with acoustically rigid inner and outer surfaces. Inside the duct there is a uniform axial mean flow of density ρ 2 , Mach M 2 and speed of sound c 2 . In the outer region the flow is also uniform and axial, with density ρ 1 , Mach M 1 and speed of sound c 1 . There is no shear layer between the two flows, instead they are separated by a vortex sheet. For this problem, termed "Munt problem", the analytical solution has been found by Munt (1977) and it has been subsequently generalized for annular ducts and lined walls [Demir & Rienstra (2006); Gabard & Astley (2006); Rienstra (1984)]. Since the analytical solution is available only for points at great distance from the cylinder exit, it is not possible to compare directly the solution of LEE obtained using the DGM with the analytical one. Instead LEE are solved only in a small computational domain, the near-field domain, and then the far-field solution is evaluated for the near-field one using the three-dimensional integral formulation of the wave equation proposed by Ffowcs Williams and Hawkings [Iob et al. (2010)]. The far-field results are then compared with the analytical solution. The LEE computational domain extends for z ∈ [−2.5 m; 5.5 m] and for r ∈ [0.0 m; 3.9 m] and is surrounded by vertical and horizontal PML layers with a thickness of 0.7 m. This domain is discretized using a uniform structured grid with about 5, 000 rectangular elements with p = 3. The discretization and the boundary conditions for the near-field domain are shown in figure 8(b). The first case considered is the "no-flow" condition: the mean flow is assumed to be at rest both inside and outside of the duct, with a speed of sound equal to c 1 = c 2 = 340.17 m/s and a density equal to ρ 1 = ρ 2 = 1.225 Kg/m 3 . In figure (9(a)) the instantaneous near-field pressure field for a plain wave, mode (0, 1), at a frequency of 956 Hz is reported, whereas figure (9(b)) displays the corresponding Sound Pressure Level (SPL) directivity pattern in the far-field. The directivity is evaluated on an arc having the center defined at the center of the duct exit section and radius equal to r = 46 m. It can be seen that the agreement between the numerical and the analytical solutions is very good.
To study the effect of the mean flow, another flow condition is analyzed: in this condition there is a mean flow velocity inside the duct, with Mach number equal to 0.447, whereas outside the duct the fluid is at rest. In figure 10(a) and 10(b), instantaneous pressure field and directivity for the mode (9, 1), frequency 866 Hz are presented. Also in this case there is a good agreement with the analytical solution. Computational times for all the conditions is about 2 minutes on an AMD Athlon 64 X2 Dual Core Processor and calculations require less than 1Gb of RAM.

Engine exhaust propagation and radiation
To study a more realistic test case, the propagation of duct modes inside a turbofan engine exhaust is considered. Numerical simulations on this geometry are carried out for two flow conditions. In the first condition, termed condition A, the mean flow is at rest with a speed of sound equal to c ∞ = 340.17 m/s and a total density equal to ρ t ∞ = 1.225 Kg/m 3 . In the second one, condition B, the external flow is at rest, c ∞ = 340.17 m/s and ρ ∞ = 1.225 Kg/m 3 , whereas there is a mean velocity both inside the by-pass duct and inside the core duct. Flow properties at duct's inlet are the following: c fan = 353.15 m/s, ρ t fan = 1.327 Kg/m 3 , and M fan = 0.35 for the by-pass duct and c turb = 508.75 m/s, ρ t turb = 0.598 Kg/m 3 , and M turb = 0.29 for the core duct. Like for the "Munt" problem, LEE are solved only in the near-field domain, and then the far-field solution is evaluated using the Ffowcs Williams and Hawkings formulation. The  figure 11). For all the presented cases the incoming wave enters from the by-pass duct and consists in the duct mode (4, 1) at a frequency of 7981.25 Hz, which corresponds to a dimensionless frequency of kr fan = 16.22 where r fan = 0.11 m is the radius of the by-pass duct at the exit plane. Elements are discretized with interpolation polynomials of degree p = 3, therefore the discretization has about 60000 nodes. Near-field instantaneous pressure perturbations for the flow condition A are shown in figure 13(a), whereas the corresponding SPL directivity pattern is shown in figure 13(b).The far-field directivity is evaluated on an arc of radius r = 12 m having the center placed at the center of the by-pass duct exit section. Same results are presented for the flow condition B in figure 13(a) and figure 13(b). It can be shown that the presence of the flow field, in particular the presence of the shear layers between core flow, by-pass flow and external region, deflects the sound waves in the upward direction. When the flow is at rest the maximum far-field SPL level is located at a directivity angle of 21 • , whereas this maximum moves at an angle of 55 • when there is a mean flow inside the by-pass and the core duct.

Conclusions
A numerical model based on the solution of the LEE in the frequency domain with a discontinuous Galerkin method has proved to be a valuable tool for the analysis of the propagation of sound waves. With the LEE it is possible to study acoustic wave propagation in the presence of rotational mean flows. The frequency domain approach suppress the Kelvin-Helmholtz instability waves which pollute LEE solutions in time-domain calculations. The static condensation technique greatly reduces the memory requirements of the DGM and make feasible the frequency domain approach. Moreover each single calculation, limited to a single frequency, is well suited to the exhaust noise radiation problem where the incoming wave can be treated as a superposition of the elementary duct modes. The model is well suited for design and optimization processes. The model has been successfully validated with the analytical solution of the Munt problem. In the case of realistic configurations, the numerical results reproduce the main expected features.

A. Duct mode theory
According to the duct mode theory the acoustic field inside a duct may be expressed by means of a series expansion in a particular family of solutions, called duct modes. Duct modes are interesting mathematically because they form a complete basis by which any solution can be represented. Nonetheless, modes are also interesting from a more physical point of view because they are solutions in their own right, not just mathematical building blocks. The simple structure of the modes can be therefore useful to understand the usually complicated behavior of the total field. Consider an infinite duct with an arbitrary cross section with inside a uniform mean flow with Mach number equal to M 0 and speed of sound equal to c 0 . The duct is aligned with the z-axis and its cross section, A, is on the (x, y)-plane ( figure 14). Assuming a harmonic time dependency, the acoustic field can be written as where ω is the angular frequency and the Fourier modesp andv satisfy continuity and where α 2 is the corresponding eigenvalue. If the duct cross section is annular and the boundary condition is uniform everywhere, the solution of the eigenvalue problem is relatively simple and may be obtained by separation of variables. Considering an annular duct with inner radius equal to h and outer radius equal to a, the eigensolutions consist of combinations of exponential and Bessel functions. In cylindrical coordinates this yields ψ n (r, θ) = U mμ (r) e Imθ , where m = 0, ±1, ±2, . . . is the azimuthal index of the mode and μ = 1, 2, . . . is the radial one. Positive values of m correspond with clockwise rotating modes and negative m with counter-clockwise rotating modes. The radial function U mμ is obtained from where U ′′ mμ and U ′ mμ are, respectively, the second and the first order derivative of U mμ (r). Solving Eq (74), the radial shape of the modes can be evaluated as where J m and Y m denote the Bessel functions of the first and second kind of order m, whereas J ′ m and Y ′ m are their derivatives. Eigenvalues α mμ = α n are evaluated via boundary conditions, imposing that U ′ mμ (a) = U ′ mμ (h) = 0, the value of α mμ corresponds to the μ-th zero of the equation Although technically speaking {α 2 mμ } are the eigenvalues of minus the cross-sectional Laplace operator, it is common practice to refer to α 2 mμ as the radial eigenvalue or radial modal wave number, to m as the circumferential eigenvalue or circumferential wave number, and to k ± n = k ± mμ as the corresponding axial eigenvalue or axial wave number. The radial and axial wave numbers satisfy For an annular duct mode the general solution may therefore be written aŝ where C ± mμ is the amplitude of the mode, U mμ r, α mμ is given by Eq. (75), the radial wave number α mμ is the μ-th zero of Eq. (76) and the axial wave number k ± mμ is given by Eq. (78).

A.1 Cut-off frequency
The z-axis pressure dependence of a duct mode takes two completely different forms depending on whether the driving frequency is below or above a critical value. The reason of this radical change in transmission properties as the driving frequency is gradually swept through this critical value is due to the fact that the nature of the z-variation of the pressure field in a duct behaves, see Eq. (71), like the term e −ik mμ z . When k mμ is purely imaginary, the amplitude of the pressure fluctuations are subjected to an exponential decay (because Im k mμ < 0). Usually the rate of decay is large enough to reduce the pressure intensity to a negligible value in a short distance compared to the duct radius. When k mμ is real, true wave motion propagates in the duct. From the definition of k mμ , Eq. (78), it is easy to see that k mμ is purely imaginary if 1 − M 2 0 α 2 mμ > ω 2 /c 2 0 . The frequency f c which tells apart this two different behaviors is called cut-off frequency Cut-off frequency depends on the geometry of the duct, on the mode of the incident wave and also on the mean flow properties. For any ω there are always a (finite) μ = μ 0 and m = m 0 beyond which 1 − M 2 0 α 2 mμ > ω 2 /c 2 0 , so that k mμ is purely imaginary, and the mode decays exponentially in z. Therefore there is always a finite number of modes with real k mμ . Since they are the only modes that propagate, they are called cut-on. The remaining infinite number of modes, with imaginary k mμ , are evanescent and therefore called cut-off.