Florian Monteghetti
email: florian.monteghetti@onera.fr

Denis Matignon
email: denis.matignon@isae.fr

Estelle Piot
email: estelle.piot@onera.fr

Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized Euler equations

Keywords:

Time-domain impedance boundary conditions (TDIBCs) can be enforced using the impedance, the admittance, or the scattering operator. This article demonstrates the computational advantage of the last, even for nonlinear TDIBCs, with the linearized Euler equations. This is achieved by a systematic semi-discrete energy analysis of the weak enforcement of a generic nonlinear TDIBC in a discontinuous Galerkin finite element method. In particular, the analysis highlights that the sole definition of a discrete model is not enough to fully define a TDIBC. To support the analysis, an elementary physical nonlinear scattering operator is derived and its computational properties are investigated in an impedance tube. Then, the derivation of time-delayed broadband TDIBCs from physical reflection coefficient models is carried out for single degree of freedom acoustical liners. A high-order discretization of the derived time-local formulation, which consists in composing a set of ordinary differential equations with a transport equation, is applied to two flow ducts.

Introduction

Across various industries, stringent regulations on sound pollution have led to a growing usage of sound absorbing materials, the design of which requires suitable computational tools. An accurate computation of the acoustic field inside an industrial sound absorbing material is costly by contemporary computing standard, and is thus the realm of specific studies that focus on the material in isolation. Practical computations of sound absorption are typically done by abstracting the geometrical features of the material using a so-called impedance boundary condition (IBC). The computation of sound absorption using an IBC is one aspect of computational aeroacoustics (CAA) [START_REF] Tam | Computational Aeroacoustics: A Wave Number Approach[END_REF]. Broadly, this paper deals with the discretization of IBCs in the time domain for the linearized Euler equations. This introduction is split into two parts: Sec. 1.1 covers impedance modeling and its applications, while Sec. 1.2 focuses on the time-domain discretization of IBCs and presents the article objectives and outline.

Modeling and applications of impedance boundary conditions

Many sound absorbing materials are variations of the acoustical (Helmholtz) resonator. Assumed locally reacting, they can be modeled using a (continuous) single-input single-output time-invariant operator known as the impedance operator [START_REF] Morse | Theoretical Acoustics[END_REF]§ 6.3] [3,Chap. 10]. At low incident sound pressure level (SPL), the impedance is linear and thus reduces to a timedomain convolution [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF]§ III.3]; it is a hereditary (i.e. non time-local) operator whose kernel exhibits a long memory behavior linked to viscous and thermal effects within both the neck and the cavity. (See [5, § III] and references therein for common linear models.) Above a SPL threshold, a turbulent dissipation mechanism that occurs at the resonator neck, known as vortex shedding, results in the addition of a nonlinearity [6, § 4] [7, § I] [8, § 2] [9, § 5]. The influence of a subsonic grazing flow can be modeled through an additional parametric dependence that does not change the mathematical nature of the impedance operator [10, § 3] [11, § 7].

The effect of a locally reacting absorbing material on the sound field is accurately described, after a homogenization distance [START_REF] Laurens | Effective conditions for the reflection of an acoustic wave by low-porosity perforated plates[END_REF], by using its impedance operator as a boundary condition. The effectiveness of this modeling and the advent of acoustically treated jet engines has led to IBCs becoming a staple part of aeroacoustics [START_REF]Noise Control, NASA[END_REF]Chaps. 13 & 14]. Early works focused on the prediction of sound absorption in a duct with flow: Cremer [START_REF] Cremer | Theory regarding the attenuation of sound transmitted by air in a rectangular duct with an absorbing wall, and the maximum attenuation constant produced during this process[END_REF] derived an approximation of the optimal impedance and Pridmore-Brown [START_REF] Pridmore-Brown | Sound propagation in a fluid flowing through an attenuating duct[END_REF] established his celebrated equation, widely used to compute duct modes. The identification of unstable surface modes by Tester [START_REF] Tester | The propagation and attenuation of sound in lined ducts containing uniform or 'plug' flow[END_REF] lead to a wealth of investigations focused on the hydrodynamic stability of a base flow with an IBC [START_REF] Khamis | Viscous effects on the acoustics and stability of a shear layer over an impedance wall[END_REF]. Inverse methodologies have been developed to identify the IBC and provide an alternative to more intrusive measurement techniques [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF]. In all of these studies, the IBC is linear and the Cauchy problem is formulated in the frequency domain. Although less popular than their time-harmonic counterparts, time-domain IBCs (TDIBCs) have been used in wave propagation problems including duct aeroacoustics [START_REF] Özyörük | Time-domain numerical simulation of a flow-impedance tube[END_REF][START_REF] Bin | Broadband impedance boundary conditions for the simulation of sound propagation in the time domain[END_REF][START_REF] Liu | Stability analysis and design of time-domain acoustic impedance boundary conditions for lined duct with mean flow[END_REF][START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF], room acoustics [START_REF] Botteldooren | Finite-difference time-domain simulation of low-frequency room acoustic problems[END_REF], as well as outdoor sound propagation [START_REF] Cotté | Time-domain simulations of sound propagation in a stratified atmosphere over an impedance ground[END_REF]. Richter et al. [START_REF] Richter | Liner Impedance Modeling in the Time Domain with Flow[END_REF][START_REF] Richter | A review of time-domain impedance modelling and applications[END_REF] and Troian et al. [START_REF] Troian | Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow[END_REF] identified an IBC in the time domain. Gabard & Brambley [START_REF] Gabard | A full discrete dispersion analysis of time-domain simulations of acoustic liners with flow[END_REF] used a time-domain formulation to investigate the (in)stability of the Ingard-Myers boundary condition.

The use of a time-domain IBC may merely be more convenient than a time-harmonic one, for instance when a broadband source is considered. However, it is indispensable when the harmonic problem cannot be formulated, such as in the presence of time-dependent domains, moving sources, or nonlinearities. Typically, the nonlinearity comes from the impedance operator, as mentioned above, or from the considered partial differential equation (PDE). Nonlinear PDEs are of particular interest when hydrodynamic phenomena cannot be neglected, such as close to a supersonic fan tip where shocks occur [29, § 5.3], or for flow control. Scalo et al. [START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF] and Olivetti et al. [START_REF] Olivetti | Direct numerical simulation of turbulent flow with an impedance condition[END_REF] performed a Navier-Stokes simulation (large eddy simulation and direct numerical simulation, respectively) with a linear TDIBC to investigate the interaction between a turbulent boundary layer and an impedance wall.

It is worth noting that IBCs have other uses and denominations. In simulations of combustion chambers, IBCs are used to truncate a part of the chamber [START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF][START_REF] Douasbin | Delayed-time domain impedance boundary conditions (d-tdibc)[END_REF] or model injectors [START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF][START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF] for instance. In mathematical control, TDIBCs (also known as Dirichlet to Neumann maps) are commonly used to stabilize the wave equation since they modify the underlying semigroup generator and can yield asymptotic or even exponential stability. Both finite-dimensional [START_REF] Abbas | The multidimensional wave equation with generalized acoustic boundary conditions I: Strong stability[END_REF] and infinite-dimensional TDIBCs, involving for example fractional derivative [START_REF] Grabowski | Stabilization of wave equation using standard/fractional derivative in boundary damping[END_REF] and time delay [START_REF] Wang | Wave equation stabilization by delays equal to even multiples of the wave propagation time[END_REF], have been studied. Surface or generalized IBCs [START_REF] Yuferev | Surface Impedance Boundary Conditions: A Comprehensive Approach[END_REF] are employed in electromagnetism to model non-perfect electric conductors. The first and most widely used model is that of Leontovich, which is a fractional operator identical to that encountered in the acoustic impedance of perforated plates [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF]. See [START_REF] Chun | High-order accurate thin layer approximations for time-domain electromagnetics. Part I: General metal backed coatings[END_REF][START_REF] Hiptmair | Fast convolution quadrature based impedance boundary conditions[END_REF] for two applications in the time domain.

Time-domain discretization of impedance boundary conditions

A discrete TDIBC consists of three components. First, the discrete impedance model, i.e. the finite-dimensional operator (in the sense of systems theory [START_REF] Brockett | Finite Dimensional Linear Systems[END_REF][START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF]) that one wishes to apply at the boundary. Second, the numerical algorithm used to evaluate the said operator; in the case of linear continuous time-invariant operators (i.e. of linear TDIBCs) this amounts to computing a time-domain convolution. Practically these first two elements go hand in hand since the expression of the convolution kernel dictates how the convolution can be efficiently computed. The third and last component is the (semi-)discrete formulation, i.e. how the TDIBC is enforced at the (semi-)discrete level.

The first two components are well-documented for linear TDIBCs. First, broadband discrete impedance models can be built from the so-called oscillatory-diffusive representation of physical models [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF]. Second, a convenient way to compute a causal convolution is to recast it into a set of ordinary differential equations (ODEs) (or, in the sense of systems theory, to realize the corresponding LTI operator [41, § 17] [42, § 1.2]), a method sometimes named "auxiliary differential equations" [START_REF] Bin | Broadband impedance boundary conditions for the simulation of sound propagation in the time domain[END_REF][START_REF] Dragna | A generalized recursive convolution method for time-domain propagation in porous media[END_REF], "canonical form implementation" [START_REF] Liu | Stability analysis and design of time-domain acoustic impedance boundary conditions for lined duct with mean flow[END_REF][START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF], or "state-space model" [START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF] in the acoustics or fluid dynamics literature. In this way, the computation of a (in)finite-dimensional hereditary operator is reduced to the integration of a (in)finite set of time-local differential equations, which can be done with arbitrary accuracy. If the impedance operator models a reflection phenomenon (e.g. a wave reflection in the cavity of a Helmholtz resonator), its realization involves not only ODEs but also time-delayed ODEs [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF]. Alternatively, when a realization of the impedance operator is beyond reach, more general algorithms can be used, such as convolution quadrature or multistep methods [START_REF] Lubich | Runge-Kutta methods for parabolic equations and convolution quadrature[END_REF][START_REF] Hiptmair | Fast convolution quadrature based impedance boundary conditions[END_REF][START_REF] Sauter | Convolution quadrature for the wave equation with impedance boundary conditions[END_REF]. The key feature of these methods is that they only require the knowledge of the Laplace transform of the convolution kernel, which is always known in the case of linear impedance models (impedance modeling, be it theoretical or empirical, is often carried out in the frequency domain).

The third component of a TDIBC is seldom discussed. For instance, there does not seem to be a consensus on whether a TDIBC should be based on the impedance z [START_REF] Özyörük | Time-domain numerical simulation of a flow-impedance tube[END_REF][START_REF] Rienstra | Impedance models in time domain, including the extended Helmholtz resonator model[END_REF][START_REF] Bin | Broadband impedance boundary conditions for the simulation of sound propagation in the time domain[END_REF][START_REF] Gabard | A full discrete dispersion analysis of time-domain simulations of acoustic liners with flow[END_REF][START_REF] Olivetti | Direct numerical simulation of turbulent flow with an impedance condition[END_REF], the admittance y [START_REF] Liu | Stability analysis and design of time-domain acoustic impedance boundary conditions for lined duct with mean flow[END_REF][START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF], or the reflection coefficient β [START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF][START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF]. This topic has been mentioned by Gabard & Brambley [START_REF] Gabard | A full discrete dispersion analysis of time-domain simulations of acoustic liners with flow[END_REF] who warned that some instabilities reported in the literature may be due to an unsuitable implementation and showed the benefit of a characteristic-based implementation in their study, which is echoed in works that focus on large eddy simulations and direct numerical simulations [START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF][START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF][START_REF] Douasbin | Delayed-time domain impedance boundary conditions (d-tdibc)[END_REF].

Another scarcely mentioned aspect is the impact of the TDIBC on the maximum admissible time step. Although this aspect is of secondary concern for stability studies, it is of crucial importance for large-scale applications like those that involve hydrodynamics or for inverse methods that need to explore the impedance parameter space. To the best of the authors' knowledge, the only known result is that for a proportional impedance (i.e. z(t) ∝ δ(t)), using the reflection coefficient yields a CFL stability condition independent of the impedance coefficient [48, § 3.3] [49, 2.3].

The objective of this paper is to contribute to addressing these gaps by establishing the computational advantage of the reflection coefficient, even for nonlinear TDIBCs, with the linearized Euler equations. This is achieved by a systematic investigation of the weak enforcement of a generic nonlinear TDIBC in a discontinuous Galerkin finite element method. To support the analysis, an elementary physical nonlinear scattering operator is derived and its computational properties are investigated in an impedance tube. Then, for application to two flow ducts, the derivation of broadband TDIBCs from physical reflection coefficient models is carried out for single degree of freedom acoustical liners. The derived formulation is purely time-local and consists in composing a set of ODEs with a transport equation.

This paper is organized as follows. Sec. 2 gathers preliminaries including a rigorous statement of impedance admissibility conditions. Two of three components of a TDIBC, namely the discrete model and the time discretization algorithm, are covered in Sec. 3, which derives both linear and nonlinear discrete models from the analysis of physical impedance models. The third component, namely the (semi-)discrete formulation, is addressed independently in Sec. 4 with an investigation of the weak enforcement of a generic nonlinear impedance boundary condition within a discontinuous Galerkin finite element method. The last two sections gather numerical validations and applications of the presented analysis. Sec. 5 validates the proposed linear and nonlinear TDIBCs by comparison with the impedance tube analytical solution. Sec. 6 presents an application to two flow ducts.

Preliminaries on the continuous formulation

The purpose of this section is to introduce the physical problem considered in this paper, as well as to gather definitions and notations for later use (such as the analysis of Sec. 4). The linearized Euler equations are briefly recalled in Sec. 2.1. Sec. 2.2 gives the three equivalent formulations of an IBC and rigorously states the admissibility conditions. Lastly, Sec. 2.3 presents the assumption made on the base flow and discusses its mathematical and physical implications.

Linearized Euler equations

In this work the hydrodynamic field is split between a steady base flow and an unsteady perturbation. The perturbations of pressure p, velocity u, and density ρ are governed by the homentropic linearized Euler equations (LEEs), defined on

(0, ∞) × with ⊂ R n a bounded open subset, ⎧ ⎪ ⎨ ⎪ ⎩ ∂ t u + (u 0 • ∇)u + z -1 0 c 0 ∇ p + (u • ∇)u 0 + z -1 0 c -1 0 p(u 0 • ∇)u 0 = 0 ∂ t p + (u 0 • ∇)p + z 0 c 0 ∇ • u + γ p∇ • u 0 = 0 p = c 2 0 ρ. ( 1
)
The specific heat ratio is denoted γ > 1. Quantities related to the base flow are designated by the subscript "0": u 0 is the base flow velocity, ρ 0 the base flow density, c 0 the speed of sound, and z 0 := ρ 0 c 0 the characteristic impedance of the propagation medium. All quantities are dimensional. Column vectors are denoted in bold and the symbol (u • ∇) denotes the convective derivative defined as (u

• ∇) f := i u i ∂ i f and (u • ∇) f := ∇ f • u = ( j u j ∂ j f i ) i
when applied to scalar-valued and vector-valued functions, respectively. The symbol ":=" is used to denote a definition. The LEEs (1) entail hypotheses on both the base flow and the perturbations. The base flow is a perfect gas, so that c 2 0 = γ rT 0 = γ p 0 /ρ 0 , and can solve either the steady Euler or Navier-Stokes equations. The perturbations are "small", inviscid, and homentropic. The rather strong homentropicity assumption enables to replace the energy equation by the algebraic relation p = c 2 0 ρ and implies that c 0 must be constant in . As a result, the LEEs (1) do not include an entropy mode, but only hydrodynamic and acoustic ones.

In spite of its simplicity, this model is commonly used in duct aeroacoustics, see Sec. 6. For the detailed derivation of the LEEs and additional physical insights, the reader is referred to [25, Chap. 2] and references therein.

Formulations and admissibility conditions of impedance boundary conditions

On a subset z of the boundary := ∂ , an IBC is applied, defined as

p(t, x) = Z[u(•, x) • n(x)](t) (t ∈ (0, ∞), x ∈ z ), ( 2
)
where n is the outward unit normal and Z is the (single-input single-output) impedance operator. If the impedance operator is linear continuous time-invariant (LTI), the IBC (2) reduces to

p(t, x) = [z u(•, x) • n(x)] (t), ( 3
)
where z is the impedance kernel and " " denotes the time-domain convolution of distributions [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF]Chap. III]. Since the LEEs (1) are linear, the LTI IBC (3) can also be formulated in the Laplace domain, provided that z is causal, as

p(s, x) = ẑ(s) û(s, x) • n(x) ( (s) > c), (4)
where the Laplace transform is defined as f (s) = ∞ 0 f (t)e -st dt for a locally integrable function with finite exponential growth, i.e. f ∈ L

1 loc ([0, ∞)) with | f (t)| ≤ Me ct .
If f is a causal distribution, as encountered with physical impedance models, the distributional extension f (s) = f , e -st is employed [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF]Chap. 6] [50,Chap. 2]. As long as the problem is linear, [START_REF] Kinsler | Fundamentals of Acoustics[END_REF] and ( 4) are equivalent. The Laplace formulation is particularly useful for stability studies (using a Briggs-Bers analysis for instance) or when deriving representations of physical impedance models, as is done in Sec. 3.2. Formally s = jω (with j 2 = -1) yields the Fourier transform but such a mere formal substitution is error-prone, especially when dealing with admissibility conditions, which justifies the present use of the Laplace transform; a discussion on this topic is provided in Appendix A. In the remainder of this paper, "impedance" can designate the single-input single-output operator Z, the convolution kernel z, or its Laplace transform ẑ. The impedance can have a spatial dependency (i.e. Z(x, •), z(x, t), or ẑ(x, s)), although it is not explicitly written for the sake of clarity. In this study, the IBC ( 2) is meant to model a passive sound absorbing material so that it obeys the admissibility conditions given below.

Definition 1 (Admissibility conditions).

Let Z be a single-input single-output operator. It is said to be an admissible impedance operator if it enjoys the following properties: causality; reality (real-valued inputs are mapped to real-valued outputs); passivity, i.e. for every smooth and compactly supported input u ∈ C ∞ 0 (R) and for every time instant t > 0,

⎛ ⎝ t -∞ Z(u)(τ )u(τ ) dτ ⎞ ⎠ ≥ 0, (5)
where the overline denotes the complex conjugate.

Physically, the left-hand side of ( 5) is the energy supplied to the system over (-∞, t). The passivity condition [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF], which applies to every passive system encountered in physics, means that the system does not produce energy [51, p. 301] [52, p. 12]. If Z is LTI, then it is interesting to note that this definition could be restricted to reality and passivity. Indeed, for real LTI systems, passivity implies causality [51, Lemma 10.3] [50, Note 8]. A consequence of this result is that a real anticausal LTI system cannot be passive: intuitively, for such a system the supplied energy over (-∞, t) can be made arbitrarily negative by modifying the input u in the future. A key fact in practice is that the admissibility of an LTI impedance can be readily read on its Laplace transform ẑ, as summarized in Proposition 3 below, which could be interpreted as a rigorous statement of [START_REF] Rienstra | Impedance models in time domain, including the extended Helmholtz resonator model[END_REF]Thm. 1].

Definition 2 (Positive-real function

). Let C + 0 := {s ∈ C | (s) > 0} be the open right half-plane. A function f : C + 0 → C is positive-real if f is analytic in C + 0 , f (s) ∈ R for s ∈ (0, ∞), and [ f (s)] ≥ 0 for (s) > 0.
Proposition 3. An LTI operator Z is an admissible impedance operator if and only if its Laplace transform ẑ(s) is a positive-real function.

Proof. This is a standard fact of Systems Theory. See [52, § 2.11] for the case where the kernel z ∈ L

(ẑ(s) = 1 / √ s); the fractional derivative z(t) = H(t) / √ πt δ (t) (ẑ(s) = √ s);
the time-delay system z(t) = aδ(t) + bδ(tτ ) with a ≥ b ≥ 0 and τ ≥ 0 (ẑ(s) = a + be -sτ ). The following impedance kernels are not admissible: ẑ(s) = √ s /j, which fails the reality condition ẑ(s) ∈ R for s ∈ (0, ∞); ẑ(s) = a + be -sτ with a < b and τ ≥ 0, which is not passive; z(t) = δ(t + τ ) with τ ≥ 0 (ẑ(s) = e sτ ), which is not causal; z(t) = e t H(t) (ẑ(s) = (s -1) -1 ), which is not passive (ẑ(s) is not analytic in C + 0 since it admits 1 as a pole).

As mentioned in the introduction and covered further below, it can be beneficial to consider other forms of the IBC [START_REF] Morse | Theoretical Acoustics[END_REF].

When possible, it may be written using the admittance operator Y as

u(t, x) • n(x) = Y[p(•, x)](t). ( 6
)
If the impedance operator Z is LTI and admissible, then the admittance operator Y exists, is unique, and is also LTI and admissible so that (6) reduces to u(t, x) • n(x) = [y p(•, x)](t), [START_REF] Cummings | Transient and multiple frequency sound transmission through perforated plates at high amplitude[END_REF] where the admittance kernel y is the convolution inverse of z, i.e. y z = δ or ŷ(s)ẑ(s) = 1. The admissibility conditions on Y are strictly identical to that of Z, so that Proposition 3 with "admittance" substituted for "impedance" holds true.

The third and last considered formulation of the IBC (2) is the so-called scattering formulation

p(t, x) -u(t, x) • n(x) = B[ p(•, x) + u(•, x) • n(x)](t), (8)
where p := z -1 0 p and B is known as the scattering operator. The expression (8) has an elementary physical interpretation based on the characteristics of the LEEs (1): the incident wave p + u • n is reflected back as pu • n. Admissibility conditions can also be formulated on the scattering operator, see the definition below.

Definition 6 (Admissibility conditions).

Let B be a single-input single-output operator. It is said to be an admissible scattering operator if it enjoys the following properties: causality; reality (real-valued inputs are mapped to real-valued outputs); passivity, i.e. for every smooth and compactly supported input v ∈ C ∞ 0 (R) and for every time instant t > 0,

t -∞ |B(v)(τ )| 2 dτ ≤ t -∞ |v(τ )| 2 dτ . ( 9
)
The equivalence between the passivity properties ( 5) and ( 9) follows from the identity 4 pu = ( p + u) 2 -( pu) 2 . The potential computational benefit of B over Z and Y can be guessed by comparing [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] and [START_REF] Zhang | Numerical investigation and modelling of acoustically excited flow through a circular orifice backed by a hexagonal cavity[END_REF] and is demonstrated on an elementary example in Sec. 5. When possible, the scattering operator B can be deduced from the (nonlinear) impedance operator Z through

B := (z -1 0 Z -I) • (z -1 0 Z + I) -1 = I -2(z -1 0 Z + I) -1 , ( 10
)
where I denotes the identity operator. An example of physical nonlinear scattering operator derived using this expression is given in Sec. 3.1. If Z is LTI and admissible, then B can be uniquely defined and is an admissible LTI scattering operator with kernel β, whose Laplace transform β(s), known as the reflection coefficient in acoustics, is given by

β(s) = z -1 0 ẑ(s) -1 z -1 0 ẑ(s) + 1 = 1 - 2 z -1 0 ẑ(s) + 1 ( (s) > 0). (11)
Intuitively, z = ∞, y = 0, and B = I (kernel β = δ) yield a hard wall while z = 0, y = ∞, and B = -I (kernel β = -δ) yield a pressure-release wall. Similarly to the impedance case, the admissibility of an LTI scattering operator can be conveniently characterized using its Laplace transform. The result, given in Proposition 8, is identical to that of Proposition 3 with "bounded-real" substituted for "positive-real".

Definition 7 (Bounded-real function). A function

f : C + 0 → C is bounded-real if f is analytic in C + 0 , f (s) ∈ R for s ∈ (0, ∞), and | f (s)| ≤ 1 for (s) > 0.
Proposition 8. An LTI operator B is an admissible scattering operator if and only if its Laplace transform β(s) is a bounded-real function.

Proof. See the references quoted in the proof of Proposition 3. Note that the conformal map from {s ∈ C | (s) > 0} onto {s ∈ C | |s| < 1} defined by s → (s-1) /(s+1) transforms a positive-real function into a bounded-real one. [45, Eq. 1].

Assumption on the base flow and well-posedness

In preparation for the analysis of Sec. 4, it is convenient to rewrite the LEEs (1) as a Friedrichs-symmetric system

∂ t v(t, x) + Av(t, x) = 0 (t ∈ (0, ∞), x ∈ ), ( 12
)
where v := (u , p) ∈ R n+1 is the perturbation vector, homogeneous to a velocity, and the spatial operator A is defined as Av := A(∇)v + B v with (I n denotes the n × n identity matrix)

A(n) = (u 0 • n)I n c 0 n c 0 n u 0 • n , B = ∇u 0 1 c 0 ∇u 0 • u 0 0 n γ ∇ • u 0 .
The boundary matrix A(n) is symmetric, with eigenvalues u 0 • n of multiplicity (n -1) and u 0 • n ± c 0 |n| of multiplicity 2.

(Due to the homentropicity assumption, the entropy mode is absent.) The absolute value | A(n)|, the positive part A(n) ⊕ , and the negative part A(n) are defined using the diagonal form of A(n) through

|A(n)| := P (n) | (n)| P (n) -1 , 2 A(n) ⊕ := |A(n)| + A(n), 2 A(n) := |A(n)| -A(n), where | (n)| := diag([| i |] i∈J1,n+1K
(Av, v) L 2 ( ) = 1 2 (C(u 0 )v, v) L 2 ( ) + 1 2 (A(n)v, v) L 2 (∂ ) , ( 13
)
where the symmetric amplification matrix C (u 0 ) and the L 2 scalar products are given by

C (u 0 ) = ∇u 0 + ∇ u 0 -(∇ • u 0 )I n 1 c 0 ∇u 0 • u 0 1 c 0 u 0 • ∇ u 0 (2γ -1)∇ • u 0 and (v, w) L 2 ( ) := (v, w) R n+1 dx, (v, w) L 2 (∂ ) := ∂ (v, w) R n+1 dx.
In this paper, the energy analysis is carried out with the standard acoustic energy defined as [2, § 6.2] [3, § 5.6]

v(t) 2 L 2 ( ) := (v, v) L 2 ( ) = | p(t, x)| 2 dx + |u(t, x)| 2 dx, (14)
so that the continuous energy balance, which expresses the evolution of acoustic energy in the domain , reads

1 2 d dt v(t) 2 L 2 ( ) = -(Av, v) L 2 ( ) = - 1 2 (C(u 0 )v, v) L 2 ( ) - 1 2 (A(n)v, v) L 2 (∂ ) , ( 15
)
where the right-hand side is given by

(C(u 0 )v, v) R n+1 = (∇ • u 0 ) (2γ -1) p2 -|u| 2 + 2u • ∇u 0 • u + 1 c 0 p u 0 (A(n)v, v) R n+1 = (u 0 • n) p2 + |u| 2 + 2c 0 p(u • n). ( 16
)
The expression of (A(n)v, v) R n+1 suggests the following assumption.

Assumption 10. The base flow u 0 obeys (at least) a slip condition u 0 • n = 0 at the impedance boundary z .

Thanks to this assumption, an admissible impedance yields

t 0 (A(n)v, v) L 2 ( z )
dτ ≥ 0 and there is energy dissipation at the impedance boundary z . Therefore, in the absence of other boundary sources, the presence of an instability is linked to the spectrum of the amplification matrix C (u 0 ). Intuitively, this a priori energy estimate yields uniqueness of the solution, and is the stepping stone to obtain well-posedness in e -κt C((0, ∞); H 1 ( ) n+1 ) with κ finite [START_REF] Benzoni-Gavage | Multi-dimensional Hyperbolic Partial Differential Equations: First-order Systems and Applications[END_REF]. It also follows from the energy estimate that κ ≥ min x∈ λ min (u 0 ), where λ min (u 0 ) is the minimum eigenvalue of C (u 0 ). A sufficient condition for stability is that min x∈ λ min (u 0 ) ≥ 0, which is a stringent condition on the base flow u 0 (x). This assumption is typically not verified in applications: it holds for a constant flow but fails for a Poiseuille flow for instance.

Remark 11. Although Assumption 10 is mathematically sufficient for (1), (3) to be well-posed, it is interesting to note that, physically, a no-slip condition u 0 = 0 is required at the impedance wall. Indeed, if u 0 = 0 the effect of the hydrodynamic boundary layer (refraction of sound waves [START_REF] Pridmore-Brown | Sound propagation in a fluid flowing through an attenuating duct[END_REF]) must be explicitly modeled, which leads to a non-locally reacting impedance, out of the scope of this paper. For example, the standard Ingard-Myers boundary condition reads u n = Y (p, ∇ p), and has been shown by Brambley using a Briggs-Bers analysis to prevent well-posedness, in the sense that there is no finite value for κ [START_REF] Brambley | Fundamental problems with the model of uniform flow over acoustic linings[END_REF]. A corrected impedance that accounts for both the acoustic and hydrodynamic boundary layers has been derived by Khamis & Brambley [17]. The computational interest of such non-local boundary conditions is that the boundary layer need not be discretized, see [START_REF] Gabard | A full discrete dispersion analysis of time-domain simulations of acoustic liners with flow[END_REF] for a use in the time domain.

Remark 12. If u 0 is uniform, the acoustical energy ( 14) satisfies a conservation law, namely [START_REF] Pridmore-Brown | Sound propagation in a fluid flowing through an attenuating duct[END_REF]. A wider class of base flows can be covered using the aeroacoustical energy introduced by Cantrell and Hart, which reads [55, Eq. ( 16)]

v 2 CH := p 2 L 2 ( ) + u 2 L 2 ( ) + 2( p, M 0 • u) L 2 ( ) = u + pM 0 2 dx + 1 -|M 0 | 2 | p| 2 dx,
and defines a norm if and only if M 0 L ∞ ( ) < 1 where M 0 := u 0/c 0 . Under Assumption 10, using

• CH instead of • L 2 ( ) in (15) yields the boundary term (A CH (n)v, v) R n+1 := 2c 0 p + M 0, • u (u • n),
where the subscript " " denotes the tangential component. The boundary term ( 16) is recovered if u 0 obeys a no-slip condition.

Physical reflection coefficient models in the time domain

As mentioned in the introduction, a discrete TDIBC consists of three components: a discrete time-domain impedance model, an algorithm to compute said discrete model, and a coupling method with the considered PDE. This section tackles the first two components, while the last one is covered independently in Sec. 4, so that Secs. 3 and 4 can be read independently. This section is purely focused on the impedance models and does not rely on (u, p) satisfying a particular equation in . The broad purpose of this section is to derive time-local discrete models from the analysis of physical models, for later use in Secs. 5 and 6. The presented analysis focuses on models suitable for acoustical liners, recalled in Sec. 3.1, which include an example of nonlinear scattering operator. The last two sections focus on the two linear models. Sec. 3.2 gives the oscillatory-diffusive representation of their reflection coefficients, which yield a formulation with an infinite number of delayed ODEs. Its practical discretization is then tackled in Sec. 3.3, where the time delay is recast using a transport equation, thus yielding a purely time-local formulation.

Physical models for acoustical liners

The numerical application presented in Sec. 6 involves two so-called single degree of freedom acoustical liners: a ceramic tubular (CT) liner made from a ceramic tubular core (channel of length l c , diameter d c , and porosity σ c ∈ (0, 1]) and a rigid backplate; a micro-perforated (MP) liner made from a honeycomb core (cell of length l c , diameter d c , and porosity σ c ∈ (0, 1]) sandwiched between a perforated facesheet (thickness l p , hole diameter d p , and porosity σ p ∈ (0, 1]) and a rigid backplate. CT liners are mostly used in academic and benchmark experiments, while MP liners are widely used in industrial applications. The standard linear impedance model for these liners reads

ẑphys (s) = σ -1 p ẑp (s) + σ -1 c ẑc (s) ( 17
)
where ẑp (resp. ẑc ) is the impedance of a perforation (resp. cavity). The CT liner model ẑCT is obtained with ẑp = 0, while σ c = 1 yields the MP liner model ẑMP . (See [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] and references therein for a wider panel of linear impedance models.)

The expression [START_REF] Khamis | Viscous effects on the acoustics and stability of a shear layer over an impedance wall[END_REF] separates the contribution of the perforated plate from that of the cavity, a feature that assumes the conservation of u • n across the perforation, which is verified as long as the perforation thickness l p is much shorter than the considered wavelengths. The chosen perforation model is given by a fractional polynomial

ẑp z 0 (s) = a 0 + a1 /2 √ s + a 1 s, ( 18
)
where each of the non-negative coefficients has a physical interpretation: a 0 models frequency-independent losses, a1 /2 frequency-dependent losses coming from visco-thermal effects, and a 1 is known as the mass reactance and does not incur any loss. Note that if a1 /2 = 0, ẑp is multivalued: in this work the cut is always defined as (-∞, 0] and √ • denotes the branch that coincides with the real-valued square root on (0, ∞)

a 0 = 3l p ν c 0 ( d p/2) 2 , a1 /2 = 2l p √ ν c 0 d p/2 , a 1 = l p c 0 , ( 19
)
where ν denotes the kinematic viscosity. For the cavity, the adopted model is a monodimensional wave equation with a fractional wavenumber k c (s)

ẑc z 0 (s) = coth(jk c (s)l c ), jk c (s)l c = b 0 + b1 /2 √ s + b 1 s. ( 20
)
The non-negative coefficients b have a physical interpretation similar to those of ( 18): b 0 models frequency-independent losses, b1 /2 frequency-dependent losses, and b 1 is half the back-and-forth traveling time in the cavity. The high-frequency approximation of a wavenumber derived by Bruneau gives [5, Eq. 26]

b 0 = 0, b1 /2 = √ ν c 0 d c/2 γ -1 √ Pr + 1 l c , b 1 = l c c 0 , ( 21
)
where Pr is the Prandtl number. From Proposition 3, to check the admissibility of these linear models it is sufficient to check that ẑ(s) is a positive-real function, which is straightforward.

By contrast with the CT liner, the MP liner is sensitive to the incident sound pressure level: above a given threshold, nonlinear effects occur in the perforation. The prevailing nonlinear term for a perforated plate with porosity σ p is based on the incompressible Bernoulli equation, and has been derived by Cummings [7, Eq. 1]

Z nl z 0 (u • n) = C nl c 0 |u • n|u • n, C nl = 1 -C 2 c 2σ 2 p C 2 c , ( 22
)
where the nondimensional coefficient C c ∈ (0, 1] is known as the contraction coefficient. Evidence for the relevance of this model has been provided by Singh & Rienstra [START_REF] Rienstra | Nonlinear asymptotic impedance model for a Helmholtz resonator of finite depth[END_REF]. Since the model is nonlinear, it is naturally expressed in the time domain. It is straightforward to check that this single-input single-output operator is real-valued, causal and passive (since (|u|u)u ≥ 0, the passivity condition ( 5) is verified): therefore, it defines an admissible impedance model. By contrast with linear models, expressing the admittance and scattering operators of a nonlinear impedance model is more intricate. Let us give an elementary algebraic example by studying the impedance model

Z z 0 (u • n) = a 0 u • n + C nl c 0 |u • n|u • n. (23)
A somewhat tedious but straightforward computation shows that the admittance and scattering operators are given by

Y(p) = 2 p a 0 + a 2 0 + 4 C nl c 0 | p| (24)
B( p + u • n) = β 0 2( p + u • n) 1 + 1 + 4C nl (1+a 0 ) 2 | p+u•n| c 0 + 1 c 0 C nl (1 + a 0 ) 2 4| p + u • n|( p + u • n) 1 + 1 + 4C nl (1+a 0 ) 2 | p+u•n| c 0 2 , ( 25
)
where β 0 = (a 0 -1)/(a 0 + 1) is the reflection coefficient associated with a 0 . To the best of the authors' knowledge, the nonlinear admittance and scattering operators [START_REF] Cotté | Time-domain simulations of sound propagation in a stratified atmosphere over an impedance ground[END_REF], [START_REF] Richter | Liner Impedance Modeling in the Time Domain with Flow[END_REF] have never been considered. The expression [START_REF] Richter | Liner Impedance Modeling in the Time Domain with Flow[END_REF] is used in Sec. 5 to illustrate the computational interest of nonlinear scattering operators.

Remark 13. In [START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF], a nonlinear model is used to model injectors of a combustion chamber in a LES simulation. The corre-

sponding scattering operator is B(v)(t) = χ a (v)(t -τ )
where τ > 0 is a time delay and χ a is an algebraic function of the incoming characteristic [START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF]Eq. (2.15)].

The effect of a grazing flow on the impedance of a sound absorbing material has been studied both experimentally [START_REF] Kirby | The impedance of perforated plates subjected to grazing gas flow and backed by porous media[END_REF] and numerically [START_REF] Zhang | Numerical investigation of a honeycomb liner grazed by laminar and turbulent boundary layers[END_REF] using direct numerical simulations. For a perforated plate backed by an air-filled cavity, a grazing flow is usually associated with a resistance increase that degrades the material performance. Kirby and Cummings proposed an experimental correlation between the impedance and the wall friction velocity, namely [10, Eqs. ( 12)-( 13)]. This result suggests that the grazing flow can be empirically modeled as an additional parametric dependence, which preserves the locally-reacting nature of the IBC.

Oscillatory-diffusive representation of physical reflection coefficient models

In [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF], it is shown that the time-domain convolution z phys u can be expressed without any approximation using an infinite number of ODEs, some of which may be delayed. This result is obtained by deriving the so-called oscillatory-diffusive representation of the physical model z phys using complex calculus. However, as shown in Sec. 4, it is computationally advantageous to enforce the IBC through the scattering formulation (8) rather than the impedance or admittance formulations (2), [START_REF] Melling | The acoustic impedance of perforates at medium and high sound pressure levels[END_REF]. This section shows that, fortunately, the reflection coefficient βphys does also enjoy an oscillatory-diffusive repre- sentation, given by ( 29). (As a matter of fact, this is also the case for the admittance ŷphys .) Since the computations are similar to that carried out in [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF], only the key steps are provided; the reader is referred to [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF][START_REF] Lombard | Diffusive approximation of a time-fractional Burger's equation in nonlinear acoustics[END_REF] and references therein for background on oscillatory-diffusive representations.

The starting point of the analysis is the identity coth(s) = 1 + 2e -2s (1e -2s ) -1 that enables to rewrite the reflection coefficient (( 11), ( 17)) as

βphys (s) = 1 + ĥ1 (s) + e -sτ ĥ2 (s), ( 26
)
where the time delay τ := 2b 1 > 0 is the cavity back-and-forth traveling time and the functions ĥ1 and ĥ2 , which induce deviations from the rigid wall βphys (s) = 1, are given by

ĥ1 (s) = - 2 R(s) , ĥ2 (s) = -ĥ1 (s)e -2(b 0 +b 1 /2 √ s) , R(s) = 1 + 1 σ c + 1 σ p ẑp (s) + -1 + 1 σ c - 1 σ p ẑp (s) e -2jk c (s)l c .
The interest of the apparently gratuitous expression ( 26) is that both ĥ1 and ĥ2 admit an oscillatory-diffusive representa- tion (( 27), ( 28)). The analytical expression of these two functions is used in the discretization methodology presented in Sec. 3.3.1. The representation of ĥi (i ∈ {1, 2}) is derived by inverting the Laplace transform using the residue theorem [56, Chap. VII] on a "keyhole-shaped" Bromwich contour C R, . Among the sufficient hypotheses needed to carry out the computations, ĥi must decay uniformly on {|s| = R} as R → ∞ and admit finite residues at every pole and branch points (the only possible branch point of ĥi is 0). If a1 /2 or a 1 are positive (which is the case for an MP liner), the oscillatory-diffusive representation of ĥi is

ĥi (s) = n∈I r i,n s -s n + ∞ 0 μ i (ξ ) s + ξ dξ, (27)
where (s n ) n with n ∈ I ⊂ Z is the sequence of poles of ĥi , each one solving R(s n ) = 0, r i,n is the associated residue r i,n :=

Res( ĥi , s n ), and μ i is the so-called diffusive weight (or measure), linked to the jump of ĥi across its cut (-∞, 0] by μ i (ξ ) :=

(j2π ) -1 [ ĥi (ξ e -jπ ) -ĥi (ξ e +jπ )]
. These quantities can be readily computed numerically (the poles s n are not always known analytically), which enables to verify [START_REF] Troian | Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow[END_REF]. The first (resp. second) term on the right-hand side of ( 27) is the oscillatory (resp.

diffusive) part of ĥi . The oscillatory part is associated with the resonances, by contrast with the diffusive part that stems from visco-thermal losses. Indeed, without the multivalued fractional term √ s (physically, if ν = 0, see [START_REF] Özyörük | Time-domain numerical simulation of a flow-impedance tube[END_REF]) the diffusive part vanishes, i.e. μ i = 0. The admissibility conditions given by Proposition 8 imply that μ(ξ) ∈ R and that if s n is a pole of ĥi with residue r i,n , then s n is a pole of ĥi with residue r i,n (hermitian symmetry). The desired representation of βMP is then obtained by combining [START_REF] Richter | A review of time-domain impedance modelling and applications[END_REF] and [START_REF] Troian | Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow[END_REF].

However, for the CT liner ẑp (s) = 0 so that ĥ1 (and ĥ2 if b1 /2 = 0) fails the decay condition on {|s| = R} as R → ∞. By considering instead ĥ1 (s) /s, the following representation can be derived, assuming σ c = 1,

ĥ1 (s) = C + n∈I r 1,n s -s n + ∞ 0 μ 1 (ξ ) s + ξ dξ, C - 1 1 + σ -1 c , ( 28
)
which differs from ( 27) by a real constant C . The representation of βCT is then obtained by combining ( 26), ( 27) with i = 2, and (28). Overall, for both liners, the physical reflection coefficient has the causal representation

β phys (t) = β 1 δ(t) + n∈I r 1,n e s n t + + ∞ 0 μ 1 (ξ )e -ξt + dξ + n∈I r 2,n e s n (t-τ ) + + ∞ 0 μ 2 (ξ )e -ξ(t-τ ) + dξ, (29)
where e st + := e st H(t) denotes the causal exponential function and β 1 ∈ R. The importance of (29) lies in the fact that it provides a structural information on β phys , namely that it "reduces" to an infinite number of first-order systems with delay

β phys v(t) = β 1 v(t) + n∈I r 1,n ϕ(t, s n ) + ∞ 0 μ 1 (ξ )ϕ(t, -ξ)dξ + n∈I r 2,n ϕ(t -τ , s n ) + ∞ 0 μ 2 (ξ )ϕ(t -τ , -ξ)dξ, ( 30
)
where ϕ(t, s) := (e st + v)(t) and v is a shorthand for " p + u • n". The computational interest of (30) stems from the ability to compute ϕ through the first-order ODEs ∂ t ϕ(t, s) = sϕ(t, s) + v(t) with null initial condition ϕ(t = 0, s) = 0. As a result, ( 30) is a time-local (but infinite-dimensional) representation of the hereditary convolution operator B(v) = β phys v. Theoretically, this representation enjoys an energy balance, useful for stability analyses (see Sec. 4 and Appendix C). The above analysis informs the discrete model proposed in the next section.

Discrete reflection coefficient model

The identity [START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF] naturally suggests to define the discrete reflection coefficient model β as

β(s) = β∞ + N s n=1 r1,n s -sn + N ξ k=1 μ1,k s + ξk + e -s τ ⎛ ⎝ N s n=1 r2,n s -sn + N ξ k=1 μ2,k s + ξk ⎞ ⎠ β v(t) = β∞ v(t) + N s n=1 r1,n ϕ(t, sn ) + N ξ k=1 μ1,k ϕ(t, -ξk ) + N s n=1 r2,n ϕ(t -τ , sn ) + N ξ k=1 μ2,k ϕ(t -τ , -ξk ), ( 31
)
where ϕ(t, •) follows the first-order ODE defined above. It involves N ϕ := N s + N ξ additional variables, where N ξ (resp. N s ) variables come from the diffusive (resp. oscillatory) part of β phys . Note that β∞ = β(+∞) ∈ [-1, 1] from Proposition 8. In summary, the analysis of the physical model ( 17) has led to the discrete model [START_REF] Olivetti | Direct numerical simulation of turbulent flow with an impedance condition[END_REF], which requires the computation and delay of N ϕ ODEs. Without the delay τ = 0, the derived model β can be interpreted as a so-called multipole model, postulated in [START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF] for instance [5, § I]. In preparation for the application of (31) presented in Secs. 5 and 6, Sec. 3.3.1 discusses the computation of the poles (s n , ξk ) n,k and weights (r i,n , μi,k ) n,k , while Sec. 3.3.2 covers the time delay discretization.

Poles and weights computation

The discrete model β is fully determined by the constant β∞ , N ϕ poles (s n , ξk ) n,k , 2N ϕ weights (r i,n , μi,k ) n,k , and time delay τ . These parameters should be such that β is a satisfactory representation of the considered sound absorbing material, typically known by its physical characteristics and possibly some experimental data ( βexp (j2π f m )) m . Criteria for a satisfactory approximation include: β is a bounded real function (see Proposition 8); N ϕ is as low as possible; ω → β(jω) has a physical behavior at frequencies not covered by the experimental data (usually low and high frequencies); the maximum frequency (2π ) -1 max n,k (|s n |, ξk ) n,k is consistent with the stability region of the time discretization scheme (otherwise the IBC could reduce the timestep of an explicit integration). In Sec. 6, the three-step methodology given below is followed. The principle of the approach is to use as much as possible the information obtained on the physical model in Sec. 3.2 to ease the optimization process.

1. (Physical model) Compute the coefficients a and b from the liner dimensions using the models ( 19), [START_REF] Liu | Stability analysis and design of time-domain acoustic impedance boundary conditions for lined duct with mean flow[END_REF]. In particular, this provides a value for the constant β 1 and the time delay τ . These coefficients can then serve as an initial point for a nonlinear least squares optimization on ω → ẑphys (jω) -ẑexp (jω) . This optimization can be interpreted as computing the various physical impedance corrections [START_REF] Melling | The acoustic impedance of perforates at medium and high sound pressure levels[END_REF][START_REF] Cummings | Transient and multiple frequency sound transmission through perforated plates at high amplitude[END_REF].

Inputs: material characteristics (l p , d p , σ c , etc.), experimental data (ẑ exp (j2π f m )) m .

Outputs: coefficients (σ -1

p a 0 , σ -1 p a1 /2 , σ -1 p a 1 , σ c , b 0 , b1 /2 , b 1 ), time delay τ = 2b 1 , constant β 1 = 1 -(1 + σ -1 c ) -1 . 2.
(Discrete model -Linear least squares) Choose N ξ diffusive poles along the cut (-∞, 0] (typically, a logarithmic repartition is satisfactory). Compute N s oscillatory poles sn by solving R(s n ) = 0. Compute the weights (r n,i ) n and ( μk,i 27) holds (add the constant if (28) holds). This is an overdetermined linear least squares optimization that is solved instantaneously by pseudo-inverse. See [5, § IV.D,VI.A] for additional details. Note that the time delay τ has no role whatsoever during this step. Inputs: ĥi (from output of step 1), number of poles N s , diffusive poles (ξ k ) k∈J1,N ξ K .

) k by minimizing ω → N ξ k=1 μk,i jω+ ξk + N s n=1 rn,i jω-s n -ĥi (jω) if (
Outputs: oscillatory poles

(s n ) n∈J1,N s K , weights (r n,i ) n∈J1,N s K and ( μk,i ) k∈J1,N ξ K (i ∈ {1, 2}).
3. (Discrete model -Nonlinear least squares) Compute new poles and weights with a nonlinear least squares optimization on ω → β(jω) -βexp (jω) with initial poles and weights given by step 2, β∞ = β 1 and τ = τ from step 1. The constant β∞ and delay τ can also be optimized along the poles and weights.

Inputs: (s n , ξk ) n,k , (r n,i , μk,i ) n,k , τ = τ , β∞ = β 1 , ( βexp (j2π f m )) m . Outputs: (s n , ξk ) n,k , (r n,i , μk,i ) n,k , (optional) τ , (optional) β∞ .
When the physical model is satisfactory, as is the case for the MP and CT liners considered herein, the first two steps deliver a discrete model β that may be sufficient for some engineering applications. If step 3 is used, its role is to adjust the poles and weights to improve the fit against experimental data. The first two steps can then be interpreted as using a physical model to help finding an initial guess for the poles and weights. When faced with a mismatch between the physical model and experimental data, tuning the former may be required, but this is not a computationally intensive task since there are few parameters and they have a physical interpretation. The third stage can be performed using any nonlinear least squares method: in Sec. 6, we simply relied on the trust region optimization method [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF] implemented in MATLAB ® lsqnonlin, whose execution takes a few seconds on a contemporary computer, but the more tailored vector fitting algorithm [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF] can also be used, as in [START_REF] Troian | Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow[END_REF]. However, note that step 3 comes with the following caveat. There is a trade-off between the fit quality against the experimental data, which is usually narrowband, and the broadband behavior of β. This trade-off is especially acute when using impedance identified with base flow, as these can be associated with significant uncertainties that ideally should be accounted for in the optimization process. This is exemplified by the question of whether to optimize the time delay τ . Physically, it is linked to the anti-resonant frequencies f n given by f n = n /τ for (17) with a lossless cavity b 0 = b1 /2 = 0 (i.e. a canonical Helmholtz resonator). Therefore, although τ = τ may enable a better fit, it may not be worth altering the broadband behavior. Eventually, knowledge about both the model and experimental data is helpful to inform the optimization process.

Time delay discretization

The discretization of the time delay is done independently of the choice of weights and poles described in Sec. 3.3.1. Theoretically, the delay can be recast into ϕ(t -τ , s) = ψ(t, s, -l τ ) where the additional function ψ(•, s, •) obeys the following transport equation on (-l τ , 0) with

l τ = c τ τ > 0 ∂ t ψ(t, s, θ) = c τ ∂ θ ψ(t, s, θ) (t ∈ (0, ∞), θ ∈ (-l τ , 0)) ψ(t, s, 0) = ϕ(t, s), ψ(0, s, θ) = 0. ( 32
)
This device is commonly used in theoretical and numerical studies of delay differential equations, see the references in [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF] and [62, § 1]. Although not needed herein, several delays τi can also be tackled by defining τ := max i τi . The monodimen- sional PDE (32) is discretized using a discontinuous Galerkin method with N K elements and N p ≥ 2 nodes per element (order N p ). The discretization accuracy is measured by the number of points per wavelength

PPW τ ( f max ) := N ψ τ f max , ( 33
)
where N ψ := N K N p is the total number of nodes and f max is the maximum frequency of interest. Target values for PPW τ as a function of N p can be found in [START_REF] Hu | An analysis of the discontinuous Galerkin method for wave propagation problems[END_REF]Tab. 1], which shows that it is advantageous to choose N K = 1 so that N p is the sole discretization parameter. The impact of PPW τ is illustrated in Sec. 6. Herein, the main argument in favor of delaying through a transport equation, compared to continuous Runge-Kutta methods [START_REF] Zennaro | Natural continuous extensions of Runge-Kutta methods[END_REF], is the meaningfulness of (33) for wave propagation problems, which makes the discretization straightforward to set up a priori. In summary, the discrete impedance model ( 31) is computed through

β v(t) = β∞ v(t) + N s n=1 r1,n ϕ(t, sn ) + N ξ k=1 μ1,k ϕ(t, -ξk ) + N s n=1 r2,n ψ(t, sn , -l τ ) + N ξ k=1 μ2,k ψ(t, -ξk , -l τ ), ( 34
)
which implies the use of (N ψ + 1)N ϕ additional variables. The two parameters to control are max k,n (ξ k , |s n |), which is to be chosen in accordance with the spatial discretization scheme, and PPW τ ( f max ), for which values are given in [START_REF] Hu | An analysis of the discontinuous Galerkin method for wave propagation problems[END_REF]Tab. 1].

This covers the first two components of the TDIBC. The next section deals with the third, necessary for the numerical applications of Secs. 5 and 6.

Remark 14. The discrete model [START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF] has been derived from a mathematical analysis that emphasizes the distinct components of the model ( 26): an oscillatory-diffusive part that models both oscillatory and diffusion phenomena (the latter of parabolic nature), and a hyperbolic part associated with the time delay that models a wave reflection phenomenon. In [32, § 4.1.1], a first-order finite-difference discretization of the monodimensional LEEs is used to truncate a combustion chamber. An alternative discretization strategy is to perform a rational approximation of the time delay, see [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF] and references therein. For a pure time delay, i.e. β(s) ∝ e -τ s , a multipole approximation method is presented in [START_REF] Douasbin | Delayed-time domain impedance boundary conditions (d-tdibc)[END_REF].

Numerical fluxes for nonlinear time-domain impedance boundary conditions

This section deals with the third component of a TDIBC, namely its (semi-)discrete formulation, independently of Sec. 3. It analyzes the weak enforcement of an admissible IBC within a discontinuous Galerkin (DG) finite element method, but the numerical flux formalism is employed to enable a straightforward transition to other methods popular in fluid mechanics. For the analysis, the IBC is only assumed admissible and need not be given by one of the models covered in Sec. 3, so that both sections are independent. The analysis shows the computational interest of a numerical flux based on the scattering operator B [START_REF] Kirby | The impedance of perforated plates subjected to grazing gas flow and backed by porous media[END_REF], the so-called B-flux ((39), ( 63)), over fluxes based on the impedance and admittance, namely the Z-flux ((39), ( 59)) and the Y-flux ((39), [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF]). This justifies a posteriori the interest in B displayed in Sec. 3 and prepare the numerical applications of Secs. 5 and 6. The section starts with a brief introduction of the standard DG method in Sec. 4.1. Sec. 4.2 then defines admissibility conditions for an impedance numerical flux, upon which the analyses of Secs. 4.3 and 4.4 rely.

Discontinuous Galerkin method

Let (T h ) h be a quasi-uniform sequence of meshes indexed by h := max T ∈T h h T , where h T denotes the diameter of the element T ∈ T h . For simplicity, each mesh T h is assumed to be simplicial, geometrically conformal, and shape-regular. Let P k n Fig. 1. Notations for an interior face F ∈ F i h .

be the space of polynomials of n variables and total degree at most k. The approximation space is taken as V h := P k n (T h ) n+1 , where P k n (T h ) is the broken polynomial space defined as

P k n (T h ) := {v ∈ L 2 ( ) | ∀T ∈ T h , v |T ∈ P k n (T )}.
The domain is assumed to be a polyhedron, so that it can be exactly covered by each mesh. See [65, Chap. 1] and [66, Chap. 1] for definitions of the mentioned properties (more detailed citations are available in Appendix B). The semi-discrete formulation of (12) reads: find

v h ∈ C 1 ([0, ∞), V h ) such that ∂ t v h + A h v h = 0, ( 35
)
where the spatial discretization is embodied by the finite-dimensional operator A h :

V h → V h defined by ∀w h ∈ V h , (A h v h , w h ) L 2 ( ) := T ∈T h (Av h , w h ) L 2 (T ) + ((A(n)v h ) * -A(n)v h , w h ) L 2 (∂ T ) ,
where the quantity (A(n)v h ) * is the so-called numerical flux function, uniquely defined at each face. Intuitively, this definition can be viewed as resulting from two integration by parts. At an interior face F ∈ F i h we use the upwind flux

(A(n F )v) * := A(n F ) ⊕ v |T 1 -A(n F ) v |T 2 = A(n F ){v} + 1 2 |A(n F )|JvK where {v} := (v |T 1 + v |T 2 )/2 is the face average and JvK := v |T 1 -v |T 2 is
the face jump, see Fig. 1 for the face-specific notations. This leads to

(A h v h , w h ) L 2 ( ) = T ∈T h (Av h , w h ) L 2 (T ) + ((A(n)v h ) * -A(n)v h , w h ) L 2 (∂ ) - F ∈F i h (A(n F )Jv h K, {w h }) L 2 (F ) + 1 2 F ∈F i h (|A(n F )|Jv h K, Jw h K) L 2 (F ) . ( 36
)
At a boundary face F ∈ F b h , the IBC is weakly enforced through a numerical flux introduced in Sec. 4.2 and analyzed in Secs. 4.3 and 4.4. In the applications of Secs. 5 and 6, the numerical fluxes (A(n

F )v) * s := A(n F ) ⊕ v |T 1 -A(n F ) v s (source v s ) and (A(n F )v) * out := A(n F ) ⊕ v |T 1 (
approximate non-reflecting outlet) are also used. The bilinear form [START_REF] Grabowski | Stabilization of wave equation using standard/fractional derivative in boundary damping[END_REF] is standard and fits within the framework proposed by Ern & Guermond [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] for the DG discretization of Friedrichs systems, and notations have been kept as close as possible to that employed in this seminal work. The reader interested in an error analysis with the standard boundary condition " A(n)v = M(n)v" is referred to [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] and [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]Chap. 3 & 7]. The well-posedness of the semidiscrete formulation [START_REF] Abbas | The multidimensional wave equation with generalized acoustic boundary conditions I: Strong stability[END_REF] follows from the Cauchy-Lipschitz-Picard theorem [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Thm. 7.3] since it is a finite-dimensional ODE.

Using the identity (13) on each element T enables to get

(A h v h , v h ) L 2 ( ) = 1 2 (C(u 0 )v h , v h ) L 2 ( ) + 1 2 |v h | 2 upw + (A(n)v h ) * - 1 2 A(n)v h , v h L 2 (∂ ) , ( 37
)
where the upwind seminorm is defined as

|v h | 2 upw := F ∈F i h (|A(n F )|Jv h K, Jv h K) L 2 (F ) .
This leads to the following semi-discrete energy balance

1 2 d dt v h 2 L 2 ( ) = -(A h v h , v h ) L 2 ( ) = - 1 2 (C(u 0 )v h , v h ) L 2 ( ) - 1 2 |v h | 2 upw -(A(n)v h ) * - 1 2 (A(n)v h ), v h L 2 (∂ ) , ( 38
)
which is to be contrasted with its continuous counterpart [START_REF] Pridmore-Brown | Sound propagation in a fluid flowing through an attenuating duct[END_REF]. The term |v h | 2 upw is the energy dissipation due to the use of an upwind flux (it would be null had a centered flux been used) while the boundary term

((A(n)v h ) * -1 2 A(n)v h , v h ) L 2 (∂ )
includes the contribution of the weakly enforced IBC, studied below.

Admissibility conditions for an impedance numerical flux

On a part of the boundary := ∂ denoted z , the IBC is weakly enforced through a centered flux

(A(n)v) * := 1 2 A(n)(v + v g ), ( 39
)
where v g is the so-called ghost state that needs to be suitably defined as a function of v and Z, Y, or B. A systematic derivation of the possible expressions for A(n)v g is carried out in Sec. 4.4 with an energy analysis formalized through the so-called admissibility and continuity conditions defined below. Recall that the spatial discretization is said to be consistent

if for v(t) ∈ V the exact solution, (A h v, w h ) L 2 ( ) = (Av, w h ) L 2 ( ) for all w h ∈ V h .
Remark 15. The ghost state v g is a function known in closed-form and in this paper the numerical flux ( 39) is enforced without adding computational nodes outside of : the terminology "ghost state" is chosen for its intuitive nature, and does not imply that "ghost nodes" or "ghost cells" are used whatsoever.

Definition 16 (Admissibility). The impedance numerical flux [START_REF] Chun | High-order accurate thin layer approximations for time-domain electromagnetics. Part I: General metal backed coatings[END_REF], uniquely determined by the expression of A(n)v g , is said to be admissible if it is both consistent and passive.

• (Consistency) Let v(t) ∈ V be the exact solution (in particular, it does obey the IBC). The consistency condition reads

(A(n)v) * = A(n)v or equivalently A(n)v g = A(n)v. ( 40
)
• (Passivity) Let v h (t) ∈ V h (in particular, it does not obey the IBC). The passivity condition reads

∀t > 0, t 0 (A(n)v h ) * - 1 2 A(n)v h , v h L 2 ( z ) dτ = 1 2 t 0 (A(n)v g h , v h ) L 2 ( z ) dτ ≥ 0. (41)
The passivity condition is to be understood in light of the semi-discrete energy balance [START_REF] Yuferev | Surface Impedance Boundary Conditions: A Comprehensive Approach[END_REF], see Example 24 for a numerical illustration of its necessity. In addition to these two admissibility conditions, the two following continuity properties are also computationally desirable.

• (Hard-wall continuity) As "Z → ∞" (or "Y → 0", or "B → I"), a hard wall (u • n = 0) is recovered without singularity, which can be written formally as

lim Z→∞ A(n)v g = M hw (n)v, M hw (n) := c 0 2ζ n n n -n 0 ( 42
)
and ζ is an arbitrary non-negative parameter. • (Pressure-release continuity) As "Z → 0" (or "Y → ∞", or "B → -I"), a pressure-release boundary ( p = 0) is recovered without singularity, i.e.

lim Z→0 A(n)v g = M pr (n)v, M pr (n) := c 0 0 n,n -n n 2ζ ( 43
)
and ζ is an arbitrary non-negative parameter.

It is straightforward to verify that A(n)v g := M hw/pr (n)v respects the consistency [START_REF] Hiptmair | Fast convolution quadrature based impedance boundary conditions[END_REF] and passivity [START_REF] Brockett | Finite Dimensional Linear Systems[END_REF] . The ability to recover both the hard-wall and pressure-release cases is of particular interest when performing an inverse method on the IBC where the parameter space needs to be explored.

Remark 17. Note that at this stage no assumption is made regarding the time-domain discretization of the IBC: Definition 16 is purely semi-discrete.

The application of the energy analysis to generic nonlinear IBCs, which leads to the derivation of the Z, Y, and B fluxes, is done in Sec. 4.4. For the sake of clarity, the elementary case of a proportional impedance Z(u • n) ∝ u • n is first fully worked out in Sec. 4.3. In spite of its simplicity, this example provides an intuitive understanding of the computational advantage of the scattering operator B.

Weak enforcement of proportional impedance boundary conditions

This section focuses on the computational properties of numerical fluxes for the so-called proportional impedance z -1 0 z(t) = a 0 δ(t) with a 0 > 0 ( p = a 0 u • n). (From Proposition 3 it is admissible.) To weakly enforce the IBC (2) using the numerical flux [START_REF] Chun | High-order accurate thin layer approximations for time-domain electromagnetics. Part I: General metal backed coatings[END_REF], the ghost state v g is sought as linearly dependent upon v v g :=

α 3 n ⊗ n α 4 n α 1 n α 2 v, A(n)v g = c 0 pg n u g • n = M 0 (n)v, ( 44
)
where

M 0 (n) = c 0 α 1 n ⊗ n α 2 n α 3 n α 4 . ( 45
)
To obtain an admissible flux from this generic expression, two of the four degrees of freedom (DoF) must be removed, as summarized in the proposition below.

Proposition 18.

The numerical flux function given by ( 39), ( 44) is admissible if and only if

α 1 = (1 -α 2 )a 0 and α 4 = (1 -α 3 )/a 0 , with α 2 ≤ 1, α 3 ≤ 1, and (α 2 + α 3 ) 2 ≤ 4(1 -α 2 )(1 -α 3 ).
Proof. The proof is elementary and consists in using the admissibility conditions. (Consistency) Let v = u , p be the exact solution. The consistency condition (40) reads

α 1 n ⊗ n (α 2 -1)n (α 3 -1)n α 4 v = 0.
The fact that v obeys the IBC p = a 0 u • n readily leads to α 1 = (1α 2 )a 0 and α 4 = (1α 3 )/a 0 , so that the two DoF left are α 2 and α 3 . (Passivity

) Let v h ∈ V h . Since (M 0 (n)v h , v h ) R n+1 = ( M0 vh , vh ) R 2 with vh = [ u h • n ph ] and M0 := α 1 α 2 +α 3 2 α 2 +α 3 2 α 4 , ( 46
)
the condition (41) holds if and only if M0 is positive semidefinite. The inequalities on α 2 and α 3 follow from the application of the following elementary lemma that is frequently used in the remainder of this paper. Remark 20. In the case of a proportional impedance, the consistency (40) and passivity [START_REF] Brockett | Finite Dimensional Linear Systems[END_REF] conditions are the (DG1) and (DG2) properties stated in [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] and needed for the error analysis.

Therefore, after examination of the admissibility conditions there are two remaining DoF, α 2 and α 3 , so that

v g := α 3 n ⊗ n 1-α 3 a 0 n (1 -α 2 )a 0 n α 2 v, M 0 (n) = c 0 (1 -α 2 )a 0 n ⊗ n α 2 n α 3 n 1-α 3 a 0 . ( 47
)
Further constraints can be obtained by considering the continuity conditions ( 42) and ( 43): lim

a 0 →0 α 2 = -1, lim a 0 →0 α 3 = 1, lim a 0 →∞ α 2 = 1, lim a 0 →∞ α 3 = -1, (48)
lim

a 0 →0 (1 -α 2 )a 0 = 0, lim a 0 →0 (1 -α 3 ) 1 a 0 ≥ 0, lim a 0 →∞ (1 -α 2 )a 0 ≥ 0, lim a 0 →∞ (1 -α 3 ) 1 a 0 = 0. ( 49
)
The line [START_REF] Delorme | Computational aeroacoustics applications based on a discontinuous Galerkin method[END_REF] suggests defining α := α 2 with α 3 = -α, which leads to the α-flux

v g := -αn ⊗ n (1+α) a 0 n (1 -α)a 0 n α v, M 0 (n) = c 0 (1 -α)a 0 n ⊗ n αn -αn 1+α a 0 . ( 50
)
Any flux of this form is admissible, as long as α ∈ [-1, 1]. To respect the remaining continuity conditions [START_REF] Ventribout | Contrôle des perturbations aéroacoustiques par impédances de parois: application à un modèle de matériaux poreux[END_REF], one can choose α = β 0 := (a 0 -1)/(a 0 + 1), which yields the β 0 -flux

v g := -β 0 n ⊗ n (1 -β 0 )n (1 + β 0 )n β 0 v, M 0 (n) = c 0 (1 + β 0 )n ⊗ n β 0 n -β 0 n 1 -β 0 , (51)
an apparent computational interest of which is the boundedness of its components with respect to a 0 . In summary, application of the admissibility and continuity conditions leads to [START_REF] Zemanian | Distribution Theory and Transform Analysis[END_REF], an expression that, to the best of the authors' knowledge, was first proposed by Ventribout [49, § 1.3.2] with a view on application to optimal control. Further insights into the benefit of choosing α = β 0 can be obtained by deriving a CFL stability condition.

For simplicity, let us consider the explicit Euler scheme

v n+1 h -v n h + t Ãh v n h + tA {0} h v n h = 0, (52)
with constant time step t and CFL number defined as CFL := tc 0 h -1 . To highlight the IBC contribution, the decomposition

A h v h = Ãh v h +A {0} h v h is used, where A {0} h : V h → V h
is the boundary term of A h given by, assuming that the IBC is applied on the whole of (i.e. z = ∂ ), [START_REF] Benzoni-Gavage | Multi-dimensional Hyperbolic Partial Differential Equations: First-order Systems and Applications[END_REF] so that

(A {0} h v h , w h ) L 2 ( ) = 1 2 (M 0 (n)v h -A(n)v h , w h ) L 2 (∂ ) ,
( Ãh v h , w h ) L 2 ( ) := T ∈T h (Av h , w h ) L 2 (T ) - F ∈F i h (A(n F )Jv h K, {w h }) L 2 (F ) + 1 2 F ∈F i h (|A(n F )|Jv h K, Jw h K) L 2 (F ) . ( 54
)
Using a discrete energy method (see Appendix C), the following sufficient stability condition can be derived. Proposition 21. Assume ∇u 0 = 0 (uniform base flow) and V h = P 0 n (T h ) n+1 (finite volume discretization). A sufficient L 2 stability condition for [START_REF] Lozano | Dissipative Systems Analysis and Control: Theory and Applications[END_REF]

is CFL ≤ 2 -1 C -2 DG 1 + |u 0 | c 0 -1 and α 1 + α 4 ≥ C 2 tr CFL(α 2 1 + (α 3 -1) 2 + α 2 4 + (α 2 -1) 2 ) (α 2 + α 3 ) 2 ≤ 4 α 1 -C 2 tr CFL(α 2 1 + (α 3 -1) 2 ) α 4 -C 2 tr CFL(α 2 4 + (α 2 -1) 2 ) , ( 55
)
where the positive constants C DG and C tr , defined in Lemmas 31 and 33, are non-dimensional and do not depend upon the initial data or the impedance.

Remark 22.

Recall that, to be admissible, the matrix M 0 (n) given by (45) must obey Proposition 18. In particular, the passivity condition (41) requires the conditions α 1 + α 4 ≥ 0 and (α 2 + α 3 ) 2 ≤ 4α 1 α 4 that are less stringent than [START_REF] Cantrell | Interaction between sound and flow in acoustic cavities: mass, momentum, and energy considerations[END_REF].

The interest of Proposition 21 lies in the condition (55) that gives the influence of the numerical flux on the CFL number.

Let us highlight three particular cases of practical interest by considering the α-flux ((39), ( 50)) with α ∈ [-1, 1], which has been derived above.

• (1-flux) If α = 1 (α 1 = 0), then pg = p and the flux does not control u h • n at the impedance boundary since

(M 0 (n)v h , v h ) R n+1 = 2c 0 a -1 0 p2 h . A stability condition, namely CFL ≤ C -2 tr α 4 /(1 + α 2 4
), can be obtained from (55) if and only if α 2 = -1 and α 3 = 1. But since α 2 = α = 1, the scheme is not provably stable with the proposed energy analysis.

• (-1-flux) If α = -1 (α 4 = 0), then u g = (u • n)n and the numerical flux does not control ph since (M 0 (n)v h , v h ) R n+1 = 2c 0 a 0 (u h • n) 2 .
Similarly, a stability condition cannot be obtained from Proposition 21 since that would require α 2 = 1 and α 3 = -1, values at odds with the definition of the α-flux with α = -1.

• (β 0 -flux) The α-flux satisfies (M 0 (n)v h , v h ) R n+1 = c 0 a 0 (1 -α)(u h • n) 2 + c 0 a -1 0 (1 + α) p2
h and thus control both u h • n and ph if α ∈ (-1, 1). Under this assumption, the stability condition is given by [START_REF] Cantrell | Interaction between sound and flow in acoustic cavities: mass, momentum, and energy considerations[END_REF] with

α 1 = (1 -α)a 0 , α 4 = (1 + α)a -1 0 ,
and α 2 = α = -α 3 . The β 0 -flux ((39), ( 51)

) yields CFL ≤ 1 C 2 tr min 1 1 + β 2 0 , 1 2(1 + β 0 ) , 1 2(1 -β 0 ) , ( 56
)
which shows that the CFL number has a positive upper bound with respect to a 0 , namely CFL ≤ 4 In conclusion, among the admissible α-fluxes (( 39), ( 50)) with α ∈ [-1, 1], the β 0 -flux (( 39), ( 51)) is the optimal choice since its components are bounded with respect to a 0 ∈ [0, ∞], it controls both ph and u h • n, and delivers the CFL stability condition [START_REF] Gamelin | Complex Analysis[END_REF]. Extension of this result to nonlinear scattering operators B is done next in Sec. 4.4.

Weak enforcement of nonlinear impedance boundary conditions

Let us now consider the generic, possibly nonlinear, IBC under the three forms defined in Sec. 2.2

p = a 0 u • n + a Q Q(u • n) (a) u • n = a 0 p + a Q Q( p) (b) p -u • n = B( p + u • n) (c), ( 57
)
where a 0 > 0, a Q ≥ 0, Q is an admissible impedance operator (see Definition 1), and B is an admissible scattering operator (see Definition 6). Following the method of Sec. 4.3, Secs. 4.4.1 and 4.4.2 derive the admissible numerical flux for each form ((57a), (57b), (57c)), namely the Z-flux (( 39), ( 59)), Y-flux (( 39), ( 61)), and the B-flux (( 39), ( 63)), and establish the superiority of the B-flux.

Remark 23. For the analysis presented below, the IBCs ( 57) are only assumed to be admissible. In particular, they need not be given by one of the models covered in Sec. 3, so that Secs. 3 and 4.4 are independent.

Numerical flux based on the impedance or the admittance

To weakly enforce the IBC (57a) (impedance operator z -1

0 Z(u • n) = a 0 u • n + a Q Q(u • n))
, the ghost state is sought as

A(n)v g = M 0 (n)v + m Q (n)Q(u • n), ( 58
)
where M 0 (n) is given by [START_REF] Sauter | Convolution quadrature for the wave equation with impedance boundary conditions[END_REF] and

m Q (n) = c 0 [γ 1 n , γ 2 ]
. Proposition 25 shows that the only admissible expression is

A(n)v g = c 0 2a 0 n n -n n 0 v + c 0 2a Q n 0 Q(u • n), ( 59
)
where there are no remaining DoF. The Z-flux (( 39), ( 59)) fulfills the pressure-release continuity condition [START_REF] Dragna | A generalized recursive convolution method for time-domain propagation in porous media[END_REF] but not the hard-wall one [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF] since there is a singularity as a i → ∞. Moreover, note that 2 so that there is no control of p2 h at the impedance boundary, a phenomenon already encountered in Sec. 4.3.

(M 0 (n)v h , v h ) R n+1 = 2c 0 a 0 (u h • n)
Example 24. One may think that the result of Proposition 25 is unnecessarily stringent, i.e. that the admissibility conditions proposed in Definition 16 are too constraining, but it can be verified numerically that it is not so. Let us consider the acoustical cavity

= (0, 1) 2 with c 0 = 1, z 0 = 1, and impedance Z(u • n) = a 0 u • n + a Q u • n applied to the whole of ∂ .
Two dispersion relations 1D (k i , ẑ(jω), ω) can be derived for eigenfunctions with separated variables p(x) = p1 (x 1 ) p2 (x 2 )

with pi (x i ) = (A i e jk i x i + B i e -jk i x i ) (note that it does not provide all the eigenvalues if ẑ(s) = +∞). Fig. 2 plots the eigenvalues computed with a sixth-order DG method (see Appendix D) and flux

A(n)v g = c 0 (1 -α)a 0 n ⊗ n αn -αn 1+α a 0 v + c 0 (1 -α)a Q n -(1 + α) a Q a 0 Q(u • n), ( 60
)
with reduces to the α-flux (( 39), ( 50)) if a Q = 0. The value α = -1 yields stable eigenvalues that match the exact ones until a cut-off frequency. However, for α = β 0 , eigenfunctions have a less pronounced decay and can even be unstable. Hence the impedance model Z is passive, but passivity is lost at the semi-discrete level. This conforms with the fact that the flux (( 39), ( 60))

is consistent for any α but passive if and only if α = -1 from Proposition 25.

Proposition 25.

The numerical flux function (( 39), ( 58)) is admissible if and only if the ghost state is given by [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF].

Proof. Consistency and passivity conditions are checked. (Consistency) The condition [START_REF] Hiptmair | Fast convolution quadrature based impedance boundary conditions[END_REF] yields

α 1 = (1 -α 2 )a 0 , α 4 = (1 -α 3 )/a 0 , γ 1 = (1 -α 2 )a 1 , and γ 2 = -(1 -α 3 ) a Q a 0 . (Passivity) Let v h ∈ V h . We have (A(n)v g h , v h ) R n+1 = (M 0 (n)v h , v h ) R n+1 + c 0 (1 -α 2 )a Q Q(u h • n)u h • n -c 0 (1 -α 3 ) a Q a 0 Q(u h • n) ph .
The passivity condition (41) holds true for any v h and Q if and only if

(M 0 (n)v h , v h ) R n+1 ≥ 0, (1 -α 2 )a Q ≥ 0, and α 3 = 1.
Lemma 19 gives α 2 = -1 and enables to conclude. 2

The admittance case (57b) is identical, as an application of the admissibility conditions shows that the only admissible ghost state of the form

A(n)v g = M 0 (n)v + m Q (n)Q( p) is given by A(n)v g = c 0 0 n,n n -n 2a 0 v + c 0 0 n 2a Q Q( p), (61)
which obeys the hard-wall continuity condition [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF] but not the pressure-release one [START_REF] Dragna | A generalized recursive convolution method for time-domain propagation in porous media[END_REF]. Note that

(M 0 (n)v h , v h ) R n+1 = 2c 0 a 0 p2
h so that there is no control of (u h • n) 2 at the impedance boundary. In view of the expressions ( 59) and ( 61), the fluxes based on impedance and admittance have symmetrical continuity properties such that the Z (resp. Y)-flux should be preferred to the Y (resp. Z)-flux when the impedance is close to a pressure-release wall (resp. hard wall). However, the next section shows these two fluxes are trumped by a flux based on the scattering operator B.

Numerical flux based on the scattering operator

To weakly enforce the IBC under its scattering formulation (57c), the ghost state is sought as

A(n)v g = M 0 (n)v + m B (n)B( p + u • n), ( 62
)
where M 0 (n) is given by ( 45) and m B (n

) = c 0 [γ 1 n , γ 2 ]
. Proposition 27 below shows that the only admissible ghost state is given by

A(n)v g = c 0 n n 0 n 0 n 1 v + c 0 n -1 B( p + u • n). ( 63
)
Note that with this flux there is control of both u h • n and ph at the boundary, and the boundary dissipation term is

1 2 t 0 (A(n)v g h , v h ) L 2 ( z ) dτ = c 0 4 t 0 ph + u h • n 2 L 2 ( z ) -B( ph + u h • n) 2 L 2 ( z ) dτ + c 0 4 t 0 ph -u h • n -B( ph + u h • n) 2 L 2 ( z ) dτ , ( 64
)
where the two terms on the right-hand side have a clear interpretation. The first one is the dissipation associated with the scattering operator B, see the passivity condition [START_REF] Zhang | Numerical investigation and modelling of acoustically excited flow through a circular orifice backed by a hexagonal cavity[END_REF], and the second one can be interpreted as a penalization term for the non-respect of the IBC.

Remark 26. The B-flux (( 39), ( 63)) derived from the energy analysis carried out in the proof of Proposition 27 reads

(A(n)v) * = c 0 2 n 1 ( p + u • n) + c 0 2 -1 B( p + u • n).
It is interesting to note that this flux is exactly that derived using flux vector splitting [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]Chap. 8], which for a linear hyperbolic PDE with constant coefficients amounts to exactly solving the monodimensional Riemann problem at the boundary. Indeed, given Assumption 10 we have at the impedance boundary z

A(n)v = A(n) ⊕ v -A(n) v = c 0 2 n n n n 1 v + c 0 2 -n n n n -1 v = c 0 2 n 1 ( p + u • n) + c 0 2 n -1 ( p -u • n).
Therefore, this flux has a direct interpretation based on characteristics, making it natural to use with a hyperbolic law. Proposition 27. The numerical flux (( 39), ( 62)) is admissible if and only if the ghost state is given by [START_REF] Hu | An analysis of the discontinuous Galerkin method for wave propagation problems[END_REF].

Proof. (Consistency) Condition ( 40) leads to α 1 = γ 1 = 1α 2 and α 4 = -γ 2 = 1α 3 so that there are only two DoF left, α 2 and α 3 . (Passivity) Let v h ∈ V h . Recall that the passivity property of the scattering operator B reads [START_REF] Zhang | Numerical investigation and modelling of acoustically excited flow through a circular orifice backed by a hexagonal cavity[END_REF], which makes the study of passivity more intricate than in the previous sections. The generic expressions of M 0 (n) and m B (n

) give 1 c 0 (A(n)v g h , v h ) R n+1 = 1 c 0 (M 0 (n)v h , v h ) R n+1 + 1 c 0 (m B (n)B( ph + u h • n), v h ) R n+1 = (1 -α 2 )(u h • n) 2 + (1 -α 3 ) p2 h + (α 2 + α 3 ) ph u h • n -B( ph + u h • n) (1 -α 3 ) ph -(1 -α 2 )u h • n . Using the identity ab = 1 2 a 2 + 1 2 b 2 -1 2 (a -b) 2 on B( ph + u h • n) (1 -α 3 ) ph -(1 -α 2 )u h • n and collecting the terms yields 1 c 0 (A(n)v g h , v h ) R n+1 = 1 2 ( ph + u h • n) 2 - 1 2 B( ph + u h • n) 2 + 1 2 (1 -α 3 ) ph -(1 -α 2 )u h • n -B( ph + u h • n) 2 - 1 2 (α 2 u h • n -α 3 ph ) 2 .
Since this quantity must be non-negative for every admissible scattering operator B and every v h ∈ V h , the passivity condition ( 41) is achieved if and only if α 2 = α 3 = 0. 2

In summary, even for a generic nonlinear scattering operator B, the B-flux (( 39), ( 63)) keeps the properties mentioned when studying the proportional case, such as the control of both u h • n and ph at the boundary. Note that this conclusion has been reached without considering how B is computed, i.e. the analysis is so far independent of Sec. 3. However, to conclude this section, it is insightful to derive a CFL stability condition, as in Sec. 4.3. For this, let us consider the following reflection coefficient, inspired by the representation of physical models derived in Sec. 3.2,

β(s) = ξ max ξ min 1 s + ξ dμ(ξ ), ( 65
)
where μ is a measure on (ξ min , ξ max ) with 0 ≤ ξ min < ξ max < ∞. A sufficient condition for β to be bounded-real is

μ L 1 μ ξ 2 L 1 ≤ 1, see Appendix C.
The corresponding semi-discrete formulation results from a coupling between the additional variables ϕ h and the acoustic field v (space and time variables are explicitly stated to avoid ambiguity)

∂ t ϕ h (t, x, -ξ) = -ξ ϕ h (t, x, -ξ) + ph (t, x) + u h (t, x) • n(x) (x ∈ ∂ , ξ ∈ [ξ min , ξ max ]) ∂ t v h (t) = -Ãh v h (t) -A {0} h v h (t) -A {β} h β ( ph + u h • n)(t). ( 66
)
The spatial discretization is given by ( 53), ( 54) with an additional boundary contribution A {β} h in the decomposition of A h compared to that of Sec. 4.3 given by, assuming that z = ∂ , (A

{β} h v h , w h ) L 2 ( ) := 1 2 (m B (n)v h , w h ) L 2 (∂ ) = - c 0 2 (v h , (w p h -w u h • n)) L 2 (∂ ) . ( 67
)
The result, stated in Proposition 28, is a natural extension of Proposition 21 and more specifically of [START_REF] Gamelin | Complex Analysis[END_REF]. It shows that, with an explicit time integration, the maximum frequency (2π ) -1 ξ max of the scattering operator stiffens the system; in practice, ξ max must be chosen consistently with the dissipation and dispersion properties of the spatial discretization, as advocated in Sec. 3.3.1. Numerical applications based on the scattering operator are considered in Secs. 5 and 6.

Proposition 28. Assume ∇u 0 = 0 (uniform base flow), V h = P 0 n (T h ) n+1 (finite volume discretization), and β bounded-real such that

μ L 1 μ ξ 2 L 1 ≤ 1. A sufficient L 2 stability condition for (66) is CFL ≤ min 1 2C 2 DG 1 + |u 0 | c 0 -1 , 1 2C 2 tr , c 0 h ξ max , ( 68
)
where the positive constants C DG and C tr , defined in Lemmas 32 and 33, are non-dimensional and do not depend upon the initial data or the impedance.

Numerical validation with nonlinear impedance tube

This section validates the proposed Z-flux ((39), ( 59)) and B-flux ((39), ( 63)), and illustrates the computational advantage of the latter. It is enough to consider the so-called impedance tube since the IBC (2) is monodimensional (i.e. only the normal component u • n is needed at the boundary, which stems from the locally reacting hypothesis). Let us briefly recall the impedance tube exact solution, expressed in the notations of Fig. 3. At x = 0 the incoming characteristic is given by ( p + u)(t, 0) = φ s (t) where φ s is causal (i.e. φ s (t) = 0 for t < 0). The nonlinear IBC ( 8) is enforced at x = L x . Assuming null initial conditions p(0, x) = u(0, x) = 0, the exact solution is

p(t, x) = 1 2 φ s t - x c 0 + 1 2 B(φ s ) t - 2L x -x c 0 , u(t, x) = 1 2 φ s t - x c 0 - 1 2 B(φ s ) t - 2L x -x c 0 . ( 69
)
Note that the solution, derived with an elementary application of the method of characteristics, is naturally expressed with the nonlinear scattering operator B. If the IBC is linear, then B(v) = β v and the analytical solution can also be derived using the Laplace transform. As detailed in Sec. 3, the computation of B(φ s ) may involve delayed ODEs.

The LEEs (1) are discretized using a fourth-order DG method (see Appendix D) and the RKF84 eight-stage fourth-order 2N-storage Runge-Kutta method [71, Tab. A.9]. The mesh of 240 triangles (2400 nodes) is shown in Fig. 3. The source is imposed using the numerical flux (A(n F )v) * s (see Sec. 4.1) and the IBC is weakly enforced with the Z-flux ((39), ( 59)) or the B-flux ((39), ( 63)). Fig. 4 shows the exact and computed solutions at the inlet x = 0 for four impedance models and a Gaussian source centered at 2 kHz. For the sake of clarity, the impedance tube is chosen long enough to avoid overlapping between the incident and reflected waves. For each of the four models the exact solution matches the computed one, which validates the proposed numerical flux functions.

The top right graph covers the linear TDIBC βA given in Table 2 and plotted in Fig. 4 that models the MP liner studied in Sec. 6. The computation of βA v through (34) involves a transport equation whose N ψ -point discretization is characterized by its number of points per wavelength [START_REF] Douasbin | Delayed-time domain impedance boundary conditions (d-tdibc)[END_REF]: here, N ψ = 4 is sufficient. The bottom left graph covers the non-passive scattering operator B = 3I, outside of the scope of the analysis presented in this paper. However, the B-flux does enable to compute this case, by contrast to the Z-flux with which no stable computations could be obtained, even with a reduced Fig. 6. Flow duct considered in Sec. 6 (aspect ratio of the GIT). Mesh of 106 triangles used in Fig. 14. Pressure measurements are taken on the lower wall y = 0. See Table 1 for geometrical dimensions corresponding to the GIT and GFIT.

Table 1

Geometrical dimensions of the GIT and GFIT. The acoustical cut-off frequency is computed through f cut-off = c 0 (2L y ) -1 , with c 0 = √ γ rT 344.32 m.s -1 (T = 295 K, γ = 1.4, and r = 287.058 J.kg -1 .K -1 ). CFL number and different meshes. Finally, the right shows the harmonic distortion, genuinely nonlinear feature, yielded by the algebraic nonlinear model [START_REF] Botteldooren | Finite-difference time-domain simulation of low-frequency room acoustic problems[END_REF] with C nl = 1. As expected, both the Z and B fluxes lead to the same solution.

L x1 (mm) L x2 (mm) L x (mm) L y (mm) f cut-off (
In Fig. 4, the CFL number has been kept at 0.85, its maximum value with a hard wall, for all cases. However, the analysis of Sec. 4 has shown that the B-flux is to be favored over the Z-flux. Fig. 5 illustrates this point with the nonlinear model used in Fig. 4. With the Z-flux, the maximum CFL number that leads to a stable computation decreases as the incident SPL increases, since the resistance increases. By contrast, the B-flux enables to keep the same CFL number, which can be understood from the passivity condition (9) that implies that B is a contraction. The value C nl = 1 has been chosen arbitrarily, to highlight the nonlinear effect; as the nonlinear coefficient C nl (resp. contraction coefficient C c ) goes to infinity (resp. zero), the SPL value above which the CFL drops with the Z-flux goes to zero. Although this example is elementary, it illustrates the advantage of the nonlinear scattering formulation. A more advanced but linear application is presented in the next section.

Numerical application in linear duct aeroacoustics

In this section, numerical simulations of two flow ducts are compared with experimental data, a summary of which is given in Sec. 6.1. Sec. 6.2 describes the employed numerical methodology. Comparison with experimental data is done in the last two sections, Secs. 6.3 and 6.4.

Experimental methodology and data

The study focuses on two ducts designed by the National Aeronautics and Space Administration (NASA), namely the Grazing Incidence Tube (GIT) and Grazing Flow Impedance Tube (GFIT). A short summary of the experiments reported in [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF] and [START_REF] Primus | ONERA-NASA cooperative effort on liner impedance eduction[END_REF] is provided below.

Experimental setup. The GIT and GFIT share a similar geometry, described in Fig. 6 and Table 1. The acoustical source is placed upstream of the entry plane x = 0, where a reference microphone is positioned on the lower wall. The source is monochromatic, with a frequency chosen below the lowest cut-off frequency so that only plane waves propagate in the duct. A near-anechoic termination is placed at the exit plane x = L. The sound absorbing material sample is mounted on the top wall between L x1 and L x 2 , while microphones are located on the lower wall. A Pitot probe system enables to measure the base flow u 0 at given cross-sections. See [18, Fig. 2] and [START_REF] Primus | ONERA-NASA cooperative effort on liner impedance eduction[END_REF]Fig. 5] for visualizations of the ducts.

Impedance identification methodology. At each source frequency f s , pressure measurements are taken on the lower wall y = 0 and used to identify the impedance value ẑid (j2π f s ). The inverse problem relies on the 2D convected Helmholtz equation (∂ t + u 0 ∂ x ) 2 pc 2 0 p = 0 derived from the LEEs (1) by taking a uniform base flow u 0 = c 0 M avg e x , where M avg is the measured average Mach number. Since the base flow is assumed uniform, the IBC used for the identification is not (3), as considered in this paper, but the standard Ingard-Myers boundary condition with uniform impedance and straight boundary: [18, Eq. 2]

-ẑ id (s)∂ n p = s c 0 p + 2M avg ∂ x p + c 0 M 2 avg s ∂ 2 x p.
Although in principle more accurate results could be obtained by using corrected versions of the Ingard-Myers boundary condition (see Remark 11), the identified impedance ẑid has proven satisfactory. Note that in addition to ẑid another impedance is identified, namely the exit impedance enforced at the outlet x = L: this impedance is close to 1 since the exit is nearly anechoic in the experiments.

Experimental data. The material considered in the GFIT experiment is a single degree of freedom micro-perforated liner The experimental data, partially reported in [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF], covers an incident SPL of 130 dB, average Mach numbers in {0,0.079,0.172,0.255,0.335,0.400}, and frequencies ranging from 0.5 kHz to 3 kHz by steps of 0.1 kHz. The 0.5 kHz data is excluded since it is less reliable [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF].

Numerical methodology

Discretization. The discretization of the LEEs ( 1) is identical to that of Sec. 5. In all of the computations, the CFL number is kept at 0.85, its maximum value with a hard wall (i.e. the TDIBC does not reduce the time step). At the entry plane x = 0, a plane wave source is imposed using the flux (A(n F )v) * s , while a non-reflecting boundary condition, exact only for plane waves, is imposed at x = L x using (A(n F )v) * out (see Sec. 4.1). Note that no exit impedance is considered, by contrast with the impedance identification methodology used in [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF][START_REF] Primus | ONERA-NASA cooperative effort on liner impedance eduction[END_REF]. The TDIBC is given by [START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF] and is weakly enforced with the B-flux 39), ( 63)) derived in Sec. 4.4.

((
Base flow. The chosen base flow u 0 = c 0 M 0 (y)e x is the hyperbolic velocity profile [17, Eq. 2.4a]

M 0 (r) = M c tanh(δ -1 (1 -|r|)) + M c 1 -tanh(δ -1 ) δ -1 (1 + tanh(δ -1 )) + 1 + |r| (1 -|r|) M c = M avg δ ln(cosh(δ -1 )) + [1 -tanh(δ -1 )][(6δ) -1 (1 + tanh(δ -1 )) + 2 /3] -1 , (70)
where r ∈ [-1, 1]. This velocity profile has two parameters: the average Mach number M avg and the nondimensional boundary layer thickness δ ∈ (0, 1]. For δ = 1, the velocity profile is almost identical to the Poiseuille profile M 0 (r) = 3 2 M avg (1 -r 2 ).

More accurate alternatives to this velocity profile include using a turbulent eddy viscosity model where M 0 is known analytically [START_REF] Marx | Effect of turbulent eddy viscosity on the unstable surface mode above an acoustic liner[END_REF], or, when possible, an interpolation of the experimental flow. Using these other velocity profiles does not significantly alter the presented results, but imply the use of a more refined mesh, hence why ( 70) is chosen herein.

In both experiments, the identified impedance values exhibit a dependency on the grazing base flow. In Secs. 6.3 and 6.4, to account for this dependency, we follow a simple approach inspired by the experimental study [START_REF] Kirby | The impedance of perforated plates subjected to grazing gas flow and backed by porous media[END_REF], whose result has been recalled in Sec. 3.1. The grazing flow is considered as an additional parameter, so that the physical quantities found in the acoustical model [START_REF] Khamis | Viscous effects on the acoustics and stability of a shear layer over an impedance wall[END_REF], such as l p or l c for instance, are tweaked when a base flow is present. This empirical approach has been found to be sufficient to match the experimental results.

Post-processing. A polychromatic source p s (t) = f s ∈I sin(2π f s t) is used, and each frequency is then separated at the postprocessing stage using a sixth-order recursive band-pass filter designed and applied using the MATLAB ® Signal Processing Toolbox ™ (functions designfilt and filter). The filtered pressure signals are then used to compute root-mean-square (RMS) values p RMS at each frequency and microphone locations along the bottom wall. The decibel (dB) values are computed with p dB = 20 log 10 (p RMS ) + C , where the constant C is chosen so that p dB (x = 0) matches the experimental value at the reference microphone. (Since the numerical scheme is linear, decibel values are indeed defined up to an additive constant.) The simulation is performed over 80 periods of the lowest frequency and convergence of the RMS value is checked.

Grazing Flow Impedance Tube (GFIT)

Following the three-step methodology described in Sec. 3.3.1 discrete models β are built to match the values identified in the GFIT at M avg = 0 and M avg = 0.271, see Table 2.

M avg = 0 case. The physical model [START_REF] Khamis | Viscous effects on the acoustics and stability of a shear layer over an impedance wall[END_REF] with coefficients given by ( 19), [START_REF] Liu | Stability analysis and design of time-domain acoustic impedance boundary conditions for lined duct with mean flow[END_REF] provides a fair initial point for the nonlinear least-squares optimization of step 1, the output of which is plotted in Fig. 7. The optimized model accounts for the resistance increase at 0.4 kHz and exhibits a high frequency behavior close to that of the non-optimized model (not plotted), with anti-resonances located around the approximate values f n = nc 0/2l c . The output of step 2, i.e. the discrete oscillatory-diffusive representation of the physical model obtained using a linear least-squares optimization, is therefore satisfactory. In step 3, it is used as an initial guess to build the final discrete model βA (the time delay τ is not optimized) shown in Fig. 7. A large number of poles, namely three pairs of oscillatory poles s n and two high frequency diffusive poles ξ k , have been ) Six proportional-integralderivative TDIBCs βPID , each matching the identified impedance at one frequency only. ( ) MP liner in the GFIT at M avg = 0 and 120 dB [START_REF] Primus | ONERA-NASA cooperative effort on liner impedance eduction[END_REF]. to validate this model above 2.6 kHz. This sensitivity is also exhibited by step 3. If the time delay τ is optimized, it is greatly reduced (i.e. τ τ ) thus modifying the high frequency behavior. To illustrate this point, Fig. 9 shows two discrete models.

The model βC is obtained by keeping l c to its physical value in the initial guess of step 1, choosing only one pair of oscillatory poles s n in step 2, and optimizing on τ in step 3. It provides an adequate approximation of the experimental data but a poor high frequency behavior, linked to its negligible delay τ τ /54. The model βB is obtained by using the physical model shown in Fig. 9 during step 1, choosing three pairs of oscillatory poles s n and one high frequency diffusive pole ξ k in step 2, and keeping τ constant in step 3.

Fig. 10 plots the computed SPL distributions. The dotted curves enable to check the relevance of the identified impedance values. The strongest disparities are obtained below 1 kHz where the measurements suggest the presence of a longitudinal resonance not modeled with our non-reflecting outlet, also noticeable in the no flow case, see Fig. 8. At these low frequen- cies, the fidelity of the SPLs obtained with βB and βC is therefore bound to be limited. Although the TDIBC B is significantly more expensive than C , both lead to similar SPLs with no significant discrepancies. The time delay of βC is negligible on the considered frequency range so that the minimal value N ψ = 2 is already too large, as shown by the value of PPW τ .

Grazing Incidence Tube (GIT)

M avg = 0 case.

Step 1 delivers an adequate set of coefficients for the physical model [START_REF] Khamis | Viscous effects on the acoustics and stability of a shear layer over an impedance wall[END_REF], see Fig. 11. The model βD is obtained by choosing the first two pairs of oscillatory poles s n in step 2 and optimizing the time delay τ in step 3 (although it can be kept constant as well). The slight decay of the physical model at high frequency is linked to a small but non-null diffusive part that can be captured by adding diffusive poles ξ k to βD [5, § VI.B]. However, this is not needed here in view of the computed SPLs shown in Fig. 12. Compared to the MP liner considered in Sec. 6.3, both a higher delay and a higher maximum frequency lead to a sensible increase in N ψ . The discrete model βE , plotted in Fig. 13, is obtained by using only 2 oscillatory poles s n and excluding the data points in/at [0.7, 0.8], [1, 1.4], 1.8, and [2.1, 2.3] kHz. Fig. 14 shows that βE compares favorably to the experimental pressure measurements. The increased resistance at 1 kHz compared with the identified value reduces the instability, although not enough to fit the experiment. The discrepancy at 3 kHz, since it is also obtained with the identified value, could be explained by the presence of higher-order modes in the experiment, which can be assessed using the cut-off frequency of a hard walled duct with uniform flow, given by c 0 (2L y ) -1 1 -M 2 avg = 3.01 kHz. The discrete model βF is obtained by excluding the data points in/at [0.7, 1.2], 1.8, and [2.1, 2.3] kHz. Compared to βE , the addition of one pair of oscillatory poles s n and two high frequency diffusive poles ξ k leads to an overfit that significantly alter the high frequency behavior, see Fig. 13. However, the corresponding SPLs are satisfactory, especially at 3 kHz. Since the delay of both models is close, they share a common value of N ψ .

Conclusion and perspectives

The semi-discrete energy analysis presented in Sec. 4 has shown the computational advantage of using the scattering operator over the impedance and admittance operators when solving the linearized Euler equations. In the analysis, the IBC has only been assumed to be admissible in the sense of Definitions 1 and 6. The derived numerical flux, the so-called B-flux ((39), ( 63)), has been interpreted in Remark 26 as that obtained from solving the monodimensional Riemann problem at the impedance boundary, which is standard in computational fluid dynamics. In practice, if the scattering operator cannot be computed efficiently or is not available, then the impedance (resp. admittance) is to be preferred a priori if the material is close to a pressure-release wall (resp. hard wall). A consistent but not passive numerical flux has been exhibited in Example 24, which justifies the need for such a careful analysis of the numerical flux and highlights the fact that the sole definition of the discrete model is not enough to fully define a TDIBC. An example of nonlinear scattering operator, derived from a physical perforation impedance model in Sec. 3, has been shown in Sec. 5 to deliver a stability CFL number independent of the incident sound pressure level, by contrast with a formulation based on the impedance operator. The analytical solution [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF], sufficient for validating of a nonlinear TDIBC, has been recalled.

Within the linear realm, the oscillatory-diffusive representation of a physical reflection coefficient model, applicable to acoustical liners, has been given in Sec. 3 and has led to the definition of the time-delayed broadband discrete model [START_REF] Olivetti | Direct numerical simulation of turbulent flow with an impedance condition[END_REF] that bears some similarities with discrete models postulated in the literature. In order to perform time integration with a Runge-Kutta method, a purely time-local formulation [START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF], based on a realization of the time delay through a transport equation, has been applied to two flow ducts in Sec. 6.

More marginally, Sec. 2 and Appendix A have recalled the interest of using the Laplace transform over the Fourier transform to deal with the admissibility conditions of impedance models.

The limitations of this work suggest ideas for future studies; five of them are listed below, the first three being numerical. First, a similar energy analysis could be applied to Ingard-Myers type IBCs in order to gain theoretical insights into their discretizations. Second, the presented nonlinear numerical application could be extended to a differential impedance operator Z. In such a case, the computational challenge lies in the efficient computation of B, which must avoid computing (z -1 0 Z + I) -1 since it may involve a nonlinear stiff ODE. Third, the present study can be seen as a first step towards the analysis of the enforcement of a TDIBC with the Euler or Navier-Stokes equations. Challenges include the filtering of the perturbation as well as the treatment of the diffusion term at the impedance boundary. The last two suggestions are physical. Fourth, the present work shows that TDIBCs benefit from physical knowledge about the impedance. This calls for additional investigations into the suitability of available models in the presence of grazing flow and broadband sources, as well as the inclusion of uncertainties during both identification and simulation. Last, the physical relevance of using a nonlinear impedance model with the LEEs could be investigated, by comparison with experimental data gathered in an impedance tube at high SPL.

Proof. Let w

h := ((w u h ) , w p h ) ∈ P k n (T h ) n+1 . The definition of A {0} h yields (A {0} h v n h , w h ) L 2 ( ) := 1 2 (M 0 (n)v h -A(n)v h , w h ) L 2 (∂ ) = c 0 2 ∂ ⎛ ⎜ ⎜ ⎝ α 1 u h • n (α 2 -1) p (α 3 -1)u h • n α 4 p ⎞ ⎟ ⎟ ⎠ • ⎛ ⎜ ⎜ ⎜ ⎝ w u h • n w u h • n w p h w p h ⎞ ⎟ ⎟ ⎟ ⎠ dx.
The Cauchy-Schwarz inequality then gives (A

{0} h v n h , w h ) L 2 ( ) ≤ c 1/2 0 √ 2 (N(n)v h , v h ) 1/2 L 2 (∂ ) w h L 2 (∂ )
∈ P k n (T h ) n+1 , A {0} h v h + A {β} h β( ph + u h • n) L 2 ( ) ≤ C tr 2 -1 /2 c 0 h -1 /2 ph -u h • n -β( ph + u h • n) L 2 (∂ ) ,
where C tr > 0 is defined in Lemma 30.

Proof. Let w

h := ((w u h ) , w p h ) ∈ P k n (T h ) n+1 . The definition of A {0} h + A {β} h
and the Cauchy-Schwarz inequality yield

(A {0} h v h , w h ) L 2 ( ) + (A {β} h β( ph + u h • n), w h ) L 2 ( ) := 1 2 (M 0 (n)v h -A(n)v h , w h ) L 2 (∂ ) + 1 2 (m B (n)β( ph + u h • n), w h ) L 2 (∂ ) = c 0 2 ( ph -u h • n -β( ph + u h • n), w p h -w u h • n) L 2 (∂ ) ≤ c 0 √ 2 ph -u h • n -β( ph + u h • n) L 2 (∂ ) w h L 2 (∂ ) ,
and the conclusion follows from Lemma 30. 2

The non-boundary terms of A h , gathered in Ãh , are covered by the last lemma below.

Lemma 33. For any v

h ∈ P k n (T h ) n+1 , Ãh v h L 2 ( ) ≤ C DG h -1 /2 |u 0 | + c 0 1/2 L ∞ ( ) |v h | upw + A(∇ h )v h + B v h L 2 ( ) ,
where

C DG := ( √ 2 + 1 √ 2 )C tr > 0 with C tr defined in Lemma 30.
Proof. This is a standard estimate that relies on the Cauchy-Schwarz and discrete trace inequalities, see [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]Chap. 3 & 7] and [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF]. Let w h ∈ P k n (T h ) n+1 with k ∈ N * . Each of the three terms in ( 54) is estimated separately. The first term readily yields

T ∈T h (Av h , w h ) L 2 (T ) = (A(∇ h )v h + B v h , w h ) L 2 ( ) ≤ A(∇ h )v h + B v h L 2 ( ) w h L 2 ( ) .
For the last two terms, we use the inequality The proof follows the energy method, in the spirit of the analysis carried out at the continuous and semi-discrete levels in Secs. 2.3 and 4.3. A comprehensive analysis of the energy method to derive CFL-type stability condition for coercive problems can be found in [START_REF] Levy | From semidiscrete to fully discrete: stability of Runge-Kutta schemes by the energy method[END_REF]. Note that the stability conditions derived in [START_REF] Levy | From semidiscrete to fully discrete: stability of Runge-Kutta schemes by the energy method[END_REF] cannot be directly applied herein since A h is not coercive in the sense of [START_REF] Levy | From semidiscrete to fully discrete: stability of Runge-Kutta schemes by the energy method[END_REF]Eq. 3.2], see [START_REF] Wang | Wave equation stabilization by delays equal to even multiples of the wave propagation time[END_REF]. See also [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]Chap. 6], [66, § 7.5.2], and [START_REF] Cohen | Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations[END_REF]Chap. 7].

Given the expression of the continuous energy [START_REF] Cremer | Theory regarding the attenuation of sound transmitted by air in a rectangular duct with an absorbing wall, and the maximum attenuation constant produced during this process[END_REF], we define the discrete energy as E n h := v n h 2 L 2 ( ) . By taking the scalar product of (52) with v n h and using the identity ab = 1 2 a 2 + 1 2 b 2 -1 2 (ab) 2 , one obtains the standard discrete energy balance [75, Eq. 3.4] 1

2 E n+1 h = 1 2 E n h + 1 2 v n+1 h -v n h 2 L 2 ( ) -t(A h v n h , v n h
). The expression of (A h v h , v h ) L 2 ( ) given by (37) yields 1 2

E n+1 h = 1 2 E n h + 1 2 v n+1 h -v n h 2 L 2 ( ) - t 2 (C(u 0 )v n h , v n h ) L 2 ( ) - t 2 |v n h | 2 upw - t 2 (M 0 (n)v n h , v n h ) L 2 (∂ ) ,
where each term has a clear interpretation:

|)h -1 |v h | 2 upw - t 2 (I(n)v h , v h ) L 2 (∂ ) ,
where I(n) := 1 2 (M 0 (n) + M 0 (n) ) - .

The conclusion then follows from Lemma 19. 2

The proof of Proposition 28 is similar, but relies on an energy balance for the additional variables ϕ h , given herein.

(Although energy balances associated with diffusive representation of positive-real functions are available in the literature, this is not the case for bounded-real functions.) Let β be given by [START_REF] Ern | Theory and Practice of Finite Elements[END_REF] and denote b = β a for concision. The realization of β is, as seen in Sec. The energy of this realization can be defined as from which we deduce that a sufficient condition for β(s) to be bounded-real (i.e. for β to be an admissible scattering operator) is that the measure μ satisfies μ L 1 μ ξ 2 L 1 ≤ 1. Using (C.3), we can then prove Proposition 28.

Proof. (Proposition 28) The proof relies on an energy analysis similar to that of Proposition 21 with one additional technicality: the energy is not reduced to v h 2 L 2 ( ) but must include a contribution E β from the diffusive variable ϕ h defined as (C.2). For the sake of clarity, we first write down the semi-discrete energy balance before considering its discretization.

(Semi-discrete energy balance) Using the energy balance associated with E β (C.3), the semi-discrete energy balances reads

1 2 d dt E β ≤ 1 2 ( ph + u h • n) 2 -β( ph + u h • n) 2 - μ L 1 2 ξ max ξ min (-ξ ϕ h + ph + u h • n) 2 dμ ξ 2 1 2 d dt v h 2 L 2 ( ) = -1 2 (C(u 0 )v h , v h ) L 2 ( ) -1 2 |v h | 2 upw -1 2 (M 0 (n)v h , v h ) L 2 (∂ ) -(A {β} h v h , v h ) L 2 (
) . The identity [START_REF] Zennaro | Natural continuous extensions of Runge-Kutta methods[END_REF] enables to rewrite the second line to explicit the energy exchange between the diffusive variables ϕ h and the acoustic field v :

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 2 d dt E β ≤ 1 2 ( ph + u h • n) 2 -β( ph + u h • n) 2 - μ L 1 2 ξ max ξ min (-ξ ϕ h + ph + u h • n) 2 dμ ξ 2 1 2 d dt v h 2 L 2 ( ) = -1 2 (C(u 0 )v h , v h ) L 2 ( ) -1 2 |v h | 2 upw -c 0 4 ph + u h • n 2 L 2 (∂ ) -β( ph + u h • n) 2 L 2 (∂ ) -c 0 4 ph -u h • n -β( ph + u h • n) 2 L 2 (∂ ) .
The semi-discrete energy balance is obtained by integrating the first inequality over ∂ , multiplying it by c 0/2 and summing it with the second one 1 2 d dt

E h ≤ - 1 2 (C(u 0 )v h , v h ) L 2 ( ) - 1 2 |v h | 2 upw - c 0 4 μ L 1 ξ max ξ min -ξ ϕ h + ph + u h • n 2 L 2 (∂ ) dμ ξ 2 - c 0 4 ph -u h • n -β( ph + u h • n) 2 L 2 (∂ ) , (C.4)
where the extended energy is defined as

E h := v h 2 L 2 ( ) + c 0 2 E β .
The third and fourth terms of the right-hand side express the dissipation that occurs at the impedance boundary (here, ∂ ).

(Discrete energy balance) The explicit Euler discretization of [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF] 1 2

E n+1 h ≤ 1 2 E n h - t 2 (C(u 0 )v n h , v n h ) L 2 ( ) - t 2 |v n h | 2 upw - c 0 4 t μ L 1 ξ max ξ min -ξ ϕ n h + pn h + u n h • n 2 L 2 (∂ ) dμ ξ 2 - c 0 4 t pn h -u n h • n -β( pn h + u n h • n) 2 L 2 (∂ ) + 1 2 v n+1 h -v n h 2 L 2 ( ) + c 0 2 μ L 1 ξ max ξ min 1 2 ϕ n+1 h -ϕ n h 2 L 2 (∂ ) dμ ξ .
The task is now to estimate the two anti-dissipative terms so that they can be provably controlled by the three dissipative ones. Let us bound the anti-dissipation on v h with Lemmas 32 and 33:

1

2 v n+1 h -v n h 2 L 2 ( ) = t 2 2 Ãh v h + A {0} h v h + A {β} h β( ph + u h • n) 2 L 2 ( ) ≤ t 2 Ãh v h 2 L 2 ( ) + t 2 A {0} h v h + A {β} h β( ph + u h • n) 2 L 2 ( ) ≤ C 2 DG t 2 (c 0 + |u 0 |)h -1 |v h | 2 upw + C 2 tr 2 t 2 c 2 0 h -1 ph -u h • n -β( ph + u h • n) 2 L 2 (∂ ) .
Using Lemma 32 and the identity

  [56, § I.4]. Values for the three coefficients can be obtained through various theoretical or empirical models; the high-frequency approximation of a model derived by Crandall, based on the axisymmetric Stokes equations, yields[START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] Eq. 7]

1 C 2 tr a 0 1 + a 2 0 ,

 1200 -1 C -2 tr . By contrast, consider the 0-flux (α 2 = α 3 = 0) which yields CFL ≤ so that the CFL number decreases to 0 as a 0 → ∞ or a 0 → 0. This recovers a result of Ventribout [49, § 2.3].

Fig. 2 . 2 1 + k 2 2

 222 Fig. 2.Eigenvalues λ h ∈ σ (-A h ) for the acoustical cavity = (0, 1) 2 with c 0 = 1, z 0 = 1, and impedanceZ(u • n) = a 0 u • n + a Q u • n on ∂ . DG6 (N K = 8triangles, 504 DoF). ( ): Z-flux ((39), (59)). ( ): flux ((39), (60)) with α = β 0 . ( ): λ = jω = ±jc 0 k 2 1 + k 2 2 with two monodimensional dispersion relations

Fig. 3 .Fig. 4 . 2 p ref z0 10 SPL 20 (

 3421020 Fig. 3. Impedance tube used in Sec. 5: L x = 1 m, L y = Lx /100, and mesh with 240 triangles.

Fig. 5 .

 5 Fig. 5. Maximum CFL number against SPL for the impedance tube of Fig. 4. (

(

  MP) made from a honeycomb core (thickness l c = 38.1 mm, cell diameter d c = 9.5 mm) sandwiched between a perforated facesheet (thickness l p = 0.8 mm, hole diameter d p = 0.3 mm, and porosity σ p = 5 %) and a rigid backplate. Impedance identifications have been done at an incident SPL of 120 dB (chosen to minimize nonlinear effects), average Mach numbers in {0,0.180,0.271}, and frequencies ranging from 0.4 kHz to 2.6 kHz by steps of 0.2 kHz. The GIT experiment focused on a ceramic tubular liner (CT57) made from a ceramic tubular core (l c = 85.6 mm, d c = 0.6 mm, σ c = 57 %) and a rigid backplate.

Fig. 8 .

 8 Fig. 8. RMS values of acoustic pressure p RMS on the lower wall y = 0 (see Fig. 6). DG4-RKF84-B-flux on 110 triangles (1100 nodes), CFL = 0.85. (

Fig. 9 . 1 c= 2 .

 912 Fig. 9. Impedance ẑ and reflection coefficient β. ( ) ẑid /z 0 . MP liner in the GFIT at M avg = 0.271 and 120 dB [72]. ( ) ẑMP /z 0 (17) with optimized coefficientsσ -1 p a 0 = 2.220 × 10 -14 , σ -1 p a1 /2 = 2.220 × 10 -14 s 1/2 , σ -1 p a 1 = 2.220 × 10 -14 s, σ -1 c = 2.116, b 0 = 2.220 × 10 -14, b1 /2 = 7.67 × 10 -3 s 1/2 , and

Fig. 10 .

 10 Fig. 10. RMS values of acoustic pressure p RMS on the lower wall y = 0. DG4-RKF84-B-flux on 326 triangles (3260 nodes), CFL = 0.85. Base flow (70) with M avg = 0.271 and δ = 0.2. ( ) TDIBC βB with N ψ = 3 so that PPW τ (2.6 kHz) = 6.79. ( ) TDIBC βC with N ψ = 2 so that PPW τ (2.6 kHz) = 248. (

M

  avg = 0.4 case. The identified impedance values display a very low resistance within [0.7, 1.1] kHz and the shape of | βid | suggests the presence of noise in the data, see the bottom left of Fig.13, as one may expect due to the high value of M avg . The output of step 1 is sensitive to the chosen experimental points since the physical model (17) cannot fit both the low and high resistance regions, namely [0.7, 1.1] kHz and[2.3, 3] kHz. However, step 1 is not sensitive to the initial guess provided that the chosen data points cover the anti-resonance. Additionally, the 1 kHz value is best removed since it leads to an instability with the LEEs as shown in the dotted curve at the top right of Fig.14, and investigated in [74, § VI.B.2]. (Including the 1 kHz point in the optimization process has been found to systematically lead to an unstable pole, i.e. (s n ) > 0 or ξ k < 0.) Due to the shortcomings of the physical model, step 3 has a tendency to overfit the experimental data, so that we do not optimize on the delay τ during step 3.

Fig. 11 .Fig. 12 .

 1112 Fig.11. Impedance ẑ and reflection coefficient β. ( ) ẑid /z 0 . CT57 liner in the GIT at M avg = 0 and 130 dB[START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF]. (

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Impedance ẑ and reflection coefficient β. ( ) ẑid /z 0 . CT57 liner in the GIT at M avg = 0.4 and 130 dB [18]. (

α 1 -C 2 tr CFL(α 2 1 + (α 3 - 1 ) 2 ) α 2 +α 3 2 α 2 +α 3 2 α 4 -C 2 tr CFL(α 2 4+

 12231222422 N(n) + N(n) ) so that (I(n)v h , v h ) R n+1 = ( Ǐ(n) vh , vh ) R 2 with vh := (u h • n, ph ) and Ǐ(n) := c 0 (α 2 -1) 2 )

3 . 2 ,

 32 ∂ t ϕ(t, -ξ) = -ξ ϕ(t, -ξ) + a(t), ϕ(0, -ξ) = 0 b(t) = ξ max ξ min ϕ(t, -ξ) dμ(ξ ), (C.1)

2 L 1 a 2 - μ L 1 2 ξ

 222 max ξ min (-ξ ϕ h + a) 2 dμ ξ 2 , (C.3)

2

  Remark 9 (IBC for wave equation). Without base flow, i.e. u 0 = 0, both p and u obey the wave equation, and the IBC (2) reduces to ∂ t p = Z (-c 0 z -1 0 ∂ n p). Hence, a hard (resp. pressure-release) wall corresponds to a Neumann (resp. Dirichlet) boundary condition on p, while a proportional impedance Z ∝ I yields a Robin boundary condition. Studies focused on the wave equation typically use the linear admittance formulation ∂ n p = -c -1

	0 z 0 y ∂ t p [37, Eq. 1.1] [36, Eqs. 5, 6] [35, Eq. 2]

  ). The symmetry of A(n) implies, by definition, that the operator ∂ t + A(∇) is a

	Friedrichs-symmetric operator; equivalently, (12) is a symmetric hyperbolic system [53, Def. 1.2]. Using an integration by
	parts, this symmetry property yields

  and the conclusion follows from Lemma 30. 2

	Lemma 32. For any v h

  follows from the fact that the real symmetric matrix A is diagonalizable with A = P P , where:= diag(u 0 • n + c 0 , u 0 • nc 0 , u 0 • n, • • • , u 0 • n) and P is an orthogonal matrix (recall that |A| := P ||P ). (This inequality yields the (DG8)condition in[START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF].) Using (B.2), the Cauchy-Schwarz inequality, and Lemma 30 giveF )Jv h K, {w h }) L 2 (F ) ≤ √ 2C tr h -1 /2 max |v h | upw w h L 2 ( ) , F )|Jv h K, Jw h K) L 2 (F ) ≤ |v h | upw w h L 2 ( ) .Remark 34. Assuming a uniform base flow (∇u 0 ≡ 0) and a finite volume discretization(k = 0), A(∇ h )v h + B v h L 2 ( ) = 0,and the estimate of Lemma 33 is reduced toÃh v h L 2 ( ) ≤ C DG h -1/2 (|u 0 | + c 0 ) 1/2 |v h | upw ,which is used in Sec. 4 to derive a sufficient stability condition using an energy method.

	Appendix C. Proofs of stability conditions
	Proof. (Proposition 21)
							1 /2
	h F ∈F i	h (A(n F ∈F i	ρ(|A(n F )|)
	where ρ(|A(n F )|) denotes the spectral radius. Similarly,
	1 2	h F ∈F i	(|A(n 1 √ 2	h F ∈F i C tr h -1 /2 max	ρ(|A(n F )|)	1 /2

(A(n F )v h , w h ) L 2 (F ) ≤ (|A(n F )|v h , v h ) 1 /2 L 2 (F ) (|A(n F )|w h , w h ) 1 /2 L 2 (F ) , (B.2) which

The final estimate follows from max

F ∈F i h ρ(|A(n F )|) ≤ max F ∈F i h (|u 0 | + c 0 ) ≤ |u 0 | + c 0 L ∞ ( ) . 2

  ( ) is the anti-dissipation due to the explicit nature of the time-marching scheme;-(C(u 0 )v n h , v n h ) L 2 ( ) is the contribution of a non-constant base flow u 0 ; -|v n h | 2upw is the dissipation due to the upwind flux (would be null with a centered flux);-(M 0 (n)v n h , v n h ) L 2 (∂ )is the dissipation due to the IBC.The key step of the proof is to estimate the anti-dissipation term1 the energy production that occurs at the impedance boundary is bounded by A ( ) . The estimate of this term is provided by Lemma 31, which does not require the finite volume hypothesis. Lemma 33 provides a bound for Ãh v n h L 2 ( ) that does need the finite volume hypothesis. Using both Lemmas 31 and 33, the discrete energy balance becomes

							1 2 v n+1 h	-v n h	2 L 2 2 v n+1 h	-v n h	2 L 2 ( ) . The decomposition of A h readily
	yields the estimate
	1 2	v n+1 h	-v n h	2 L 2 ( ) =	t 2 2	A h v n h	2 L 2 ( ) ≤ t 2 Ãh v n h	2 L 2 ( ) + t 2 A h v n {0} h	2 L 2 ( ) ,
	{0} h v n h L 2 1 2 2 E n+1 h ≤ 1 2 1 -t min

where x∈ λ min (u 0 ) E n h -t 2 1 -2C 2 DG t(c 0 + |u 0

  leads to the identity (ab)b = 1 2 a 2 -1 2 b 2 -1 2 (ab)2 enables to rewrite this with the discrete extended energy E n

								ξ max
	(v n+1 h	-v n h , v n h ) L 2 ( ) +	c 0 2	μ L 1		(ϕ n+1 h	-ϕ n h , ϕ n h ) L 2 (∂ )	dμ ξ
							ξ min
										ξ max
	≤ -	t 2	(C(u 0 )v n h , v n h ) L 2 ( ) -	t 2	|v n h | 2 upw -	c 0 4	t μ L 1	-ξ ϕ n h + pn h + u n h • n 2 L 2 (∂ )	dμ ξ 2
										ξ min
		-	c 0 4	t pn
	c 0 2 μ L 1 ξ max ξ min ϕ n h	2 L 2 (∂ )	dμ ξ :

hu n h • nβ( pn h + u n h • n) 2 L 2 (∂ )

,

and h := v n h 2 L 2 ( ) +

Acknowledgements

This research has been financially supported by the French ministry of defense (Direction Générale de l'Armement) and ONERA (the French Aerospace Lab). The authors are grateful to Dr. Michael Jones (NASA) for providing the CT57 and MP experimental data, Dr. Sébastien Pernet (ONERA) for his insights into the analysis of Discontinuous Galerkin methods, and Dr. Hugues Deniau (ONERA) for invaluable exchanges on computational fluid dynamics.

Table 2

TDIBCs (31) obtained in Sec. [START_REF] Melling | The acoustic impedance of perforates at medium and high sound pressure levels[END_REF]. Only poles in the upper half-plane {s | (s) ≥ 0} are given (the full set is obtained by complex conjugation).

(s n , -ξk ) n,k (rad.s -1 ) r1,n , μ1,k r2,n , μ2,k βA β∞ = 1 τ = 2.074709e-04 s -9.425004e+04 ) ẑMP /z 0 (17) with optimized coefficients σ -1

kHz, and τ = 2.075 × 10 -4 s. (This curve assumes that the time delay term e -2π j f τ is perfectly approximated.) chosen to build a broadband approximation. The corresponding SPLs along the lower wall are given in Fig. 8 for two choices of N ψ to illustrate the impact of the delay discretization: a value of N ψ = 4 is sufficient here. The agreement with experimental pressure measurements is satisfactory, the largest error occurring at 1.4 kHz. However, based on the plot of | βid | given at the bottom left of Fig. 7, this point appears to be an outlier. To confirm this, the SPLs obtained with the identified impedance values are also given. They are computed using six proportional-integral-derivative (PID) impedance models ẑPID,i (s) = d 0,i s -1 + d 1,i + d 2,i s, implemented as [START_REF] Tudisco | Application of the time-domain impedance boundary condition to large-eddy simulation of combustion instability in a shear-coaxial high pressure combustor[END_REF] with τ = 0, each chosen so that ẑPID,i (j2π f i ) = ẑid (j2π f i ). The plot shows that at 1.4 kHz the tuned PID model also exhibits a significant error.

M avg = 0.271 case. The identified values exhibit a strong resistance increase compared to the no flow case. Here, step 1 requires some care: the optimization of the physical model coefficients ( 17) is strongly dependent on the initial point, in stark contrast to the other three cases considered in this paper. A contributing factor to this sensitivity is the lack of anti-resonance in the experimental data, which stops at 2.6 kHz. Fig. 9 plots an optimized model that exhibits two anti-resonances, obtained by doubling l c in the initial guess (( 17), ( 19), ( 21)). Additional experimental data would be needed

Appendix A. Difference between Fourier and Laplace transforms for impedance models and admissibility conditions

Many works related to impedance boundary conditions use the Fourier transform, at least formally. The objective of this appendix is to briefly present why the Laplace transform should be preferred to the Fourier transform when dealing with causal functions. The Fourier transform is defined over the set of tempered distributions S (R) (i.e. generalized functions with at most polynomial growth at infinity) [4, Chap. V] [51, Chap. 7] [50, § 1.8]. This implies that the use of the Fourier transform is not always possible, as it does not apply to an exponentially growing function like t → H(t)e t for instance. The Laplace transform is defined over the set of causal distributions D + (R) with at most finite exponential growth at infinity

Fourier transform of physical impedance models

The Laplace transform ẑ of a physical impedance model is a holomorphic function of the open right half-plane { (s) > 0}, see e.g. the examples of Sec. 3.1. Since it is a positive-real function, the Fourier transform F (z) exists and belongs to S (R) [50, § 3.5]. However, the mere substitution of jω (ω ∈ R) for s in the expression of ẑ(s) does not yield the Fourier transform, i.e. F (z)(ω) = ẑ(jω). Mathematically, this is due to the fact that the limit lim σ →0 ẑ(σ + iω) may only exist in the space of tempered distribution S (R).

The impedance model ẑ(s) = 1 /s, which arises in the modeling of ground layers [5, Eq. 24], is a textbook example of this phenomenon. The function ω → 1 /jω exhibits a non-integrable singularity at ω = 0. As a result it is not a distribution and thus cannot be a Fourier transform. Taking the limit in S (R) yields the standard expression

where pv denotes the Cauchy principal value. It is important to note that F -1 pv 1 jω (t) = 1 2 sign(t) is not causal: the term "π δ" is needed for causality. Similarly, consider the impedance of a lossless cavity (see Sec.

,

where the infinite sum converges in S (R). Note that the π /b 1 -periodic distribution F -1 [pv (coth(jb 1 ω))] is not causal, as illustrated by its Fourier series expansion pv (coth(jb

c n e j2nb 1 ω , with c 0 = 0, c n = -1, and c -n = -c n for n ≥ 1.

Admissibility conditions

In Sec. 2.2, admissibility conditions have been formulated using the Laplace transforms ẑ, ŷ, and β, see Propositions 3 and 8. These conditions are straightforward to verify in practice using the expressions of physical models. However, formulating the same admissibility conditions using the Fourier transform is more involved. Let us illustrate this point by considering the causality condition for a tempered distribution z ∈ S (R). with a sufficiently large real part. This is consistent with the fact, recalled in Sec. 2.2, that a real (continuous) LTI system cannot be both anticausal and passive. Note that a formal use of the Fourier transform can lead to a mistake, since this kernel satisfies both ẑ(jω) = ẑ(-jω) and (ẑ(jω)) ≥ 0 for any ω ∈ R but is not passive.

Appendix B. Discontinuous Galerkin estimates

This appendix gathers the few technical results needed for the stability analysis presented in Sec. 4. The background material for the analysis can be found in [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF], [START_REF] Ern | Theory and Practice of Finite Elements[END_REF], and [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF]. As stated in Sec. The key result needed is the standard trace inequality, which is recalled below.

Lemma 30 (Discrete trace inequality). Let (T h

) h be a quasi-uniform sequence of simplicial, geometrically conformal, and shape-regular meshes. For any dimension n ∈ N * and degree k ∈ N * , there is a constant C tr > 0 such that ∀h > 0, ∀T ∈ T h , for any face F T ,

The constant C tr only depends upon n, k, and the mesh regularity parameter ρ 1 := inf h inf T ∈T h r T h T with r T the radius of the largest ball inscribed in T .

Proof. Let h > 0 and T ∈ T h . Since T h is a simplicial, geometrically conformal, and shape-regular mesh, we have [66, Eq. (1.39)] h

, where C tr,s only depends upon ρ 1 , n, and k. By definition, the quasi-uniformity of T h gives a constant C > 0 such that Ch ≤ min T ∈T h h T . 2

The estimation of the impedance boundary terms in A h ((53), ( 67)) relies on the two following lemmas. Note that these lemmas imply a departure from the framework developed in [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF], since the condition (DG3a) does not hold.

Lemma 31. For any v

and C tr > 0 is defined in Lemma 30.

Appendix D. Implementation details

The implemented numerical scheme is a standard DG method. It is nodal in the sense that on each of the N K triangles the solution v h is computed through its components in the Lagrange basis associated with the N p points described in [69, § 6.1]. For a N-th order DG, N p = N(N + 1)/2. Let us now assume that the linear TDIBC is given by

which corresponds to the TDIBC [START_REF] Olivetti | Direct numerical simulation of turbulent flow with an impedance condition[END_REF], where w is a shorthand for " ph + u h • n". The global assembly is performed in three steps.

(1) Spatial discretization. The global DG formulation reads

where v := (v h (x k i )) i,k is the discrete acoustic field (3N K N p elements), and M, K , and F are the standard mass, stiffness, and flux DG matrices. The IBC manifests itself through the (rectangular) observation matrix C z that are associated with the DoF that belong to the impedance boundary z . The operators Q i are applied to each component of the vector C z v of length N z . Note that, at this stage, the actual definition of the operators Q i does not matter.

(2) State-space realization. By construction, the operators Q i have a state-space realization with state vector ϕ :

where v now denotes the extended state

and the matrices are obtained by concatenation (the same notations M, K , and F s are used for the sake of concision). For instance,

is of length N ϕ N z . This is a finite-dimensional delay differential equation that can be advanced in time with a variety of methods.

(3) Time delay computation. Following Sec. where ψ k j is vector of length N ψ . Note that (D.4) can be written as (D.2), i.e. it is a state-space realization of the time delay. Combining (D.3) with (D.4) leads to the final global formulation M v + K v = F s v s , where the vector v, of length 3N K N p + (N ψ + 1)N ϕ N z , now includes the acoustic field [v h (x k i )] i,k as well as the additional variables (ϕ k j ) k, j and (ψ k j ) k, j . The addition of a nonlinear term Q nl to the discrete model B does not change the assembly process described above, which yields

In the case of (25), Q nl (w) can be directly computed. If it exists, Q nl (w) can be replaced by its nonlinear state-space realization.