domain: phys spi chim

text:
Ilyes Moufid

Rémi Roncen

Denis Matignon

Estelle Piot

Time-domain simulation of the acoustic nonlinear response of acoustic liners at high sound pressure level

This paper focuses on the numerical modeling of the acoustic response of perforated-plate liners at high sound pressure levels in the time domain. In order to do so, a porous-based description of the perforated plate is used to represent the visco-thermal processes occurring inside the perforated plate. It is done through the use of the equivalent fluid model (EFM) containing two irrational transfer functions described herein by a generic model covering the Johnson Champoux Allard Pride Lafarge model (JCAPL), the JCAL and JCA models. The nonlinear phenomena occurring at high sound pressure levels are taken into account by using the Forchheimer's correction in the time-domain EFM, introducing a quadratic nonlinearity in the equations. The formulation of the nonlinear EFM equations in the time domain leads to an augmented system for which a proof of stability is given. From the nonlinear EFM, an approximate model is built for numerical simulations with a multipole approximation of the transfer functions. Sufficient stability conditions are provided for the nonlinear multipolebased approximate EFM. A numerical scheme using a discontinuous Galerkin method is developed to validate the model against experiments with perforated-plate liners.

Introduction

The mitigation of aircraft engine noise is currently done through the use of acoustic liners [START_REF] Michael | Broadband and Low-Frequency Acoustic Liner Investigations at NASA and ONERA[END_REF]. The most classical geometry is a single-degree-of-freedom liner (SDOF), which consists in a honeycomb core placed between a perforated plate and a rigid backplate. The perforated plate is the main source of the acoustic-energy dissipation which comes from visco-thermal friction occurring in the vicinity of the perforation of the plate, and from vortex shedding appearing at high sound pressure level (SPL) [START_REF] Ingård | Acoustic Circulation Effects and the Nonlinear Impedance of Orifices[END_REF]Labate, 1950, Melling, 1973]. The depth of the cavities characterizes the range of frequencies where the sound attenuation is maximal. These liners are known to be very efficient in a narrow range of frequencies, tailored to meet the noise signature of the treated system [START_REF] Girvin | Aircraft noise-abatement and mitigation strategies[END_REF]. Moreover, the acoustic behavior of liners depends nonlinearily on the SPL, even at moderate SPL values where the laws of linear acoustics can still be considered valid, far from the liner.

In order to find materials which produce sound attenuation on a larger range of frequencies, more involved geometries of sound absorbing materials are studied experimentally and in numerical simulations. In this paper, the focus is on time-domain numerical simulations, which are more suited to evaluate broadband signals, to account for moving acoustic sources, or to consider nonlinear phenomena stemming from the liner's response.

Basic examples of other studied absorbent materials are double degree of freedom liners [START_REF] Gautam | On the acoustic performance of double degree of freedom Helmholtz resonator based acoustic liners[END_REF], perforated plates backed by a porous medium [START_REF] Tayong | Sound absorption of a micro-perforated plate backed by a porous material under high sound excitation: measurement and prediction[END_REF][START_REF] Peng | Sound absorption of a porous material with a perforated facing at high sound pressure levels[END_REF] or porous materials alone [START_REF] Cao | Porous materials for sound absorption[END_REF]. The latter is usually a two-phase medium with a solid skeleton filled by a fluid (air herein) by means of pores of relatively small size. Consequently a wide variety of media can be considered as a porous medium, including perforated plates [START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF]. The approach used to model the sound absorption of material in the time domain generally depends on the type of material: surface impedance Z as a boundary condition for liners, volume modeling for porous media.

On the one hand, numerical investigations used to simulate the acoustic response of liners under broadband noise are classically based on time-domain impedance boundary conditions (TDIBC) [Tam andAuriault, 1996, Fung and[START_REF] Fung | Broadband Time-Domain Impedance Models[END_REF]. This approach relies on the assumption that the material is locally reacting, meaning that the transmitted waves propagate only in the normal direction within the material. The numerical impedance approximations currently used in the literature are known as multipole models (MM) [START_REF] Reymen | Time-Domain Impedance Formulation Based on Recursive Convolution[END_REF][START_REF] Li | Improved Multipole Broadband Time-Domain Impedance Boundary Condition[END_REF]. They are defined by a discrete sum of elementary first or second-order low-pass systems, also called IIR filter in digital signal processing [START_REF] Alexander | Chapter 2 -Fundamental DSP Concepts[END_REF], and have been shown to have two important benefits: straightforward verification of admissibility conditions [Rienstra, 2006, Dragna and[START_REF] Dragna | Physically admissible impedance models for time-domain computations of outdoor sound propagation[END_REF] and versatility. An MM involves convolution products in the time domain, which can be readily computed by means of additional first-order ordinary differential equations (ODE) using the recursive convolution technique [START_REF] Luebbers | A frequency-dependent finitedifference time-domain formulation for dispersive materials[END_REF]. This technique was already present in previous works [START_REF] Cockburn | Numerical Resolution of Maxwell's Equations in Polarisable Media at Radio and Lower Frequencies[END_REF][START_REF] Carcione | Viscoacoustic wave propagation simulation in the earth[END_REF] and was more recently called the auxiliary differential equation (ADE) method [START_REF] Dragna | A generalized recursive convolution method for timedomain propagation in porous media[END_REF]. It results in an augmented system which can be solved by classical time-integration schemes, without requiring storage of previous time-step solutions.

On the other hand, the equivalent fluid model (EFM) [Lafarge, 2009] is generally used to describe the macroscopic acoustic behavior of rigid porous media for aeronautic applications1 [START_REF] Attenborough | Outdoor ground impedance models[END_REF]. In this framework, porous media are equivalent, in terms of acoustic behavior, to that of a fluid with a frequency-dependent effective density and a frequency-dependent effective compressibility. The latter two are irrational transfer functions defined by means of the dynamic tortuosity α and the normalized dynamic compressibility β, respectively. Both functions are known in the low (LF) and high (HF) frequency limits when additional assumptions are made on the pore shape [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF][START_REF] Philibert | à suivre et des formules à employer dans les questions de distribution d'eau, etc[END_REF][START_REF] Linton | Theory of dynamic permeability and tortuosity in fluidsaturated porous media[END_REF][START_REF] Lafarge | Propagation du son dans les matériaux poreux à structure rigide saturés par un fluide viscothermique: Définition de paramètres géométriques, analogie électromagnétique, temps de relaxation[END_REF]. However, no exact definition can be retrieved at all frequencies except in the simplest cases. Several works following different approaches have resulted in various semi-phenomenological models predicting α and β on the whole frequency range, such as the Höroshenkov model [Horoshenkov et al., 2019[START_REF] Kirill | Erratum: A three-parameter analytical model for the acoustical properties of porous media[END_REF], the Wilson model [START_REF] Wilson | Relaxation-matched modeling of propagation through porous media, including fractal pore structure[END_REF][START_REF] Wilson | Simple, relaxational models for the acoustical properties of porous media[END_REF], or the Johnson Champoux Allard Pride Lafarge (JCAPL) model [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF][START_REF] Linton | Theory of dynamic permeability and tortuosity in fluidsaturated porous media[END_REF][START_REF] Lafarge | Propagation du son dans les matériaux poreux à structure rigide saturés par un fluide viscothermique: Définition de paramètres géométriques, analogie électromagnétique, temps de relaxation[END_REF][START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF][START_REF] Pride | Drag forces of porous-medium acoustics[END_REF]. As for the impedance of liners, they are usually approximated by MM in numerical simulations [START_REF] Zhao | An equivalent fluid model based finitedifference time-domain algorithm for sound propagation in porous material with rigid frame[END_REF][START_REF] Alomar | Time-domain simulations of sound propagation in a flow duct with extended-reacting liners[END_REF], and a recent study [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] showed that real-parameter multipole approximations of α and β are well suited for the aforementioned models related to conventional porous material, as opposed to metamaterial which can have negative material properties on some frequency ranges [START_REF] Craster | Acoustic Metamaterials: Negative Refraction, Imaging, Lensing and Cloaking[END_REF]Guenneau, 2013, Bellis and[START_REF] Bellis | Simulating transient wave phenomena in acoustic metamaterials using auxiliary fields[END_REF].

In that respect, a liner is modeled in time-domain numerical simulations by TDIBC either by a direct approximation of its related surface impedance Z computed by means of experimental data or using a porous-based definition of the perforated plate [START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF]. A volume-based modeling of perforated plate using the EFM in the linear regime (i.e. at low SPL) induces a higher numerical cost compared to a TDIBC, due to the meshing of the material. Hence, there was no motivation to explore this volume-based approach for classical locally-reacting liners.

At high SPL, the impedance of absorbent materials was suggested to have an additional term depending on the square of the normal velocity:

Z nl (u • n) = ρ 0 C l nl |u • n| u • n, with n the outward unit normal to the surface of the liner, u the particle velocity, ρ 0 the density and a nonlinear coefficient C l nl . The latter has several definitions, depending on a vena contracta coefficient [START_REF] Cummings | Transient and multiple frequency sound transmission through perforated plates at high amplitude[END_REF] or the discharge coefficient [START_REF] Melling | The acoustic impendance of perforates at medium and high sound pressure levels[END_REF], and all definitions depending on the porosity ϕ of the plate.

In porous media, it is known that the Darcy law [START_REF] Philibert | à suivre et des formules à employer dans les questions de distribution d'eau, etc[END_REF] characterizing the static resistivity σ 0 of the medium has to be replaced by the total resistivity σ at high Reynolds number (defined according to the microgeometry of porous media) because the porous-medium resistivity varies with the amplitude of the velocity [START_REF] Mcintosh | Nonlinear wave propagation through rigid porous materials. I: Nonlinear parametrization and numerical solutions[END_REF]Lambert, 1990, Aurégan and[START_REF] Aurégan | Measurement of the nonlinear behavior of acoustical rigid porous materials[END_REF]. Two regimes can be distinguished: a low-Reynolds regime where the increase of the resistivity is quadratic in the amplitude of the velocity, and a high-Reynolds regime with a linear dependency following the Forchheimer's law [START_REF] Forchheimer | Wasserbewegung durch Boden[END_REF]:

σ = σ 0 (1 + C p nl U ) ,
where U = ∥u∥ represents the amplitude of the velocity u and C p nl is related to the Forchheimer's coefficient. This nonlinear coefficient C p nl was first measured [START_REF] Umnova | Behavior of rigid-porous layers at high levels of continuous acoustic excitation: Theory and experiment[END_REF], then deduced by a method similar to matched asymptotic expansions [START_REF] Wang | Sound absorption of porous metals at high sound pressure levels[END_REF], and more recently straightforwardly defined with an analytical expression depending on the porosity and the discharge coefficient [START_REF] Laly | Acoustical modeling of micro-perforated panel at high sound pressure levels using equivalent fluid approach[END_REF]. The porous-based description of perforated plates coupled with the Forchheimer's correction of the resistivity was used by [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF] to define a time-domain admittance boundary condition (TDABC) modeling the acoustic behavior of liners by a surface impedance. This approach showed a good prediction from numerical simulations of the acoustic behavior of liners under high SPL.

Another approach not yet explored to model the nonlinear acoustic behavior of perforated-plate liners is to use a porous-based acoustic description of perforated plate with an appropriate time-domain EFM, i.e. a volume-based modeling of perforated plates. It would be beneficial to the study of extended-reacting liner and modeling more complex sound absorption materials with porous media coupled with perforated plates for example. Linear timedomain EFM giving good prediction of the acoustic behavior of porous media [START_REF] Alomar | Time-domain simulations of sound propagation in a flow duct with extended-reacting liners[END_REF][START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF], the objective of the present paper is to extend this approach to predict the acoustic behavior of perforated plate liners at high SPL in the Forchheimer regime. To that end, the starting point for this study is the time-domain EFM presented in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] based on oscillatory-diffusive (OD) representation [START_REF] Hélie | Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses[END_REF]Matignon, 2006, Monteghetti et al., 2016] used to describe the two irrational transfer functions α and β in the time domain.

In the present work, a brief description of the linear EFM based on OD representations of α and β [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] is first given in Section 2. It is followed by a description of the proposed nonlinear EFM with the investigation of its well-posedness in Section 3. Then in Section 4, a multipole approximation of α and β is built to approximate the nonlinear EFM using the ADE method, as done in [START_REF] Blanc | Wave simulation in 2D heterogeneous transversely isotropic porous media with fractional attenuation: a Cartesian grid approach[END_REF][START_REF] Alomar | Extended-reacting liners in time-domain simulations for broadband attenuation with flow[END_REF][START_REF] Xie | A Discontinuous Galerkin Method for Wave Propagation in Orthotropic Poroelastic Media with Memory Terms[END_REF]. Next, the validity and the efficiency of the proposed approach is demonstrated numerically in Section 6, 7 and 8, with a Discontinuous Galerkin (DG) finite element method combined with a Runge-Kutta (RK) time scheme. Finally, a conclusion is drawn in Section 9.

Generic equivalent fluid model

The frequency-dependent equations controlling the acoustic behavior of a rigid porous material are expressed in the Laplace domain:

ρ 0 α(s) s u = -∇ ∇ ∇ p , χ 0 β(s) s p = -∇ ∇ ∇ • u, (1a) (1b)
where h is the one-sided Laplace transform of a locally integrable function h with finite exponential growth, i.e. h ∈

L 1 loc ([0, ∞)) with |h(t)| < A e σ h t , defined as h(s) := ∞ 0 h(t) e -st
dt in the right half-plane ℜ(s) > σ h , with σ h the convergence abcissa. The macroscopic fluid velocity u(t, x) and the macroscopic acoustic pressure p(t, x) are obtained by averaging the microscopic velocity field and microscopic pressure field over a representative elementary volume [Lafarge, 2009], ρ 0 is the fluid density and χ 0 is the adiabatic compressibility.

Different models describing α and β can be found in the literature. Herein, the generic model (2) introduced in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] is used, covering the JCAPL, JCAL and JCA models.

α g (s) = α ∞     1 + M s + N 1 + s L -1 s     , (2a)
β g (s) = γ - γ -1 1 + M ′ s + N ′ 1 + s L ′ -1 s , ( 2b
)
with the two physical parameters: the high-frequency limit of the tortuosity α ∞ and the heat capacity ratio γ.

The model parameters M , N , L, M ′ , N ′ and L ′ are defined in Table 1 and depend on the kinematic viscosity ν, the porosity ϕ, the static viscous permeability k 0 (σ 0 = ρ 0 ν/k 0 is the static resistivity), the characteristic viscous length Λ, the static viscous tortuosity α 0 , the static thermal permeability k ′ 0 , the characteristic thermal length Λ ′ , the static thermal tortuosity α ′ 0 and the Prandtl number Pr. Note that the second term in the brackets in (2a), with M = ϕσ 0 /ρ 0 α ∞ , is related to Darcy's law. [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] describing the JCA, JCAL and JCAPL models.

Model

L M N L ′ M ′ N ′ JCA νϕ 2 Λ 2 /4k 2 0 α 2 ∞ νϕ/k 0 α ∞ νϕ/k 0 α ∞ 16ν/Λ ′ 2 Pr 8ν/Λ ′ 2 Pr 8ν/Λ ′ 2 Pr JCAL νϕ 2 Λ 2 /4k 2 0 α 2 ∞ νϕ/k 0 α ∞ νϕ/k 0 α ∞ νϕ 2 Λ ′ 2 /4k ′ 0 2 Pr νϕ/k ′ 0 Pr νϕ/k ′ 0 Pr JCAPL ν/Λ 2 α0 α∞ -1 2 νϕ/k 0 α ∞ 2ν/Λ 2 α0 α∞ -1 ν/Λ ′ 2 (α ′ 0 -1) 2 Pr νϕ/k ′ 0 Pr 2ν/Λ ′ 2 (α ′ 0 -1) Pr
It was shown in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] that the dynamic tortuosity (2a) admits the OD representation

α g (s) = α ∞ 1 + M s + N ∞ L µ g (ξ) s + ξ dξ , (3a) µ g (ξ) := 1 π ξ L -1 ξ , ( 3b
)
where µ g is called the diffusive weight of the diffusive part ∞ L µ g (ξ)/(s + ξ) dξ of the transfer function α g , and the normalized dynamic compressibility (2a) admits the OD representation

β g (s) = 1 + (γ -1) r 0 s -s 0 + N ′ ∞ L ′ ν g (ξ) s + ξ dξ , (4a)
ν g (ξ) := 1 π ξ ξ L ′ -1 (ξ -M ′ + N ′ ) 2 + N ′2 ξ L ′ -1 , ( 4b
)
where ν g is the diffusive weight of the diffusive part ∞ L ′ ν g (ξ)/(s + ξ) dξ of the transfer function β g and an isolated term r 0 /(ss 0 ) exists if and only if M ′ -N ′ -L ′ < 0. In this latter case,

s 0 := L ′ λ 2 0 -1 < 0, (5)
is a negative pole of β g with an associated positive residue

r 0 := 2L ′ λ 0 M ′ + N ′ (λ 0 -1) 2L ′ λ 0 + N ′ > 0, (6)
where the parameter λ 0 is defined as

λ 0 := -N ′ + N ′ 2 -4L ′ (M ′ -N ′ -L ′ ) 2L ′ ∈ (0, 1). ( 7
)
On the contrary; if M ′ -N ′ -L ′ > 0, β g has no pole and r 0 is taken null. Before applying the inverse Laplace transform on the linear EFM (1) in order to have a time-domain EFM, two diffusive variables are introduced [START_REF] Haddar | Theoretical and numerical analysis of the Webster Lokshin model[END_REF]Matignon, 2008, Lombard andMatignon, 2016]:

φ φ φ(ξ) := u/(s + ξ), (8a)
ψ(ξ) := p/(s + ξ). (8b)
A rearrangement of the latter two definitions gives, respectively, s φ φ φ(ξ) = -ξ φ φ φ(ξ)+ u, and s ψ(ξ) = ξ ψ(ξ)+ p. Hence, the convolution products appearing in the time-domain EFM, when expressed with the generic model ( 2) can be computed by means of additional ordinary differential equations (ODE):

                         ∂ t u(t, x) + M u(t, x) + N ∞ L µ g (ξ) ∂ t φ φ φ(ξ) dξ = - 1 ρ 0 α ∞ ∇ ∇ ∇p(t, x), ∂ t p(t, x) + (γ -1) r 0 ∂ t ψ(-s 0 ) + N ′ ∞ L ′ ν g (ξ) ∂ t ψ(ξ) dξ = - 1 χ 0 ∇ ∇ ∇ • u(t, x), ∂ t φ φ φ(ξ) = -ξ φ φ φ(ξ) + u(t, x), ∂ t ψ(ξ) = -ξ ψ(ξ) + p(t, x). ( 9a
) (9b) (9c) (9d)
The initial conditions of the diffusive variables are set non-null [START_REF] Lombard | Diffusive Approximation of a Time-Fractional Burger's Equation in Nonlinear Acoustics[END_REF]:

   ϕ ϕ ϕ(ξ; 0, x) = u(0, x)/ξ, ψ(ξ; 0, x) = p(0, x)/ξ, (10a) (10b)
in order to keep a finite problem at t = 0. The proof of the well-posedness of the time-domain EFM (9) with the initial conditions (10) was given in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF]. In the latter, a positive-definite global energy functional

E lin (t) := E m (t) + E diff (t), (11)
was defined by summing the classical mechanical energy E m with an energy E diff defined for the diffusive variables. The global energy E l was proven to decrease over time for classical porous media when there is no external contribution. More details about E l are gathered in Appendix A.

This Section 2 presented the linear time-domain EFM based on OD representations of α and β for the JCA, JCAL and JCAPL models. The two irrational transfer functions expressed with the Wilson model and the Horoshenkov model also admit an OD representation, given in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF]. Therefore, the same methodology can be used to verify the well-posedness of the EFM expressed with these models.

The next section focuses on the extension of the model ( 9) to take into account the nonlinear phenomena at high SPL by using the Forchheimer's correction.

Nonlinear generic equivalent fluid model

In the linear EFM, the resistivity σ 0 considered in the dynamic tortuosity α is assumed to be constant. As explained in the introduction, the total resistivity of a porous material is dependent on the velocity magnitude, with either a linear [START_REF] Irmay | On the theoretical derivation of Darcy and Forchheimer formulas[END_REF][START_REF] Beavers | Non-Darcy Flow Through Fibrous Porous Media[END_REF][START_REF] Joseph | Nonlinear equation governing flow in a saturated porous medium[END_REF][START_REF] Herbert | Attenuation of intense sinusoidal waves in air-saturated, bulk porous materials[END_REF] or a quadratic [START_REF] Rasoloarijaona | Nonlinear seepage flow through a rigid porous medium[END_REF][START_REF] Wodie | Correction non linéaire de la loi de Darcy[END_REF][START_REF] Firdaouss | Nonlinear corrections to Darcy's law at low Reynolds numbers[END_REF]] relation depending on the velocity regime. In this section, the Forchheimer's correction of the resistivity is written in the high Reynolds regime:

σ = σ 0 (1 + C f ∥u∥) . (12)
With the total resistivity σ, the linear EFM leads to a nonlinear EFM given in Section 3.1. An analysis of the energy of the proposed nonlinear model is then carried out in Section 3.2 to compare it to that of the linear EFM.

Analytical nonlinear model

The EFM is studied with the dynamic variables α and β expressed by the generic model (2), where the static resistivity σ 0 is included in the parameters M , N and L; except for the JCAPL model, where only M depends on σ 0 = η/k 0 (see Table 1). Thus, a correction of σ 0 into σ would bring different modifications depending on the model used. Moreover, for the JCAL and JCA models, the introduction of σ in L includes a term containing the velocity magnitude in the irrational part of the dynamic tortuosity. Such a correction [START_REF] Umnova | Behavior of rigid-porous layers at high levels of continuous acoustic excitation: Theory and experiment[END_REF]] leads to a time-domain expression of the EFM not straightforwardly usable for theoretical and numerical analysis, due to a velocity-dependent diffusive weight µ g in (9a).

A focus on the asymptotic behavior of α [Johnson et al., 1987, Avellaneda and[START_REF] Avellaneda | Rigorous link between fluid permeability, electrical conductivity, and relaxation times for transport in porous media[END_REF] α

LF (jω) = ϕσ 0 ρ 0 ν ν jω + α 0 + O ω→0 (jω), (13a)
α HF (jω) = α ∞ + 2α ∞ Λ ν jω + O ω→∞ 1 jω , (13b)
shows an absence of dependence on σ 0 at HF, frequencies at which the irrational term 1/ √ s prevails. Thus, if the resistivity σ 0 is present in the square root in α g (2a), it is there only artificially, introduced by the relaxation function used to link the LF behavior to the HF behavior of α in the JCA and JCAL models. The independence of the irrational part of α with respect to the resistivity in the JCAPL model is thus an advantage, not present in the JCA and JCAL models, when one wishes to apply the Forchheimer's correction (12) to take into account the total resistivity. [START_REF] Turo | Influence of Forchheimer's nonlinearity and transient effects on pulse propagation in air saturated rigid granular materials[END_REF]] used a dynamic tortuosity expression similar to the JCAL model while circumventing the dependence of the irrational term on σ 0 by using the expression

α(s) = α ∞ + ϕσ 0 ρ 0 s + 2α ∞ Λ η ρ 0 s , ( 14
)
built by summing the contributions of the LF (second term) to the contributions of the HF (first and third terms), instead of linking them by a relaxation function. Thus, the application of the Forchheimer's correction is done by the single correction of the second term of ( 14) related to Darcy's law. Note that this second term is none other than the second term of the generic dynamic tortuosity (2a) where M = ϕσ 0 /ρ 0 α ∞ . Knowing that the Forchheimer's correction is applied to Darcy's law, expressed through the term M/s in the generic model (the second term of (9a) in the time domain), we choose to correct the static resistivity σ 0 by the total resistivity σ only in M , keeping the diffusive part unchanged. This correction leads to the following nonlinear EFM

                         ∂ t u + M u + M C f u ∥u∥ + N ∞ L µ g (ξ) ∂ t φ φ φ(ξ) dξ = - 1 ρ 0 α ∞ ∇ ∇ ∇p, ∂ t p + (γ -1) r 0 ∂ t ψ(-s 0 ) + N ′ ∞ L ′ ν g (ξ) ∂ t ψ(ξ) dξ = - 1 χ 0 ∇ ∇ ∇ • u, ∂ t φ φ φ(ξ) = -ξ φ φ φ(ξ) + u, ∂ t ψ(ξ) = -ξ ψ(ξ) + p, (15a) (15b)
(15c) (15d)
When the JCAPL model is used, system (15) is the one obtained in a direct way by replacing each static resistivity σ 0 with the total resistivity σ. For the JCA and JCAL models, system (15) does not take into account the Forchheimer correction that would be present in the diffusive part of (15a).

Well-posedness of the time-domain nonlinear EFM

Before the verification of a correct representation of the nonlinear effects by the model ( 15), a first question can be asked about the modification of the properties of the EFM when adding the term M C f u ∥u∥. In the linear case, the EFM has been shown to be stable for different expressions of α and β [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF]. Here, we investigate if this property is preserved, knowing that the coefficient C f is necessarily positive when the porous material is a passive medium.

By defining exactly the same energy functional E nl as the one for the linear EFM (9), but this time for the nonlinear EFM (15), leads to

d dt E nl (t) = d dt E lin (t) -ρ 0 α ∞ M C f Ω ∥u∥ 3 dx. ( 16
)
The result (16) makes possible to extend the stability of the linear EFM to the nonlinear case, based on the known behavior of E lin . It is synthesized in the next theorem.

Theorem 3.1. Given a bounded domain Ω, with no contribution at the boundary ∂Ω (i.e., p = 0, or u • n = 0 on ∂Ω), the nonlinear EFM (15) with the dynamic variables described by the generic model (2) is stable for positive model parameters.

Indeed, noting that for the same initial conditions, E(0) = E nl (0) with E ′ nl (t) ≤ E ′ (t) ≤ 0, the decay of the energy of the nonlinear EFM is straightforward, as soon as C f is positive.

To conclude, the use of the Forchheimer correction in the time-domain EFM does not change its well-posedness from the linear EFM as long as C f is positive. Moreover, the theoretical result given in Theorem 3.1 could be extended to a quadratic correction (low Reynolds regime), or to a more general correction under more general conditions as explained in Appendix B. Now that the analytical nonlinear model is given and proven stable, the next section focuses on developing an approximation more suitable for numerical solvers that keeps this stability property.

Multipole approximation of the nonlinear equivalent fluid model

In order to perform numerical simulations, an approximate model has to be built. To do so, the multipole-based approximations

α mm (s) = c ∞ 1 + c -1 s + K k=1 r k s -s k , (17a)
β mm (s) = c ′ ∞   1 + c ′ -1 s + K ′ k=1 r ′ k s -s ′ k   , (17b)
of α and β are used with real parameters r k , r ′ k , s k and s ′ k , following the same approach as [START_REF] Alomar | Time-domain simulations of sound propagation in a flow duct with extended-reacting liners[END_REF][START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF]. The coefficients c ∞ , c -1 , c ′ ∞ and c ′ -1 are related to the asymptotic behavior of α and β. Thus,

it immediately comes that c ∞ = α ∞ , c -1 = νϕ/α ∞ k 0 = M , c ′ ∞ = 1 and c ′ -1 = 0.
The parameters are not replaced by their values in this section, in order to keep a symmetrical expression of the equations and make the reading easier.

The approximations ( 17) can be obtained by using an optimization approach, as the vector-fitting (VF) algorithm [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF][START_REF] Gustavsen | Improving the pole relocating properties of vector fitting[END_REF][START_REF] Deschrijver | Macromodeling of Multiport Systems Using a Fast Implementation of the Vector Fitting Method[END_REF] which gives the optimal parameters of the multipole approximations α mm and β mm to best fit α and β for a given range of frequencies. A direct approach can also be used with a discretization of the diffusive part of the OD representations (3a) and (4a) based on quadrature formula [START_REF] Monteghetti | Time-local discretization of fractional and related diffusive operators using Gaussian quadrature with applications[END_REF]:

∞ L µ g (ξ) s + ξ dξ ≈ K k=1 r k s -s k . ( 18
)
Regardless of the approach used, the diffusive part in the time-domain equations ( 15a) and ( 15b) are approximated as follows:

∞ L µ g (ξ) ∂ t φ φ φ(ξ) dξ = ∞ L µ g (ξ) [-ξ φ φ φ(ξ) + u] dξ (19a) ≈ K k=1 [r k s k φ φ φ(-s k ) + r k u], (19b)
with an additional ODE for each evaluation of the diffusive variables (φ φ φ(-s k ) ≡ φ φ φ k and ψ(-s ′ k ) ≡ ψ k ). Moreover, Forchheimer's correction is applied to the coefficient c -1 , being the term associated with Darcy's law, in an analogous way to the analytical nonlinear EFM. It leads to the approximate nonlinear EFM:

                                       ∂ t u + 1 ρ 0 c ∞ ∇ ∇ ∇p + c -1 (1 + C f ∥u∥) u + K k=1 r k u + K k=1 r k s k φ φ φ k = 0, ∂ t p + 1 χ 0 c ′ ∞ ∇ ∇ ∇ • u +   c ′ -1 + K ′ k=1 r k   p + K ′ k=1 r ′ k s ′ k ψ k = 0, ∂ t φ φ φ k = s k φ φ φ k + u (∀k ∈ [[1, K ′ ]]), ∂ t ψ k = s ′ k ψ k + p (∀k ∈ [[1, K ′ ]]), φ φ φ k (0, x) = -u(0, x)/s k (∀k ∈ [[1, K ′ ]]), ψ k (0, x) = -p(0, x)/s ′ k (∀k ∈ [[1, K ′ ]]), (20a)
(20b) (20c) (20d) (20e) (20f)
Following [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF], a global energy functional E a nl can be defined for the approximate EFM (20), representing the discretization of E nl . The approximate energy E a nl , detailed in Appendix A, leads to sufficient conditions that keeps the approximate nonlinear EFM well-posed when describing classical porous media. These conditions are the same as the ones given for the linear approximate EFM in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF], summarized in the following property.

Property 4.1. Given a bounded domain Ω, with no contribution at the boundary ∂Ω (i.e., p = 0, or u • n = 0 on ∂Ω), the augmented energy E a nl of the approximate nonlinear EFM (20), with the dynamic variables described by a multipole model (17), is positive definite and decreasing, if the following two conditions are satisfied:

-the poles (s k , s ′ k ) k are negative real numbers, (C 1 ) -the weights (r k , r ′ k ) k are positive real numbers. (C 2 )
Note that the first condition (C 1 ) is a necessary one, while the second condition (C 2 ) on the weight sign is a sufficient one. Now that conditions ensuring the stability of the approximate nonlinear EFM are known, numerical simulations can be done after developing the numerical scheme, when proper multipole approximations of α and β are used.

Numerical approximation of the quadratic term in velocity

The numerical scheme used in this work is based on a fourth-order Runge-Kutta scheme [START_REF] Toulorge | Optimal Runge-Kutta schemes for discontinuous Galerkin space discretizations applied to wave propagation problems[END_REF] for the discretization in time and the space discretization is handled by a Discontinuous Galerkin (DG) scheme, a method well suited for acoustic problems [START_REF] Cohen | Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations[END_REF]. The RK-DG scheme is briefly detailed in Appendix C. In this section, a special attention to the approximation of the quadratic term in velocity in the DG solver is required because, it introduces a nonlinear term in the equations. Therefore, the methods used to discretize this term are presented in this section.

To understand how the nonlinear velocity term is accounted for in the numerical scheme, the variational formulation (50) used in the DG method is rewritten by omitting the time and space variables

0 = d dt Ti q h λ i j dΩ + Ti F(q h ) • ∇λ i j dΩ - ∂Ti F * (q e h , q i h ) • n i λ i j dσ + Ti b(q h ) λ i j dΩ, (21)
where q h = (u h v h p h ) is the solution vector with u h = (u h v h ), the λ i j are the j Lagrange polynomials defined on each triangle T i , F is associated to the terms coupling the two equations in u and p, and F * is the numerical flux.

The additional nonlinear term that must be calculated in the variational formulation is

∀j ∈ 1, d , Ti ∥u h ∥ u h λ i j dΩ, (22)
with

u h (t, x) = d k=1 u i,k h (t)λ i k (x), (23)
on a cell T i . The term ( 22) is taken into account in the last term of the right-hand side of the equality ( 21), and more precisely in the definition of b. The term u h λ i j in ( 22) is treated classically by calculating the products λ i j λ i k , but the presence of ∥u h ∥ prevents an exact calculation of ( 22) for any polynomial u h .

First strategy (quadrature).

A first method to approximate this integral is to use a quadrature formula (ω k , x k ) k , with ω k the quadrature coefficients and x k the quadrature nodes, so as to have

Ti ∥u h ∥ u h λ i j dΩ ≈ k ω k ∥u h (x k )∥ u h (x k ) λ i j (x k ). ( 24
)
On a cell T i where one of the components of u h changes sign, or if ∥u h ∥ is not a polynomial, the quadrature rule will give an approximation of the integral ( 22). In the case where the components of u h do not change sign and ∥u h ∥ can be expressed by a polynomial of degree p (in particular when the velocity field is unidirectional, along x or y), an exact approximation of the integral can be made with a quadrature rule by adapting the number of quadrature points to the degree of the polynomial ∥u h ∥ λ i j . In this latter case, the polynomial resulting from the product ∥u h ∥ u h λ i j is of the order of three times that of u h , and requires a larger number of quadrature nodes than the degrees of freedom for the DG. It is therefore necessary to interpolate u h on a larger set of quadrature nodes in order to keep an exact computation of (24) when the term in the integral is a polynomial. This first approach has the disadvantage of requiring a large number of quadrature points and an interpolation of the solution on these points at each iteration.

Second strategy (local).

A second method, less expensive, but without theoretical results on its accuracy, is to use the approximation

Ti ∥u h ∥ u h λ i j dΩ = Ti d k=1 u i,k h λ i k d k=1 u i,k h λ i k λ i j dΩ, (25a)
≈ Ti d k=1 u i,k h u i,k h λ i k λ i j dΩ, (25b)
≈ d k=1 u i,k h u i,k h Ti λ i k λ i j dΩ, ( 25c
)
where the velocity amplitude at each degree of freedom is directly multiplied by the value of the velocity at the associated degree of freedom. This gives a local approximation of the velocity amplitude contribution, in the sense that the velocity amplitude is taken to be spatially independent, implying a local contribution at each point x i j only. Therefore, the approximation ( 25) is solved numerically like the other terms in Ti b(q h ) λ i j dΩ, except that for (25), each component of the velocity vector is multiplied by its amplitude at each iteration due to the time dependency of the velocity amplitude.

The first method is called the quadrature approach (quad), while the second is called the local approach (loc). For a verification of the correct implementation of these methods and an analysis of their error, a numerical case where an analytical solution can be obtained is studied in the next section.

Study of a one-dimensional one-equation case

This section addresses the study of a one-dimensional numerical case in space where the theoretical solution is known, in order to obtain first results on the approximation error of the nonlinear term. The one-dimensional reduced case considered herein is

∂ t u + ∂ x u + C f |u|u = 0, u(0, x) = u 0 (x), (26a) (26b)
the theoretical solution of which is given by

u th (t, x) = u 0 (x -t) 1 + C f t|u 0 (x -t)| . ( 27
)
The function u th is obtained by using the method of characteristics on system (26). Note that the solution u th has a constant derivative u ′ th along the characteristics passing through a point x k where u 0 is equal to zero. This property of the u th function is illustrated in Figure 1 where the tangents at the points of sign change are shown in black dotted line for an initial solution u 0 (x) = xe -10x 2 and a nonlinear coefficient C f = 5. One can observed that the slope remains unchanged at these x k points. Particular attention will be paid to the sign change zones in the simulations to check the behavior of the numerical solution u when the latter has a different sign on the same mesh cell.

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 x -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15

u(t, x) t = 0 t = 1 t = 2 y = x -b 0 Figure 1: Solution of (26) at different times for u0(x) = xe -10x 2
in blue line with C f = 5. The slopes y = x + b0 in black dotted line are drawn for b0 equal to 0, 1 and 2.

System (1) is a 1D sub-case of the EFM, where the visco-thermal effects are not taken into account. Hence, this reduced system is solved by the numerical scheme implemented for the EFM, with u ≡ p, the addition of the term |p|p in the second equation of the EFM and all the coefficients taken as a unit except for C f , giving an artificial system of two identical equations

     ∂ t u + ∂ x p + C f |u|u = 0, ∂ t p + ∂ x u + C f |p|p = 0, u(0, x) = p(0, x) = u 0 (x). (28a) (28b) (28c)
Note that the nonlinear term |p|p of ( 28b) is added to keep the symmetry of the equations, allowing us to solve problem (26) through system (28). This artificial addition implies a higher numerical computation cost, but it allows us to verify the implementation and approximation of the nonlinear term directly on the RK-DG 2D code presented in Appendix C. Two different problems are solved, one where the initial solution

u pos (x) = exp -(x -x 0 ) 2 2W 2 0 ( 29
)
is positive, and a second with a mexican-hat initial solution (also known as a Ricker wavelet)

u mex (x) = 1 - (x -x 0 ) 2 W 2 0 exp -(x -x 0 ) 2 2W 2 0 , (30)
which changes in sign on the domain of definition. The two cases are shown in Figure 2 with W 0 = 0.08 and x 0 = 0.5. According to what was discussed in Section 5, the quadrature method is expected to give a better approximation of the nonlinear term in the first case than in the second case, where a larger error can be expected wherever the sign changes.

0.0 0.5 1.0 1.5 2.0 2.5 3.0 x 0.0 0.2 0.4 0.6 0.8 1.0 u(t, x) t = 0 t = 1 t = 2 (a) u0 = upos 0.0 0.5 1.0 1.5 2.0 2.5 3.0 x -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 u(t, x) t = 0 t = 1 t = 2 (b) u0 = umex
The numerical simulation based on ( 28) is carried out on the mesh given in Figure 3, with a DG solver of order varying from 2 to 5, a time step ∆t = 10 -3 and a final time t f = 2. A calculation of the error made on the numerical solution with respect to the theoretical solution is presented in Figure 4 for different values of C f . The error is first calculated without any nonlinear effects (C f = 0), giving a reference case, and then for two values of C f . The solid lines correspond to the error obtained by the local approach and the dashed lines to the quadrature approach. Relative error Relative error In Figure 4a, the difference between the theoretical solution associated with the case u 0 = u pos and the numerical solution calculated with this same initial condition is very similar whatever the approach used (local or quadrature), with an error decreasing by one decade per degree of polynomial used in the DG. This is therefore in agreement with what was expected and with the accuracy of the numerical scheme. In the second case where the initial condition is u mex , the order of convergence is lost as soon as the nonlinear coefficient is taken non-null, as shown in Figure 4b. The error nevertheless remains decreasing, with a very small difference between the two approaches (of the order of 1% of the error).

C f =0 (loc) C f =5 (loc) C f =20 (loc) C f =0 (quad) C f =5 (quad) C f =20 (quad) 4 • 10 -(2+k) 4 • 10 -(2+k)
C f =0 (loc) C f =5 (loc) C f =20 (loc) C f =0 (quad) C f =5 (quad) C f =20 (quad) 4 • 10 -(2+k) 4 • 10 -(2+k) (b) u0 = umex
To better understand the loss of convergence observed in Figure 4, the error is calculated at each point of the mesh and displayed in Figure 5a for t = 2. In Figure 5b, it is clear that the error is mainly located in the vicinity of the sign changes of the solution. The slope of the solution at these points nevertheless remains close to the theoretical one, which we know to be unchanged over time.

In conclusion, a clear loss of accuracy comes from the approximation of the nonlinear term in the numerical scheme. However, with a fine mesh and a DG solver of sufficient order (typically higher than 3), the error from the approximation of the nonlinear term is relatively small. Moreover, it should be noted that the error of the additional term is very local, occurring only at the sign changes of the solution. For this study, sufficiently fine meshes will be used to limit the error during numerical simulations, in which the local approach is used (the quadrature one giving similar results but with a greater time cost).

Numerical simulations of impedance-tube experiments

This section focuses on reproducing numerically the experiments carried out by [START_REF] Howerton | Application of Swept Sine Excitation for Acoustic Impedance Eduction[END_REF] and [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF] to compare the results of the approximate nonlinear EFM with the experimental results recently published in the literature. It contains a description of the studied case in Section 7.1, and a presentation of the numerical results together with their analysis in Section 7.2.

Description of the studied case

The two materials selected from the work of [START_REF] Howerton | Application of Swept Sine Excitation for Acoustic Impedance Eduction[END_REF] and [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF] are acoustic liners, composed of a perforated plate and cavities. The goal is to model the acoustic behavior of these liners using the approximate nonlinear EFM (20) through a porous modeling of the perforated plates. Numerically, this leads to adopt a volume description of the liners, different from the TDIBC usually performed in the literature [START_REF] Fung | Broadband Time-Domain Impedance Models[END_REF][START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF][START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF]. A first Section 7.1.1 describes the materials at the centre of this study, followed by a second Section 7.1.2 giving a description of the configuration and the experiment to be modelled. A last subsection 7.1.3 introduces the multipole approximations of the dynamic variables α and β used for the resolution of the approximate nonlinear EFM.

Porous-based modeling of acoustic liners

The geometric properties of the studied liners are given through the thickness e p of their perforated plate, the perforation radius r p , the porosity ϕ p , and the length l c of their cavities, all represented in Figure 6a. The values of these properties are given in Table 2, for the macro-perforated plate of Billard [2021, Table 3.1], referred to here as BMP1, and the perforated plate studied in [Howerton et al., 2019, Tab. 1] referred to as GE01, keeping the same name used in their study. Table 2: Geometric properties of perforated plates studied in [Billard, 2021, Tab. 3.1] and [Howerton et al., 2019, Tab. 1].

Plaque

ϕ p r p (mm) e p (mm) l c (mm) BMP1 0.1 0.8 1.5 16 GE01 0.087 0.5 0.635 37.465

Perforated plates can be modelled by a porous medium by relating the geometric properties of the plate to the intrinsic properties of a porous material by [Atalla andSgard, 2007, Allard and[START_REF] Allard | Propagation of Sound in Porous Media[END_REF]

α ∞ = 1, Λ = Λ ′ = r p , k 0 = k ′ 0 = ϕ p r 2 p 8 . (31)
By doing so, two porous media are defined with physical parameter given in Table 3 to have acoustic behavior equivalent to the two plates in Table 2.

In the Table 3, the static tortuosity α ∞ is set equal to 1 by definition, but it must be corrected by the relation (32) to take into account the acoustic radiation at the exit of perforations. The static tortuosity to be used is therefore

α ∞,e = 1 + 2ε e e p , (32)
with the collar correction ε e = 0.85 r p (1 -1.14 ϕ p ) [START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF]. Given the intrinsic parameters of the porous media, the dynamic tortuosity α and the dynamic compressibility β can be expressed by the JCA and JCAL model through the generic model. If values of α 0 and α ′ 0 had been known, the use of the JCAPL model could also have been considered. Finally, as the objective is to work at high SPL (between 120dB and 150dB), it is necessary to define the nonlinear coefficient C f to be used. It is done by identification from the expression of the total resistivity used in Laly et al. [2018, Eq. (20)] and Diab et al. [2022, Eq. (4)]:

ϕ p α ∞ α ∞,e Λ (m) Λ ′ (m) k 0 (m 2 ) k ′ 0 (m 2 ) σ 0 (N• s•m -4 ) BMP1 0
σ = σ 0 + d ρ 0 (1 -ϕ 2 p ) πϕ p e p C 2 D ∥u∥ , ( 33
)
where it is recalled that C D is the discharge coefficient, taken equal to 0.76 in [START_REF] Motsinger | Design and performance of duct acoustic treatment: aeroacoustics of flight vehicles[END_REF], and d is a constant value set equal to 1.6 in [START_REF] Laly | Acoustical modeling of micro-perforated panel at high sound pressure levels using equivalent fluid approach[END_REF]. Consequently, the coefficient C f to be used in the approximate nonlinear EFM (20) is

C f = d σ 0 ρ 0 (1 -ϕ 2 p ) πϕ p e p C 2 D , (34)
equal to 3.17 for BMP1 and to 2.93 for GE01 (values rounded to the hundredth).

In the aeronautical field, an incident SPL varying between 120 dB and 150 dB (which is equivalent to an acoustic velocity amplitude U between 0.06m.s -1 and 2.15m.s -1 ) can be characterized as a high sound level regime. An initial calculation of the Reynolds number (Re = ρ 0 U Λ/ϕη) associated with the flow within the pores of the perforated facesheets based on their geometric properties and the range of the studied SPL shows that Re is in the Forchheimer's regime. Indeed, a porous medium having a porosity in the order of 10% gives Re ≈ 4 and Re ≈ 146 at 120dB and 150dB, respectively, with ρ 0 = 1.226 kg • m 3 , η = 1.803 • 10 -5 kg • m -1 • s -1 and a characteristic length of order Λ ≈ 10 -4 m. Moreover, it has been observed that the resistivity of a porous medium subjected to a flow tends towards a maximum threshold value at very high Reynolds [START_REF] Barree | Beyond Factors: A Complete Model for Darcy, Forchheimer, and Trans-Forchheimer Flow in Porous Media[END_REF], suggesting that the Forchheimer law has a limited domain of validity at high Reynolds. [START_REF] Huang | Applicability of the Forchheimer Equation for Non-Darcy Flow in Porous Media[END_REF] evaluated a domain of validity of the Forchheimer regime for a pore-scale Reynolds number between Re ≈ 4 and Re ≈ 180. It is thus relevant to consider the nonlinear effects for the two plates given in Table 2 for high sound levels (120 ≥ SPL ≥ 150).

Description of the experiment

The experiments to be simulated are based on the impedance tube configuration modelled herein in 2D by the domains Ω 1 = [0, 0.1675] × [0, 0.012] and Ω 2 = [0, 0.1881] × [0, 0.012] whose meshes are presented in Figure 7. Note that all distance are given in meters, unless otherwise noted. 20) is solved. Furthermore, each liner contains a single closed cavity, and the cavity walls are represented by a rigid wall (except for the cavity surface in contact with the porous plate). A non-reflective wall condition is imposed on the left edge at x = 0. Finally, at the interface between the domains, the continuity conditions (35) are imposed.

ϕ air u air • n x = ϕ por u por • n x , (35a)
p air = p por , (35b)
In the study of [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF], the impedance is calculated for frequencies from 500 Hz to 5000 Hz, while Howerton et al. [2019] gives impedance measurements between 400 Hz and 3000 Hz. In order not to make two meshes each specific to a case, the meshes have been adapted to be sufficiently accurate at f = 5 kHz. To do this, the mesh size h is taken equal to 12 mm and the degree of the polynomials used in the DG solver is fixed at 3, giving more than 10 points per wavelength (λ min = c p /f max ≈ 50 mm). Moreover, the height of the duct being l y = 12 mm, the cut-off frequency is

f c = c 0 2l y ≈ 14.22 kHz (36)
Therefore, only the plane mode propagates in the tube as long as the sound signals has a frequency content lower than f c . [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF] conducted the impedance tube experiment with two types of sound source: a swept-sine signals and a white-noise signal. Depending on whether the source is monochromatic or broadband, the experimental results differ. In the study of [START_REF] Howerton | Application of Swept Sine Excitation for Acoustic Impedance Eduction[END_REF], only monochromatic signals are used. To compare the numerical results with their experimental results, the study is restricted to harmonic waves, thus having only one frequency in their spectrum. These waves are generated on the left edge of the domain (at x = 0) by the pressure and velocity fields:

p(t, (x = 0, y)) = A + sin(2πf t), (37a)
u(t, (x = 0, y)) = p(t, (x = 0, y))/(ρ 0 c 0 ), (37b)

with A + the amplitude of the incident wave adjusted to to obtain the target SPL at the liner surface and f the frequency of the generated harmonic wave. The pressure field and velocity field are taken to be zero in the entire domain at t=0.

In order to verify the correct modeling of the acoustic behavior of the liners with the approach based on a porous medium model, simulations will also be carried out in the linear regime with a broadband sound signal similar to the Mexican hat wavelet (30) :

p(0, (x = 0, y)) = A + p mex (t) (t < 10 -3 ms), (38a)
u(0, (x = 0, y)) = A + ρ 0 c 0 p mex (t) (t < 10 -3 ms), (38b)
with

p mex (t) = 1 - c 2 0 W 2 0 (t -t 0 ) 2 exp - 1 2 c 2 0 W 2 0 (t -t 0 ) 2 , ( 39
)
and null incoming flux conditions p(t, (x = 0, y)) and u(t, (x = 0, y)) for t > 10 -3 ms. The parameter values are defined so that the frequency content of the sound signal ranges from 500Hz to 5kHz: t 0 = 0.5 • 10 -3 s and W 0 = 3.35 • 10 -2 m (c 0 being the speed of sound in air).

Multipole approximations of dynamic variables

The JCAL model is used to express the transfer functions α and β over all frequencies of interest. The multipole approximations of these functions are obtained with the VF algorithm, applied on the frequency intervals I 1 =[500Hz, 5000Hz] for the BMP1 material, and I 2 =[400Hz, 3000Hz] for the GE01 material. Moreover the VF is only applied on the diffusive part of α and β, i.e. α/α ∞ -M/s for the dynamic tortuosity and β -1 for the dynamic compressibility. The number of parameters is chosen to have a relative error between the JCAL model and the MM under 1%. Therefore, the multipole approximations for the BMP1 case and the GE01 case are defined with 3 parameter pairs each. The values of the coefficients are gathered in Tables 4 and5.

Numerical results compared with experimental data

The two materials BMP1 and GE01 are modelled numerically by a 2D domain detailed in Section 7.1.2, with a description of the acoustic behavior of the perforated plates using the multipole approximations presented in Section 7.1.3. Before studying how they behave at high sound levels, a first study in Section 7.2.1 is carried out at low sound levels on the BMP1 liner, i.e. in the linear regime with a coefficient C f = 0. Following this, a Section 7.2.2 presents the results obtained when nonlinear effects are taken into account for materials under monochromatic sound signal.

Linear regime

At low noise levels, the equations solved in the perforated plate are those of the approximate linear EFM, i.e. system (9) with C f = 0. In the normal-incidence case under the assumption than only plan waves propagate, the impedance of the studied liners can be calculated with the theoretical impedance:

z(ω) := z eq (ω) ϕ p -jϕ p z cav (ω) cot(e p k eq (ω)) + z eq (ω) ϕ p z cav (ω) -j cot(e p k eq (ω))z eq (ω) , ( 40
)
with z cav (ω) = -j cot(l c ω/c 0 ) the cavity impedance, z eq (ω) = α(ω)/β(ω) the equivalent characteristic impedance and k eq (ω) = ω α(ω) β(ω) the equivalent wavenumber. It is thus possible to compare the impedance obtained numerically with the theoretical (linear) impedance model ( 40).

Numerical calculation of the impedance of the BMP1 liner is done thanks to simulations carried out with a broadband signal (38) where the amplitude A + is normalized, equal to 1. The choice of this amplitude does not change the numerical impedance due to the fact that the simulations are done in the linear regime (with the approximate linear EFM). Simulations are carried out over a relatively long time (t = 20ms), compared to the time interval over which the sound signal is generated (1ms). The velocity u and the pressure p are extracted at each iteration at the surface of the liner, in order to be able to compute the impedance

Z(x) := p(x) u(x) , (41a)
z(x) := 1 ρ 0 c 0 Z(x). ( 41b
)
after the use of the fast Fourier transform (FFT) on u and p. The theoretical impedance (40) can then be compared with the numerically calculated impedance, both presented in Figure 8. In Figure 8, the impedance from the numerical simulation is visibly close to the theoretical impedance with a relative error of the order of 1% at intermediate frequencies, rising to about 10% at the extremes of the interval [500Hz,5kHz]. The increasing discrepancy between the two curves at HF is mainly due to the accuracy of the multipole approximations. The addition of a pair of parameters (r k , s k ) in the VF approximation significantly reduces the relative error. It should be noted that one of the resonant frequencies f r of the BMP1 liner is around 2500Hz, visible on Figure 8b where f r is the value for which ℑ(Z) cancels.

The previous experiment is simulated again, but this time using monochromatic signals of the type (37) as the source term generated on the left edge of the domain. Thus, about fifteen simulations are launched, each with a different harmonic wave of frequency within [500Hz, 5000Hz] and amplitude A + = 1 (i.e. ∼ 91dB). The final time of each simulation is again taken equal to 20ms to have a sufficiently large steady state time window on which to apply the FFT. The impedance obtained by simulation is presented in Figure 9 compared with the impedance measured experimentally by Billard for the BMP1 material. [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF]. Vertical black dashed lines delimit the frequency range of interest.

A first observation of Figure 9 shows that with harmonic sound signals, the theoretical linear impedance is again well approximated. Nevertheless, we can see that the experimental results of [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF] obtained with harmonic signals at 130dB and 150dB are very different from those predicted by the theoretical linear impedance (40) and from the numerical results calculated with the approximate linear EFM. Therefore, nonlinear effects are taken into account (C f ̸ = 0) to see if the model proposed in Section 4 can predict the experimental results.

Nonlinear regime

In this section, the equations solved in the porous medium are those of the system (20) with the nonlinear coefficient C f defined by equation ( 34). The numerical simulations to be carried out are identical to those presented in the previous section 7.2.1, but with harmonic waves of amplitudes A + to be determined.

In the experiments of [START_REF] Howerton | Application of Swept Sine Excitation for Acoustic Impedance Eduction[END_REF] and [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF], the sound level is measured at the surface of the liner, which implies that the experimental results given at a certain SPL are not necessarily obtained with an incident pressure of the same SPL. It is therefore necessary to first obtain the amplitudes of the incident waves that produce the desired sound levels at the surface of the perforated plate. To do so, a first approximation of the amplitude can be obtained by using the reflection coefficient R := p -/p + , where p + is the pressure of amplitude A + of the wave propagating to the right and p -is the pressure of amplitude A -of the wave propagating to the left. It is known that the pressure at a point x p on the surface of the liner (the origin of the axis is taken at the surface of the liner: x p = 0) is directly given by the sum of the incident and reflected amplitudes.

p(t, x = x p ) = A + + A - = A + + A + R = A + (1 + R)
Hence, using the order of magnitude of the expected impedance z = (1+R)/(1-R) (e.g. via experimental data) and the target SPL gives an approximation of the amplitudes A + of the incident waves to use. An iterative adjustment by simulation is then made to refine the value of the incident sound wave amplitudes to converge towards the target SPL.

For the simulations associated with the BMP1 liner case, the incident wave SPL, denoted SPL i , and the SPL extracted from the numerical simulations at the liner surface, denoted SPL w , are given in Figure 10. The values of the incident wave amplitudes used to obtain the SPLs in Figure 10 are given in Table 6. Thus, in order to keep a constant SPL w at the surface of the studied liner, an important modification of the sound level of the incident wave is necessary. Moreover, we recall that f r ≈ 2500 Hz is one of the resonant frequencies of BMP1 in the vicinity of which the attenuation of the sound level is the highest (over the studied frequency range). This attenuation is to be correlated with the strong SPL i that one must have in the vicinity of f r to reach the target SPL w . On the other hand, as soon as the frequency moves away from the resonant frequency, the SPL i to be imposed decreases, becoming lower than the SPL w , up to a difference close to 6dB for frequencies where the liner reflection is almost total. This 6db difference was expected because the acoustic behaviour of a liner far from its resonance frequency is similar to that of a rigid wall where A + = A -.

Once the source terms to impose are known for each frequency, simulations can be performed by solving the LEE in the tube and cavity, and solving the approximate nonlinear EFM in the perforated plate. These simulations led to the numerical impedance of the BMP1 liner shown in Figure 11 using the pressure field and velocity field extracted on the surface of the material as performed in Section 7.2. The absorption coefficient computed from the impedance is displayed in Figure 12. One can see that the impedance pattern measured experimentally by Billard is well captured by the present numerical approach. The amplitude of the resistance in the vicinity of the resonance is very similar between the experimental and the numerical results, but a slight deviation of the frequency at which the maximum is reached can be observed. In this study, the static tortuosity α ∞ , which is one of the physical properties impacting the value of the resonant frequency, has a the constant value. Therefore, the discrepancy might come from the dependence of the static tortuosity on the velocity amplitude [START_REF] Mcintosh | Nonlinear wave propagation through rigid porous materials. I: Nonlinear parametrization and numerical solutions[END_REF]Lambert, 1990, Maa, 1998] omitted in this work. Furthermore, a larger difference between the experimental and numerical results is visible when moving away from the resonant frequency, in particular at LFs at 130dB. This may be due to a physical phenomenon not taken into account in the nonlinear EFM, or to a bias in the experimental results. Despite this, the results of Figures 11 and12 remain very satisfactory. It should be noted that the numerical simulations are not carried out at all frequencies where experimental data are available.

An additional point to highlight is the comparison of the frequency content of the reflected sound wave with that of the monochromatic incident wave p + . Figure 13 presents a harmonic analysis for two monochromatic waves, one at f 1 = 2000Hz and the other at f 2 = 3000Hz, with and without nonlinear effects taken into account in the EFM.

(ω)| p -(C f = 0) p -(C f = 0) p + p -(C f = 0) p -(C f = 0) p + (a) Incident wave of frequency f 1 = 2000Hz.
(ω)| p -(C f = 0) p -(C f = 0) p + p -(C f = 0) p -(C f = 0) p + (b) Incident wave of frequency f 2 = 3000Hz.
Figure 13: Frequency content of the incident pressure p+ and the reflected pressure p-for two simulations at 150dB (with and without consideration of nonlinear effects). A vertical black dashed line delimits the maximum frequency of the frequency range of interest.

Figure 13 shows that the reflected wave contains additional harmonics in the nonlinear regime, modifying the spectrum of the total pressure at the liner surface. The increase in the modulus of the pressure for frequencies above 5000Hz, and in particular the presence of a small peak around 11.25kHz in the linear regime (which is also found in the nonlinear regime) is most probably related to mesh which is no longer fine enough for this frequency range, but which nevertheless makes it possible to discern the additional harmonics.

Thus, in the case where the sound source has a broadband frequency content, the waves generated at a certain frequency will produce harmonics at higher frequencies and change the sound amplitude distribution. This can lead to a greater or lesser contribution of nonlinear effects for a broadband incident wave compared to a monochromatic incident signal, depending on the frequency and phase of the generated harmonics. In addition, the impedance will be different depending on the frequency spectrum distribution of the incident sound signal, due to a superposition of the harmonics generated by the LF waves on the HF waves.

All the steps taken to obtain the impedance of the BMP1 material are repeated for the GE01 material. The simulations with harmonic signals lead to the impedance presented in Figure 14, calculated for several frequencies. The real and imaginary parts of the impedance obtained experimentally by [START_REF] Howerton | Application of Swept Sine Excitation for Acoustic Impedance Eduction[END_REF] are given at 120dB and 140dB and compared to the impedance calculated numerically with the approximate nonlinear EFM. The theoretical linear impedance (40) is displayed, showing very little change in reactance (imaginary part of the impedance) at high SPL as expected. The resistance obtained from the simulations again shows a maximum near the resonant frequency (f r ∈ [1650Hz, 1700Hz]), however the experimental data appears to show a peak in the resonance at a frequency slightly below f r . The maximum resistance amplitude is more or less in agreement with the experimental measurements, with a more noticeable deviation at 140dB. Despite this, we again observe a similar overall impedance pattern between the experimental measurements and the numerical calculations.

In conclusion, the approximate nonlinear EFM correctly represents the acoustic behavior of liners subjected to a high SPL, when studied with a monochromatic sound source. The increase in resonance resistance is well reproduced with the value of the nonlinear coefficient (34). A discrepancy on the frequency at which the resistance peak occurs is however visible, and may be due to the modeling of the static tortuosity which has been taken here independent of the sound level.

Simulations in grazing incidence tube

In practice, acoustic liners in aircraft engines used for noise reduction are subject to grazing flow. The objective of this section is to model the acoustic behavior of acoustic liner studied by [START_REF] Michael | Evaluation of Spanwise Variable Impedance Liners with Three-Dimensional Aeroacoustics Propagation Codes[END_REF] and [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF] in a grazing incidence configuration, however without mean flow. To do so, a porous-based modeling of the perforated plate is used with a similar approach to that of Section 7, leading to a volume modeling of the liners (where cavities are meshed). A first Section 8.1 presents the materials and the numerical modeling of the experiment, followed by a Section 8.2 presenting the numerical results.

Numerical modeling of the experiment

The experiments to be numerically reproduced are based on those presented in [Jones et al., 2017, Sec. II] and [Diab et al., 2022, Sec. 5]. In these experiments, the studied materials are SDOF acoustic liners. Sound waves propagates with a grazing incidence along the liners, as schematized in Figure 15. In this study, monochromatic signals are generated upstream of the duct (x = 0m) and propagate mainly from the left to the right. Experimental pressure measurements made on the opposite edge of the material y w ) make possible to study their acoustic properties, by looking for example at the attenuation of the sound level downstream of the duct (x = 0.812m) for different incident wave frequencies. In order to perform the simulations reproducing the grazing-flow impedance-tube experiments, details about the liners and the duct dimensions are given in Section 8.1.1. Then, a porous-based approximation of the acoustic behavior of the two studied materials is performed in Section 8.1.2, followed in Section 8.1.3 by the numerical modeling of the domain representing the configuration in Figure 15.

Description of the liners and the configuration

The geometric properties of the liners selected from the literature are given in table 7. The liner from [Diab et al., 2022, Tab. 1] is named DPP1, referring to the PP#1 absorber in their article, and the second liner is the SVC1 from [Jones et al., 2017, Sec. II.B].

Table 7: Geometric properties of the DPP1 and SVC1 liners studied in [START_REF] Michael | Evaluation of Spanwise Variable Impedance Liners with Three-Dimensional Aeroacoustics Propagation Codes[END_REF][START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF] In [START_REF] Michael | Evaluation of Spanwise Variable Impedance Liners with Three-Dimensional Aeroacoustics Propagation Codes[END_REF], the 50.8 mm-wide and 406.4 mm-long SVC1 liner consists of a 4 × 30 array of separated square chambers with 22 perforations, each of size 15.748 mm × 15.748 mm, and the cavities of l c = 76.2 mm in depth. Note that the porosity of its perforated plate directly calculated from the geometry of the liner is

ϕ p = ( 4 × 30 × 22 × πr 2 p ) / (50.8 × 406.4 × 10 -6 ), ≈ 10.892 %, (42)
which is slightly different from the porosity given in the article [START_REF] Michael | Evaluation of Spanwise Variable Impedance Liners with Three-Dimensional Aeroacoustics Propagation Codes[END_REF]: ϕ p = 11%. In a first step, the porosity in Table 7 is used.

Porous-based modeling of the liners

As in Section 7.1, a porous-based description of the SVC1 and DPP1 perforated plate is adopted using the relations (31) for Λ, Λ ′ , k 0 , k ′ 0 and (32) for the static tortuosity. Moreover, the nonlinear coefficient, defined by equation ( 34), is equal to 2.32 for the SVC1 liner and 0.47 for the DPP1 liner.

The VF algorithm is used to compute the multipole approximation of α and β described by the JCAL model for the SVC1 and DPP1 liners,respectively,on [400 Hz,3000 Hz] and [200 Hz,3000 Hz]. The VF is applied on the diffusive part of α and β, leading to the 3-parameter multipole approximations detailed in Table 8 and9. Knowing the equivalent porous media describing the acoustics of the DPP1 and SVC1 liners and their thickness e p given in Table 7, a mesh of the domain can be built. Contrary to the normal incidence configuration, the locally reacting property of the acoustic liners must now be taken into account. To do this, two different meshes have been made and are presented in Figure 16 for the SVC1 liner (similar meshes were done for the DPP1 liner). These two meshes are of the same dimensions, both composed of a blue mesh duct with a thin porous medium in brown/red positioned on the low edge and backed by cavities in green. The latter are meshed in two different ways: the cavity zone is meshed in one block in the regular SVC1 reg mesh, while in the SVC1 cav mesh, 20 cavities are meshed and separated by rigid walls with no thickness. Thus, in this second case, the mesh forces the wave propagation in the cavities to be along y. For the first SVC1 reg mesh, the propagation along y in the cavities is imposed by cancelling the contribution in x in the LEE. To do so, only the y component of the waves at the interfaces of the mesh cells in the green cavity is taken into account in the DG solver. Note that the mesh is slightly refined on the upstream and downstream front of the liner to improve the modeling of the impedance discontinuity zones.

In both studies Jones et al. [2017], [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF], the impedance is calculated for frequencies up to 3000Hz. The meshes have been adapted to be sufficiently accurate at f = 3kHz. To do so, the mesh size h is taken to be at most 3cm and the degree of the polynomials used in the DG solver is set to 2, giving more than 10 points per wavelength. The time step is ∆t = 5 • 10 -7 , , selected as a function of the smallest cell sizes of the mesh and the multipole approximations. The total time simulation is at least 0.03 s.

Numerical results

This section is split in two parts: Section 8.2.1 focuses on the SVC1 liner with a sensitivity analysis of physical parameters, and Section 8.2.2 reproducing the numerical simulations done in [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF] for the DPP1 liner.

SVC1 liner

Simulations are done on both meshes (with cavities or with a forcing of the locally-reacting property in the equations). More specifically, three simulations at different frequencies are performed with a SPL maintained at 140dB upstream of the duct (x = 0). To do so, an incident acoustic wave is generated at x = 0 with the appropriate amplitude: A + = 278 (resp. 268) at 600 Hz, A + = 180 (resp. 180) at 1000 Hz and A + = 335 (resp. 332) at 1400 Hz for the regular mesh SVC1 reg (resp. cavity mesh SVC1 cav ). The pressure extracted at the duct wall (y w ) opposite to the liner is compared to the measured SPL from [START_REF] Michael | Evaluation of Spanwise Variable Impedance Liners with Three-Dimensional Aeroacoustics Propagation Codes[END_REF] in Figure 17. The main acoustic behaviour seen on the experimental study is numerically reproduced in spite of visible deviations, especially at frequencies far from the resonant frequency (∼1000 Hz). Moreover, thanks to a non-reflecting outlet condition at the end of the duct, there is no oscillations of the SPL at its end in the simulations, compared to the experiments which could cause the difference in results, especially at low frequencies (600Hz) where the liner's attenuation is low, and where the anechoïcity is more difficult to achieve experimentally. Furthermore, a full homogeneous approximation of the liner may neglect some physical effects than can occur between two adjacent cavities that are separated, as for the SVC1 liner. Finally in terms en meshes, both lead to very similar results for this liner. Among the three simulations, the noise attenuation is maximal at 1000 Hz, which is close to a resonant frequency of the SVC1 liner. At this frequency, the numerical simulation gives a SPL at the end of the tube equal to 85dB instead of 75dB according to the experimental data. This discrepancy may be due to the definition of the equivalent porous medium not being precise enough (e.g. by neglecting some physical phenomena) or to the sensibility of the modeling at the resonance. To illustrate that, modified equivalent porous medium are defined with a ±5% deviation on the radius r p . It implies new multipole approximations of α and β, a new value of the porosity ϕ p calculated by using equation ( 42 The numerical simulations done with a radius 5% lower or 5% higher than r p = 0.5207 mm given in Table 7 shows a significant impact on the SPL at the end of the tube with a ±2 ∼3dB difference. Therefore, very small manufacturing defects on the liners can considerably change the acoustic properties of the liner, especially around the resonant frequency. Moreover, given the SPL upstream of the liner impedance discontinuity, we can surmise that the impedance realised in practice by the SVC1 liner is higher than that used in the simulation. This can be deduced from the amplitude of the SPL variation upstream of the liner, which is lower in the experimental case, indicating a lower impedance discontinuity, and thus a more resistive liner. Having a finite size cavity network in a SDOF sample imposes a cavity porosity lower than 1, which artificially increases the overall liner resistance. To model this effect, one would have to represent the cavity network as a porous material of equivalent properties, which can be done with the same formalism as the one used to represent a perforated plate [START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF].

DPP1 liner

Unlike the SVC1 liner, the DPP1 liner was studied using a numerical approach presented in [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF]. Consequently, the comparison is made with numerical results. The latter are obtained at 1600 Hz and for 3 different incident SPL inc (120dB, 130dB and 140dB). In Figure 19, the data from [Diab et al., 2022, Fig. 24] (labeled Diab) are compared to the numerical results (labeled Num) obtained with the present approach. Figure 19: Normalized SPL along the opposite wall to the DPP1 liner obtained by the present numerical approach and obtained by [Diab et al., 2022, Fig. 24] for an incident signal at 1600 Hz and different incident SPL. Vertical black dashed lines delimit the liner position.

Simulations based on the approximate nonlinear EFM show the same relative increase of the SPL at the end of the duct as the incident SPL increase. Moreover, the SPL along the duct at SPL inc = 140 dB obtained by the porous-based approach is very close to the numerical data of [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF]. However, at lower SPL inc , differences between the simulations appears with a more significant decrease of the SPL along the duct wall in the present simulations compared to Diab's results for SPL inc = 120 dB. It might come from the fact that in the approach presented herein, the static tortuosity α ∞ is taken independent of the velocity. A modification of α ∞ changes the resonant frequency, hence the SPL attenuation. As in Section 8.2.1, a slight difference around the resonant frequency can lead to important discrepancies. The latter are very small at 140 dB, which could come from the fact that at this SPL, the static tortuosity α ∞ has a similar values in both simulations, or it can simply be a particular case giving rise to very similar results.

Conclusion

A nonlinear EFM based on diffusive variables α and β described by a generic model and the Forchheimer's correction was developed. Following [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF], a proof of the stability of the proposed nonlinear model was given as soon as the nonlinear coefficient C f is positive. Then, an approximate nonlinear model was built using multipole approximation, for which sufficient conditions on the multipole parameters are given to keep a stable approximate model.

In the DG solver, two approaches were used to discretize the additional nonlinear term coming from the Forchheimer's correction. Both approaches led to similar results, showing a loss of precision of the DG solver, which can be compensated by a fine mesh and a sufficiently high DG order. However, it would be interesting to investigate more precisely how to improve the approximation of this term in the numerical scheme.

The approximate nonlinear model was used to reproduce the experiments done by [START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF] and [START_REF] Howerton | Application of Swept Sine Excitation for Acoustic Impedance Eduction[END_REF] in normal incidence, showing numerical results very close to the experimental data, especially in terms of resistance amplitude. It highlighted that the Forchheimer's correction mainly accounts for the nonlinear phenomena at high SPL. Small discrepancies were observed on the frequency for which the resistivity peak is reached. This could come from a static tortuosity assumed to be independent to the velocity in the present work. It would be interesting to see how to account for the velocity correction in the proposed nonlinear EFM.

A last section focused on grazing incidence tube with two liners studied in [START_REF] Michael | Evaluation of Spanwise Variable Impedance Liners with Three-Dimensional Aeroacoustics Propagation Codes[END_REF][START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF]. The approximate nonlinear EFM gave an attenuation of the SPL along the acoustic duct with similar variations to the experimental and numerical data available in the literature. However, some visible differences could be seen when comparing the results. The sensitivity of the SVC1 acoustic behavior we tried to reproduce numerically at the resonant frequency could explain the differences for the SVC1 liner, and a different definition of the static tortuosity in [START_REF] Diab | Nonlinear broadband time-domain admittance boundary condition for duct acoustics. Application to perforated plate liners[END_REF] might be the main reason of the differences in SPL decrease, because it changes the resonant frequency of the material, hence the impedance at the chosen frequency. This last point about the grazing incidence tube experiments requires further study for a better understanding.

A theoretical study at low velocity magnitude, where this type of correction is necessary, leads to the same theoretical results given in Theorem 3.1. It can even be extended to a more general case, where the total resistivity is defined as σ = σ 0 1 + Φ ∥u∥ .

(47)

with Φ : R + → R + of C 1 (e.g. quadratic at low velocity amplitude and linear at high velocity amplitude as in Figure 20). Starting from the nonlinear EFM corrected with (47), the method used to obtain the energy balance variations ( 16) leads in this case to

d dt E nl (t) = d dt E(t) -ρ 0 α ∞ M Ω Φ ∥u∥ ∥u∥ 2 dx, (48)
allowing to conclude on the stability of the more general nonlinear model, since Φ > 0.

x Φ(x) Φ(0)

Appendix C Numerical scheme

The code used in this article is based on the one used in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] which solves the linearized Euler equations and the linear EFM in 2D, using a Discontinuous Galerkin method detailed in Section C.1 and a Runge-Kutta method detailed in Section C.2

C.1 Space discretization

Given a domain Ω ⊂ R 2 for which there is a partition (T h ) where h denotes the maximum diameter of the elements in the partition. The approximation space is taken as V h := {v, | ∀T ∈ T h , v |T ∈ P k (T )} where P k (T ) is the space of polynomials of degree at most k. A basis {λ i j ∈ P k (T i ), j = 1}, with d = (k + 1)(k + 2)/2, is defined for each element T i ∈ T h . Thus, a scalar function q is approximated on a cell T i by q h (t, x) := d j=1 q i,j h (t)λ i j (x).

(49)

The solution vector q h (t, x) = (u h v h p h ) is then defined based on (49). Applying the DG method [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF], it is possible to write: 0 = d dt Ti q h (t, x) λ i j dΩ + Ti F(q h (t, x)) • ∇λ i j dΩ -∂Ti F * (q e h (t, x), q i h (t, x)) • n i , λ i j dσ + Ti b(q h (t, x)) λ i j dΩ,

where F(q h ) = (A x q h , A y q h ), b(q h ) = Bq h , n i = (n i x , n i y ) is the outgoing unit normal vector with respect to the edge ∂ t T i and F * is the numerical flux. The solution q h on the edge of a cell T i is denoted q i h or q e h depending on whether it is, respectively, the internal or the external value of T i that is taken into account. In the code, the numerical flux used is the flux vector splitting F * (q i h , q e h ) = A + q i h + A -q e h ,

where the incoming and outgoing flows are separated into two, respectively associated with A + and A -. The latter are matrices containing respectively the positive and negative eigenvalues of A = A x n x + A y n y . Note that this flow

Figure 2 :

 2 Figure2: Solutions of the nonlinear system (28) at different times, for the two initial conditions (29) and (30), with W0 = 8 • 10 -2 and x0 = 0.5.

Figure 3 :

 3 Figure 3: Rectangular mesh [0,3]×[0,0.05] of 248 triangles.

Figure 4 :

 4 Figure 4: Relative mean square error of the numerical solutions at t = 2 for the initial conditions (29) and (30), with W0 = 8 • 10 -2 , x0 = 0.5 and different values of the nonlinear coefficient C f .

  Theoretical and numerical solutions with theoretical slopes at sign change points.

Figure 5 :

 5 Figure 5: Solutions at t = 2 for the initial condition (30), with W0 = 8 • 10 -2 , x0 = 0.5 and C f = 20.

  Geometry of a single-degree-of-freedom acoustic liner. Perforated plate modeled by a porous medium.

Figure 6 :

 6 Figure 6: Modeling of a perforated plate in a liner by a porous medium.

Figure 7 :

 7 Figure 7: Meshes associated with the impedance tube for liners BMP1 and GE01 (see Table2) with the blue domain (left) representing the air domain, the green domain (right) the cavity also containing air and the brown domain (middle) associated with the porous medium.

Figure 8 :

 8 Figure 8: Theoretical and numerical surface impedance of the BMP1 material with a broadband sound source. Vertical black dashed lines delimit the frequency range of interest.

Figure 9 :

 9 Figure9: Real part of the theoretical surface impedance of the BMP1 material, calculated numerically with monochromatic signals and measured experimentally by[START_REF] Billard | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range[END_REF]. Vertical black dashed lines delimit the frequency range of interest.

Figure 10 :

 10 Figure10: Sound levels at the surface of the BMP1 liner as a function of the incident wave frequency, calculated from numerical simulations with monochromatic signals.

Figure 11 :Figure 12 :

 1112 Figure11: Real part of the impedance of the BMP1 material given by the theoretical linear impedance, and calculated from experimental and numerical results at high SPL. Vertical black dashed lines delimit the frequency range of interest.

Figure 14 :

 14 Figure14: Impedance of the GE01 material given by the theoretical linear model (in blue line), and calculated from experimental (green diamond and purple square) and numerical results (red dots) at high SPL. Vertical black dashed lines delimit the frequency range of interest.

Figure 15 :

 15 Figure 15: Configuration of the grazing flow impedance tube with an acoustic liner positioned on its bottom edge. Dimensions are given in meters.

Figure 16 :

 16 Figure 16: SVC meshes of the studied domain representing the grazing-flow impedance tube with the liner located on the low edge of the duct.

Figure 17 :

 17 Figure 17: SPL along the opposite wall to the SVC1 liner for different sound source frequencies. Vertical black dashed lines delimit the liner position.

Figure 18 :

 18 Figure18: SPL along the opposite wall to the SVC1 liner and modified SVC1 liners with an incident signal at 1000 Hz. Vertical black dashed lines delimit the liner position.

Figure 20 :

 20 Figure 20: Example of a function Φ : R + → R + of C 1 , having a quadratic behavior for small x and linear for large x.

Table 1 :

 1 Model parameters of the generic model from

Table 3 :

 3 Values of the intrinsic parameters of the porous media modeling the perforated plates given in Table2(rounded to the hundredth).

	Paramètres

Table 4 :

 4 Parameters of the multipole approximations of α and β described by the JCAL model for the BMP1 liner and computed with the VF algorithm.

		c ∞	c -1	r 1	r 2	r 3	s 1	s 2	s 3
	α mm 1.58 116.32 275.07 511.66 3 494.56 -1 493.83 -9 864.41 -91 761.40
	β mm	1	0	137.43 232.39 1 603.64	-973.34	-9 102.34 -86 455.56

Table 5 :

 5 Parameters of the multipole approximations of α and β described by the JCAL model for the GE01 liner and computed with the VF algorithm.

		c ∞	c -1	r 1	r 2	r 3	s 1	s 2	s 3
	α mm 1.89 249.09 330.04 667.19 4 762.25 -1 388.11 -7 562.82 -66 587.35
	β mm	1	0	199.86 281.72 2 211.19	-762.40	-7 288.13 -65 096.33

Table 6 :

 6 Amplitudes in Pa of the incident harmonic waves as a function of the frequency for a target sound level SPLw at the surface of the BMP1 liner at 130dB and 150dB.

	Frequency (Hz) 500	900 1300 1700 2000 2200 2400 2500
	SPL w = 130dB	45	47	52	65	98	156	224	225
	SPL w = 150dB 454	476	539	713	925 1034 1079 1078
	Frequency (Hz) 2600 2800 3000 3400 3800 4200 4600 5000
	SPL w = 130dB 200	136	99	70	60	55	52	50
	SPL w = 150dB 1063 1002 919	753	641	577	539	515

.

  Linerϕ p (%) r p (mm) e p (mm) l c (mm) x 1 (mm) x 2 (mm) y w (mm)

	SVC1	11	0.5207	0.8636	76.2	203.2	609.6	63.5
	DPP1	1.5	0.25	1.0	10	203	609	50.8

Table 8 :

 8 Parameters of the multipole approximations of α and β described by the JCAL model for the SVC1 liner and computed with the VF algorithm.

		c ∞	c -1	r 1	r 2	r 3	s 1	s 2	s 3
	α mm 1.32 328.98 279.60 642.78 4 761.05 -1 749.40 -8 299.13 -71 934.73
	β mm	1	0	189.33 272.23 2 110.28	-744.24	-7 175.68 -64 260.67

Table 9 :

 9 Parameters of the multipole approximations of α and β described by the JCAL model for the DPP1 liner and computed with the VF algorithm.

		c ∞	c -1	r 1	r 2	r 3	s 1	s 2	s 3
	α mm 1.18 1 591.22 352.94 1 371.20	12 50	-4 381.59 -12 373.69 -110 504.85
	β mm	1	0	634.33	490.80	5 496.65 -2 059.66 -12 190.22 -101 169.15
	8.1.3 Numerical configuration

As opposed to the Biot model[Biot, 1956a,b] used for example in biomechanics and geophysics for poroelastic materials where the solid phase is in motion.

Appendix A Energy analysis

The energy balance of the EFM with the generic model defined in [START_REF] Moufid | Energy analysis and discretization of the timedomain equivalent fluid model for wave propagation in rigid porous media[END_REF] is briefly introduced in this Appendix. First, the classical mechanical energy is defined and composed of the kinetic energy and the potential energy:

Then, an energy is defined for the diffusive variables:

where

Note that the energy E ψ0 is similar to E ψ but for a fixed pole ξ = -s 0 with a positive residue r 0 , and all the model parameters (M , N , L, M , N ′ , L ′ ) are positive. The energy balance:

is computed based on the definitions ( 43) and ( 44), and its derivative is found to be

when there is no contribution at the boundary ∂Ω (i.e., p = 0, or u • n = 0 on ∂Ω) and the functions are defined in the appropriate functional spaces [START_REF] Haddar | Theoretical and numerical analysis of the Webster Lokshin model[END_REF]. For the nonlinear EFM, we make use of the same energy functional E nl = E lin . The only difference appears in the computation of the derivatives, where an additional negative term can be found for d dt E nl (see equation ( 16)).

Appendix B General resistivity correction

In section 3, the dependence of the resistivity on the velocity amplitude is taken to be linear, in agreement with the Forchheimer's correction. However, the resistivity dependence is quadratic at low Reynolds, leading to

is an exact solution of the one-dimensional Riemann problem with constant coefficients. At the boundary of the domain ∂Ω, the imposed flux is the centered flux

where q e h represents a ghost state defined with q i h . Finally, the semi-discrete equation reads:

with q h the unknown and S the source term.

C.2 Time discretization

The inversion of the mass matrix M in is direct thanks to its diagonal structure per block. Indeed, each block is of the order of a few tens of rows, thus small enough to calculate the inverse of each small block by a direct method. Thus, the spatial discretization ( 53) can be rewritten

with L h the semi-discrete operator, E = M -1 K and G = M -1 S. The Runge-Kutta method used for the spatial discretization is the fourth-order RKF84 method given by [START_REF] Toulorge | Optimal Runge-Kutta schemes for discontinuous Galerkin space discretizations applied to wave propagation problems[END_REF] and which was notably used for the results given in [Monteghetti et al., 2018]. This method is shown to be very efficient when used in combination with a DG method for wave propagation problems in the works of [START_REF] Toulorge | Optimal Runge-Kutta schemes for discontinuous Galerkin space discretizations applied to wave propagation problems[END_REF].

Let {t n } N n=0 be a partition of [0, T ] ⊂ R + , ∆t = t n+1t n be the time step and q n h be the approximate solution at time t n . The steps of the RKF84 algorithm are

dq (i) = A i dq (i-1) + L h t n + c i ∆t, q (i-1) (55b) q (i) = q (i-1) + B i dq (i) , for i = 1 . . . 8, (55c)

where the coefficients A i , B i and c i are given by Toulorge and Desmet [2012, Table A.9].