Benoit Nennig

Emmanuel Perrey-Debain

Mabrouk Ben Tahar

A mode matching method for modelling dissipative silencers lined with poroelastic materials and containing mean flow

Keywords: numbers: 43.20.Mv, 43.50.Gf, 43.28.Py

A mode matching method for predicting the transmission loss of a cylindrical shaped dissipative silencer partially filled with a poroelastic foam is developed. The model takes into account the solid phase elasticity of the sound absorbing material, the mounting conditions of the foam and the presence of an uniform mean flow in the central airway. The novelty of the proposed approach lies in the fact that guided modes of the silencer have a composite nature containing both compressional and shear waves as opposed to classical mode matching methods in which only acoustic pressure waves are present. Results presented demonstrate good agreement with finite element calculations provided a sufficient number of modes are retained. In practice, it is found that the time for computing the transmission loss over a large frequency range takes a few minutes on a personal computer. This makes the present method a reliable tool for tackling dissipative silencers lined with poroelastic materials.

I. INTRODUCTION

Dissipative silencers are becoming indispensable elements of modern exhaust systems. In many cases of practical interest, they consist of an expansion chamber filled with an absorbent material that dissipates the acoustic energy into heat. Their frequent use in the automotive industry, in heating, ventilation and air-conditioning (HVAC) applications, or in the aircraft industry necessitates the development of fast and accurate numerical methods in order to analyze the effectiveness of these devices under different operating conditions such as the possible presence of mean flow and the frequency range of interest. There are now numerous models for computing sound attenuation through dissipative silencers each having advantages and drawbacks according to the configuration in hand. Although this is not the place for a complete survey, we can refer to a very recent paper of Kirby [START_REF] Kirby | A comparison between analytic and numerical methods for modelling automotive dissipative silencers with mean flow[END_REF] and the references therein. All these techniques have in common that they have been developed with the assumptions that the lining is either locally reacting, represented by its normal impedance, or bulk reacting. In this latter, the absorbing material is assumed to have either zero or infinite stiffness, only a compression wave is allowed to propagate and the material is considered to be an effective fluid characterized by a complex wave number and density [START_REF] Beranek | Acoustical properties of homogeneous, isotropic rigid tiles and flexible blankets[END_REF][START_REF] Delany | Acoustical properties of fibrous absorbent materials[END_REF] . Fibrous media such as glass fiber or glass wool, for instance, are usually modelized as bulk reacting materials.

However, if the solid structure of the absorbent material has a finite stiffness which is the case for a wide range of poroelastic materials such as foams, it is known that three type of waves (two longitudinal and one transa) Electronic address: benoit.nennig@utc.fr verse) are allowed to propagate through the medium [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range. II. higher frequency range[END_REF] . This makes the behavior of foams more complicated and somewhat less predictable than that of glass fiber. One major reason for this is that the acoustical behavior of foam treatments having finite dimensions is known to be very sensitive to the boundary conditions that exits at its surfaces [START_REF] Song | Investigation of the vibrational modes of edge-constrained fibrous samples placed in a standing wave tube[END_REF][START_REF] Pilon | Behavioral criterion quantifying the edge-constrained effects on foams in the standing wave tube[END_REF] . A proper analysis of these foams is more demanding and the Biot theory of poroelasticity must then be considered [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range. II. higher frequency range[END_REF][START_REF] Allard | Propagation of Sound in Porous Media: Modeling Sound Absorbing Materials[END_REF] . In this context, there has been considerable effort into devising finite element formulations based on the Biot model and we can refer to Refs. 5, 8, 9 for the most recent advances in this matter. In a recent communication, Nennig et al. [START_REF] Nennig | A 3d numerical method for studying poroelastic liners with mean flow[END_REF] have presented a full 3D Finite Element (FE) model for solving the Biot's equations using the mixed displacement-pressure formulation of Atalla [START_REF] Atalla | Enhanced weak integral formulation for the mixed (u, p) poroelastic equations[END_REF] . The model was coupled with the Galbrun equations for the acoustic displacement in the airway in order to take into account the presence of a sheared mean flow in the central duct. Though the method is capable of modeling arbitrary shaped silencers, the 3D FE model is still limited by the computer resources and this can have a negative impact when, for instance, some efficient optimizations of the silencers (geometry, material properties and mounting conditions) are needed.

When dissipative silencers contain an axially invariant cross section with a uniform mean flow in the airway, it is judicious to take advantage of the separability of the wave equations. In this context, semi-analytical methods such as the very popular Mode Matching Methods (MMM) are usually considered as the ultimate tool. MMM have been successfully applied in many fields when sudden transitions are present in the waveguide. The discontinuity may be of geometric-type, due to a change in the medium properties or an abrupt change in boundary conditions. The method first requires computing the eigenmode basis in each segment of the wave guide. By expanding the wave field in this basis, a matching procedure must then Mode matching method for poroelastic silencers 1 be applied in order to respect the dynamic and kinematical continuity condition at the discontinuity interface between two consecutive segments. This method is attractive because it leads to an immediate reduction in the dimension of the problem and can be much quicker than a full FE model. This is why the method is now widely used within the engineering community.

In this work we propose a new MMM for dealing with expansion chambers partially filled with an poroelastic material in which two longitudinal and one transverse waves are present. One crucial step of the method is to find with sufficient accuracy the eigenvalues and the transverse eigenfunctions of the cylindrical chamber. When the silencer contains an arbitrary cross-section, the eigenmode basis is usually numerically computed via finite element eigenvalue analysis over the crosssection [START_REF] Kirby | Transmission loss predictions for dissipative silencers of arbitrary cross section in the presence of mean flow[END_REF][START_REF] Albelda | A transversal substructuring mode matching method applied to the acoustic analysis of dissipative mufflers[END_REF] . Besides the size of the associated FE matrix which is expected to grow at least quadratically with respect to a typical wavenumber; the advantage of such an approach is that the approximated eigenvalues are obtained via standard algebraic solvers. In the present paper, we shall restrict ourselves to silencers with circular cross-section. In this case, eigenmodes admit closed-form solutions in terms of Bessel functions of integer orders and the key issue is to then solve the resulting governing eigenvalue equation. This is not a trivial task though there are numerous available numerical techniques for this purpose [START_REF] Kirby | Analytic mode matching for a circular dissipative silencer containing mean flow and a perforated pipe[END_REF] . It is remarkable that these techniques largely employed in the field of duct acoustics have in common the requirement of initial approximations for the eigenvalues to start the algorithm. Generally the approach relies on finding asymptotic approximations that are exact or at least known with sufficient accuracy and then track the roots in the complex plane with respect to some parameter which is either the frequency, the Mach number of the gas flow or the admittance in the context of a locally reacting liner. The method can be time consuming and more importantly all these root finding techniques suffer from the problem of missing roots [START_REF] Kirby | A point collocation approach to modeling large dissipative silencers[END_REF] . The complicated frequency dependence of Biot's coefficients renders asymptotic approximations in the low frequency limit very tedious if not intractable. For this reason, we favored another root-finding technique which does not rely on any initial guess. The method based on the Argument Principle has been previously employed in a somewhat different context [START_REF] Brazier-Smith | On the determination of the roots of dispersion equations by the use of winding number integrals[END_REF][START_REF] Cristini | Calcul des zéros d'une fonction analytique avec points de branchements (Roots calculation of an analytic function with branch points)[END_REF][START_REF] Chen | Efficient and accurate numerical analysis of multilayer planar optical waveguides in lossy anisotropic media[END_REF] and this is presented in this paper.

This work begins by deriving the eigenequation for the modal axial wavenumbers. Each root of the equation corresponds to a mode of the dissipative silencer. This mode has a composite nature containing both pressure and displacement type waves as opposed to modes of bulk reacting lined duct where only acoustic pressure waves are present. This will be commented on with respect to their physical interpretation. In particular, we show that different families can be clearly distinguished. The matching procedure at the inlet and outlet planes of the silencer is then explained in detail. We finish the paper by presenting results showing the impact of elastic resonances of the foam on the Transmission Loss for silencers with typical dimensions encountered in the au-tomotive industry.

II. PROBLEM STATEMENT

A. Constitutive equations and continuity conditions

The silencer considered here consists of a cylindrical chamber duct of length L surrounded by a soundabsorbing poroelastic material (region ii in Figure 1). In the airway, a mean flow of speed U 0 is present. The inlet and outlet pipes (regions i and ii) are identical, each having a circular cross section with rigid walls. In the airflow domain (0 < r < r 1 ) the acoustic pressure p obeys the convected wave equation ∆p -

1 c 2 0 D 2 p Dt 2 = 0, (1)
where D Dt ≡ ∂ ∂t + U 0 ∂ ∂z stands for the material derivative along the mean flow. Here, t is time, c 0 denotes the sound speed and equations are written in the usual cylindrical coordinate system (r, θ, z). It is convenient for the analysis to introduce the particle perturbation displacement w within the fluid. Assuming an homogeneous flow of density ρ 0 , this is given via the linearized momentum equation as

ρ 0 D 2 w Dt 2 = -∇p. (2)
Though these relations are commonly used within the assumption of uniform flow [START_REF] Astley | A finite element scheme for attenuation in ducts lined with porous material: comparison with experiment[END_REF] , the displacement perturbation formalism remains valid for more general flow fields [START_REF] Galbrun | Propagation d'une onde sonore dans l'atmosphère terrestre et théorie des zones de[END_REF] .

In the expansion chamber, the wave propagation in the sound-absorbing media (r 1 < r < r 2 ) is described via Biot's model [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range. II. higher frequency range[END_REF] . This latter is grounded on the superposition of a fluid phase and a solid phase which are coupled together and respectively described by the fluid phase displacement U and the solid phase displacement u. For time-harmonic representation (e -iωt ), we have the following coupled system 4

∇ • σ s + ω 2 (ρ 11 u + ρ 12 U) = 0, ( 3
) ∇ • σ f + ω 2 (ρ 12 u + ρ 22 U) = 0. ( 4
)
Solid and fluid phase stress tensors are given by

σ s = (A ∇ • u + Q ∇ • U) I + 2N ε s , (5)
σ f = (Q ∇ • u + R ∇ • U) I, (6)
Mode matching method for poroelastic silencers 2

where ε s = 1/2(∇u + (∇u) T ) is the usual strain tensor and I is the identity matrix. The total stress tensor σ t is, by definition, the sum of σ f and σ s . Biot's coefficients A, N, Q, R are related to the material properties (given in Table IV) by the Allard-Johnson model. Their expressions can be found in the literature or in the reference textbook 7 as well as the other quantities introduced in this section. These quantities are all complex and frequency-dependent, A and N correspond to the Lamé coefficients, R is the effective bulk modulus of the fluid phase and Q indicates the coupling of the two phases volumic dilatation. The imaginary part of A and N includes the structural damping and, in Q and R this part includes the thermal dissipation. The imaginary parts of the effective density coefficients ρ 11 , ρ 22 and ρ 12 take into account viscous damping. At the air-porous interface (denoted by Γ in Figure 1), the following compatibility conditions must be verified

σ t n = -p n , (7a)
p p = p , ( 7b
) φU • n + (1 -φ)u • n = w • n . ( 7c
)
Here φ is the porosity and the pore pressure p p is obtained from the fluid phase tensor as -3φp p = trσ f . These conditions are classical in the non-flow case [START_REF] Debergue | Boundary conditions for the weak formulation of the mixed (u, p) poroelasticity problem[END_REF][START_REF] Kang | Finite element modeling of isotropic elastic porous materials coupled with acoustical finite elements[END_REF][START_REF] Lovera | Boundary conditions for a fluid-saturated porous solid[END_REF] though the last one is sometime expressed in terms of the normal volume velocity. The advantage of keeping the displacements as primitive variables is that (7c) remains valid when flow is present in the airway. On the hard surface (bold solid line Γ r ), the poroelastic layer is assumed to be clamped, that is the solid phase displacement as well as the normal fluid phase displacement are set to zero

u = 0 and U • n = 0. (8)
In conditions (7) and (8), it is understood that n denotes the normal unit vector at the boundary of the porous domain which can then be either along the axial or the radial direction. For the sake of illustration, we restricted our attention to the 'clamped' conditions (8) but other mounting configurations where, for instance, the foam is sliding may be treated by simply changing Eq. ( 8) accordingly [START_REF] Nennig | Influence of solid phase elasticity in poroelastic liners submitted to grazing flows[END_REF][START_REF] Kang | Sound propagation in circular ducts lined with noise control foams[END_REF] . The presence of a perforated sheet at the air-porous interface can also be taken into account.

In this case, the condition (7c) remains unchanged but (7b) and (7a) are modified with the addition of the pressure jump between both sides of the perforate plate.

B. Helmholtz decomposition

Because the previous wave equations are written in terms of the physical variables, these are not appropriate for a modal analysis. Instead, we shall rewrite each quantity in a potential form using the Helmholtz decomposition. In the porous media, both displacement fields are then written as

u = ∇ϕ + ∇ ∧ ψ and U = ∇χ + ∇ ∧ Θ. (9)
After equations decoupling, we have

ϕ = ϕ 1 + ϕ 2 , χ = µ 1 ϕ 1 + µ 2 ϕ 2 , (10)
where

µ i = P k 2 i -ω 2 ρ 11 ω 2 ρ 12 -Qk 2 i , i = 1, 2, (11)
are the waves amplitude ratios between the two phases in the porous material. Similarly, the vector potential Θ is simply obtained as

Θ = µ 3 ψ with µ 3 = ρ 12 /ρ 22 . (12)
Under this form, each potential ϕ 1 , ϕ 2 , ψ fulfills the Helmholtz equation with the associated complex wave numbers

k 2 1 = ω 2 2(P R -Q 2 ) (P ρ 22 + Rρ 11 -2Qρ 12 + √ D), (13)
k 2 2 = ω 2 2(P R -Q 2 ) (P ρ 22 + Rρ 11 -2Qρ 12 - √ D), ( 14
)
k 2 3 = ω 2 N ρ 11 ρ 22 -ρ 2 12 ρ 22 . (15)
Here, P = A + 2N , D stands for the discriminant of a quadratic equation and D = (P

ρ 22 + Rρ 11 -2Qρ 12 ) 2 - 4(P R -Q 2 )(ρ 11 ρ 22 -ρ 2 12
). Physically, there are two compressional waves associated with ϕ 1 , ϕ 2 and one rotational (shear) wave associated with ψ. They all propagate in the two phases and their relative contribution are given by µ i . In elastodynamics, one compressional and one rotational waves are already present, but for poroelatic materials, the presence of two phases, adds a fluid-born compressional wave.

Similarly, we can treat the displacement w in the fluid domain using a similar decomposition. Now, since the flow is uniform, we can look for purely acoustic modes by putting

w = ∇ϕ 0 . (16)
Thus the rotational terms which are known to be associated with the hydrodynamic modes are discarded in the present analysis. Note that ( 16) is valid within the assumption that the mean flow is uniform as acoustic and hydrodynamic modes are found to be decoupled in this particular case.

III. MODE CALCULATION A. Eigenvalue equation

In this section, the silencer chamber is assumed first to be infinite in length. The eigenvalue analysis is performed by assuming that all perturbative quantities (call it X) have the same dependence in both z and t, so we write

X(r, θ, z, t) = X(r, θ) e i(βz-ωt) (17)
Mode matching method for poroelastic silencers 3

The vectorial nature of the Helmholtz equation for the potential ψ renders the θ-dependence a bit more tricky as its three components in the cylindrical vector basis do not behave with the same azimuthal parity. Following a similar work on hollow elastic cylinders as in Ref. 25, we find that the general form for the potentials splits into two independent families of modal solutions. Symmetric modes in the poroelastic media are given explicitely by

φ1 = [A 1 J m (α 1 r) + B 1 Y m (α 1 r)] cos mθ, (18a) φ2 = [A 2 J m (α 2 r) + B 2 Y m (α 2 r)] cos mθ, (18b) ψr = [A 3 J m+1 (α 3 r) + B 3 Y m+1 (α 3 r)] sin mθ, (18c) ψθ = -[A 3 J m+1 (α 3 r) + B 3 Y m+1 (α 3 r)] cos mθ, (18d) ψz = [A ′ 3 J m (α 3 r) + B ′ 3 Y m (α 3 r)] sin mθ, (18e)
where J m and Y m are respectively the first kind and the second kind Bessel functions. Antisymmetric modes are found by simply swapping the trigonometric functions cos mθ and sin mθ and by removing the minus sign for the azimuthal component ψθ . Note the radial and azimuthal components are dependent via their coefficients. This is the gauge condition borrowed from Ref. 25 stemming from the non-unicity of the Helmholtz decomposition. Radial wave numbers α i are obtained from the dispersion relation as

α i = k 2 i -β 2 , i = 1, 2, 3. (19)
where the choice of the branch cut for the square root operator is irrelevant here as we are not dealing with radiation conditions (which would then imply a careful examination of the behaviour of each physical quantity as r goes to infinity). The different branch cuts may, at most, alter the sign of the transverse wavenumber, which has no physical significance as they are not associated with other modes. If a precise definition of the square root can sometimes be judicious for reasons of symmetry [START_REF] Rienstra | Contributions to the theory of sound propagation in ducts with bulk-reacting lining[END_REF] , this can not be followed in the present study due to the presence of the mean flow.

Similarly, the acoustic displacement scalar potential is given by φ0 = A 0 J m (α 0 r) cos mθ (or sin mθ) (20) and from Eq. ( 16) the acoustic pressure is recovered as p = ρ 0 c 2 0 Ω 2 φ0 . Here, the tranverse wavenumber is obtained from the associated dispersion relation and we take

α 0 = Ω 2 -β 2 , (21)
where Ω = k 0 -M β, M = U/c 0 is the Mach number and k 0 = ω/c 0 denotes the usual acoustic wavenumber. The last step is to express the components of the total stress tensor σ t , the displacements u, U, w as well as the fluid pressures p and p p in terms of the displacement potentials and their derivatives. After substitution in the boundary conditions (7) and ( 8), we find that the modal vector containing the wave potential amplitudes

V = (A 0 , A 1 , B 1 , A 2 , B 2 , A 3 , B 3 , A ′ 3 , B ′ 3 )
T must be a non-trivial solution of the 9×9 algebraic system

M(α 0 , α 1 , α 2 , α 3 , β)V = 0. (22)
In this work, we shall focus on purely axisymmetric solutions (m = 0). In this particular case, both axial and radial components of the vector potential are discarded and the 7 × 7 matrix coefficients are given explicitly in Appendix A. At this point, we shall observe that the matrix M depends solely on the axial wavenumber β once the dispersion relations have been inverted via (19) and (21). Thus, finding the modes amounts to finding the zeros of the complex function f (β) = det M(β), i.e. we have the following eigenvalue equation in β,

f (β) = 0. (23)

B. Root finding

The accurate computation of the roots of Eq. ( 23) represents a crucial stage as it provides the incident and reflected axial wavenumber as well as the associated eigenmode transverse profiles obtained from (22). There are numerous available numerical techniques for this purpose. Although this is not the place for a complete survey we can cite the Newton-Raphson method, Muller's method, the Secant method (see discussion on this topic in Ref. 13) and the Nelder-Mead simplex method [START_REF] Alonso | Eigenvalue solution for the convected wave equation in a circular soft wall duct[END_REF] . All these techniques have in common that they all require initial approximations for the zeros to start the algorithm. To avoid any missing roots, application of the Argument Principle is commonly favored [START_REF] Lawrie | Mode-matching without rootfinding: Application to a dissipative silencer[END_REF] . This fundamental result, based on the Cauchy's Theorem, allows to compute the number of zeros (including its multiplicity) of f from the following integral relation

S 0 = 1 2πi C f ′ (β) f (β) dβ = N β k=1 n k , (24)
where N β is the number of zeros and n k the corresponding multiplicity of the k-th zero lying in the interior of the closed curve C. Note f is chosen to be analytic for the result to hold but the presence of poles may be included by modifying the formula accordingly. In the present work, none of the root-finding techniques cited above will be used. Instead, we shall employ an analytic method that does not require knowledge of initial guesses. This is a classical procedure based on the generalization on the previous relation to any monomial β n , that is 29

S n = 1 2πi C f ′ (β) f (β) β n dβ = N β k=1 n k β n k , (25)
where β k denote the position of the k-th zero. To simplify further the analysis, function f is assumed to have only simple zeros so the integral ( 24) is exactly the number of zeros contained in the interior of C. Note this is not a stringent assumption as the occurrence of multiple zeros is highly improbable (in a somewhat different context, it is known that multiplicity may exist for specific material Mode matching method for poroelastic silencers 4

properties set for which simple zeros are merging [START_REF] Nilsson | The propagation of sound in cylindrical ducts with mean flow and bulk-reacting lining. I. modes in an infinite duct[END_REF] ). We now introduce the associated polynomial Π for the interior of C, that is a polynomial having the same zeros as f , i.e.

Π(β) = N β k=1 (β -β k ) = N β n=0 C n β n . (26)
Of course the zeros are not known yet, but the polynomial coefficients C n can be efficiently computed from the following recursive algorithm 16,17

C n =        1, if n = S 0 1 n -S 0 S0-n j=1 S j C n+j , if n = S 0 -1, . . . , 0. (27)
Here S n (n = 0, 1, 2, . . . , N β ), denote the values of the integral (25). Once Π is known, finding its zeros simply requires computing the eigenvalues of the companion matrix [START_REF] Edelman | Polynomial roots from companion matrix eigenvalues[END_REF] which is a very fast procedure for moderate size matrices. In fact, the most time consuming operation is the numerical evaluation of Bessel functions with complex arguments at some quadrature points when computing integrals (25). Since these are integrals of an analytic periodic function over a complete period, the trapezoidal rule is the optimal quadrature rule. Let C be defined from a regular function γ(s) over the interval s ∈ [0, 2π], then the q-point trapezoidal rule approximation to S n is given by

S n ≈ 1 iq q-1 i=0 γ n f d f ds (i/q), (28)
where we put f (s) = f (γ(s)). Though the s-derivative of f may be obtained formally via symbolic software, it is far more efficient to evaluate the derivative using high order central finite difference scheme, that is

d f ds (i/q) ≈ 1 i/q J j=-J a |j| f i + j q . ( 29
)
The procedure is extremely fast since the discrete values of f on the regular grid are already calculated. Here we used either the 5 th or 9 th -order scheme. The corresponding coefficients are displayed in Table I. The choice of the curve C must depend upon the region of the complex plane where eigenvalues are searched. Since eigenvalues are expected to be symmetric with respect to the origin (at least in absence of mean flow), choosing the circular path γ = ae is appears to be a good compromise.

To avoid possible round-off errors, it is then preferable to factorize the term a n and exclude it for the trapezoidal summation of Eq. ( 25). Now, as most of the zeros corresponding to very evanescent fields in the axial direction are located near the imaginary axis, elliptical contour integration may also be used. In this case we take γ = a cos s + ib sin s with b > a.

If f has many zeros in the search area, then the map from the S n 's to the polynomial coefficients computed from ( 27) is likely to be ill-conditioned. Moreover, the location of the zeros of the associated polynomial are also very sensitive to its coefficients. The reasons for this are discussed in Ref. 29 but the nature of the problem can be identified from (26). It suffices to observe that C N β /C 0 = N β k=1 β k can be extremely large as some of the β k 's are located far from the origin. The consequences on the roots location may be significant. In addition our numerical tests have shown that the lowest attenuation modes are the most affected by this problem, leading to a discernable frequency shift on the Transmission Loss. One way to overcome this is to use high precision arithmetic 32 which is not convenient and timeconsuming. A simple trick to avoid this is to reduce the search area so that the number of zeros is kept to a reasonable value [START_REF] Delves | A numerical method for locating the zeros of an analytic function[END_REF] . The solution employed here is to split the integration path into concentric circular (or elliptical) rings. In order to illustrate the procedure, we take the silencer A, which dimensions are given in Table III, filled with a XFM foam (see Table IV). We consider computing roots located in the interior of a circle of radius a = 180. In the first approach, we take one circular ring. We find that the search region contains 20 roots. By increasing the number of quadrature points q, the numerical value of a specific root chosen to be associated with one of the lowest propagation mode and for which error was maximum, is given in Table II. This shows poor convergence due to round-off errors. In this example, increasing the number of points does not yield a better solution. In the second approach, we take two circular rings. The search region is then divided into two subregions each containing 10 roots. The improvement is clearly shown as 2500 quadrature points on each ring provides up to 12 digits of accuracy. Here the number of points is chosen so that the numerical cost remains comparable for both approaches.

In practice, we found that very good accuracy was obtained if (i) there is no more than 10,15 roots within each search area in order to avoid conditioning problems and possible round-off errors, (ii) a sufficient number of quadrature points on each ring must be taken. In the work presented in this paper, the root finding algorithm is entirely automated. For a given frequency, the first step is to identify the largest search area (or the radius of the largest circle) in the complex plane containing approximately 2K modes (K is the number of modes in each direction of propagation and is prescribed by the user). This is done heuristically, by first considering that the chamber is empty, that is we consider classical acoustic modes in a rigid cylinder of radius r 2 . In this case the relation between the number of roots and the radius is an easy task since zeros of Bessel functions are already Mode matching method for poroelastic silencers 5 known. In the presence of a poroelastic material, there are other types of modes (see discussion in section III.C) and it was observed that the number of modes in each search area is approximately 4 times larger. Given this, the largest radius can be found. Now, since roots are distributed rather regularly in the complex plane, the largest search area is then divided using regularly spaced concentric circles with the condition that each search area contains, say 10 up to 15 roots. Finally, the number of quadrature points on each circular path points is found on the ground that between 500 and 1000 points for a circle of radius a = 100 should provide a sufficient accuracy (see Table II). For a given radius, the number of points is chosen accordingly as to grow linearly with the radius.

C. Eigenmodes

Without flow, the roots distribution in the complex plane is remarkable as shown in Figure 2. Concentric rings used in our calculation have also been displayed for the sake of illustration. Two families of roots can be identified.

The first set comprises roots located near and along the imaginary axis. These are essentially associated with fluid modes which are also present when considering equivalent fluid models [START_REF] Nilsson | The propagation of sound in cylindrical ducts with mean flow and bulk-reacting lining. I. modes in an infinite duct[END_REF] . Here, the fluid embedded in the foam pore is strongly coupled with the one in the airway whereas the solid frame displacement is nearly zero. This is illustrated in Figures 34where the pressure and displacement modes profiles are plotted for three roots identifiable with the corresponding marker in the complex plane of Figure 2 (note that in all examples described in this section, results are obtained with silencer A with a XFM foam). The eigenmode profiles are computed once the root has been found by simply inverting the one-rank deficient system (22) and applying the Helmholtz decomposition ( 9) and ( 16).

In the second set of modes, the skeleton motion is strongly coupled with the airway and their presence in the modal series is essential in order to respect the clamped boundary conditions on the hard wall of the silencer, this fact is illustrated in section IV.C. The fact that the associated wave numbers do not lie on the real or imaginary axis is due to the shear modulus and not to the dissipative properties of the material (see for instance Ref. 35 in the case of the elastic cylinder). In the low frequency limit, Biot's coefficients tend to be real quantities and the non-dissipative case is recovered. In this regime, eigenvalues occur in pairs with the same negative conjugate wave number and are associated with the same eigenmode. Thus their sum yields standing wave vanishing when z → ±∞. It is worth noting that standing waves do not carry energy, so the decrease in amplitude do not correspond to an energy dissipation. In our case, because of the dissipation, the zeros are not negative conjugate as illustrated on Figure 5 especially when the frequency increases (i.e. when the dissipation increases). As well, the eigenmodes are not identical : the solid phase displacement remains almost unchanged but the fluid phase pressure is different as illustrated in Figure 3 and4. One may also observe in Figure 5 that, for both family of modes, some zeros may cross the imaginary axis as the frequency increases. This behavior which is attributed to the dissipation has already been observed by Nilsson 30 for bulk reacting liner.

The evolution of the roots location with respect to the Mach number is conveniently displayed in Figure 6. We can see that eigenvalues of the second family are not affected by the presence of the mean flow whereas roots of the first family (the fluid modes) are very sensitive. Now, the presence of another mode of a different nature can also be identified. This mode, only present with gas flow has deep connection with the stability of the vortex sheet (of infinitely small thickness) that exits in the vicinity of the deformable boundary Γ. A quick inspection of the radial profiles shows that displacement and pressure fields tend to be confined to within a small neighborhood of the boundary and decrease exponentially away from it. In the context of a locally reacting liner, this mode is sometimes called 'surface wave mode' [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF][START_REF] Brambley | Classification of aeroacoustically relevant surface modes in cylindrical lined ducts[END_REF] . This surface mode has been kept in the matching procedure (see next section) but its presence was found to have no noticeable effect in our applications.

Mode matching method for poroelastic silencers 6 [START_REF] Debergue | Boundary conditions for the weak formulation of the mixed (u, p) poroelasticity problem[END_REF][START_REF] Panneton | Modélisation numérique tridimensionnel par éléments finis des milieux poroélastiques[END_REF] . With the resistivity σ, the tortuosity α inf , the viscous and thermal characteristic lengths Λ and Λ ′ , the poisson coefficient ν and the effective skeleton density ρ1. The effective skeleton density ρ1 = (1φ)ρs, where ρs is the density of the material of the frame. •. For r < r1, acoustic pressure p and for r1 < r < r2 pore pressure pp. The amplitude max. is scaled to one.

Foam φ σ (kNm -4 s) α inf Λ (µm) Λ ′ (µm) ρ1 (kgm -3 ) N (kPa) ν FM4 0.

IV. MODE MATCHING A. Matching conditions

The classical set of matching conditions for bulk material [START_REF] Kirby | Simplified techniques for predicting the transmission loss of a circular dissipative silencer[END_REF][START_REF] Selamet | Analytical approach for sound attenuation in perforated dissipative silencers[END_REF] is extended to poroelastic material. For a brief nomenclature, we note w = (w r , w z ), u = (u r , u z ) and U = (U r , U z ), the vectorial displacements in the cylindrical basis. In the airway (r ≤ r 1 ) we have the continuity of pressure as well as the normal acoustic displacement Mode matching method for poroelastic silencers 7 over the inlet plane Γ 1 :

p I = p II , ( 30a
) w I z = w II z , (30b)
where the superscript (i=i,ii) indicates the associated region in which the physical quantity of interest is defined.

On the lateral hard wall (r 1 ≤ r ≤ r 2 ), the kinematic condition yields

u II r = u II z = U II z = 0, (31)
and similar conditions hold on the exit plane Γ 2 . The condition (31) expresses the fact that the foam is clamped on the lateral side of the silencer chamber. If continuity conditions (30) are standard in the no-flow case, the transmission problem becomes far less trivial when an uniform flow is present. In this latter case, the pressure and displacement fields exhibit some singular behavior at the plane of discontinuity which must then be taken into account in the matching procedure. The difficulty is partially solved and discussed in a recent paper from Gabard & Astley 40 in the case of a locally reacting liner with impedance discontinuity. In the present scenario involving a transition with a poroelastic material, it is expected that the physical mechanisms that are taking place are a hard problem deserving a proper analysis. Therefore we will leave this point for further studies.

B. Numerical implementation

In regions i=i, ii, iii, each quantity is expanded via its truncated modal decomposition. For instance, the acoustic pressure in the airway is expressed as (the time dependence is omitted for clarity)

p i = K k=1 e iβ i,+ k z Φ i,+ p,k (r)A i,+ k + e iβ i,- k z Φ i,- p,k (r)A i,- k . (32)
In the above summation, the origin of the duct axis coordinate is arbitrary and in this work, we chose to fix it at the inlet plane Γ 1 . In the dissipative silencer (i=ii), none of the eigenvalues lies on the real axis and we took the convention that superscript ± stands for the sign of the imaginary part of the axial wave number. If Im β > 0 the mode is said to be right running otherwise the mode is said to be left running. In the inlet and outlet pipes some modes are propagative, right (+) and left (-) running modes are then defined in the usual way as follows

β i,± = -k 0 M ± k 2 0 -(1 -M 2 )(α i 0 ) 2 1 -M 2 , (33)
where i= i or iii and transverse wavenumbers α i 0 are obtained from the rigid wall condition: J 1 (α i 0 r 1 ) = 0. Coefficients A i,± k are the modal amplitudes, Φ i,± p,k (r) are the corresponding radial eigenfunctions and subscript p indicates the physical quantity (here the pressure) it is associated with. A similar decomposition holds for the displacements w r , w z , u r , u z , U r and U z . For the sake of simplicity, the number of modes K of the truncated modal series is taken identical in each region, but this might be changed if necessary [START_REF] Kirby | Analytic mode matching for a circular dissipative silencer containing mean flow and a perforated pipe[END_REF] . More importantly, the correspondence k ↔ β II,+ k is given by ordering the imaginary part of the eigenvalues is ascending order so that lowest attenuation modes are always included in the series. This is also done for the leftrunning modes. In regions i and iii, the same ordering holds once purely propagative modes have been taken into account. Now, to clarify the analysis further, we introduce the line vector Φ i,± p containing the radial eigenfunctions

Φ i,± p = Φ i,± p,1 , . . . , Φ i,± p,K , (34)
and similarly for the other variables. With this definition, we have at the inlet plane

p II = Φ II,+ p A II,+ + Φ II,- p A II,-, (35)
and

p I = Φ I,+ p A I,+ + Φ I,- p A I,-. (36)
Here A i,± is the column vector containing the modal amplitudes. At the exit plane (z = L), the same formula holds with the addition of the phase factors.

Matching conditions in the airway are imposed in a weighted sense, i.e. we proceed by choosing a weighting function and then integrating over the cross-sectional area of the inlet (or outlet) pipe:

r1 0 p I Ψ p r dr = r1 0 p II Ψ p r dr, ( 37a
) r1 0 w I z Ψ wz r dr = r1 0 w II z Ψ wz r dr. ( 37b
)
Clamped conditions on the hard wall of the silencer must be imposed with more care. Indeed, through early numerical experiments, it was found that imposing (31) separately for each quantity is likely to produce illconditioned systems. It is therefore preferable to relax Mode matching method for poroelastic silencers 8 the constraints by summing the three conditions giving the weighted form

r2 r1 u II r Ψ ur + u II z Ψ uz + U II z Ψ Uz r dr = 0. (38)
Here, weighting functions Ψ p , Ψ wz , Ψ ur , Ψ uz and Ψ Uz are chosen among the radial eigenmodes associated with each physical variable. We will discuss the most appropriate choice further on. Now, after modal substitution we obtain a classical scattering system with unknowns A I,-, A II,± and A III,+ whereas the known incident field given by the value of A I,+ and A III,-are placed on the right-hand side. The inversion of such a system is unfortunately subject to round-off errors due to the presence of strongly evanescent waves. Because of these limitations, we consider an alternative procedure that was proposed by Cummings [START_REF] Cummings | Sound attenuation of a finite length dissipative flow duct silencer with internal mean flow in the absorbent[END_REF] . The idea is to treat each matching condition separately by considering all incoming waves as known quantities. The system is then solved for the outgoing waves. This yields, using the same formalism as in Refs. 40, 41, the following coupled scattering system

X 1 A I,- A II,+ = Y 1 A I,+ A II,-, (39a)
X 2 E X A III,+ A II,-= Y 2 E Y A III,- A II,+ . (39b)
Here, E X and E Y are diagonal matrices containing the phase factors at the exit plane z = L,

E X = diag e iβ III,+ 1 L . . . e iβ III,+ K L , e iβ II,- 1 L . . . e iβ II,- K L , (40)
E Y = diag e iβ III,- 1 L . . . e iβ III,- K L , e iβ II,+ 1 L . . . e iβ II,+ K L . (41)
The scattering matrices have the following form

X 1 =   P I ,--P II ,+ W I ,--W II,+ 0 -ŴII,+   , Y 1 =   -P I ,+ P II ,- -W I,+ W II,- 0 ŴII,-   , (42)
X 2 =   P III ,+ -P II,- W III ,+ -W II ,- 0 -ŴII,-   , Y 2 =   -P III ,-P II ,+ -W III,-W II ,+ 0 ŴII,+   , (43)
where block matrices P i,± , W i,± and ŴII,± are computed from the overlap integrals

P i,± = r1 0 (Ψ p ) † Φ i,± p r dr, for i = I,II,III, (44)
W i,± = r1 0 (Ψ wz ) † Φ i,± wz r dr, for i = I,II,III, (45)
and

ŴII,± = r2 r1 (Ψ ur ) † Φ II,± ur r dr + r2 r1 (Ψ uz ) † Φ II,± uz r dr + r2 r1 (Ψ Uz ) † Φ II,± Uz r dr. (46)
Here Ψ p , Ψ wz , Ψ ur , Ψ uz and Ψ Uz are the line vectors containing the corresponding weighting functions and symbol † means that we take the transpose conjugate of these vectors. As pointed out by Gabard & Astley [START_REF] Gabard | A computational modematching approach for sound propagation in threedimensional ducts with flow[END_REF] , the choice of these functions can have a significant impact on the performance of the method. In order to enhance the conditioning of the scattering matrices X i (i = 1, 2) it is preferable to use appropriate weighting functions in order to maximize their diagonal terms. Thus the eigenfunction of the each diagonal block of X i (i = 1, 2) are chosen as weighting functions. For example rigid wall modes are used for the matrices P i,± , i.e. we take Ψ p = Φ I,- p and Ψ p = Φ III,+ p (these are in fact the same in this study since the inlet and outlet pipes are identical) respectively in matrices X 1 and X 2 . Modes in the expansion chamber are used for the matrices W i,± with Ψ wz = Φ II,+ wz and Ψ wz = Φ II,- wz respectively for matrices X 1 and X 2 . The same procedure is followed for the other matrices.

Finally equations corresponding to the continuity of displacements are all multiplied by the factor (ωρ 0 c 0 ) 2 in order to be of comparable magnitude with the ones corresponding with the pressure. Overlap integrals along the radial coordinate are performed using classical Gaussian quadratures though analytical form could possibly be found as in Ref. 13. In practice, between 20 and 50 points were found sufficient in our applications shown in the next section.

The system is then solved iteratively. At the first iteration, A II,-is fixed at zero and the system (39a) is inverted. This produces an initial value for the rightgoing waves A II,+ which is then used at the exit plane system (39b) and so on. The pseudo-inverse of the scattering matrix is performed using the pinv command in Matlab. Note that the factorization of the diagonal matrix E X is one of the main advantages of this iterative process. Indeed its condition number is very large and can exceed 10 40 due to the combination of rightgoing and leftgoing 'evanescent' waves. By inverting E X analytically, effects of possible round-off errors are therefore minimized.

C. Results and validation

In order to investigate the robustness of the proposed methodology, we present here some results with MMM in comparison with FE calculations carried out by the authors and presented in Refs. 10, 23. The silencers configurations from Table III have the same dimensions as in Ref. 34 and poroelastic material properties can be found in Table IV. In the airway, the density and the celerity of sound are ρ 0 = 1.213 kgm -3 and c 0 = 342.35 m/s. As the frequency range of interest is considered below the first cut-on mode of the rigid pipes, only the plane wave propagates in these regions and the Transmission Loss for the silencer is simply given by (we take |A I (47)

In Figure 7, the convergence of the MMM is presented for the silencer A by taking successively 5, 7, 13 and 15 modes in each propagation direction. The XFM polymer foam has been chosen to illustrate the effects of the Mode matching method for poroelastic silencers 9

strong coupling between fluid and solid phases that exists in the poroelastic material [START_REF] Nennig | Influence of solid phase elasticity in poroelastic liners submitted to grazing flows[END_REF] . We can see that 15 modes yields good agreements with FE results. The significant number of required modes, especially near the elastic resonances identifiable at the peaks, is attributed to the clamped conditions on the exit planes Γ 1 and Γ 2 .

In particular, modes of the second family associated with the solid phase are unavoidable in the matching procedure whereas the highly attenuated modes of the first family (fluid modes) may be ignored. In the low frequency range (below 500 Hz) and in the no flow case, it was found that only the lowest attenuation fluid mode could be kept in the modal series without deteriorating the results. In Table V is shown a convergence analysis (in K) for the TL (dB) for three frequencies 500, 900 and 1000 Hz. The frequency 900 Hz with the XFM foam corresponds to a resonance peak and the convergence is slower in this case. At 500 Hz, and more generally below the first resonances, the TL converges rapidly to one or two decimal places and this is in line with the TL predictions of Kirby 1 using Analytical Mode Matching.

It is worth noting that if Eq. ( 38) is added to Eq. (37b), as it is generally done when dealing with bulk reacting materials [START_REF] Kirby | Simplified techniques for predicting the transmission loss of a circular dissipative silencer[END_REF] , the convergence is not as good and 15 modes in the series was found to produce up to 3 dB error between 700 and 1300 Hz. This is particularly relevant for foams presenting strong solid phase resonances. However, for soft foams like FM4 (Table IV), accurate results can still be obtained with this 'condensed' formulation.

Calculations on silencer B which presents a strong area ratio (r 2 /r 1 ) 2 ≈ 7.5 also compare very favorably with FE results as illustrated in Figure 8. Here again, resonance effects due to the skeleton elasticity are clearly visible below 1000 Hz.

The presence of the mean flow slightly affects the quality of the results especially with the RGW2 wool (see Figure 9) though the error does not exceed 1 dB. In particular, discernable discrepancies are visible in the low frequency regime (say below 1000 Hz) and these differences can not be explained on the ground that the FE model is underdiscretized with respect a typical wavelength. The physical reason for the nature of the solution to this difficult problem is not clear and the authors are reluctant to speculate in this paper as to its cause. We shall just mention, in passing, two remarks: (i) no particular care has been taken for the edge conditions at the inlet and outlet planes, (ii) hydrodynamic modes which can be generated and discretized in the FE model are deliberately discarded in the present modal approach.

Figure 10 shows the TL predictions for different Mach number: M = -0.3, -0.1, 0, 0.1 and 0.3. The mean flow modifies the apparent acoustic wave length in the airway and changes the number of oscillations in the absorbing material. The flow is found, in general, to reduce the attenuation for the downstream propagation and to increase it for the upstream propagation. However near the resonance peaks the behavior may be different due to some spatial coincidence. Note that comparisons with FE results were found to be in good agreement in all cases (for the sake of illustration, a comparison is shown in Figure 11 for M = 0.2). The behavior of the RGW2 wool, which corresponds to a more common material in dissipative silencers for HVAC applications, is depicted in Figures 9 and12. We also present a comparison with FE calculations 23 using the limp model 2 ; that is, an equivalent fluid model which considers infinitely soft frame and takes into account only skeleton inertial effects. For rigid frame materials, this model gives similar results as the rigid frame model [START_REF] Panneton | Comments on the limp frame equivalent fluid model for porous media[END_REF] . These results indicate that properly modeling the skeleton can lead to significant improvements of few dBs on the Transmission Loss even for wool-like materials. The result is not too surprising as it was already reported that an edge constrained-fibrous sample placed in a standing wave tube can exhibit shearing resonances [START_REF] Song | Investigation of the vibrational modes of edge-constrained fibrous samples placed in a standing wave tube[END_REF] . Thus, a precise description of the mounting conditions can have serious consequences on the effectiveness of the acoustic treatment.

Before we finish this section, we shall comment on the computational efficiency of the MMM. All previous results were obtained using four integration circles each containing 500 to 1000 quadrature points. In Matlab implementation, it takes approximately 5 seconds to find the roots for a given frequency. When 15 modes are included in the modal series, about 20 seconds are needed for the overall matching procedure on a pentium IV personal computer. Note that the CPU time can be largely reduced (say by a factor of 5 to 10) with compiled languages such as Fortran or C and this renders the present method an efficient one. The authors are aware that the root finding algorithm presented here was favored because it does not suffer from missing roots. The method is certainly slower (in terms of number of function evaluations) than other techniques such as the popular Newton Raphson method and we think that the computational time could be further reduced by judiciously exploiting the benefit of the two methods.

Mode matching method for poroelastic silencers 10

V. CONCLUSIONS AND PROSPECTS

In this paper, we have presented a new Mode Matching Method (MMM) for acoustic propagation through circular dissipative silencers partially filled with a poroelastic material. The model, based on the Helmholtz decomposition of the Biot's equations, takes into account the solid phase elasticity of the sound absorbing material as well as the mounting conditions of the foam. In the noflow case, comparisons with a full FE model show good agreements even at the resonance peaks corresponding to elastic resonances of the foam which is assumed to be clamped on the hard wall of the silencer. It is also shown that wool materials which are usually described via fluid equivalent models can, when the elastic phase is included, exhibit some resonance effects clearly visible on the Transmission Loss curves.

When an uniform mean flow is present in the airway, the problem remains separable and the modal analysis still holds. However, a mode of a different nature is also present. This mode has connection with the stability of the flow in the vicinity of the permeable surface separating the gas flow and the foam. It was found, however, that the inclusion of this mode in the modal series has no noticeable effects regarding the numerical exam-Mode matching method for poroelastic silencers 11 ples shown in this work. Comparisons with FE results showed small but discernible discrepancies. In all cases, it is believed that this is is not a major issue for engineering applications, as the error on the TL does not exceed 1 dB. In practice, the present method is shown to be very efficient as the time for computing the TL over a large frequency range takes a few minutes on a personal computer. Thus, this new MMM remains very competitive compared with classical mode matching techniques in which the sound absorbing material is assumed to be bulk reacting. Now, there is still room for some improvement of the method. As computational time is essentially spent in computing Bessel functions with complex arguments, the use of asymptotic expansion could possibly be utilized for sufficiently large arguments. Another route of particular interest to us is to investigate bi-orthogonality relations [START_REF] Scandrett | Bi-orthogonality relationships involving porous media[END_REF] that exist for poroelastic media in order to extend it to our configurations. It is hoped that this could then yield a better understanding of the matching procedure at the interface leading to better conditioned scattering matrices. Work is on-going to include the effect a perforate plate at the air-porous interface so that more realistic silencers could be studied. In principle, following previous studies [START_REF] Nilsson | The propagation of sound in cylindrical ducts with mean flow and bulk-reacting lining. I. modes in an infinite duct[END_REF][START_REF] Kirby | Simplified techniques for predicting the transmission loss of a circular dissipative silencer[END_REF] , this is not a major task once the transmission properties of the plate are known.

FIG. 1 .

 1 FIG.1. Geometry of the silencer.

1 FIG. 4 .FIG. 5 .

 145 FIG. 4. Modulus of modal displacement profile related to markers shown in Fig. 2 at 400 Hz for silencer A filled with XFM foam:⋄, ,•. For r < r1, acoustic axial displacement wz and for r1 < r < r2 solid phase axial displacement uz. The amplitude max. is scaled to one.

FIG. 6 .

 6 FIG.6. Trajectory of the roots in the complex β -plane at 400Hz with respect to the Mach number ranging from 0 (•) to 0.8 ( * ) for XFM foam with silencer A.

1 ,

 1

FIG. 7 .

 7 FIG. 7. TL for silencer A without flow and XFM foam, 78,000 dofs 3D FE calculation 10 , MMM with 5,7 and 13 modes, • MMM with 15 modes.

FIG. 8 .FIG. 9 .FIG. 10 .

 8910 FIG. 8. TL for silencer B without flow and XFM foam, 100,000 dofs 3D FE calculation 10 , • MMM with 16 modes.

FIG. 12 .

 12 FIG. 12. TL for silencer A with a M = 0 flow and RGW2 wool, 78,000 dofs 3D FE calculation 10 , • MMM with 16 modes, 14,000 dofs Limp model FE calculation 2,23 .

TABLE I .

 I Finite difference scheme coefficients.

		a0	a1	a2	a3	a4
	5 th order	0	2/3	-1/12	-	-
	9 th order	0	4/5	-1/5	4/105	-1/280

TABLE II .

 II Numerical computation of a specific root at 1730 Hz. With one circular ring (left) and two circular rings (right). Bold character indicates the correct digits.

		One ring		Two rings
	q	β	q	β
	500	18.1310333520 + 17.2414041458 i	250	17.1619265538 + 15.2060360731 i
	1000	17.6094544451 + 14.9503659291 i	500	17.6293456921 + 15.0650726700 i
	5000	17.6308169754 + 15.0352734983 i	2500	17.6294474640 + 15.0651513163 i
	20000	17.6362164376 + 15.0289704442 i	10000	17.6294474640 + 15.0651513164 i

TABLE III .

 III Silencer dimensions 34 in meter.

	Silencer	L	r1	r2
	A	0.315	0.037	0.0762
	B	0.330	0.037	0.1016

TABLE IV .

 IV Materials properties used in numerical tests

  Roots in the complex β -plane at 400 Hz for XFM foam with silencer A.

			99			65			1.98	37	121	16	18(1 -0.1i)	0.3
	XFM		0.98			13.5			1.7	80	160	30	200(1 -0.05i)	0.35
	RGW2		0.99			9			1	192	384	16.3	200(1 -0.1i)	0
	600
	400
	200
	β
	Im	0
	-200
	-400
	-600
		-800	-600	-400	-200	0	200	400	600	800
						Re β
	0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 FIG. 2. r 0 0.01 0.02 0.03 0.04 1	0.05	0.06	0.07	0.08
	FIG. 3. Modulus of modal pressure profile related to markers
	shown in Fig. 2 at 400 Hz for silencer A filled with XFM foam:
	⋄,	,

TABLE V .

 V Convergence of TL predictions (in dB) with respect to the number of modes K. FIG. 11. TL for silencer A with a M = 0.2 flow and XFM foam, 78,000 dofs 3D FE calculation 10 , • MMM with 15 modes.

	K			500 Hz	900 Hz		1000 Hz
	1			14.5333	17.8914		23.2962
	2			13.9663	16.6979		15.5509
	5			15.5231	30.2493		26.2362
	7			15.7968	27.0677		25.0583
	10			15.8620	25.9387		24.6308
	13			15.9510	26.2165		24.8859
	17			15.9862	25.5172		24.3866
	22			15.9857	25.1448		24.2645
		60
		50
		40
	TL (dB)	30
		20
		10
		0	0	500	1000	1500	2000	2500
					f (Hz)

Acknowledgments

The authors would like to thank Dr. G. Gabard as well as the reviewers for their helpful discussions and comments.

APPENDIX A: MATRIX M COEFFICIENTS

The V vector is given by

and the matrix M coefficients are

Mode matching method for poroelastic silencers 12

M 7,5 = -J 1 (α 1 r 1 )α 1 (φ µ 1 -φ + 1) (A48)

M 7,7 = J 1 (α 0 r 1 )α 0 (A50)

with