<figure><img src="image_1.png" /></figure>

HAL
open science

# Discontinuous Galerkin Implementation of the Extended Helmholtz Resonator Model in Time Domain

Nicolas Chevaugeon, Jean-François Remacle, Xavier Gallez

► To cite this version:

Nicolas Chevaugeon, Jean-François Remacle, Xavier Gallez. Discontinuous Galerkin Implementation of the Extended Helmholtz Resonator Model in Time Domain. 12th AIAA/CEAS Aeroacoustics Conference (27th AIAA Aeroacoustics Conference), May 2006, Cambridge, United States. 10.2514/6.2006-2569. hal-03920882

HAL Id: hal-03920882

https://hal.science/hal-03920882

Submitted on 3 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Discontinuous Galerkin Implementation Of The Extended Helmholtz Resonator Model In Time Domain

Nicolas Chevaugeon, Jean-François Remacle, Université Catholique de Louvain – www.gce.ucl.ac.be
1, place du levant, batiment Vinci, B1348 Louvain-la-Neuve, Belgium

Xavier Gallez, Free Field Technologies SA – www.fft.be
Axis Park Louvain-la-Neuve 1 rue Emile Francqui B-1435 Mont-Saint-Guibert Belgium

## Summary

For being the major source of noise pollution in aero-engines, fan forward and rearward noise is of primary importance in the aeronautics industry. The computation of aircraft engine noise is usually treated in frequency domain. Relevant quantities like acoustic power spectrum or sound directivity can be easily post processed out of such simulation.

Considering that a constant number of grid points per wavelength is required, the size of the linear systems is growing like $f^3$, if $f$ is the frequency. Because the systems to solve are stiff, especially in the presence of a mean flow, frequency domain solvers cannot handle very large problems. Typically, if $D$ is the diameter of the engine nacelle and $\lambda$ is the wavelength of the signal, frequency domain solvers are limited to a ratio $k_r = D/\lambda$ of around 30. The limited efficiency and scalability of direct solvers do not make them good candidates for addressing problems where the Helmholtz number $k_r$ is as high as 50, especially on distributed memory computers.

A possible way to increase the range of frequencies of the calculations is to switch to time domain. It is still possible to compute both power spectrum and sound directivity but not as directly as it was done in the frequency domain. So, a time domain approach has to be proven to be at least as efficient as the frequency domain one.

In the framework of the Messiaen project (European collaborative project under the Sixth Framework Programme), an efficient high order discontinuous Galerkin method (DGM) has been developed that solves the linearized Euler equations. The DG method$^{1-4}$ is a popular scheme for the resolution of hyperbolic conservation laws. The properties of the quadrature free implementation of the method$^5$ allows to obtain an efficiency that is close to the peak efficiency of the processor. Moreover, good scalability properties are obtained in parallel. This method is now implemented in an industrial framework (Actran DGM). It has been demonstrated that such a method is a good alternative to the frequency domain at high Helmholtz number.$^6$

The aim of our work is to participate to a more general research that aim is to find ways to reduce the noise of aircraft engines. One technology that enables noise reduction the use of acoustic liners. Liners notably reduce the noise power spectrum and allow to change the sound directivity.

Liners are usually modeled in the frequency domain as a frequency dependent wall impedance. This model has to be translated in the time domain where all the frequencies are present, even when we try to solve for a fixed frequency input data. Therefore, an impedance model that matches the design impedance and that is valid for a range of frequency is needed. In$^7$ S.W. Rienstra has developed a modified Helmholtz resonator model in time-domain. Starting from the quadrature free RK-DGM framework, we provide a description of the data structures and algorithms that are required to implement this time domain impedance model. We will discuss the conversion from the direct expression of the relation between acoustic pressure and normal velocity to boundary conditions implemented in the discontinuous Galerkin code for aeroacoustic problems. First, we will detail the non-flow case. Then, we will extend the model when a mean flow is present. The
latest problem is, from far, the most difficult one. We will show what has already been conjectured, i.e. that the Ingard-Myers condition⁸,⁹ that is usually used in the flow case is unstable in the time domain. In this work, we propose a way to circumvent this instability. Here, the model is itself unstable and the usual numerical tricks that prevent numerical instabilities do not apply. Some examples are provided that show the new approach to be stable and accurate. Comparisons with a frequency domain code are provided.

# I. introduction

In the frequency domain, soft walls (liners) are modeled by a complex impedance $Z(\omega)$ that depends on the pulsation $\omega$. This impedance defines a relation between the pressure $\hat{p}$ and normal velocity $\hat{v}$:

$$ \hat{p}(\omega) = Z(\omega)\hat{v}(\omega). \quad (1) $$

Equation (1) is valid when there is no mean-flow. In the presence of a mean-flow, the relation becomes:⁸,⁹

$$ (\hat{\mathbf{v}} \cdot \mathbf{n}) = (i\omega + \mathbf{v}_0 \cdot \nabla - \mathbf{n} \cdot (\mathbf{n} \cdot \nabla \mathbf{v}_0)) \frac{\hat{p}}{i\omega Z} \quad (2) $$

where $\mathbf{n}$ is the outward normal and $\mathbf{v}_0$ is the mean-flow velocity. Since we are interested in the relation on wall, the no-penetration of the mean-flow condition($\mathbf{v}_0 \cdot \mathbf{n} = 0$), allows to recast (2) to

$$ (\hat{\mathbf{v}} \cdot \mathbf{n}) = (i\omega + \mathbf{v}_0 \cdot \nabla + \mathbf{v}_0 \cdot (\mathbf{n} \cdot \nabla \mathbf{n})) \frac{\hat{p}}{i\omega Z} \quad (3) $$

From equation (1) in frequency domain, one can cast the equation in time domain by mean of an inverse Fourier transform to get, for the no mean-flow case:

$$ p(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} z(t-\tau)v(\tau)d\tau \quad (4) $$

and

$$ \frac{\partial}{\partial t} p + \mathbf{v}_0 \cdot \nabla p + \mathbf{v}_0 \cdot (\mathbf{n} \cdot \nabla \mathbf{n}) p = \frac{1}{2\pi} \int_{-\infty}^{+\infty} z(t-\tau) \frac{\partial}{\partial t} v(\tau) d\tau \quad (5) $$

in presence of a mean-flow, where $p$ is the time domain pressure, $v$ is the time domain normal velocity and $z(\tau)$ denotes the inverse Fourier transform of $Z(\omega)$.

When one wants to solve acoustic problem in time domain, one needs a model for $Z(\omega)$, since detailed knowledge of $Z(\omega)$ is not available. What is usually given is the impedance for one or more design frequency $\omega_0$ such as $Z_0 = Z(\omega_0)$.

The challenge is to develop a model for which $Z(\omega)$ has an explicitly known inverse Fourier transform, and which is physically possible. In S.W. Rienstra paper,⁷ conditions on the expression of the impedance function are derived and applied on different example. Mass-spring-damper, rational function of $\omega$, basic Helmholtz resonator and finally the so-called Extended Helmholtz Resonator were tested.

The impedance function for the Extended Helmholtz Resonator model can be written has follow:

$$ Z(\omega) = R + i\omega m - i\beta \cot\left(\frac{1}{2}\omega\nu\Delta t - i\frac{1}{2}\varepsilon\right) \quad (6) $$

Where $\varepsilon$ is the damping in the cavity's fluid, $\beta$ is the cavity reactance, $\nu\Delta t = 2L/c_0$ represent the depth of the cavity, $R$ is the face-sheet resistance and $\omega m$ is the face-sheet mass reactance. This model is proved to insure reality condition, is passive and causal, if we have:

$$ R \ge 0 \quad (7) $$

$$ m \ge 0 \quad (8) $$

$$ \nu\Delta t \ge 0 \quad (9) $$

$$ \beta > 0 \quad (10) $$

$$ \epsilon > 0 \quad (11) $$
The inverse Fourier transform $z(t)$ of $Z(\omega)$ can be written explicitly as:

$$ \frac{z(t)}{2\pi} = (R + \beta)\delta(t) + m\delta'(t) + 2\beta \sum_{n=1}^{\infty} e^{-\epsilon n} \delta(t - n\nu\Delta t) \quad (12) $$

which, in the no-flow case, gives an equation for the boundary condition in terms of $p(t)$, $v(t)$, $v'(t)$ and values of $v(t)$ in the past:

$$ p(t) = (R + \beta)v(t) + mv'(t) + 2\beta \sum_{n=1}^{\infty} e^{-\epsilon n}v(t - n\nu\Delta t) \quad (13) $$

An alternative derivation give the following expression :

$$ p(t) - e^{-\epsilon}p(t - \nu\Delta) = (R + \beta)v(t) - (R - \beta)e^{-\epsilon}v(t - \nu\Delta) + m(v'(t) - v'(t - \Delta t)) \quad (14) $$

Other forms of the time domain pressure, normal velocity relation are available. In the first part of this study, we will consider the implementation of relation (14). If the impedance is given for one design frequency only, the model is over parameterized and some parameters can be chosen freely. For any of these choices, the response in frequency domain to a forced harmonic input will be the same, while the transitory response will be different. In particular, it's always possible to choose $m = 0$, and avoid the computation of the inertia term in the equation. This is an advantage in term of computation time of course. Nevertheless, we will consider this term for several reason. First, the model with mass fully represents a Helmholtz resonator in time domain. Analytical solution of the response of such Helmholtz resonator in time domain are available.¹⁰ Then, the treatment of the mass term will be mandatory in the flow case anyway.

## II. Liner boundary condition : the no-flow case

In order to discuss the implementation of the boundary condition, we need first to briefly remind the construction of the discontinuous Galerkin discretization. We start with a conservation law of the form:

$$ \partial_t \mathbf{u} + \partial_x \mathbf{F} + \partial_y \mathbf{G} + \partial_z \mathbf{H} = \mathbf{r} \quad (15) $$

where $\mathbf{u}$ is a set of Eulerian variables, $\mathbf{F}$, $\mathbf{G}$ and $\mathbf{H}$ are convective fluxes and $\mathbf{r}$ is a source term. The following variational form

$$ \int_{\Omega} \partial_t \mathbf{u} \hat{\mathbf{u}} \, dv = \int_{\Omega} \mathbf{F} \partial_x \hat{\mathbf{u}} \, dv + \int_{\Omega} \mathbf{G} \partial_y \hat{\mathbf{u}} \, dv + \int_{\Omega} \mathbf{H} \partial_z \hat{\mathbf{u}} \, dv \\ - \int_{\partial \Omega} \mathbf{f} \hat{\mathbf{u}} \, ds + \int_{\Omega} \mathbf{r} \hat{\mathbf{u}} \, dv, \quad \forall \hat{\mathbf{u}} \quad (16) $$

is obtained by multiplying by an arbitrary test function $\hat{\mathbf{u}}$, and integrating by part over the domain $\Omega$. In (16), $\partial\Omega$ is the boundary of $\Omega$ and $\mathbf{f}$ is the flux in the direction $\mathbf{n}$ normal to $\Omega$. The domain is then divided in a set of elements $e$, where the field variables are approximated by $\mathbf{u} = \mathbf{u}_k^e \phi_k$, where $\phi_k$ are the polynomial shape functions and $\mathbf{u}_k^e$ are the coefficients of the shape function attached to the element $e$. The test function $\hat{\mathbf{u}}$ are approximated on the same space than $\mathbf{u}$, so that we get symmetric mass matrices. Note that no continuity conditions are imposed between elements. After the assumption that the fluxes are linear, and that the Jacobian of the element is constant, a linear system of equations is obtained in each element that enables to compute the derivative of the discretized field variables:

$$ \begin{aligned} \partial_t \mathbf{u}_k^e \int_e \phi_k \phi_j \, dv = & (\mathbf{F}_k^e) \int_e \phi_k \partial_x \phi_j \, dv \\ & + (\mathbf{G}_k^e) \int_e \phi_k \partial_y \phi_j \, dv \\ & + (\mathbf{H}_k^e) \int_e \phi_k \partial_z \phi_j \, dv \\ & + (\mathbf{r}_k^e) \int_e \phi_k \phi_j \, dv \\ & - \sum_{l=1}^{n_e} \bar{\mathbf{f}}(\mathbf{u}_L^{e_l}, \mathbf{u}_R^{e_l})_k \int_{e_l} \phi_k \phi_j \, ds = 0. \end{aligned} \quad (17) $$
In (17), $e_l$ is the $l^{th}$ boundary of element $e$ (e.g the $l^{th}$ edge of a triangular element for two dimensional domain or the $l^{th}$ triangular face of a tetrahedron for three dimensional problems). $\tilde{\mathbf{f}}(\mathbf{u}_{Lk}^{e_l}, \mathbf{u}_{Rk}^{e_l})$ is the numerical normal flux between the values of the field at left and right of boundary $e_l$, in the direction normal to $e_l$, pointing outside $e_l$. This flux is usually computed by solving a Riemann problem between the two states, to get the right up winding that ensures the convergence of the scheme. In the case of linearised Euler equation, the Riemann problem can be solved exactly and the solution is linear with regard to the left and the right states $\mathbf{u}_L$ and $\mathbf{u}_R$. Boundary conditions are usually inserted in the scheme by imposing the normal flux on the boundary of the domain.

In the case of aeroacoustic, we aim to solve the linearized Euler equation with mean flow in primitive variable :

$$p' = \nabla \cdot \mathbf{f}_p \qquad (18)$$

$$u' = \nabla \cdot \mathbf{f}_u \qquad (19)$$

$$v' = \nabla \cdot \mathbf{f}_v \qquad (20)$$

$$w' = \nabla \cdot \mathbf{f}_w \qquad (21)$$

where $p$ is the acoustic pressure, $u, v, w$ are the three component of the velocity. $'$ is the time derivative operator, $f_p, f_u, f_v, f_w$ are the fluxes and are linear functions of $p, u, v$ and $w$. We have

$$\mathbf{f}_p = \left[ \begin{array}{llll} \rho_0 c_0^2 u + p u_0 & \rho_0 c_0^2 v + p v_0 & \rho_0 c_0^2 w + p w_0 & \end{array} \right]^T \qquad (22)$$

$$\mathbf{f}_u = \left[ \begin{array}{llll} \frac{1}{\rho_0} p + u u_0 & u v_0 & u w_0 & \end{array} \right]^T \qquad (23)$$

$$\mathbf{f}_v = \left[ \begin{array}{llll} v u_0 & \frac{1}{\rho_0} p + v v_0 & v w_0 & \end{array} \right]^T \qquad (24)$$

$$\mathbf{f}_w = \left[ \begin{array}{llll} w u_0 & w v_0 & \frac{1}{\rho_0} p + w w_0 & \end{array} \right]^T \qquad (25)$$

where $u_0, v_0, w_0$ are the component of the mean flow velocity, $c_0$, $\rho_0$ are the meanflow sound speed and density. If the mean flow is not constant, an additional source term need to be computed, that depends on spatial derivatives of the mean flow quantity.

Following the DG scheme, we have, for each element $e$, the followings uncoupled systems of discretized ODE to solve.

$$\mathbf{M} \mathbf{p}'^e = \mathbf{K}_x^e \mathbf{F}_p^e + \mathbf{K}_y^e \mathbf{G}_p^e + \mathbf{K}_z^e \mathbf{H}_p + \sum_{l=0}^{nb} \mathbf{K}_l \mathbf{f}_p \frac{s_l}{v_e} = \mathbf{R}_p \qquad (26)$$

$$\mathbf{M} \mathbf{u}'^e = \mathbf{K}_x^e \mathbf{F}_u^e + \mathbf{K}_y^e \mathbf{G}_u^e + \mathbf{K}_z^e \mathbf{H}_u + \sum_{l=0}^{nb} \mathbf{K}_l \mathbf{f}_u \frac{s_l}{v_e} = \mathbf{R}_u \qquad (27)$$

$$\mathbf{M} \mathbf{v}'^e = \mathbf{K}_x^e \mathbf{F}_v^e + \mathbf{K}_y^e \mathbf{G}_v^e + \mathbf{K}_z^e \mathbf{H}_v + \sum_{l=0}^{nb} \mathbf{K}_l \mathbf{f}_v \frac{s_l}{v_e} = \mathbf{R}_v \qquad (28)$$

$$\mathbf{M} \mathbf{w}'^e = \mathbf{K}_x^e \mathbf{F}_w^e + \mathbf{K}_y^e \mathbf{G}_w^e + \mathbf{K}_z^e \mathbf{H}_w + \sum_{l=0}^{nb} \mathbf{K}_l \mathbf{f}_w \frac{s_l}{v_e} = \mathbf{R}_w \qquad (29)$$

In (26), $v_e$ is the volume of element $e$ and $s_l$ is the area of surface element $l$ on the boundary of $e$. $\mathbf{K}_x^e$, $\mathbf{K}_y^e$, $\mathbf{K}_z^e$ are the elementary matrices corresponding to the components of the gradient operator. $\mathbf{R}_p, \mathbf{R}_u, \mathbf{R}_p, \mathbf{R}_u$ are shortcuts for the residual of each equation. We have

$$\mathbf{K}_x^e = \frac{1}{v_e} \int_e \phi_k \partial_x \phi_j \, dv \qquad (30)$$

$$= \int_{ref} \phi_k (\partial_\xi \phi_j \partial_x \xi + \partial_\eta \phi_j \partial_x \eta + \partial_\zeta \phi_j \partial_x \zeta) d\xi d\eta d\zeta \qquad (31)$$

$$= \mathbf{K}_\xi \partial_x \xi + \mathbf{K}_\eta \partial_x \eta + \mathbf{K}_\zeta \partial_x \zeta \qquad (32)$$
where $\xi$, $\eta$, $\zeta$ are the coordinate in the reference element. $\partial_x\xi$ and term of the like are the component of the Jacobian of the mapping from the reference element to element $e$. Note that the Jacobian is supposed to be constant. This is the case for simplexes (tetrahedron, triangles). This simplification allows to write the quadrature free version of the DG scheme. With a simple rearrangement of 26, the volume term can be computed all at once for all the elements by three matrix-matrix computation since the matrices $K_\xi$, $K_\eta$ and $K_\zeta$ are the same for each element.

For each face $l$ of element $e$, the interface term can be computed as :

$$ \int_l \phi_i \phi_j \mathbf{f}_j dh \qquad (33) $$

If nodal shape function are used, the $\phi_i$ are null for each $i$ corresponding to a node other than the one on the face considered. Therefore, only the flux on the node on the face have to be computed. Provided that we know how to link boundary degree of freedom to their numbering in the full element space, we can compute the matrix

$$ \int_l \psi_k, \psi_l \mathbf{f}_j dh \qquad (34) $$

where $k$ and $l$ are indices on the space on the edge. (e.g. for order 3 triangle, we have 10 degrees of freedoms per element and each edge uses 4 of them). $\psi_k$ are now the face shape functions. In the case of triangular meshes, there exists 6 possibilities to pass from edge degree of freedom numbering to element degree of freedom numbering, depending on which edge of the element is used and in which direction it is used.

Now comes the impedance boundary conditions. At time step $t$, we take the normal velocity of the fluid on the liner boundary $p_l$ and update pressure using relation (13). The normal flux of this variables are the computed and included in the part of the residual corresponding to the right equation. The normal flux $\mathbf{f}$ are computed with the new pressure as if the field was continuous on the boundary using the new pressure.

$$ v_n = n_0 u + n_1 v + n_2 w \qquad (35) $$

$$ f_p = (un_0 + vn_1 + wn_2) \qquad (36) $$

$$ f_u = n_0 P_l \qquad (37) $$

$$ f_v = n_1 P_l \qquad (38) $$

$$ f_w = n_2 P_l \qquad (39) $$

Note that we used here the assumption that the mean flow was not penetrating the liner($v_0 \cdot n = 0$). with more detail, at each boundary node we need to compute :

$$ f_u = n_0((R + \beta)v + mv' + 2\beta \sum_{n=1}^{\infty} e^{-\epsilon n}v(t - n\nu\Delta t)) \qquad (40) $$

The term in $v$ need no particular treatment. The term in $v'$, the inertia term will need special treatment since the time derivatives are not known explicitly. The last term, the memory term need some extra storage that will be detailed in the next subsection.

## A. Memory term

To compute $p_l(t)$, we need the following memory term $I(t)$, recast from 40 in a more numerically usable way:

$$ I(t) = \sum_{n=1}^{\infty} e^{-\epsilon n} v(t - n\nu\Delta t) \qquad (41) $$

By making the assumption that for $t < 0$, $v(t) = 0$, we can born the sum to a value $q$ such as $t - q\nu\Delta t > 0$ and $t - (q + 1)\nu\Delta t < 0$

The sum in equation 40 should be expressed differently in order to show that only $\nu$ values need to be stored.

The time step of the explicit time marching procedure being $\Delta t$, lets assume that the $k$ time step was computed, and we want to compute the solution for $k + 1^{th}$ time step, the solution procedure being started
at $t_0 = 0$ Lets define $v_k = v(k\Delta t)$ and $p_k = p(k\Delta t)$ and let be $r$ the rest of the division of $k$ by $\nu$. With this definition of $k$, $q$ is the integer quotient of $k/\nu$ such as $k = q\nu + r$

at time step $k$, we need $I_k = I(k\Delta t)$ to compute the flux.

$$I_k = I_{q\nu+r} \qquad (42)$$

$$= \sum_{n=1}^{\infty} e^{-\epsilon n} v(((q-n)\nu + r)\Delta t) \qquad (43)$$

$$= \sum_{n=1}^{q} e^{-\epsilon n} v_{r+(q-n)\nu} \qquad (44)$$

since for $t < 0$, $v(t)$ is supposed to be null, and then $v_{r+(q-n)\nu} = 0 \quad \forall n > q$. The expression of $I_k$ can then be rewrite as an expression of $I_{k-nu}$:

$$I_{q\nu+r} = e^{-\epsilon} (I_{(q-1)\nu+r} + v_{r+(q-1)\nu}) \qquad (45)$$

or

$$I_k = e^{-\epsilon} (I_{k-\nu} + v_{k-\nu}) \qquad (46)$$

From previous equation, it is clear that we only need to store $\nu$ values for each point where we need to compute the field : at time step $k$, we have to have access to $I_k$ which was computed at time step $k-\nu$, and we are able to compute $I_{k+\nu}$. At the beginning of time step $k$, the values of $I_{k+i}$ $0 < i < \nu$ are in storage.

Since any value of $I$ is only needed every $\nu$ time step, it is very convenient to implement storage and computation of it as a stack with $\nu$ container. At the beginning of time step $k$, the back of the stack contain $I_k$, the front contain $I_{k+\nu-1}$. The last container is then removed from the stack, the value $I_k$ is used to compute the fluxes, and together with $v_k$ to compute $I_{k+\nu}$ whose value is then put at the front of the stack. Note that the value $I_{k+\nu}$ won't be used until time step $k+\nu$ when it will be at the back of the stack.

<figure><img src="image_2.png" /><figcaption>Figure 1. state of the queue at begin (top) and end (bottom) of time step $k$</figcaption></figure>

This operation have to be done for each degree of freedom on the liner boundary, and they can be done for all this degree of freedom at once. Each container of the stack will in practice contain all the $I$ term for each degree of freedom. For example, for a mesh of tetrahedron, all the triangle classified on boundary faces attached to liner need to be marked as such. For each of this triangle, the field is described by the restriction on this triangle of the field on the tetrahedron to which this triangle is attached to (only one tetrahedron here, since the triangle is on the boundary). In terms of lagrange or hierarchic shape functions basis, this restriction can be computed by getting from the shape functions and coefficients of the tetrahedron on the field those attached to the face we consider. Such operator already exist in the code since that's exactly what we do to extract "left" and "right" degrees of freedoms to compute the normal flux on an interface between two elements. This said, we need in case of order $p$ shape functions on the tetrahedron to keep $(p+1) * (p+2)/2$ values of degree of freedom per Triangle on the liner, per field variable (here just $v \cdot n$) for each time step we need to keep a memory of. So to keep memory of velocity on the liner for n time step, we need an extra storage of : $\nu * n_l(p+1) * (p+2)/2$ double. where $n_l$ is the number of boundary elements on the liner. Since the number of boundary elements on the liner is typically much smaller than the number of elements to model a full engine, this extra storage is negligible.
## B. inertia term

The inertia term need to be treated carefully. Indeed, at each time step, the explicit time marching procedure consist of computing the residual for given value of the field at time t, and use it has the time derivative of the fields to update the variables. The time derivative of the field are there for not known explicitly.

A first naive implementation could approximate the time derivative of the field, since they are not known explicitly when the boundary flux is computed (In fact this flux is computed in order to evaluate this same derivative). Therefore, the problem become nonlinear calculated this way and could lead to instability. If no formal proof of the above hypotheses are provided in this paper, all the numerical experiments that we did in this way could not converge and were blowing up quite quickly. The inertia term need to be to take it into account implicitly to obtain stable solutions. Keeping the time derivatives $v'_n$ as a unknown in the computation of the flux integral, an additional mass matrix appear for each element that share one of his boundary with the discretization of the liner.

f is split into an explicit part, which depend only on u, v and w, and a part that will be treated implicitly. With, $f_u = n_0(R + \beta)(un_0 + vn_1) + 2\beta I + m(u'n_0 + v'n_1 + w'n_2)$ let's call $L_u$ the inertia term of the flux $f_u$.

$$L_u = n_0 m(u'n_0 + v'n_1 + w'n_2) \quad (47)$$

In the framework of the DG scheme, $L_u$, multiply by the test function, is integrated over the liner boundary :

$$R_l = \int_l L_u \hat{u} \, ds + \int_l L_v \hat{u} \, ds + \int_l L_w \hat{u} \, ds \quad (48)$$

$L_u$ is then approximated using the shape function by :

$$L_u = (L_u)_k \phi_k \quad (49)$$

where the k are the indices of nodes on the boundary face. Together with the shape function, and after integration, we get the discrete term associated to a face on the boundary :

$$\begin{pmatrix} \mathbf{L}_p \\ \mathbf{L}_u \\ \mathbf{L}_v \\ \mathbf{L}_w \end{pmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & n_0 n_0 \mathbf{M}_l s_l & n_0 n_1 \mathbf{M}_l s_l & n_0 n_2 \mathbf{M}_l s_l \\ 0 & n_1 n_0 \mathbf{M}_l s_l & n_1 n_1 \mathbf{M}_l s_l & n_1 n_2 \mathbf{M}_l s_l \\ 0 & n_2 n_0 \mathbf{M}_l s_l & n_2 n_1 \mathbf{M}_l s_l & n_2 n_2 \mathbf{M}_l s_l \end{bmatrix} \begin{pmatrix} \mathbf{p}'_l \\ \mathbf{u}'_l \\ \mathbf{v}'_l \\ \mathbf{w}'_l \end{pmatrix} \quad (50)$$

Where $(\mathbf{M}_l)_{i,j} = \int_{reftri} \phi_i \phi_j ds$, and reftri is the reference triangle. $s_l$ the surface of the triangle considered. Note that the matrix $\mathbf{M}_l$ is independent of the element considered. To be able to add this term to the equations to solve for each element, we have to expand the matrix $M_l$, which is defined on the $n_l$ degree of freedom of a boundary triangle, to a matrix $M_l^e$ defined on the $n_e$ degrees of freedom of the tetrahedron it is attached to. This is much like a classical assembly procedure, which is perform at the element level, to assemble the degrees of freedom associated to one boundary of the element to the degree of freedom of the whole element. The assembly procedure will depend only on the position and orientation of the face with regard to the tetrahedron it is linked to. In case the Tetrahedron case, there is 24 possible orientation/position and so as much possible $M_l^e$.

The inertia term of the normal flux of u now depend on $u'$, $w'$, $v'$ which are the unknown of the system we wanted to solve. The four decoupled linear system 26 on each element are now coupled for elements having a boundary on the liner. and we have to solve.

$$\begin{bmatrix} \mathbf{M} & 0 & 0 & 0 \\ 0 & \mathbf{M} + n_0 n_0 \mathbf{M}_l^e \frac{s_l}{v_e} & n_0 n_1 \mathbf{M}_l^e \frac{s_l}{v_e} & n_0 n_2 \mathbf{M}_l^e \frac{s_l}{v_e} \\ 0 & n_1 n_0 \mathbf{M}_l^e \frac{s_l}{v_e} & \mathbf{M} + n_1 n_1 \mathbf{M}_l^e \frac{s_l}{v_e} & n_1 n_2 \mathbf{M}_l^e \frac{s_l}{v_e} \\ 0 & n_2 n_0 \mathbf{M}_l^e \frac{s_l}{v_e} & n_2 n_2 \mathbf{M}_l^e \frac{s_l}{v_e} & \mathbf{M} + n_1 n_2 \mathbf{M}_l^e \frac{s_l}{v_e} \end{bmatrix} \begin{pmatrix} \mathbf{p}' \\ \mathbf{u}' \\ \mathbf{v}' \\ \mathbf{w}' \end{pmatrix} = \begin{pmatrix} \mathbf{R}_p \\ \mathbf{R}_u \\ \mathbf{R}_v \\ \mathbf{R}_w \end{pmatrix} \quad (51)$$

The Mass matrix is therefore not independent of the element any more, and the system is coupled between $u'$ and $v'$. For each element charing a boundary with the liner, a special mass matrix have to be computed and stored.

For all other element the matrix are the same, and therefore we can compute all the solution at once for these element and keep BLAS level 3 efficiency.¹¹ For all the element on the boundary we are stuck with
<table><caption>Table 1.</caption><thead><tr><th>Case</th><th>Δt</th><th>ν</th><th>ε</th><th>R</th><th>β</th></tr></thead><tbody><tr><td>1</td><td>0.004</td><td>190</td><td>0.87</td><td>0.0140774</td><td>1.3428</td></tr><tr><td>2</td><td>0.004</td><td>190</td><td>0.1</td><td>0.899</td><td>0.944</td></tr><tr><td>3</td><td>0.004</td><td>190</td><td>0.0001</td><td>0.999</td><td>0.939</td></tr><tr><td>4</td><td>0.004</td><td>248</td><td>0.0001</td><td>0.9908</td><td>0.025138</td></tr></tbody></table>

BLAS level 2. Usually the number of element sharing a boundary with the liner is small enough so that the impact on performance and memory size is not to huge.

## C. validation

### 1. convergence to design impedance

The aim of the first test case is to validate that the design impedance indeed correspond to the numeric impedance. On a rectangular domain, we impose an incoming plane wave coming from the left. The top and bottom boundary of the rectangle are hard wall. The right boundary is a soft wall modeled with the liner model. Pressure and normal velocity are sampled in time. The Incoming mode 0 is set at frequency $f = 1$ (plane wave). The length of the domain is a of one wave length (since $\rho_0$ and $c$ are set to one, wave length is one). We choose a design impedance such as $Z(2\pi f) = 1 + i$. We fix our time step to $\Delta t = 0.004$, sufficient to assure stability for the explicit Runge Kutta time stepping with $p$ going from 0 to 4. The Coefficient of the Extended Helmholtz resonator model can then be selected to fit design impedance.⁷ $m$ is set to 0, $\nu$ must be chosen between 125 and 249 $\epsilon$ can then be chosen between 0 and 0.868793. Once this choices are made, $R$ and $\epsilon$ are uniquely defined. We report the result for four choices of the parameters summarized in table 1. On figure 2 pressure as a function of time is plotted for case 1 and 4. Case 4 takes much more time to converge to the harmonic solution.

A convergence analysis in time is reported on table 1. The effective impedance is computed by Fourier transforming $v(t)$ and $p(t)$ on sampling on a time interval of $T = 1/f$, starting at $T$, $2T$, $3T$, $4T$, $5T$ and $100T$ for each of the 4 choice of the parameter. For case 1 to 3 we observe a fast convergence to harmonics solution. When the highest the value of the damping parameter, the faster is the convergence. Case 4 is extremes : it uses a high number $\nu$, far from the recommended values and convergence is very slow, since we need to wait $100T$ to be close to the design impedance.

<figure><img src="image_3.png" /><figcaption>Figure 2. $p(t)$ with $Z = 1 + i$, left : case 1, right case 4</figcaption></figure>

This first experiment show that the time domain Helmholtz resonator as implemented is able to model a given design impedance.
<table><caption>Table 2.</caption><thead><tr><th>Case</th><th>T</th><th>2T</th><th>3T</th><th>4T</th><th>5T</th><th>100T</th></tr></thead><tbody><tr><td>1 Z</td><td>1.533+0.321i</td><td>0.975+0.979i</td><td>0.999+1.000i</td><td>1+i</td><td>1+i</td><td>1+i</td></tr><tr><td>2 Z</td><td>2.098+0.507i</td><td>1.028+1.185i</td><td>0.956+0.984i</td><td>1.005+0.995i</td><td>1+i</td><td>1+i</td></tr><tr><td>3 Z</td><td>2.214+0.563i</td><td>1.008+1.229i</td><td>0.936+0.970i</td><td>1.011+0.992i</td><td>1+i</td><td>1+i</td></tr><tr><td>4 Z</td><td>1.024+0.1i</td><td>1.076+0.003i</td><td>1.129+0.008i</td><td>1.183+0.017i</td><td>1.512+0.153i</td><td>0.938+0.939i</td></tr></tbody></table>

## 2. transient time depend 1d case

In time domain, the response of the Dirac wave impacting an impedance wall in one dimension with the Helmholtz resonator model admit an exact solution computed by S. W. Rienstra.¹⁰ Since it's the response in time domain of a non-harmonic problem that we have here, the response for a given design impedance will change if we use different parameters that matches the design impedance. More phenomenon can appear. Therefore we don't limit ourselves to the m = 0 case.

Since the response of a Dirac impulsion can not be solved in our numerical framework, we will solve for a Gaussian pulse instead. if the time length of the pulse is small with regard to the time length of the liner (we typical used a factor of 20 in our computations) the response is similar to the one to a pulse.

Again, the one dimensional problem is solved on a two dimensional mesh in our case. A square divided in two element of order p = 3 was used. On the left side of the square, a right going wave is imposed such as $p(t) = e^{-\frac{(t-t_{ref})^2}{\sigma^2}}$ Results compared to the analytical solution are plotted on figure 1. Dcompana, where pressure versus time is plotted for different value of the model parameters, showing a wide range of possible responses.

the analytical solution was computed according to:¹⁰

$$g(t) = \delta(t) - \frac{2}{m} \sum_{n=0}^{N(t)} \exp(-\frac{R+2}{m}t_n - \varepsilon n) L_n^{(-1)}(2t_n/m) \quad (52)$$

where $N(t) = tc_0/2L$ and $L_n^{(-1)}$ is a generalized Laguerre polynomial :

$$L_n^{-1}(x) = \sum_{k=0}^{n} \frac{(-1)^k}{k!} \begin{pmatrix} n-1 \\ k-1 \end{pmatrix} x^k \quad (53)$$

For all the tested parameters set, the transient solution match very well with the analytical solution, which validate the implementation for the case m > 0.

## 3. harmonic two-dimensional case

To validate the code on a more complex problem, we show result for the following model problem : a plane wave is generated at the left of a rectangular domain has shown on figure 4. On the first half of the rectangle, hard wall boundary condition are applied, on the second half, soft wall with typical impedances are modelized using the time domain impedance model. An incoming mode is imposed on the left of the domain and an outlet characteristic boundary condition is imposed on the right of the computation domain. To avoid spurious reflection on both side, two buffer zones, where a version of the "PML" (perfectly matched layer¹²), are connected to the area of interest. Briefly, a dissipative source term, proportional to the difference between the actual field and a target field, multiplied by a damping factor is added to the residual. The damping factor vary smoothly from 0 on the boundary of the area of interest to 5 at the end of the buffer zone. The target field is 0 on the outlet side and is set to the value of the chosen duct mode on the inlet side.

For all the presented results we have a sound speed, $c_0 = 1$ a mean flow density of $\rho_0 = 1$, the plane wave is of pulsation $\omega == 2\pi$, and the length $L = 5$. The impedance of the liner is set to $Z = 1 + i$. We use polynomial order $p = 6$ for the shape functions, and the mesh is a structured 10 * 42 grid of rectangles divided in triangles.
<figure><img src="image_4.png" /><figcaption>Figure 3. Comparing DG Result versus Analytical results.</figcaption></figure>
Results are compared with those obtained by the Actran/TM code, a tool developed by Free Field Technologies. Briefly, Actran/TM¹³ computes the propagation of harmonic acoustic disturbances in non-uniform, homotropic mean flow and their radiation to the far field. Using Myers formulation of acoustic boundary conditions the effect of acoustic liner is accurately accounted for. Strong coupling with duct acoustic modes allows one to both specify incident modes and compute the amplitude of reflected ones. The algebraic system of equations resulting from the finite and infinite element discretization of the acoustic potential equation in frequency domain is solved with a parallel, out-of-core, direct solver. For the frequency domain solution, no buffer zone are needed, instead boundary condition are applied by imposing incoming duct mode.

Figure 5 show the pressure map for both Actran/TM and DG results. Figure 6 show the result along the line y = 2.5 (axis of symmetry of the duct) and y = 1 near the liner. Again, the results are very close which validate our implementation considering that the two methods are completely different.

<figure><img src="image_5.png" /><figcaption>Figure 4. DG (top) and Actran (bottom) problem setup</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Pressure map, top Actran/tm result, bottom DG results</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Comparing DG (black) Result versus actran/TM (red) pressure plot along $y = 1$ (top and $y = 2$ bottom).</figcaption></figure>
### III. mean-flow case

If a mean flow exist on top of which the linearized Euler Equation are solved, the liner boundary condition need to be corrected in order to take into account the shear-layer associated along the liner. The classical way to take into account the mean flow is to apply the Ingard-Myers$^{8,9}$ limit of a vanishing mean flow boundary layer along a lined wall. The equation for the liner, developed for the case of the extended Helmholtz resonator become:

$$p'(t) = (R + \beta)v'(t) + mv''(t) + 2\beta \sum_{n=1}^{\infty} e^{-\epsilon n} v'(t - n\nu\Delta t) + \mathbf{v}_0 \cdot \nabla p \quad (54)$$

Note that we have neglected here the curvature term $\mathbf{v}_0 \cdot (\mathbf{n} \cdot \nabla n)p$, which is valid if the curvature on the liner is not too large. The equation is now written in terms of time derivatives of $p$ and $v$ (here $v$ is the normal velocity to the wall). The approach that we used in the no-flow case which consisted in computing the explicit part of the flux and add it to the residual is no longer possible. Instead, the contribution of term depending on derivative are computed implicitly, the same way as we did in the no mean flow case for the inertia term : by coupling the system of four equations though a modified mass matrix. Since the inertia term is not needed to obtain the harmonic solution we will drop it from the development. Keeping it is possible, one just need to keep memory of the derivative of normal velocity of one preceding time step and compute the implicit contribution of the acceleration using a backward Euler. What need to be discussed is the computation of the contribution of the convective term $\mathbf{v}_0 \nabla p$

#### A. convective term

We need to compute $v_0 \cdot \nabla p$ on the liner. $p$ is defined on the three dimensional mesh as a discontinuous finite element field. A first approach could be to compute the gradient every where and then scale with the velocity. But since the convective term is needed only on the liner, we need only to compute it only on it's surface, and since $v_0 \cdot n = 0$, we only need to compute the projection of the gradient on the tangent plane to the liner surface.

In order to compute the gradient, we search the discrete field $Gx$ = and $Gy$ and $Gz$ defined on the liner surface $s$ such that on each surface element $s_e$ describing the surface of the liner, we have $G_x = G_{xk}^e \phi_k$, $G_y = G_{yk}^e \phi_k$ and $G_z = G_{zk}^e \phi_k$:

$$\mathbf{G}_{\mathbf{x}_k}^e \int_{s_e} \phi_k \phi_j - \partial_x p \phi_j ds = 0 \quad (55)$$

where $p$ is approximated using computed value on the elements at previous time step. Remember that the discretization of $p$ is itself discontinuous. Integration by part give a system of equation to solve on each surface element $s_e$:

$$\mathbf{M}_{el} \mathbf{G}_{\mathbf{x}_k}^e = \int_{s_e} p_i \phi_i \partial_\xi \phi_j d\xi d\eta + \sum_k \int_{s_{ek}} \mathrm{tr}[p_i] \phi_i \phi_j n_0 d\xi d\eta \quad (56)$$

Where $\xi$ and $\eta$ are local coordinate on liner element $s_e$, $s_{ek}$ is the $k^{th}$ boundary of liner element $s_e$ and $\mathrm{tr}$ is the trace of the discontinuous discrete $p$ on the edge. In our case we choose $\mathrm{tr}[p] = (p_L + p_R)/2$ where $p_L$ and $p_R$ are respectively the value of the discrete pressure at right and left of the $k^{th}$ boundary of $s_e$. Solving this system on each liner element yield discrete value for $G_x$. Similar systems are constructed for $G_y$ and $G_z$ to compute the other components of the gradient.

The algorithm to compute the gradient on the surface of the liner was tested as well as the fully implicit version of the boundary condition without mean flow, both giving satisfactory result, enough to move on to the computation of the a lined duct simulation with constant mean flow.

#### B. Presence of instability

To test the implementation we first did the same test case as in section 3, at the only difference that we add a constant meanflow of $u_0 = -0.5$, going from right to left. When the problem is under resolved, we are in general able to get a solution. However, when the resolution is increased, either by increasing the polynomial order or by increasing the number of elements, instability starts to appear. Figure 7 and 8 show
examples of these instabilities. We display the pressure map after 15s and 20s of simulation along with the pressure curve along the top side of the duct. A t 15s the instability starts to kick in, and jump in the pressure appear near the liner edge. At 20s the instability is fully developed at greatly pollute the solution. Attempts where made to stabilize the scheme by changing the computation of the trace of the pressure between two elements. We tried different up winding schemes to no significant effects. It seems in fact that this instability is not an artifact of the numerical scheme but is inherent to the Ingard-Myers correction to mean flow. Analytical development¹⁴⁻¹⁶ seems to show serious indications that the boundary conditions coming from the limit of a vanishing mean flow boundary layer along a lined wall allows for the presence of instabilities for any frequency. Özyörük¹⁷ reported the same type of instabilities in his numerical experiment with mean flow. It seems that the instabilities is triggered by the short wave length component that appear in time domain at the transition between hard and soft wall. Discussion about the validity of the boundary condition are clearly beyond the scope of this paper. We remind that our goal is to be able to reproduce frequency domain result in mean flow case with lined duct. In frequency domain result the instability is never triggered, and we need to be able to filter out this instability from time domain result to be able to do meaning full comparisons. In the next section we therefore propose a practical way to do just that.

### C. Smoothing operator

In order to stabilize the scheme in presence of mean flow, we made the empirical assumption that the short wave component of the pressure gradient in the convective term were responsible for the triggering of the instability. To verify this assumption we develop a smoothing procedure, valid for discontinuous Galerkin solution on unstructured meshes.

The field we are interested in are defined on a mesh. On each element of the mesh the field is defined as a lagrange polynomial of relatively high order (typical order range from 1 to 10 in our simulation), and the field is discontinuous between element. The overall strategy here is to replace the value of the approximation at each interpolation point by an approximation of the mean value of the field integrated on a disk of given radius $r_0$. The radius serve as the filtering parameter, any signal component of wave length much higher than $r_0$ will be almost not affected, short wave length will be smeared out.

Around one interpolation point, as shown in figure 9, the disk span on several element, and therefore, the function we which to integrate is discontinuous. Integration formula on a disk are therefore useless here. If an element if fully inside the circle, gauss quadrature rule can give the exact contribution to the integral for this element. For an element that intersect the boundary of the disk, we could try to use a gauss quadrature rule on this element and set the contribution to any gauss point outside the disk to zero. In this case, it's effectively trying to compute the integral of a discontinuous function inside an element, and we therefore be better of by using a

Monte-Carlo method. Another solution would be to precisely triangulate the intersection of a disk with an element and then

integrate over each sub-element. This would probably be over-kill and have a computation price we don't want to pay. The solution we propose here is to replace the integration of the discontinuous function on element that cut the circle be the integration of the field we want to integrate multiplied by a smooth blending function $b(r)$ such as the integral of $b$ over the full surface of the liner is one, value of $b$ inside is close to 1, value outside is tend to 0 when $r$ increase. Typically we choose :

$$b(r) = \frac{1}{2} - \frac{1}{\pi} atan((r - r_0)/\delta) \qquad (57)$$

with typical values of $\delta$ which represent the size of the blending zone that range between $r_0/5$ and $r_0/10$.

This procedure is general enough to treat all the cases we need to address. In case of a structured mesh on the liner, a simpler smoothing procedure could be devised.

### VALIDATION OF THE SMOOTHING PROCEDURE
To validate the above procedure and experiment the parameters, we set up the following problem : on the unit square, discretized by an unstructured mesh of 100 elements, and a polynomial order of $p = 5$ we set up by an $L_2$ projection an initial solution of the following form:

$$p = \cos(2\pi X) + \cos(10\pi X)$$
So that we have a main signal of wave length 1 superposed to an other signal of wave length 1/5. We test the smoothing procedure to see if we are able to filter out the small wave length while keeping the component of wave length 1. The plot of figure 9 reproduce the pressure value on an horizontal line at y = 0.5 crossing the unit square. On each plot, we have for different value of r₀, the smoothing radius, the initial solution, the "target" solution, e.g. the solution we wish to obtain after filtering out(p = cos(2piX)), and the solution for one and two path of the smoothing procedure. Best result are obtain by setting the radius to half the wavelength we wish to filter out and applying 2 times the filter.

## D. Two dimensional duct with mean flow : comparison with actran

The computation in section B that failed before due to the presence of an instability was completed again, but this time, the smoothing operator was applied on the gradient of p before adding its contribution to the residual. The smoothing operation permitted to get rid of the instability. The hypothesis that the short wave length component of the gradient of p trigger the instability seems to be verified numerically. Results are shown on figures 11 and 12 and compared to frequency domain results obtain with Actran/tm. We achieve very good agreement with the Actran/tm results. It's seems that at least in this particular case and all the cases that we have tested so far this simple smoothing procedure, only applied on the discretized $\nabla p$ is enough to prevent the triggering of the instability.

## E. two-dimensional nacelle with mean flow

To assess the efficiency and flexibility of the code and therefore to show that the implementation is ready to solve industrial test case, we solved the two dimensional nacelle problem depicted on figure 13, which describe the geometry and the mean-flow velocity distribution. Figure 14 show the pressure map near the fan for a lined and a hard walled case. Note that directivity and pressure maximum value are modified due to the presence of the liner.

## Conclusion

In this paper we have shown how to implement the extended Helmholtz resonator model in the context of the DGM both with and without mean flow. Instability in the presence of mean-flow were encountered in time domain, as reported by several authors. Those instabilities are usually not triggered in frequency domain, and we provided a way to filter them out to be able to obtain the same results as in frequency domain. From our numerical experiment it seemed that it was the short wave length component of the convective term in the Ingard-Meyers boundary conditions, that triggered the instabilities. Since the wave length that trigger the instabilities are smaller than the well resolved wave length of interest, we propose to filter them out of the convective term using a gaussian smoothing procedure. Other filtering techniques could be employed, in particular it is certainly possible to Fourier transform in space the convective term and filter out the high wave number in the Fourier space before projecting back to physical domain. In the chapter of interesting extensions, we will consider in the future to implement the curvature term that we decided to neglect from the Ingard-Meyers boundary conditions and investigate its influence on the results.

## Acknowledgments

This work was carried out under the European collaborative project "Messiaen" under the Sixth Framework Programme (Project Co-coordinator Jean-Louis Migeot). The financial support is greatly acknowledged. The present paper is the companion of the paper entitled "Impedance Models in Time Domain including the Extended Helmholtz Resonator Model" by Sjoerd W. Rienstra,⁷ to be presented at the present conference, The authors would like to express their special thanks to Sjoerd W. Rienstra and the partners of the Messiaen project. Support also came in the form of top computing facilities. In this connection, we would like to thank the CENAERO research center (www.cenaero.be).
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Pressure at $t = 15s$ and pressure along the liner at $Y = 5$.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Pressure at $t = 20s$ and pressure along the liner at $Y = 5$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9. sketch of the smoothing procedure</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10. Pressure in space</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11. Pressure at $t = 30s$ top: DG results, bottom: Actran results</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12. Pressure along a cut line in the duct. Red curve: Time domain DG Results. Black curve: Frequency Domain Actran TM Results top: $Y=1$. bottom: $Y=2.5$</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 13. two dimensional nacelle problem set up. geometry and Mach number</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14. pressure map near the fan left : lined wall, right hard wall</figcaption></figure>
References

¹Cockburn, B., Karniadakis, G., and Shu, C.-W., *Discontinuous Galerkin Methods*, Vol. 11 of *Lecture Notes in Computational Science and Engineering*, Springer, Berlin, 2000.

²Cockburn, B. and Shu, C.-W., "Runge-kutta discontinuous galerkin methods for convection-dominated problems," *Journal of Scientific Computing*, Vol. 16, 2001, pp. 173-261.

³Remacle, J.-F., Flaherty, J. E., and Shephard, M. S., "An Adaptive Discontinuous Galerkin Technique with an Orthogonal Basis Applied to Compressible Flow Problems," *SIAM Review*, Vol. 45, 2003, pp. 53-72.

⁴Peyret, C. and Delorme, P., "Optimal hp Discontinuous Galerkin Method Application to Computational Aeroacoustics," *11th AIAA/CEAS Aeroacoustics Conference*, 2005.

⁵Atkins, H. L. and Shu, C. W., "Quadrature-Free Implementation of Discontinuous Galerkin Method for Hyperbolic Equations," *AIAA Journal*, Vol. 36, No. 5, 1997, pp. 775-782.

⁶Chevaugeon, N., Remacle, J.-F., Gallez, X., Ploumhans, P., and Caro, S., "Efficient Discontinuous Galerkin Methods for solving acoustic problems," *11th AIAA/CEAS Aeroacoustics Conference*, 2006, AIAA-2006-2686.

⁷Rienstra, S. W., "Impedance Models in Time Domain including the Extended Helmholtz Resonator Model," *12th AIAA/CEAS Aeroacoustics Conference*, 2006, AIAA-2006-2686.

⁸Ingard, K., "Influence of Fluid Motion Past a Plane Boundary on Sound Reflection, Absorption, and Transmission." *Journal of the Acoustical Society of America*, Vol. 31(7), 1959, pp. 1035-1036.

⁹Myers, M. K., "On the acoustic boundary condition in the presence of flow." *Journal of Sound and Vibration*, Vol. 71 (3), 1980, pp. 429-434.

¹⁰Rienstra, S., "1D Reflection at an Impedance Wall," *Journal of Sound and Vibration*, Vol. 125, 1988, pp. 43-51.

¹¹Dongarra, J. J., Du Croz, J., Duff, I. S., and Hammarling, S., "A set of Level 3 Basic Linear Algebra Subprograms," *ACM Trans. Math. Soft.*, Vol. 16, 1990, pp. 1-17.

¹²Bérenger, J. P., "A Perfectly Matched Layer for the Absorption of Electromagnetic Waves," *Journal of Computational Physics*, Vol. 114, 1994, pp. 185-200.

¹³Ploumhans, P., Meerbergen, K., Knapen, T., Gallez, X., Lielens, G., and Coyette, J.-P., "Development and validation of a parallel out-of-core propagation and radiation code with validation on a turbofan application," *Proceedings of the Eighteenth International Congress on Acoustics. ICA*, 2004.

¹⁴Rienstra, S. W., "A classification of Duct Modes Based on Surface Waves," *Wave motion*, Vol. 37, No. 2, 2003, pp. 119-135.

¹⁵Rienstra, S. W. and Peake, N., "Modal Scattering at an Impedance Transition in a Lined Flow Duct," 2005, 11th AIAA/CEAS Aeroacoustics Conference.

¹⁶Koch, W. and Möhring, W., "Eigensolutions for liners in uniform mean flow ducts," *AIAA Journal*, Vol. 21, 1983.

¹⁷Özyörük, Y. and Long, L. N., "A time-domain implementation of surface acoustic impedance condition with and without flow," *2nd AIAA/CEAS Aeroacoustics Conference*, AIAA, 1996.