Unified solver for fluid dynamics and aeroacoustics in isentropic gas flows

Seediscussions,stats,andauthorprofilesforthispublicationat:https://www.researchgate.net/publication/319490349 Unifiedsolverforfluiddynamicsand aeroacousticsinisentropicgasflows ArticleinJournalofComputationalPhysics·March2017 DOI:10.1016/j.jcp.2018.02.029 CITATIONS 0 READS 89 4authors: Someoftheauthorsofthispublicationarealsoworkingontheserelatedprojects: ExtensiveUnified-domainSimulationoftheHumanVoiceViewproject Enhancedaccuracycomputationalandexperimentalframeworkforstrainlocalizationandfailure mechanisms(EACY)Viewproject ArnauPontRibas CIMNEInternationalCenterforNumericalMeth… 11PUBLICATIONS10CITATIONS SEEPROFILE RamonCodina UniversitatPolitècnicadeCatalunya 224PUBLICATIONS5,741CITATIONS SEEPROFILE JoanBaiges UniversitatPolitècnicadeCatalunya 41PUBLICATIONS292CITATIONS SEEPROFILE OriolGuasch UniversitatRamonLlull 70PUBLICATIONS506CITATIONS SEEPROFILE AllcontentfollowingthispagewasuploadedbyJoanBaigeson05September2017. Theuserhasrequestedenhancementofthedownloadedfile. Unified solver for fluid dynamics and aeroacoustics in isentropic gas flows Arnau Ponta,∗, Ramon Codinaa, Joan Baigesa, Oriol Guaschb aInternational Centre for Numerical Methods in Engineering. C/ Gran Capita, s/n, Campus Nord UPC, 08034 Barcelona, Catalonia, Spain bGTM Grup de recerca en Tecnologies M`edia, La Salle, Universitat Ramon Llull. C/Quatre Camins 30, 08022 Barcelona, Catalonia, Spain Abstract The propagation of acoustic waves is a physical phenomenon which can only be described taking into account the compressibility of the medium. The high computational cost of solving numerically the fully compressible Navier-Stokes equations, together with the poor performance of most numerical formulations for compressible flow in the low Mach number regime, has led to the necessity for more affordable numerical models for Computational Aeroacoustics. For low Mach number subsonic flows with neither shocks nor thermal coupling, both flow dynamics and wave propagation can be considered isentropic. Therefore, a joint isentropic formulation for flow and aeroacoustics can be devised which avoids the need for segregating flow and acoustic scales. Under these assumptions density and pressure fluctuations are directly proportional, and a two field velocity-pressure compressible formulation can be derived as an extension of an incompressible solver. On the other hand, the linear system of equations which arises from the proposed isentropic formulation is better conditioned than the homologous incompressible one due to the presence of a pressure perturbation term. Similarly to other compressible formulations the prescription of boundary conditions will have to deal with the backscattering of acoustic waves. In this sense, a separated imposition of boundary conditions for flow and acoustic scales which allows the evacuation of waves through Dirichlet boundaries ∗Corresponding author: apon[email protected]c.edu Preprint submitted to Journal of Computational Physics May 23, 2017 without using any tailored damping model will be presented. Keywords: Isentropic flow, Computational Aeroacoustics, Numerical Methods, Finite Elements, Lighthill Analogy 2017 MSC: 76Q05 1. Introduction The compressibility behind the acoustics in Computational Fluid Dynamics (CFD) has been widely treated for several purposes along the history of numerical methods. Towards the 70’s, the artificial compressibility method, [1], was developed with the objective of reducing the computational cost of solv-5 ing the incompressible Navier-Stokes equations in 3D domains, a research field which would also lead to projection methods, known nowadays as fractional step schemes. In this framework, the artificially added compressibility through a density or pressure perturbation term was not only a numerical artifact, but a term that could be easily associated to the acoustics of a low speed com-10 pressible flow. However, the artificial compressibility method did not aim to describe the acoustic scales of the flow, but to introduce a numerical relaxation parameter which allowed an easier fulfillment of the continuity condition. The main modification of the incompressible Navier-Stokes consisted in adding an artificial time derivative of the density or the pressure to the dimensionless con-15 tinuity equation, which improved the condition number of the final system to be solved. A similar method was later applied by [2] to the low speed compressible Navier Stokes equations, in which a time derivative of the primitive variables was added to the energy equation in order to reduce the big disparity between the flow velocity and the sound speed. The Chorin method was extended for20 both incompressible and slow compressible flows by [3] by adding similar terms to all equations in order to obtain a symmetric hyperbolic problem. In other cases such as low Mach number (M) compressible flows, the goal consisted precisely in going in the opposite direction and identifying the acoustic scales of the flow in order to remove them from the problem, [4], because they led to an25 2 ill-conditioning of the system and to the backscattering of sound waves into the computational domain. While the addition of a certain amount of compressibility has made the calculation of incompressible flows easier without taking into account the consequent acoustic field, the inclusion of compressibility in the flow formulation has been30 a drawback for calculating acoustics when dealing with low speed flows. The conservative compressible flow equations are considered the complete representation of the aeroacoustic problem because they describe directly all flow and acoustic scales without any need for modeling, which in terms of Computational Fluid Dynamics (CFD) is called Direct Numerical Simulation (DNS), [5], and35 in acoustics is referred as Direct Noise Computation (DNC), [6]. However, as stated above, this formulation performs poorly for Mach numbers tending to zero due to the huge difference between flow velocity and wave propagation speed, which causes convergence problems. In order to avoid the bad conditioning of the problem, a series of hybrid methods, which segregate the acoustics from40 the CFD, were developed. The so called acoustic analogies resolve the acoustic scales by means of an inhomogeneous wave equation where the source term that represents the aerodynamic noise comes from a previous flow calculation. The pioneer work in this field is presented in [7], which computes sound waves propagation under the hypothesis of far-field flow conditions. The method has45 been progressively extended to include diffraction by solid boundaries [8] and moving surfaces [9]. Other hybrid methods, such as the incompressible-acoustic split method presented in [10, 11] enrich the incompressible flow equations with a variable density linked to pressure perturbations. Then, the time derivative of these perturbed density is translated into isentropic fluctuations of velocity50 and pressure that are propagated using a purely acoustic compressible solver after subtracting the incompressible component of the flow field. In a similar way, some formulations propagate the near field flow information to the far field with the Linearized Euler Equations (LEE), [12, 13, 14] or with the acoustic perturbation equations [15, 16, 17], which consist in an acoustic filtering of the55 LEE source term. All these methods allow a considerable flexibility, for example 3 the use of a different discretization for each problem, as well as different flow and acoustic models. However, these models must be adapted to every case and the approximation errors need to be properly assessed. In some cases, acoustic source terms need to be modeled and might not be straightforward to implement60 in a FEM code. Moreover, the segregated calculation of the flow and acoustic components only assumes a one-way coupling from flow to acoustics, but not the other way around. The formulation proposed in this work aims for a simplification of Computational Aeroacoustics (CAA) of isentropic compressible flows and proposes a65 general framework that can be applied to any geometry, spatial discretization or flow regime below the transonic range. It consists in a compressible formulation with primitive variables without solving for the energy equation, since the flow is considered to be isentropic, which after condensing the density field becomes a system of equations in terms of the velocity and the pressure, like in70 incompressible flow solvers. As a consequence, the implementation cost is very low when one departs from an already implemented incompressible flow solver. Also, the computational cost is reduced with respect to other methodologies due to the following reasons: getting rid of the fully compressible approach and solving only for velocity and pressure, solving all scales at once without acoustic75 analogies and improving the condition number of the system the incompressible problem. The only drawback of such a compact system will be, of course, the lack of visualization of the acoustic fluctuations at the near field, where the aerodynamic scales are totally dominant and the wave propagation cannot be extracted like in [18] or [19]. As in all compressible flow models, an adequate80 equation of state needs to be chosen, in this case relating only density and pressure. Since the present paper aims at solving both aerodynamics and acoustics scales in a single calculation, the prescription of compatible and accurate boundary conditions for both components of the solution has been an important moti-85 vation for this work. From a numerical point of view, the imposition of boundary conditions can be performed as in the incompressible case, avoiding the difficul4 ties found in compressible flows. However, omitting the acoustic scales in the treatment of the external boundaries leads to undesired wave reflections which affect the accuracy and the stability of the unified solver. Therefore, a new90 method including the combined imposition of essential boundary conditions in a weak sense on the mean flow variables, [20], and a Sommerfeld boundary condition for the acoustic component of the pressure will be presented, [21]. This combination will allow the acoustic wave to leave the domain through boundaries where the flow has been prescribed a certain boundary condition.95 The paper is organized as follows: a detailed presentation of the isentropic compressible equations is shown in Section 2. The details of the aforementioned prescription of boundary conditions are presented in Section 3, and the stabilized time-discrete finite element formulation is derived in Section 4. Finally, numerical results are presented in Section 5: two cases consisting in a 2D flow100 around a cylinder (M= 0.058) and a 3D flow around an airfoil (M= 0.4) will be presented and benchmarked against the Lighthill analogy, [7], with incompressible flow and the Ffowcs Williams Hawkings (FWH) acoustic analogy, [9], with a compressible formulation respectively. This analysis will allow to validate the present method in its whole application range.105 2. Problem formulation The present work focuses in the study of the aerodynamic and acoustic behavior of an ideal gas undergoing a reversible thermodynamical process, which is a realistic hypothesis in most aeroacoustic problems without heat transfer or shocks. This initial assumption allows a drastic simplification of the compressible Navier-Stokes equations, since the energy equation does not need to be solved and the primitive variables of the problem can be used. Moreover, a general formulation can be derived for both slow and high speed isentropic flows taking into account the following equalities, [22]: p0 p=1 + γ−1 2M2γ γ−1 ,(1) 5 ρ0 ρ=1 + γ−1 2M21 γ−1 ,(2) where γis the adiabatic constant of the gas, pand ρare the total pressure and density fields including perturbations caused by the compressibility of the medium, whereas p0and ρ0are the same fields at stagnation conditions, [22]. Mis the Mach number, defined as: M:= |u| c0 ,(3) where |u|is either the modulus of the pointwise velocity (or a characteristic value of it if one wants to define a global Mach number) and c0is the speed of sound in an ideal gas defined as c0=qγRT0 M, where T0is the temperature field at stagnation, R[J/K·mol] is the universal gas constant and M[kg/mol] is the110 molar mass of the gas. From Eq. (2), the following equality between both fields can be easily obtained: p0 p=ρ0 ργ .(4) Then, deriving with respect to time both sides of Eq. (4) and using the equation of state for an ideal gas, p0=ρ0RT0 M, the next expression connecting pressure and density time derivatives can be found: ∂p ∂t =p0 ρ0 γ1 + γ−1 2M2−1∂ρ ∂t =RT0 Mγ1 + γ−1 2M2−1∂ρ ∂t .(5) The final equation of state for a low speed gas flow can then be simplified to ∂p ∂t =c2 0 ∂ρ ∂t +OM2,(6) whereas for a compressible flow the speed of sound cwill have to be computed as follows: c2=c2 01 + γ−1 2M2−1 .(7) and the following equation will hold: ∂p ∂t =c2∂ρ ∂t .(8) 6 The same procedure can be applied to the pressure gradient obtaining the same relationship with respect to the density gradient. This explicit connection between pressure and density variations will allow to greatly simplify the115 compressible Navier-Stokes equations, since the density perturbations will be expressed in terms of the pressure. It is important to highlight that the limit M→0 will lead to a problem which will be very similar to the one resulting from the artificial compressibility method and will contain the acoustic scales of the flow. This is remarkable if it is compared to other non-isentropic formula-120 tions for low Mach numbers (see for instance [23]), where density variations are linked exclusively to temperature oscillations, and as a consequence no acoustics are captured. Let us consider a computational domain Ω ⊂Rd(where d= 2,3 is the number of space dimensions) with a domain boundary Γ = ∂Ω and let (0, T) be the time interval of analysis. The isentropic compressible equations are then: ρ∂u ∂t +ρ(u· ∇)u−µ∇2u−1 3µ∇(∇ · u) + ∇p=0in Ω,(9) ∂tρ+u· ∇ρ+ρ∇ · u= 0 in Ω,(10) where uis the velocity, pthe pressure and µthe dynamic viscosity. Boundary and initial conditions need to be appended to this problem. Using Eq. (5) ρcan be expressed in terms of pand the continuity equation becomes 1 c2∂tp+1 c2u· ∇p+ρ∇ · u= 0 in Ω,(11) where c(x, t) is given by Eq. (7) and xis the spatial coordinate vector. Despite all simplifications, the previous equation still depends on the function cand125 two density dependent terms remain in the momentum equation. Calculating these two fields as implicit functions of (u, p) would increase the complexity of the new scheme with new non-linearities. In this sense, the finite element approximation will include the necessary elimination of these variables in order to obtain a problem depending only on the velocity and pressure fields.130 The next step consists in deriving the variational formulation of the previous problem. Let us denote with h·,·iωthe integral of the product of two functions 7 in the domain ω(with the subscript omitted when ω= Ω) and (·,·) the L2(Ω)- inner product. Let Vand Qbe the functional spaces where for each time tthe velocity and pressure solutions live respectively with appropriate regularity that we will not analyze here. Then, defining the velocity and pressure test functions v∈Vand ρq ∈Qthe variational formulation can be written in terms of the forms: B([u, p],[v, q]) = (ρv, ∂tu) + hρv,(u· ∇)ui+µ(∇v,∇u) (12) +1 3µ(∇ · v,∇ · u)−(∇ · v, p) (13) +1 c2q, ∂tp+1 c2q, u· ∇p+ (ρq, ∇ · u),(14) ˜ BB([u, p],[v, q]) = − hv,n·σ(u, p)iΓ,(15) where Band ˜ Bare two forms and the stress tensor is defined as σ(u, p) = −pI+µ∇u+1 3µ(∇ · u)I. The Galerkin weak form of the problem prior to applying boundary conditions can be written as follows: for all time t > 0, find u∈Vand p∈Q, with appropriate regularity in time, such that: B([u, p],[v, q]) + ˜ BB([u, p],[v, q]) = 0 (16) for all v∈Vand ρq ∈Q. Equations (2) and (7) are used to close the problem. Moreover, initial conditions need to be appended. Boundary conditions will be defined in the following section proposing a new formulation for the form ˜ BB.135 This will give rise to a decomposition of the form ˜ BB=BB−LB, with BB depending on the unknowns and LBon the boundary data, so that it can be moved to the right-hand-side of (15). 3. Imposition of boundary conditions 3.1. Mean and acoustic components140 Although the intricate prescription of boundary conditions of the fully compressible formulation is avoided in the present problem, new challenges arise 8 Figure 2: Reflection of the sound generated by vortices approaching the outflow. where α∗=αρ2c2and r0and rfare the small and big radius of the PML, respectively. Finally, the problem to be solved will be in this case B([u, p],[v, q]) + BP M L([u, p],[v, q]) + BB([u, p],[v, q]) = LB([v, q]) (24) Unlike Section 3, the importance of absorbing both hydrodynamic and acoustic scales on the outlet justifies the application of the PML to the whole variables235 (un+1 h, pn+1 h). The performance of this numerical tool will be presented in the Section 5. 4. Numerical approximation In this section we present the finite element formulation for the space approximation of the isentropic Navier-Stokes equations, including the stabilization240 15 Figure 3: A PML is attached to the original outlet of the domain Ω. terms required for obtaining a stable formulation when using P1/P1 velocitypressure elements, as well as the time discretization using finite differences. Let us consider a finite element partition of the domain Ω of size h, and use this letter as subscript to denote finite element functions and spaces. Only conforming finite element approximations will be considered in what follows.245 Let Vh⊂Vbe the finite approximation space for the discrete velocity field and let us also define Qh⊂Q, the pressure approximation space. 4.1. Time discretization Concerning the time integration, the monolithic approach for solving the incompressible Navier-Stokes equations consists in building a system with both velocity and pressure degrees of freedom, which leads to the coupled calculation of the momentum and mass equations in one single step. To approximate the first order time derivatives, a second order backward finite difference scheme 16 (BDF2) has been used. Let us partition the time interval [0, T ] into Nequal time steps of size δt := tn+1 −tnso that 0 ≡t0< t1< . . . < tn< . . . < tN≡T. Given a generic time dependent function g(t), the following notation will be used for the BDF2 approximation to the first time derivative: ∂tg|tn+1 ≈δtgn+1 := 1 δt 3 2gn+1 −2gn+1 2gn−1,(25) where gndenotes evaluation of gat time step tn. 4.2. Discrete boundary conditions250 At an arbitrary time step of the numerical simulation, the final fully discretized implicit scheme in space and time can be derived using the finite element formulation described below. Moreover, the mean flow values must be expressed according to the chosen integration scheme and the penalty parameters of the weak essential condition on ΓLmust be defined. We do this as255 follows: •As mentioned above, µp,lpcan be taken as µp=µ+|u|h, lp=h, [36]. •If the temporal window presented at Eq. (18) is defined at a discrete level as Tw=Nwδt and we use the trapezoidal rule for the integration, then the mean values can be expressed as follows:260 ¯ un+1 h=δt Tw 1 2un+1 h+ n X k=n−Nw+2 uk h+1 2un−Nw+1 h!(26) Bearing in mind the sharp initial pressure transient and the absence of a minimally developed mean flow, it is important to run several time steps (Nw) before using the present formulation in order to obtain representative mean flow variables. The same procedure is applied to pand the fluctuating components will be also expressed now on in terms of the full variables evaluated at tn+1.265 4.3. Finite element approximation For a better understanding of the derivation, the formulation will be arranged in five forms: B,BB,BP ML,LBand BS, which corresponds to the Algebraic Subgrid Scale (ASGS) stabilization terms and will be presented next. The 17 final formulation reads as follows: from known un−2 h,un−1 hand un h, compute the compressible velocity and pressure at time step tn+1,un+1 h, pn+1 h∈ Vh× Qh, such that B([uh, ph],[vh, qh]) + BP ML([uh, ph],[vh, qh]) + BB([uh, ph],[vh, qh]) +BS([uh, ph],[vh, qh]) = LB([vh, qh]),(27) for all test functions, where B([uh, ph],[vh, qh]) = ρn+1vh, δtun+1 h+hρn+1vh,un+1 h· ∇un+1 hi +µ∇vh,∇un+1 h+1 3µ∇ · vh,∇ · un+1 h−∇ · vh, pn+1 h + 1 (c2)n+1 qh, δtpn+1 h!+ 1 (c2)n+1 qh,un+1 h· ∇pn+1 h!+ρn+1qh,∇ · un+1 h. (28) As mentioned before, the condensation of ρn+1 and cn+1 is essential for keeping the complexity of the formulation low. For this reason, these two variables have been included in non-linearity iterative loop so they are updated at every time step at each Gauss point. In this way, they are evaluated with the converged unknowns of the problem at tn+1, as the implicit scheme requires: ρn+1 =ρ01 + γ−1 2 |un+1 h|2 c2 0γ−1 , c2n+1 =c2 01 + γ−1 2 |un+1 h|2 c2 0−1 .(29) Since it is understood that ρand care only evaluated at tn+1, now on they will be referred as ρand cinstead of ρn+1 and c2n+1.Next, the bilinear form BB and the linear form LBcan be easily obtained using (21): BB([uh, ph],[vh, qh]) = − hvh,n·σh¯ un+1 h,¯pn+1 hiΓL− h¯ un+1 h,n·σh(vh, qh)iΓL +βµp lp hvh,¯ un+1 hiΓL+hcρvh·n,u0 h n+1 ·niΓL∪ΓO, LB([vh, qh]) =βµp lp hvh,uLiΓL− huL,n·σh(vh, qh)iΓL,(30) where we have assumed that tO=0. Applying the definition of the mean values presented in Eq. (26) and expressing the fluctuating components in terms of the 18 problem unknowns, BBcan be rewritten as follows: BB([uh, ph],[vh, qh]) = −1 2Nw hvh,n·σhun+1 h, pn+1 hiΓL −1 2Nw hun+1 h,n·σh(vh, qh)iΓL+1−1 2Nwhcρvh·n,un+1 h·niΓL∪ΓO +1 2Nw βµp lp hvh,un+1 hiΓL+1 Nw βµp lp"n X k=n−Nw+2 hvh,uk hiΓL−1 2hvh,un−Nw+1 hiΓL# −1 Nw"n X k=n−Nw+2 hvh,n·σhuk h, pk hiΓL−1 2hvh,n·σhun−Nw+1 h, pn−Nw+1 hiΓL# −1 Nw"n X k=n−Nw+2 huk h,n·σh(vh, qh)iΓL−1 2hun−Nw+1 h,n·σh(vh, qh)iΓL# −1 Nw n X k=n−Nw+2 hcρvh·n,uk h·niΓL∪ΓO−1 2Nw hcρvh·n,un−Nw+1 h·niΓL∪ΓO. LB([vh, qh]) = βµp lp hvh,uLiΓL− huL,n·σh(vh, qh)iΓL (31) When a PML is mandatory BP ML must be included in the formulation. Using (22) the discrete bilinear form for the PML can be easily derived: BP ML([uh, ph],[vh, qh]) = α∗vh,un+1 hΩP ML +αqh, ρpn+1 hΩP ML .(32) The last step for a robust and consistent formulation consists in developing an appropriate stabilization for B. On the one hand, one can profit from all terms coming from the stabilization of the incompressible Navier-Stokes. In the present case, the Algebraic Subgrid Scale (ASGS) method for incompressible flows presented in [38] has been taken as reference. On the other hand, the two pressure terms must be included in the residual of the continuity equation for consistency and the pressure convective term must be added to the stabilization operator. Since linear P1/P1 elements will be used, the second order viscous terms have been removed from the residual and the stabilization operator of the momentum equation. As it will be shown in the next chapter, this extension of the incompressible flow stabilization also prevents the dissipation of the acoustic 19 waves: BSun+1 h, pn+1 h,[vh, qh]=X K τ1,K ρun+1 h· ∇vh+∇qh, ρδtun+1 h +ρun+1 h· ∇un+1 h−µ∇2un+1 h−1 3µ∇∇ · un+1 h+∇pn+1 hK +X K τ2,K ρ∇ · vh+1 c2un+1 h· ∇qh, ρ∇ · un+1 h+1 c2un+1 h· ∇pn+1 h+1 c2δtpn+1 hK , (33) where Kdenotes the element domain, and τ1,K and τ2,K are suitable stabilization parameters defined in each element, [39]: τ1,K =c1 ν h2+c2 |un+1 h| h−1 τ2,K =h2 c1τ1,K (34) 5. Results For a proper validation of the present formulation two different scenarios have been taken as reference. First, a 2D problem consisting in a low speed Re = 1000 flow around a cylinder has been calculated with the isentropic compressible270 equations for comparing the CFD results and the acoustic propagation to those provided by an incompressible solver and the Lighthill analogy. Second, a M = 0.4 flow around a 3D NACA 0012 airfoil has been calculated in order to evaluate the performance of the formulation against a compressible flow solver and the Ffowcs Williams & Hawkings (FWH) acoustic analogy. The main advantage275 of the isentropic compressible formulation is that it can be treated numerically like the incompressible formulation although the flow regime might not be in the incompressible range anymore. From the point of view of an end user, the only further requirement consists in introducing the three following parameters: the gas universal constant R= 8.31 J/Kmol, the molar mass, the sound propagation280 speed of the working gas and the bulk temperature. In both cases the values of air at room temperature have been considered (M= 28.97 g/mol, c0= 343 m/s and T0= 293.15 K). 20 5.1. Aerodynamic sound radiated by flow past a cylinder. M = 0.0583 The first benchmark case consists in a 2D flow around a D= 0.3 cylinder285 which allows evaluating the aeolian tones of a low Mach viscous flow, [40]. The incident velocity of 20 leads to a Reynolds and Mach numbers at the far field (away from the cylinder) of Re = 1000 and M = 0.0583 for a sound speed of c0= 343 (all units are in SI). The problem has been solved in an unstructured mesh of nearly 1 million triangular linear elements using equal interpolation for290 velocity and pressure, with a size of 3 ·10−3Dnear the cylinder surface. The case has been run up to 1.5 s with a time step δt = 1 ·10−3s, departing from an initial incompressible solution in order to ease the initial convergence of the iterative solver. In fact, the same solver does not yield convergence when the problem is computed with the incompressibility condition, even when departing295 from a converged solution. For the weak imposition of boundary conditions it has been enough taking a penalty parameter β= 1. The original case in [40] was computed with the incompressible Navier-Stokes equations and Lighthill’s analogy in the frequency domain (Helmholtz equation). Regarding the CFD, this calculation provides a shedding frequency of 15.3 Hz,300 whereas the present formulation has obtained a very similar value of 15.6 Hz. The fully developed velocity profiles are compared in Fig. 4. This very good fitting between both velocity profiles does not only assess the accuracy of the isentropic compressible equations at low Mach regimes, but illustrates the possibility of replacing the incompressible monolithic solvers when their convergence305 is not satisfactory even if the acoustics are not relevant. Moreover, it confirms the good performance of the weakly imposed inlet condition. Of course, one may think that, despite this huge benefit, the compressibility brings the big drawback of waves being reflected by the boundaries and polluting the flow solution. However, Fig. 5 shows that this inconvenience is completely resolved310 by the previously presented boundary conditions as no reflections are observed on the external boundaries. This plot also validates qualitatively the acoustic propagation at the far field. The present formulation is capable of capturing the anisotropy of the aeolian tones as well as the amplitude of the acoustic waves. 21 On the other hand, Fig. 6 aims for a quantitative validation of the phenomenon.315 All three curves show a very good fitting in the near field, but some discrepancies appear beyond the 50 m. The stripped curve shows the importance of the extension of the incompressible stabilization terms presented in Eq. (33). If the compressible terms are not taken into account in the residual, the wave propagation is not affected in the region where the incompressible flow scales are320 dominant but causes a drastic dissipation of the waves at the far-field. When the full stabilization is deployed, the resulting curve approximates the reference one more accurately. The phase error can be associated to the shedding frequency differential between both formulations, whereas the small amplitude difference can be explained by the inaccuracy of Lighthill’s analogy in the near-field re-325 gion, where the flow stagnation hypothesis is not fulfilled, and which can lead to an over-prediction of the amplitude of the acoustic signal, [16]. 5.2. Aerodynamic sound radiated by flow past an airfoil. M = 0.4 The second benchmark case consists in a 3D flow around a NACA 0012 airfoil with an angle of attack of 5◦, [41]. The flow Reynolds number based330 on the airfoil chord (d= 0.1524) is Rec= 408000 whereas the incident Mach number is M = 0.4. The problem has been solved in an unstructured mesh of nearly 20 million tetrahedral linear elements using equal interpolation for velocity and pressure, with a size of 4 ·10−4on the leading edge and 6.5·10−4 on the rest of the airfoil surface (all units are in SI). The case has been run up335 to 0.050 s with a time step δt = 10−5s, departing from an initial incompressible solution in order to ease the initial convergence of the iterative solver. For the weak imposition of boundary conditions a penalty parameter β= 125 has been taken. Unlike the previous low-speed flow, the present case generates an airjet that cannot be dissipated before reaching the outlet, for which a PML has340 been placed in this region. On the external boundaries the flow field has been prescribed separately following the presented method. The original case in [41] was computed with a compressible Large Eddy Simulation (LES) for the flow scales and the Ffowcs Williams & Hawkings (FWH) 22 Figure 4: Incompressible flow velocity (top), isentropic compressible velocity (bottom). acoustic analogy for the acoustic component, [9]. The structured finite dif-345 ference mesh has been able to capture high frequencies up to more than the Helmholtz number kd = 40, which is approximately 15 kHz. This spatial resolution, together with a much smaller time step, could not be reproduced with an unstructured mesh of tetrahedral elements under a reasonable cost taking into account the available resources. Since the object of the present work does350 not consist in assessing the performance of the solver in specific mesh typologies or in reproducing all the details of a particular problem, but in establishing a general framework for the calculation of a wide range of flows, the goal of this analysis has been restricted to the following points: the suitability of the present 23 (a) (b) Figure 5: Imaginary component of the acoustic pressure (a), isentropic compressible pressure (b). isentropic compressible formulation for reproducing the mean flow patterns at355 M=0.4, the proper propagation of the captured acoustic modes, and finally the validation of the proposed boundary conditions. The validation of the formulation in regard to the mean flow variables is perhaps the most demanding aspect of this simulation. The lack of resolution for properly capturing the boundary layer around the airfoil makes it difficult to360 evaluate the aerodynamics of the present formulation. In a first approach, the surface of the airfoil was prescribed a non-slip boundary condition, but it yielded an early separation of the boundary layer. In order to prevent this scenario, a wall-law with both buffer and logarithmic regions has been prescribed and the result in Fig. 7 has been obtained:365 Although the mean velocity field values are properly reproduced, the boundary layer still suffers an early detachment from the airfoil. In order to analyse in what extent the mesh element size, and not the formulation, was the reason for this discrepancy, the same problem has been run in a 2D section of the original domain using a much finer mesh. Fig. 7c shows that the element size around370 the wall was indeed the cause of the early boundary layer detachment. The same dependence on the mesh resolution can be found in the capturing of the acoustic modes. However, in this case the lack of accuracy can be restricted 24 [25] V. Granet, O. Vermorel, T. L´eonard, L. Gicquel, T. Poinsot, Comparison of nonreflecting outlet boundary conditions for compressible solvers on515 unstructured grids, AIAA journal 48 (10) (2010) 2348–2364. [26] K. W. Thompson, Time dependent boundary conditions for hyperbolic systems, Journal of computational physics 68 (1) (1987) 1–24. [27] T. J. Poinsot, S. Lele, Boundary conditions for direct simulations of compressible viscous flows, Journal of computational physics 101 (1) (1992)520 104–129. [28] R. Prosser, Improved boundary conditions for the direct numerical simulation of turbulent subsonic flows. i. inviscid flows, Journal of Computational Physics 207 (2) (2005) 736–768. [29] R. Prosser, Towards improved boundary conditions for the dns and les of525 turbulent subsonic flows, Journal of Computational Physics 222 (2) (2007) 469–474. [30] C. S. Yoo, H. G. Im, Characteristic boundary conditions for simulations of compressible reacting flows with multi-dimensional, viscous and reaction effects, Combustion Theory and Modelling 11 (2) (2007) 259–286.530 [31] J.-P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of computational physics 114 (2) (1994) 185–200. [32] F. Q. Hu, A perfectly matched layer absorbing boundary condition for linearized euler equations with a non-uniform mean flow, Journal of Computational Physics 208 (2) (2005) 469–492.535 [33] C. K. Tam, Z. Dong, Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow, Journal of Computational Acoustics 4 (02) (1996) 175–201. [34] C. Bogey, C. Bailly, Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases,540 Acta Acustica united with Acustica 88 (4) (2002) 463–471. 31 [35] M. Juntunen, R. Stenberg, Nitsches method for general boundary conditions, Mathematics of computation 78 (267) (2009) 1353–1374. [36] R. Codina, J. Baiges, Approximate imposition of boundary conditions in immersed boundary methods, International Journal for Numerical Methods545 in Engineering 80 (11) (2009) 1379–1405. [37] H. Takemoto, P. Mokhtari, T. Kitamura, Acoustic analysis of the vocal tract during vowel production by finite-difference time-domain method, The Journal of the Acoustical Society of America 128 (6) (2010) 3724– 3738.550 [38] R. Codina, On stabilized finite element methods for linear systems of convection–diffusion-reaction equations, Computer Methods in Applied Mechanics and Engineering 188 (1) (2000) 61–82. [39] R. Codina, J. Principe, O. Guasch, S. Badia, Time dependent subscales in the stabilized finite element approximation of incompressible flow problems,555 Computer Methods in Applied Mechanics and Engineering 196 (21) (2007) 2413–2430. [40] O. Guasch, R. Codina, An algebraic subgrid scale finite element method for the convected helmholtz equation in two dimensions with applications in aeroacoustics, Computer Methods in Applied Mechanics and Engineering560 196 (45) (2007) 4672–4689. [41] W. R. Wolf, S. K. Lele, Trailing edge noise predictions using compressible les and acoustic analogy, in: Proceedings of the 17th AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. 2784, 2011, pp. 1–25. 32 (a) (b) Figure 8: Reference wave propagation contour for kd = 2.45 (a), reference wave propagation contour for kd = 4.91 (b), calculated total pressure (c) 33 View publication statsView publication stats