A high-order stabilized finite-element model for the linearized Navier-Stokes equations

Simone Olto 1 , Hadrien Bériot, Sophie Le Bras Siemens Industry Software NV Interleuvenlaan 68, 3001 Leuven, Belgium

Hervé Denayer KU Leuven and DMMS Lab Flanders Make Celestijnenlaan 300, 3001 Leuven, Belgium

ABSTRACT A numerical approach based on high-order stabilized finite elements is proposed to solve thermo- acoustic and visco-acoustic problems accounting for mean flow e ff ects. The approach is based on the linearized Navier-Stokes equations written in conservative form in the frequency domain. A high-order finite-element method using hierarchic shape functions is applied for enhanced accuracy. A novel enrichment strategy, inspired by the extended Finite Element Method (X-FEM), is developed to model the finer scales near the walls at a reasonable computational cost. It relies on a re-orthogonalization procedure proposed to preserve both the continuity of the solution and the conditioning properties of the discrete model. The performance of the method is evaluated by performing two-dimensional simulations of acoustic waves a ff ected by visco-thermal wall losses and mean flow e ff ects while propagating in a duct. Numerical results are in good agreement with analytical models from the literature.

1. INTRODUCTION

Acoustic dissipation due to visco-thermal losses are important in many applications like musical wind instruments or turbofan liners. The sound propagation in the presence of acoustic viscous and thermal boundary layers can be modelled by the linearized Navier-Stokes equations (LNSE). They also provide a description of the convection and refraction e ff ects of acoustic waves propagating through sheared mean flows. The LNSE have been widely solved using frequency-domain finite- element methods (FEM) [1, 2]. However, they often require very fine meshes to resolve the acoustic viscous and thermal boundary layers. These constraints on the mesh resolution can lead to computational costs that can be prohibitive for industrial applications. In this work, a numerical approach based on high-order stabilized finite elements is proposed to solve the LNSE. The equations are written in conservative form in the frequency domain. This high- order finite-element method, named p -FEM, is based on the use of hierarchic shape functions and leads to enhanced accuracy as well as reduced memory and CPU-time compared to conventional

1 corresponding author: simone.olto@siemens.com

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

low-order FEM. It has been already successfully applied to solve problems using various linearized acoustic operators such as the Helmholtz equation [3], the linearized potential equation [4] or the linearized Euler equations [5]. In order to overcome the limitation due to the mesh resolution, an enrichment strategy is proposed in this study. It consists in adding to the finite-element polynomial basis analytical functions which reproduce the visco-thermal and mean flow e ff ects close to the walls. In that way, LNSE simulations can be performed on coarser meshes that do not resolve the acoustic boundary layers, leading to a reduced computational cost. Similar enrichment techniques have been successfully applied in other fields in the literature. As an example, the extended finite element method (X-FEM) considers discontinuous basis functions to better approximate cracks in structural components [6]. An enrichment strategy has also been proposed by Krank et al. as a wall-modelling approach for Reynolds-averaged Navier-Stokes simulations [7] and detached eddy simulations [8]. To the authors’ knowledge, it is the first time an enrichment strategy is developed to model visco-thermal e ff ects on sound propagation. The performance of the enriched high-order FEM approach is evaluated by performing two- dimensional simulations of a closed-end waveguide containing a quiescent medium and of the sound attenuation in duct with uniform mean flow. Enriched LNSE solutions are compared to resolved LNSE solutions and analytical results. The paper is structured as follows. Section 2 describes the LNSE formulation together with the high-order FEM and the enrichment strategy. The closed-end waveguide problem is presented in Section 3.1 and the sound attenuation in a duct in Section 3.2. Conclusions are given in Section 4.

2. FINITE ELEMENT MODEL

2.1. Governing equations We consider a perfect gas and small perturbations represented by the density ρ ′ , the velocity vector u ′ = ( u ′ , v ′ ), the pressure p ′ and the temperature T ′ around a steady mean flow of density ρ 0 , velocity u 0 = ( u 0 , v 0 ), pressure p 0 and temperature T 0 . The behaviour of these perturbations is governed by the linearized Navier-Stokes equations (LNSE). These equations can be written in two- dimensional Cartesian coordinates as:

∂ t + ∂  A i + A µ i  q

E i , j ∂ M − 1 q

∂ q

∂ x i + ∂

! = 0 (1)

∂ x i

∂ x j

where q = [ ρ ′ , ( ρ u ) ′ , ( ρ v ) ′ , ( ρ E ) ′ ] T is the vector of the unknown variables. The energy perturbation writes as:

( ρ E ) ′ =  E 0 −|| u 0 || 2 2  ρ ′ + ρ 0 c V T ′ + u 0 ( ρ u ) ′ + v 0 ( ρ v ) ′ (2)

where c v is the specific heat at constant volume and energy E 0 = c V T 0 + || u 0 || 2 2 / 2. The matrices A i , A µ i , E i , j and M − 1 are defined from the mean flow properties in Appendix. These equations are solved in the frequency domain considering an implicit time dependence for the solution vector q ( x , t ) = q ( x , ω ) e + i ω t where ω is the angular frequency. The LNSE are formulated using a Galerkin finite-element approach. For a computational domain Ω with a boundary Γ , the variational formulation writes as: Z

ψ † i ω q − ∂ ψ †

h A i + A µ i i q − ∂ ψ †

∂ x i E i , j ∂ M − 1 q

! d Ω+ Z

 ψ † F q + ψ † F LNSE  d Γ = 0 (3)

∂ x i

∂ x j

where ψ is the test function vector, † is the Hermitian operator, F = A x n x + A y n y is the normal flux, n = [ n x , n y ] T is the unit normal vector on Γ pointing out of domain Ω and F LNSE is a flux term defined as:

F LNSE = A µ i n i q + E i , j n i ∂ M − 1 q

∂ x j (4)

The conventional Galerkin finite-element formulation is known to su ff er from stability issues for convection-dominated problems [9]. The LNSE formulation is thus stabilized using a streamlined upwind Petrov Galerkin stabilization technique [10].

2.2. High-order finite element method The weak form Equation 3 is solved using a high-order finite element method, named p -FEM. The p -FEM approach is less sensitive to dispersion errors and has shown to provide substantial reductions in memory and CPU time when compared to conventional low-order FEM on a variety of time-harmonic problems including Helmholtz [3], linearized potential theory [4] and linearized Euler equations [5]. The domain Ω is partitioned into a set of non-overlapping quadrangular elements. The grid is generated using the open source mesh generator Gmsh [11]. In each element, the numerical solution is made up of high-order H 1 -conforming hierarchical polynomial shape functions [12]. In two dimensions, these functions can be classified into three categories, namely, nodal, edge, and bubble functions. In the simple situation where the grid is aligned with the Cartesian coordinates, the number of degrees of freedom contributing to the solution on any element is given by

N dof = 4 + h 2( p x − 1) + 2( p y − 1) i + ( p x − 1)( p y − 1) (5)

where p x and p y correspond to the order enforced in the x and y directions respectively. Defining direction-dependent polynomial approximations allows to better control the accuracy on distorted elements, and / or when the physics is anisotropic [13]. The three terms in Equation 5 correspond to the numbers of nodal, edge, and bubble functions, respectively. The latter are internal to the element and may thus be condensed out at assembly level to further improve the computational e ffi ciency.

2.3. Enrichment strategy Solutions of the linearized Navier-Stokes equations are challenging, with the presence of very short length scales in the thermo-viscous boundary layers. While refining the mesh is always an option, it comes at a very high computational cost. Instead, one may leverage the a priori knowledge of the solution and resort to an enrichment strategy. In this paper an innovative enrichment strategy is proposed to enrich a high-order hierarchical polynomial basis for H 1 -conforming elements with one or more continuous enrichment functions. The enrichment procedure is based on the H 1 - orthogonalization of the enrichment functions over the polynomial basis. This procedure ensures the continuity of the approximation basis over the computational domain and improves the conditioning of the resulting system matrices [14]. In a one-dimensional domain, the enrichment function f can be seen as a summation of its projection over the polynomial space and an enriched basis function

v i a v i ψ v i + X

f = X

k i a k i e ψ e k i | {z } projection

+ ψ e (6)

where ψ v i are the linear nodal basis functions and ψ e k i are the high-order edge basis functions of order k i = 2 : p where p is the edge order. ψ e is the edge enriched basis which is defined as the component of the enrichment function that cannot be represented by the polynomial basis. a v i and a k i e are coe ffi cients obtained by a projection-based interpolation approach [12]. a v i are equal to f ( v i ), while a k i e verify the following minimisation problem:

f − X

H 1 e . (7)

v i f ( v i ) ψ v i − X

k i a k i e ψ e k i

min

This minimisation problem is local, computed edge by edge. It also ensures the global continuity of the approximation. Starting from the minimisation problem, the following weak formulation can be derived to compute a k i e :

 ∂ψ e i k i ∂ t e d e = 0 (8)

 f − X

∂ ∂ t e

Z

v i f ( v i ) ψ v i − X

k i a k i e ψ e k i

e

where t e is the edge tangent. It can be noted that only integrals over the edge e are present. Using Equation 8 the coe ffi cients a k i e are computed and from Equation 6 the edge enriched basis is evaluated. When more than one enrichment function is needed, the process restarts from Equation 7: the new enriched function is projected onto the polynomial basis and the already computed enriched edge basis functions. After computing the enriched edge functions for each edge, the approach can be extended to bubble shape functions applying a similar strategy. In the current paper the enrichment strategy is used to represent the viscous and thermal acoustic boundary layer profiles which originate along no-slip and isothermal walls, respectively. While the exact solution of the LNSE problem cannot be inferred a priori in wall regions on a complex use case, analytical models may be leveraged to approximate their behaviour. In this work, three enrichment functions are introduced. Analytical viscous and thermal boundary layer profiles f v and f t in quiescent mediums [15] are used to account for wall losses. In addition, in order to account for steep gradients in the vicinity of the walls, an hyperbolic tangent profile f hyp is also considered. For a two-dimensional duct of height H , these functions are defined as a function of the wall-normal distance y as:

f v = 1 − cosh 1 + i

! , f t = 1 − cosh 1 + i

δ v y !

δ t y !

! and f hyp = tanh y

! (9)

δ v

cosh 1 + i

cosh 1 + i

H

H

δ v

δ t

2

2

where δ v = p

2 µ 0 / ( ρ 0 ω ) and δ t = p 2 λ t / ( ρ 0 ω c p ) are the viscous and thermal boundary layer thicknesses, respectively. µ 0 is the dynamic viscosity, λ t is the thermal conductivity, c p is the specific heat at constant pressure, ρ 0 the density and ω the angular frequency. After projection and re-orthogonalization, these functions are simply appended to the original high-order polynomial approximation basis. Figure 1 shows a quadrilateral element with quadratic approximation basis p x = p y = 2 and the re-orthogonalization of the hyperbolic tangent enrichment function. The enriched finite element strategy should allow to more e ffi ciently address the multiple length scales present in the solution, while not having to drastically refine the mesh in the wall regions. It is worth to point out that the enrichment procedure proposed here is completely general, as it is independent from the operator and the enrichment functions at hand.

3. APPLICATIONS

The performance of the newly developed high-order finite element formulation for the linearized Navier-Stokes equations together with the enrichment strategy is evaluated by simulating two two- dimensional test cases, namely a closed-end waveguide containing a quiescent medium and a duct in the presence of a uniform mean flow (see Figure 2).

3.1. Closed-end waveguide A two-dimensional closed-end waveguide in a quiescent medium is first considered. The waveguide is presented in Figure 2a. Its length and height are equal to L = 170 mm and H = 2 . 5 mm, respectively. The waveguide open-end boundary is represented by a plane piston

Polynomial basis Enriched basis

Figure 1: Second order hierarchic polynomial basis functions ( p x = p y = 2), together with the re- orthogonalization of an hyperbolic tangent enrichment function.

(a)

(b)

Figure 2: Representation of (a) the closed-end waveguide geometry and (b) the duct geometry.

source, corresponding to an imposed inlet velocity of V in = 1 m / s. At the closed-end, a rigid adiabatic slip wall is considered. At the top and bottom walls of the waveguide, three loss configurations are successively studied: thermal losses only, viscous losses only and visco-thermal losses. They are associated to isothermal slip walls, adiabatic no-slip walls and isothermal no-slip walls, respectively. The frequency range of interest extends from 2900 Hz to 3150 Hz. It is selected around the resonance frequency of the third axial mode of the waveguide. The ambient medium is characterized by a speed of sound c 0 = 347 m / s, a temperature T 0 = 300 K, a density ρ 0 = 1 . 2 kg.m − 3 , a specific heat ratio γ = 1 . 4, a dynamic viscosity µ 0 = 1 . 829 · 10 − 5 kg.m − 1 s − 1 and a thermal conductivity λ t = 0 . 026 W.m − 1 K − 1 . Six di ff erent cases are considered for the FEM simulations. They are defined by the LNSE strategy, the resolution of the mesh and the polynomial orders p x and p y of the solution. The parameters of the simulations are given in Table 1. Simulations are performed with and without the enrichment strategy (referred to as eLNS) on three Cartesian meshes with a fine, medium and coarse spatial discretization. The fine mesh is made of N x = 85 and N y = 120 points in directions x and y . The grid spacing is uniform in direction x while a geometric progression of ratio 0.9 from the duct centerline to the duct walls is applied in direction y , leading to an acoustic viscous boundary layer thickness δ v discretized by at least 25 elements. The medium grid has the same resolution as the fine grid in direction y whereas N x = 20 in direction x . For the coarse grid, N x = 20 and N y = 2 are used. Simulations using the enrichment strategy are carried out using the enrichment functions f hyp , f v and / or f t . They are all performed using the coarse mesh. The FEM polynomial orders p x and p y are selected based

a y plane wave mia PML > PML | |Ha a < La v R

on the resolution of the FEM mesh. Using the fine mesh, p x and p y are equal to 2 as typically found in conventional low-order FEM simulations. For the Medium mesh, p y = 2 and p x = 10 as the grid spacing is large in direction x . Finally, for the coarse mesh, the simulations are performed using p x = 10 and p y = 2 and 8.

Table 1: Parameters of the simulations for the closed-end waveguide.

Case LNSE strategy Mesh N x N y p x p y Enrichment function(s)

LNS-Fine LNS Fine 85 120 2 2 -

LNS-Medium LNS Medium 20 120 10 2 -

LNS-Coarse LNS Coarse 20 2 10 2, 8 -

eLNS-hyp eLNS Coarse 20 2 10 2, 8 f hyp eLNS-v eLNS Coarse 20 2 10 2, 8 f v eLNS-vt eLNS Coarse 20 2 10 2, 8 f v and f t

The fluctuating pressure p ′ wall computed numerically at the center of the duct closed-end (see Figure 2a) is compared to the analytical solution p ′ ref given by Bossart [16]:

√

 1 + 1 − i √ k 0

l ′ h    , (10)

p ′ ref ( L z ) = ρ 0 c 0 V in k 0 k 00

1 i sin( k 00 L z ) with k 2 00 = k 2 0

2 H

 p

l ′ v + ( γ − 1) q

where l ′ v = µ 0 / ( c 0 ρ 0 ) and l ′ h = λ t / ( c 0 ρ 0 c p ) are the characteristic length for the viscous and thermal dissipation and k 0 = ω/ c 0 . A relative error ϵ with respect to the analytical solution is computed as:

ϵ = | p ′ wall − p ′ ref | | p ′ ref | (11)

Figure 3 shows the pressure p ′ wall and the relative error ϵ for vertical polynomial order p y = 2. For the three loss configurations, the LNS-Fine solution is close to the reference solution with ϵ < 4%, indicating the validity of the LNSE implementation. Without enrichment and using a fine mesh, as expected, visco-thermal losses are well predicted. The LNS-medium results do not significantly deviate from LNS-Fine, which demonstrates that the accuracy of LNS without enrichment is driven by the mesh resolution in the wall-normal direction. This is indeed confirmed by the high error levels obtained for LNS-Coarse, especially when viscous losses are considered. Using LNS enrichment, improved FEM results are obtained compared to the LNS-Coarse solution. For the three loss configurations, accurate results are obtained when the enrichment functions include the analytical boundary profiles f v and f t for viscous and thermal losses. For eLNS-hyp, higher error levels are obtained compared to eLNS-vt. This is attributed to the enrichment function f hyp which is not an analytical boundary layer solution of the problem. Figure 4 shows the pressure p ′ wall and the relative error ϵ for vertical polynomial order p y = 8. For LNS-Coarse performed without enrichment, increasing the polynomial order p y from 2 to 8 leads to lower error levels. However, these levels remain overall higher than those of the eLNS simulations using the viscous and thermal boundary layer enrichment functions. This demonstrates that for a coarse mesh, increasing the order of the solution without using enrichment is not su ffi cient to obtain accurate results. When enrichment is used, increasing the polynomial order p y from 2 to 8 has no e ff ect on the eLNS-vt solution. It leads to improved accuracy only when the enrichment function is not an analytical boundary layer profile for the type of losses accounted in the problem.

The profiles of the normalized velocity | u ′ | obtained at f = 3050 Hz and x = 4 L / 5 for the waveguide configuration including viscous losses are represented in Figure 5 as a function of the waveguide half-height, with the wall located at y = 0. The profile obtained from LNS-Fine is used as reference here. The LNS-Coarse solution shows large discrepancies with the reference solution. This is especially true for p y = 2 where the velocity gradient near the wall is not captured. Using enrichment, the solutions eLNS-v and eLNS-vt computed using the function f v are in very good agreement with the LNS-Fine solution for both p y = 8 and p y = 2. In eLNS-hyp, a poor agreement with LNS-Fine is reported for p y = 2. Using p y = 8, the discrepancy with LNS-Fine is reduced but remains higher than in eLNS-v and eLNS-vt.

180

180

100

100

170

170

10

10

160

160

1

2800 2900 3000 3100 3200 145 150

2800 2900 3000 3100 3200 145 150

2900 2950 3000 3050 3100 3150 0.1

1

2900 2950 3000 3050 3100 3150

180

100

170

10

160

2800 2900 3000 3100 3200 145 150

2900 2950 3000 3050 3100 3150 0.5 1

(a)

(b)

Figure 3: (a) Pressure p ′ wall at the closed-end of the waveguide and (b) relative error ϵ , for vertical polynomial order p y = 2.

180

180

100

50

170

170

10

10

160

160

1

2800 2900 3000 3100 3200 145 150

2800 2900 3000 3100 3200 145 150

2900 2950 3000 3050 3100 3150 0.1

1

2900 2950 3000 3050 3100 3150

180

50

170

10

160

2800 2900 3000 3100 3200 145 150

1

2900 2950 3000 3050 3100 3150

(a)

(b)

Figure 4: (a) Pressure p ′ wall at the closed end of the waveguide and (b) relative error ϵ , for vertical polynomial order p y = 8.

The normalized velocity profile obtained in eLNS-v in Figure 5 is decomposed into its polynomial and enriched components. Results are presented for p y = 2 and 8 in Figure 6, together with the

1.5

1.5

1

1

0.5

0.5

0 0.1 0.2 0.3 0.4 0.5 0

0 0.1 0.2 0.3 0.4 0.5 0

(a)

(b)

Figure 5: Normalized velocity | u ′ | in closed-end waveguide for viscous loss configuration at x = 4 L / 5 at f = 3050 Hz for polynomial order (a) p y = 2 and (b) p y = 8.

1.2

1.4

1

1

0.8

Solution eLNS-v Polynomial part Enriched part Enrichment function

0.6

0.6

0.4

0.2

0.2

0 0.1 0.2 0.3 0.4 0.5 0

0 0.1 0.2 0.3 0.4 0.5 0

(a)

(b)

Figure 6: Decomposition of the normalized velocity | u ′ | into its polynomial and enrichment components for the FEM simulation eLNS-v using enrichment function f v and polynomial orders (a) p y = 2 and (b) p y = 8. Results for viscous loss configuration at x = 4 L / 5 and f = 3050 Hz.

analytical enrichment function f v . At p y = 2, the enriched component of the solution plays a fundamental role in the representation of the high gradient at the wall responsible for the viscous losses. This also explains the di ff erences between the LNS-Coarse and eLNS-v solutions both computed using the Coarse mesh. At p y = 8, the enriched component also contributes to the accuracy of the solution, but to a lesser extent as the higher polynomial basis better captures the small length scales near the wall compared to the simulation at p y = 2. Thanks to enrichment, LNS simulations can be performed using coarser meshes leading to a reduction of the size of the FEM problem. As an illustration, Table 2 shows the number of degrees of freedom (DoFs), which determines the size of the global system, for the six LNSE cases and for various polynomial orders p y . Using enrichment, the largest number of Dofs is obtained for eLNS-vt as two enrichment functions are considered. However, the number of Dofs in eLNS-vt is 31 times smaller than LNS-Fine for p y = 2. The number of Dofs being related to the memory requirements, the enrichment strategy allows to reduce drastically the memory needed to solve the LNSE.

3.2. Sound attenuation in a duct with uniform mean flow The enrichment strategy for LNSE is now applied to predict the viscous losses in a straight duct in the presence of a uniform mean flow. The duct height and length are H d = 2 . 5 mm and L d = 1 m, respectively (see Figure 2b). The uniform mean flow is characterized by a Mach number M 0 =

Table 2: Number of degrees of freedom for the numerical simulations of the closed-end waveguide.

Case polynomial order in the y direction

p y = 2 p y = 8

LNS-Fine 164 800 -

LNS-Medium 193 800 -

LNS-Coarse 4 000 13 700

eLNS-hyp & eLNS-v 4 600 14 200

eLNS-vt 5 200 14 800

u 0 / c 0 = 0 . 05, a speed of sound c 0 = 347 m / s, a temperature T 0 = 300 K, a density ρ 0 = 1 . 2 kg / m 3 , a specific heat ratio γ = 1 . 4 and a dynamic viscosity µ 0 = 1 . 829 · 10 − 5 kg.m − 1 s − 1 . Under these conditions, the flow is considered as laminar (with Reynolds number Re = ρ 0 u 0 H d /µ 0 < 2900) and the acoustic dissipation associated to the flow turbulence is assumed to be negligible. The horizontal duct walls are modelled by an adiabatic no-slip boundary condition, to account for viscous losses. Thermal losses are not considered here. At the duct ends, Perfectly Matched Layers (PML) [3] are used to absorb outgoing waves. Only downstream propagating waves are considered in this study. A plane wave is therefore injected at the duct inlet located at the left-hand side of the domain. The frequency of interest ranges f rom 2 kHz to 5 kHz, corresponding to shear numbers S h = H d /δ varying from 70 to 120, with δ = p

µ 0 / ( ρ 0 ω ). For S h > 10, the duct can be considered wide. Under this assumption the attenuation coe ffi cient, which is the imaginary part of the wavenumber, of the fundamental mode has been computed analytically by Weng et al. [17]. Originally derived for circular ducts [17], the analytical solution for the downstream propagating mode can be written for two-dimensional ducts in the presence of only viscous losses as follows:

α = δ √

2 H d (1 + M 0 ) 3 / 2 . (12)

In the numerical simulations, the attenuation coe ffi cient is computed from the pressure signal on the duct centerline using the matrix pencil method [18,19]. In order to quantify the losses in the duct, the normalized attenuation α/α 0 is introduced where α 0 is the attenuation coe ffi cient of the fundamental mode in a quiescent medium. Four di ff erent cases are considered for the FEM simulations. As in Section 3.1, they are defined by the LNSE strategy, the resolution of the mesh and the polynomial orders p x and p y (see Table 3).

Table 3: Parameters of the simulations for the duct with uniform mean flow.

Case LNSE strategy Mesh N x N y p x p y Enrichment function(s)

LNS-Fine LNS Fine 510 120 2 2 -

LNS-Coarse LNS Coarse 120 2 10 2, 8 -

eLNS-hyp eLNS Coarse 120 2 10 2, 8 f hyp eLNS-v eLNS Coarse 120 2 10 2, 8 f v

The normalized attenuation coe ffi cient α/α 0 obtained numerically is compared to the analytical solution in Figure 7. For p y = 2, the solutions obtained with LNS-Fine and eLNS-v overlap and

show a deviation of about 4 . 2% with the analytical solution. Thanks to enrichment, viscous losses are captured using a coarse mesh. More pronounced discrepancies with the analytical solution are obtained with eLNS-hyp, with deviation of about 20%. The poorest results are obtained for LNS- Coarse, i.e. without enrichment. Increasing the polynomial order p y from 2 to 8, no significant change is observed in eLNS-v. On the contrary, a better agreement with analytical solution is obtained in LNS-hyp when p y = 8, as also reported for the waveguide problem in Section 3.1.

1

1

0.8

0.9

0.6

0.8

0.4

0.7

2000 2500 3000 3500 4000 4500 5000 0.2

2000 2500 3000 3500 4000 4500 5000 0.6

(a)

(b)

Figure 7: Normalized attenuation of the fundamental mode obtained from the pressure perturbation on the waveguide centerline computed using a vertical polynomial order equal to (a) p y = 2 and (b) p y = 8.

The decomposition of the eLNS-v velocity profile | u ′ | into its polynomial and enriched components has also been studied. The results obtained for p y = 2 and 8 are very similar to those presented in Figure 6. For the sake of conciseness, these results are not shown here.

4. CONCLUSIONS

A new high-order stabilized finite element formulation has been proposed to solve the linearized Navier-Stokes equations in the frequency domain. This formulation can accurately predict the propagation of acoustic waves in ducts accounting for visco-thermal losses and mean flow e ff ects provided that the mesh is su ffi ciently fine to capture the acoustic boundary layers. For industrial problems, these simulations, however, are known to require high computational memory. To tackle this shortcoming, an innovative enrichment strategy has been proposed in this study. It consists in enriching the finite-element polynomial basis with analytical functions that reproduce the visco-thermal and mean flow e ff ects close to the walls. The enrichment strategy is based on a H 1 -orthogonalization of each enrichment function with respect to the high-order polynomial basis and to the other enrichment functions. The performance of the enriched strategy is examined for a closed-end waveguide containing a quiescent medium and for a duct with uniform mean flow. The enriched solutions are in good agreement with the analytical reference results. The use of the enrichment strategy leads to a significant reduction of the amount of memory required for the resolution compared to linearized Navier-Stokes simulations performed without enrichment. The analysis reveals the benefits of the enrichment strategy when low finite-element polynomial orders and coarse meshes are used. It is also worth to mention that the construction of the enrichment basis is independent from the operator (the linearized Navier-Stokes operator in this study). Therefore, the enrichment procedure could be easily extended to other types of multiscale problems.

ACKNOWLEDGEMENTS

This work is part of the Marie Skłodowska-Curie Innovative Training Network Pollution Know-How and Abatement (POLKA). We gratefully acknowledge the financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 813367.

APPENDIX: MATRICES OF THE 2D LINEARIZED NAVIER-STOKES EQUATIONS

The flux matrices A x and A y are defined as:









0 0 1 0

0 1 0 0

− u 0 2 − H 0 u 0 (3 − γ ) − γ m v 0 γ m − u 0 v 0 v 0 u 0 0

− u 0 v 0 v 0 u 0 0

A y =

A x =

− v 0 2 − H 0 − γ m u 0 v 0 (3 − γ ) γ m − ( G 0 + H 0 ) v 0 − u 0 v 0 γ m G 0 − v 0 2 γ m γ v 0

− ( G 0 + H 0 ) u 0 G 0 − u 0 2 γ m − u 0 v 0 γ m γ u 0

(13) where R is the gas constant, H 0 = ( γ − 1) ( E 0 − c V T 0 −|| u 0 || 2 2 ), G 0 = E 0 + RT 0 and γ m = γ − 1. The matrices M − 1 , A µ i and E i , j write as:





1 0 0 0





0 0 0 0

− u 0

1 ρ 0 0 0

ρ 0

0 0 0 0

M − 1 =

A µ i =

− v 0

ρ 0 0 1 ρ 0 0

0 0 0 0 u 0 τ 0 ( i , 1)

ρ 0 + v 0 τ 0 ( i , 2)

ρ 0 − τ 0 ( i , 1)

ρ 0 − τ 0 ( i , 2)

04 − E 0 −|| u 0 || 2 2 ρ 0 c V − u 0

ρ 0 0

ρ 0 c V − v 0

1 ρ 0 c V

ρ 0 c V

(14)





0 0 0 0

3 µδ i 1 δ j 1 − µδ i 2 δ j 2 2 3 µδ i 1 δ j 2 − µδ i 2 δ j 1 0

0 − 4

E i , j =

(15)

0 − µδ i 1 δ j 2 + 2

3 µδ i 2 δ j 1 − µδ i 1 δ j 1 − 4

3 µδ i 2 δ j 2 0

0 E i , j (4 , 2) E i , j (4 , 3) − λ t δ i j

! + v 0

− δ i 1 δ j 2 + 2

!# (16)

E i , j (4 , 2) = µ " u 0

− 4

3 δ i 1 δ j 1 − δ i 2 δ j 2

3 δ i 2 δ j 1

! + u 0

− δ i 2 δ j 1 + 2

!# (17)

E i , j (4 , 3) = µ " v 0

− 4

3 δ i 2 δ j 2 − δ i 1 δ j 1

3 δ i 1 δ j 2

where τ 0 = ∇ u 0 + ∇ u T 0 − 2

3 ∇· u 0 I (18)

∇ u 0 and ∇· u 0 are the gradient and the divergence of the velocity, I is the identity matrix and · T the transpose operator. δ ij is the Kronecker delta.

REFERENCES

[1] J. Gikadi. Prediction of acoustic modes in combustors using linearized Navier–Stokes equations in frequency space . PhD thesis, Technische Universität München, München, 2014.

4% t 4 x * t+ y

[2] M. Meindl, A. Albayrak, and W. Polifke. A state-space formulation of a discontinuous galerkin method for thermoacoustic stability analysis. Journal of Sound and Vibration , 481:115431, 2020. [3] H. Bériot and A. Modave. An automatic perfectly matched layer for acoustic finite element simulations in convex domains of general shape. International Journal for Numerical Methods in Engineering , 122(5):1239–1261, 2021. [4] G. Gabard, H. Bériot, A.G. Prinn, and K. Kucukcoskun. Adaptive, high-order finite-element method for convected acoustics. AIAA Journal , 56(8):3179–3191, 2018. [5] K. Hamiche, S. Le Bras, G. Gabard, and H. Beriot. Hybrid numerical model for acoustic propagation through sheared flows. Journal of Sound and Vibration , 463:114951, 2019. [6] K. Agathos, E. Chatzi, and S.P.A. Bordas. Stable 3D extended finite elements with higher order enrichment for accurate non planar fracture. Computer Methods in Applied Mechanics and Engineering , 306:19–46, 2016. [7] B. Krank and W.A. Wall. A new approach to wall modeling in LES of incompressible flow via function enrichment. Journal of Computational Physics , 316:94–116, 2016. [8] B. Krank, M. Kronbichler, and W.A. Wall. Wall modeling via function enrichment: Extension to detached-eddy simulation. Computers and Fluids , 179:718–725, 2019. [9] J. Donea and A. Huerta. Finite element methods for flow problems . John Wiley & Sons, 2003. [10] P. Rao and P. Morris. Use of finite element methods in frequency domain aeroacoustics. AIAA Journal , 44(7):1643–1652, 2006. [11] C. Geuzaine and J.F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering , 79(11):1309–1331, 2009. [12] P. Šolín, K. Segeth, and I. Doležel. Higher-order finite element methods . Chapman & Hall / CRC, 2004. [13] H. Bériot and G. Gabard. Anisotropic adaptivity of the p-FEM for time-harmonic acoustic wave propagation. Journal of Computational Physics , 378:234–256, 2019. [14] K. Agathos, E. Chatzi, and S.P.A. Bordas. A unified enrichment approach addressing blending and conditioning issues in enriched finite elements. Computer Methods in Applied Mechanics and Engineering , 349:673–700, 2019. [15] S. Moreau, H. Bailliet, J.C. Valière, R. Boucheron, and G. Poignand. Development of laser techniques for acoustic boundary layer measurements. Part II: Comparison of LDV and PIV measurements to analytical calculation. Acta Acustica united with Acustica , 95(5):805–813, 2009. [16] R. Bossart, N. Joly, and M. Bruneau. Hybrid numerical and analytical solutions for acoustic boundary problems in thermo-viscous fluids. Journal of Sound and Vibration , 263(1):69–84, 2003. [17] C. Weng and F. Bake. An analytical model for boundary layer attenuation of acoustic modes in rigid circular ducts with uniform flow. Acta Acustica united with Acustica , 102(6):1138–1141, 2016. [18] C. Weng, A. Schulz, D. Ronneberger, L. Enghardt, and F. Bake. Flow and viscous e ff ects on impedance eduction. AIAA Journal , 56(3):1118–1132, 2018. [19] T. Humbert, R. Delalande, G. Gabard, J. Golliard, and Y. Auregan. Performance of the matrix pencil algorithm in direct impedance eduction of liners: some numerical experiments. 25th AIAA / CEAS Aeroacoustics Conference , 2019.