Shaoqi Wu
email: shaoqi.wu@ec-nantes.fr

| Olivier Dazel

Grégory Legrain

A computational approach based on X-FEM for thin porous layers in acoustic problems

Keywords: X-FEM, Transfer Matrix Method, Generalized interfaces, Porous media, Fluid-Solid Interaction

The simulation of stationary acoustic fields involving thin porous layers is addressed in the present paper. Layer thickness is assumed to be relatively thin in comparison to the overall computational domain, but its non-negligible acoustic impact must be taken into consideration in the numerical model. Within the classical Finite Element Method (FEM), meshes are compatible at material interfaces and the element distortions needs to be avoided. These requirements usually lead to an excessively costly spatial discretization for these problems of interest, as it forces mesh refinement surrounding the thin layer. This paper provides a computational approach to relax this restriction. A generalized interface model derived from the plane wave Transfer Matrix Method (TMM) is established for modelling thin layers. We develop variationally consistent formulations to impose the interface conditions from models of thin layers for diverse coupling configurations, in which both acoustic fluid and poro-elastic media are taken into account. The computational domain is discretized using the eXtended Finite Element Method (X-FEM) in order to introduce strong discontinuities in elements independently of the mesh. Implementation of the proposed formulations within X-FEM is verified to be capable of providing accurate and robust solutions. The efficiency and flexibility of the present approach for multi-layer and complex geometry problems are demonstrated compared to classical interface-fitted finite element models through different simulation scenarii.

INTRODUCTION

This paper proposes an alternative computational approach for modelling thin layers in relatively large-scale vibro-acoustic problems. Thin layers are common in various physical problems such as shells [START_REF] Rubin | Cosserat Theories: Shells, Rods and Points[END_REF] , membranes, material coating, inter-phases in heterogeneous composite materials [START_REF] Bigoni | Asymptotic Models of Dilute Composites with Imperfectly Bonded Inclusions[END_REF] and resistive porous film in acoustic absorption systems [START_REF] Sugie | Absorption Characteristics of Fibrous Material Covered with Perforated Facing and Film[END_REF] . Even though the thickness of such layers is excessively small compared to the characteristic dimension of the entire domain, their effects are pronounced and cannot be ignored. However, it is prohibitively expensive and time-consuming for traditional numerical methods such as the Finite Element Method (FEM) to capture the behaviour of thin layers, for two reasons: (i) elements need to coincide without distortion on the layers interface, leading to high resolution meshes close to and in the layers, especially when their geometry is complex; (ii) Such meshes lead to large and potentially ill-conditionned discrete systems due to the amount of unknowns and mesh geometry. Therefore, finding an appropriate model and effective numerical method to tackle these issues is of prime interest.

The most employed strategy is to reduce the thickness of the layers to a limiting situation: non-thickness interface, and their physical behaviours are described through specific asymptotic models. To this end, abundant models were investigated, Hashin 4-6 , Benveniste [START_REF] Benveniste | The Effective Mechanical Behaviour of Composite Materials with Imperfect Contact between the Constituents[END_REF] , Gu et al. [START_REF] Gu | Interfacial Discontinuity Relations for Coupled Multifield Phenomena and Their Application to the Modeling of Thin Interphases as Imperfect Interfaces[END_REF] proposed a set of asymptotic models for inter-phase problems considered as a first-order spring-layers. Gurtin and co-works and Steigmann [START_REF] Steigmann | Plane Deformations of Elastic Solids with Intrinsic Boundary Elasticity[END_REF] developed a more general theory including second-order interface conditions [START_REF] Gu | Interfacial Discontinuity Relations for Coupled Multifield Phenomena and Their Application to the Modeling of Thin Interphases as Imperfect Interfaces[END_REF] . When it comes to the numerical solution of these asymptotic models within the framework of standard FEM, surface elements [START_REF] Nairn | Numerical Implementation of Imperfect Interfaces[END_REF] accounting for the additional stiffness of the interface behaviour are implemented. Double-nodes or unfitted meshes to interfaces [START_REF] Yedeg | A Nitsche-Type Method for Helmholtz Equation with an Embedded Acoustically Permeable Interface[END_REF] could be another manner to discretize the solution. More recently, Tiirats et al. [START_REF] Tiirats | How to Efficiently Apply Soft Thin Coating to Existing Finite Element Contact Model[END_REF] implemented a polynomial approximation to soft thin layer in Abaqus. However, the efficiency of these methods would decrease substantially as the geometry of the thin layers become more complicated, because standard FEM requires the boundaries of the elements to be aligned with the interfaces. Therefore, embedded discontinuity discretization techniques which involve non-smooth functions in the elements have gained a strong attention for these imperfect interface problems. These methods allow meshing to be independent of the interface geometry, leading to flexible approaches. Burman et al. [START_REF] Burman | CutFEM: Discretizing Geometry and Partial Differential Equations[END_REF] discussed the implementation of Cut-FEM for first-order discontinuity in various problems, and Han et al. [START_REF] Han | Consistent Discretization of Higher-Order Interface Models for Thin Layers and Elastic Material Surfaces, Enabled by Isogeometric Cut-Cell Methods[END_REF] employed isogeometric cut-cell methods for high-order interface models. Yvonnet et al. [START_REF] Yvonnet | An XFEM/Level Set Approach to Modelling Surface/Interface Effects and to Computing the Size-Dependent Effective Properties of Nanocomposites[END_REF] , Zhu et al. [START_REF] Zhu | Three-Dimensional Numerical Modelling by XFEM of Spring-Layer Imperfect Curved Interfaces with Applications to Linearly Elastic Composite Materials[END_REF] and Benvenuti [START_REF] Benvenuti | Variationally Consistent eXtended FE Model for 3D Planar and Curved Imperfect Interfaces[END_REF] used the eXtended Finite Element Method (X-FEM) to implement cohesive interface and spring-layer models in composites.

To prescribe weakly these interface conditions, Lagrange multipliers [START_REF] Mourad | A Bubble-Stabilized Finite Element Method for Dirichlet Constraints on Embedded Interfaces[END_REF][START_REF] Béchet | A Stable Lagrange Multiplier Space for Stiff Interface Conditions within the Extended Finite Element Method[END_REF] are frequently used. Consequently, more unknowns are introduced using this method, leading to a larger algebraic system to be solved. Besides, potential stability issues could be encountered when the multiplier space does not satisfy the inf-sup condition [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] . To avoid these difficulties, Nitsche's formulation [START_REF] Nitsche | Über Ein Variationsprinzip Zur Lösung von Dirichlet-Problemen Bei Verwendung von Teilräumen[END_REF][START_REF] Hansbo | A Finite Element Method for the Simulation of Strong and Weak Discontinuities in Solid Mechanics[END_REF][START_REF] Schillinger | The Non-Symmetric Nitsche Method for the Parameter-Free Imposition of Weak Boundary and Coupling Conditions in Immersed Finite Elements[END_REF][START_REF] Cai | A Nitsche Extended Finite Element Method for the Biharmonic Interface Problem[END_REF] (which can be considered as condensed augmented multiplier method [START_REF] Hautefeuille | Robust Imposition of Dirichlet Boundary Conditions on Embedded Surfaces[END_REF] ) has been investigated and applied to deal with Dirichlet-type or jump-type conditions. Indeed, the stabilization parameter in Nitsche's formulation needs to be defined carefully in order to ensure coercivity and stability of the solution as discussed in Wen et al. [START_REF] Jiang | A Robust Nitsche's Formulation for Interface Problems with Spline-Based Finite Elements: A Robust Nitsche's Formulation for Interface Problems with Spline-Based Finite Elements[END_REF] , which is not trivial for timeharmonic problems involving strong discontinuities. On the basis of the characteristics of the interface model and boundary integrals in the variational formulation, such discontinuity could be enforced directly as in the work of Zhu et al. [START_REF] Zhu | Three-Dimensional Numerical Modelling by XFEM of Spring-Layer Imperfect Curved Interfaces with Applications to Linearly Elastic Composite Materials[END_REF] , which is naturally consistent with the original problem and more straightforward for an implementation in the standard X-FEM.

In vibro-acoustic systems, a variety of thin layers exist such as thin elastic plates, impervious and porous screens, all of which have significant acoustical impacts. Unfortunately, there rarely exists systematic investigations in the community. Pressure drop [START_REF] Pierce | Acoustics: An Introduction to Its Physical Principles and Applications[END_REF] , which was first proposed by Pierce, is the most common model for fluid thin layers. Then, based on the asymptotic theory, a model for thin elastic and fluid layers in acoustic scattering problem was given by Bövik [START_REF] Bövik | On the modeing of thin interface layers in elastic and acoustic scattering problems[END_REF] . Wave propagation in acoustic media could also be modelled with aid of a matrix representation, known as Transfer Matrix Method (TMM) [START_REF] Folds | Transmission and Reflection of Ultrasonic Waves in Layered Media[END_REF][START_REF] Brouard | A General Method of Modelling Sound Propagation in Layered Media[END_REF] . This method is widely used for multiple layered acoustic systems, such as in [START_REF] Atalla | Modeling of Perforated Plates and Screens Using Rigid Frame Porous Models[END_REF][START_REF] Allard | Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials[END_REF] . More recently, Gaborit et al. [START_REF] Gaborit | A Simplified Model for Thin Acoustic Screens[END_REF] developed a simplified matrix representation for thin acoustic screens. However, these models are still solved in a "semi-analytical" manner and for simple geometrical configurations only. Yedeg et al. [START_REF] Yedeg | A Nitsche-Type Method for Helmholtz Equation with an Embedded Acoustically Permeable Interface[END_REF] considered permeable interfaces in homogeneous acoustic media where a pressure jump condition is imposed by a Nitsche-type method within standard FEM. In addition, the thickness of most of resistive films in sound package are typically around 0.1 to 1 mm rather than micro or nano scale as in previous mechanical problems where the volume of the thin layers can be directly neglected at the macroscopic scale. It is unclear if reduced interface models without considering the thickness is adequate for such films. As a result, it is quite natural to carry out a similar investigation to propose and validate an appropriate model for thin layers at macro scale for vibro-acoustic problems and solving it numerically within advanced mesh independent simulation techniques.

Sound absorption systems consisting of multiple porous materials are the main focus of current work. Energy loss exists in this variety of materials by means of structural, viscous and thermal effects. Two types of equations are commonly used to describe the sound behaviour of the materials: equivalent fluid model governed by Helmholtz equation with complex-valued material parameters such as JCA model [START_REF] Johnson | Theory of Dynamic Permeability and Tortuosity in Fluid-Saturated Porous Media[END_REF][START_REF] Champoux | Dynamic Tortuosity and Bulk Modulus in Air-saturated Porous Media[END_REF] or Limp model [START_REF] Katragadda | A Model for Sound Absorption by and Sound Transmission through Limp Fibrous Layers[END_REF][START_REF] Panneton | Comments on the Limp Frame Equivalent Fluid Model for Porous Media[END_REF] , and a series of coupling equations developed from Biot's thoery [START_REF] Biot | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low-Frequency Range[END_REF][START_REF] Biot | Mechanics of Deformation and Acoustic Propagation in Porous Media[END_REF] . The dissipation in equivalent fluid model is accounted by the imaginary part of material properties. It models only an airborne compressional wave in the material, assuming that the elastic frame is rigid enough or motionless beyond sufficiently high frequencies as opposed to air movement. On the contrary, Biot's theory provides a more complete description of porous materials denoted as poro-elastic materials (PEMs in the following) with fully coupled equations between fluid poro and elastic frame phases. These two governing equations are both considered in this paper for modelling porous materials.

This paper can be seen as an extension of our former work [START_REF] Wu | High-Order X-FEM for the Simulation of Sound Absorbing Poro-Elastic Materials with Coupling Interfaces[END_REF] , the application of higher-order X-FEM for acoustic analysis of poro-elastic structures. Thin porous layers are now introduced in the problem as strong discontinuities. The contribution FIGURE 1 Problem statement and computational domain definition.

of the current work has three components. First, providing a general methodology and formulation for the modelling of thin layers in time-harmonic problems, which is applicable for any acoustic media. Second, taking full advantage of X-FEM in space discretization and solution circumvents the meshing constraints for complex geometries and reduces computational costs compared to standard FE models. Third, an analytical and numerical assessment of the proposed model will be provided. This paper is organized as follows: section 2 presents the general problem and computational framework, as well as variational formulations of the target materials. Section 3 applies the proposed methodology to porous absorption systems. Interface models for porous layers based on the TMM with leading generalized variational interface operators for different coupling types are given. The structure of resulting systems is also discussed in this section. Section 4 is dedicated to offer more details on the proposed interface model, including a discussion of different transfer matrices for porous materials, a comparison with simplified model and sensitivity studies of incidence angle in TMM. In section 5, the proposed variational formulations are discretized in X-FEM space. The enrichment strategy as well as resulting discrete linear systems are provided. Section 6 reports the computational results that verify, validate and illustrate the accuracy and benefits of the proposed approach with three numerical examples.

METHODOLOGY AND FORMULATIONS

Acoustic problems in steady state (frequency domain) with a harmonic convention 𝑒 +j𝜔𝑡 are addressed. The domain of interest is divided into bulk parts with a characteristic length 𝐷 and thin layers parts with thickness 𝑑, where 𝑑∕𝐷 ≪ 1. To avoid the discretization of such narrow layers, we condense them as interfaces without thickness. All the material and geometry effects including thickness will be taken into account by appropriate interface conditions.

General formulation and interface operator

For the purpose of simplicity, we formulate the problem and equations in a bounded two-media configuration containing an embedded layer, as shown in fig. 1. The domain Ω (bounded by 𝜕Ω 1 and 𝜕Ω 2 ) consists of two subdomains Ω 1 and Ω 2 for bulk media with Ω 1 ∪ Ω 2 = Ω and Ω 1 ∩ Ω 2 = ∅. An interface Γ = 𝜕Ω 1 ∩ 𝜕Ω 2 with unit outward normal 𝐧 from Ω 1 to Ω 2 represents the condensed thin acoustic layer.

Within the finite element framework, variational formulations of such problem with specific boundary/interface conditions need to be given first. Without any loss of generality, a variational formulation can be written in primal variable 𝑢 𝑖 as the solution for each subdomain Ω 𝑖 . After combining the two subdomains, we obtain the final form as following: find

𝑢 𝑖 ∈ 𝑈 such that 2 ∑ 𝑖 [  𝒊 ( v𝑖 , 𝑢 𝑖 ) + 𝓘 𝒊 ( v𝑖 , 𝑢 𝑖 ) ] = 2 ∑ 𝑖 𝓵 𝒊 ( v𝑖 ), ∀𝑣 𝑖 ∈  ( 1
)
where v𝑖 is the test function of field 𝑢 𝑖 , which are in a dual Hilbert function space  of 𝑈 , • denotes the complex conjugation.  𝑖 ( v𝑖 , 𝑢 𝑖 ) and 𝓵( v𝑖 ) represent general bilinear and linear forms for bulk parts and boundary conditions defined on each subdomain Ω 𝑖 ∪ 𝜕Ω 𝑖 . The functional 𝓘 𝒊 ( v𝑖 , 𝑢 𝑖 ) is defined on the interface, responsible for the interface conditions.

Here, we focus on these interface terms, which will be connected to our thin layers model. The functional 𝓘 𝒊 ( v𝑖 , 𝑢 𝑖 ) for two subdomains expanded at the interface is written as:

𝓘 𝟏 ( v1 , 𝑢 1 ) + 𝓘 𝟐 ( v2 , 𝑢 2 ) = - 1 ∫ Γ v1 ũ1 d𝑆 +  2 ∫ Γ v2 ũ2 d𝑆, (2)
with ũ𝑖 given by prescribed interface conditions, which can be normally expressed as a function of the normal gradient of field 𝑢 𝑖 defined as the dual variable at the interface:

ũ𝑖 = 𝑓 (𝛁𝑢 𝑖 .𝐧), (3)
where the opposite signs in front of two above terms come from the change of normal vector at the interface Γ. Coefficients  𝑖 represent a series of constants related to the media properties. To solve eq. ( 1) within the framework of the finite element method, interface conditions need to be determined first and applied in eq. ( 2). Similar to the imposition of Neumann boundary conditions, the function of dual variables • in interface integral eq. ( 2) is commonly assigned with a prescribed function /expression. Under inspiration of Fluid-Structure Interaction (FSI), the variables on both sides of the interfaces might be expressed one from the others, resulting in coupling terms at the interaction interfaces. Therefore, we are able to find a linear relationship between primal and dual variables for subdomains in the problems of interest.

Here, this relationship is associated with specific interface conditions, which are required to represent properly the acoustic behaviour across thin layers. The function of dual variables in eq. ( 2) expressed by primal ones is written with a relationship in matrix form as:

{ ũ1 ũ2 } = [𝐀] { 𝑢 1 𝑢 2 } , ( 4
)
where [𝐀] is called the generalized admittance matrix in acoustics, since it reveals the ratio of dual and primal variables as admittance of an acoustic media. Then, substituting this equation into eq. ( 2) by replacing dual variables ũ𝑖 , the interface terms at the interface are rewritten as following:

𝓘 𝟏 ( v1 , 𝑢 1 ) + 𝓘 𝟐 ( v2 , 𝑢 2 ) = - 1 ∫ Γ v1 ( 𝐀 11 𝑢 1 + 𝐀 12 𝑢 2 ) d𝑆 +  2 ∫ Γ v2 ( 𝐀 21 𝑢 1 + 𝐀 22 𝑢 2 ) d𝑆 (5)
where the integral terms at the interface remain coupled through primal variables in each subdomain only. Once the coefficients in [𝐀] are known, the corresponding problem governed by eq. ( 1) can be solved. It is noted that the proposed generalized formulation is inherently consistent, as no extra operators such as Lagrange multiplier and supplementary conditions need to be added in the original problems.

The framework presented here can be applied to any acoustical materials such as fluid, solid and even poro-elastic materials through giving specific governing equations. The key point for the following is to determine the matrix [𝐀], whose dimension and coefficients depend on the specific problem and the adopted interface model. In this paper, an interface model based on the transfer matrix method and various specific admittance matrices are provided in section 3. In addition, the proposed interface operator is suitable for any interface model that can be rewritten in matrix form as eq. ( 4), for instance, the classical pressure drop model [START_REF] Pierce | Acoustics: An Introduction to Its Physical Principles and Applications[END_REF] . The adaption of the present framework to this model is discussed in section 4.

As mentioned in section 1, our interest is restricted to the acoustic response of porous materials that can be modelled by Helmholtz equation as equivalent fluid or Biot's equations [START_REF] Biot | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low-Frequency Range[END_REF][START_REF] Biot | Mechanics of Deformation and Acoustic Propagation in Porous Media[END_REF] as PEMs. The variational formulations of these governing equations are introduced following.

Variational formulations for porous materials

Helmholtz equation controls the acoustic behaviour in porous materials when the media are considered as equivalent fluids. The weak formulation of this governing equation for fluid pressure 𝑝 𝑓 with its test function q𝑓 consists of bilinear forms  on the bulk domain and linear form 𝓘, 𝓵 at the interfaces and boundaries as:

1 𝜔 2 ρeq ∫ Ω ∇ q𝑓 .∇𝑝 𝑓 dΩ - 1 Keq ∫ Ω q𝑓 𝑝 𝑓 dΩ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟  -∫ Γ q𝑓 𝑢 𝑡 𝑛 d𝑆 ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟ 𝓘 = ∫ 𝜕Ω q𝑓 𝑢 𝑡 𝑛 d𝑆 ⏟⏞⏞⏞⏟⏞⏞⏞⏟ 𝓵 , ( 6
)
with total displacement of fluid 𝑢 𝑡 𝑛 defined at the boundary/interface that is a function of gradient of pressure as the ũ in eq. ( 3). Coefficients ρeq and Keq denote effective density and compressibility of the fluid respectively, which are complex numbers and frequency dependent values for porous materials as given by JCA model [START_REF] Johnson | Theory of Dynamic Permeability and Tortuosity in Fluid-Saturated Porous Media[END_REF][START_REF] Champoux | Dynamic Tortuosity and Bulk Modulus in Air-saturated Porous Media[END_REF] . This equation can also be used to solve classical acoustic fluids such as air, but with real-valued material constants. Here, the coefficient 1∕ ρeq corresponds to the material's parameter  in the previous general interface operator 5.

The mixed 𝑢 𝑠 -𝑝 𝑓 Biot's equations proposed by Atalla et al. [START_REF] Atalla | A Mixed Displacement-Pressure Formulation for Poroelastic Materials[END_REF] is adopted in this work for PEMs. This formulation comprises two dynamic equations for fluid and solid phases, which are connected by a coupling term. Weak formulations are derived by multiplying by test function v𝑠 for solid displacement and q𝑓 for poro pressure and applying Green's formula:

∫ Ω σ𝑠 𝑖𝑗 ( v𝑠 )𝜖 𝑠 𝑖𝑗 (𝐮 𝑠 ) dΩ -𝜔 2 ∫ Ω ρ v𝑠 𝑖 𝐮 𝑠 𝑖 dΩ -∫ Ω γ v𝑠 𝑖 ∇𝑝 𝑓 dΩ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟  -∫ Γ v𝑠 𝑖 𝒏 𝑗 σ𝑠 𝑖𝑗 (𝐮 𝑠 ) d𝑆 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝓘 = ∫ 𝜕Ω v𝑠 𝑖 𝒏 𝑗 σ𝑠 𝑖𝑗 (𝐮 𝑠 ) d𝑆 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝓵 , ∫ Ω [ 1 𝜔 2 ρeq ∇ q𝑓 .∇𝑝 𝑓 - 1 Keq q𝑓 𝑝 𝑓 ] dΩ -∫ Ω γ∇ q𝑓 𝐮 𝑠 𝑖 dΩ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟  -∫ Γ q𝑓 𝑢 𝑡 𝑛 d𝑆 ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟ 𝓘 = ∫ 𝜕Ω q𝑓 𝑢 𝑡 𝑛 d𝑆 ⏟⏞⏞⏞⏟⏞⏞⏞⏟ 𝓵 , (7a) (7b)
where σ𝑠 represents the in-vacuo stress tensor of the solid phase and 𝜖 𝑠 is the associated strain tensor. Coefficients ρ and γ are the effective density of the solid phase and coupling parameter, respectively. The additional information concerning all these parameters is detailed in appendix A. The coefficients written in eq. ( 7) are slightly modified with respect to the original paper in order to facilitate the coupling with other media. Compared to the equivalent fluid model, three distinct waves exist in Biot's theory: one compressional airborne wave, one compressional and one shear frame-borne waves, resulting in three different wave numbers.

𝑝 𝑓 , 𝑢 𝑠 and 𝑢 𝑡 𝑛 , σ in interface functional 𝓘 are the primal and dual variables of governing equations eq. ( 6) and eq. ( 7) for porous materials. We then need to investigate expressions/relationships between these dual and primal variables in form of eq. ( 4). These relationships are related to our interface model for the thin layers and coupling configurations that are established in the next section.

APPLICATION TO SOUND ABSORPTION SYSTEMS

In this section, the methodology presented in section 2 is applied to deal with waves propagating from acoustic fluid (air) to a porous absorption system composed of thin films and bulk porous part. Both bulk media and thin films are modelled by two types of governing equations: Helmholtz and Biot's equations, it results in four different coupling configurations and therein four interface conditions. We first proceed to the interface model for thin acoustic media. Then, the associated interface conditions are derived.

As shown in fig. 2 (a,b,c,d), the interface conditions are described for a bidimensional geometry, in the incident (𝑥, 𝑦) plane. To facilitate the derivation of these conditions between sides Γ 1 and Γ 2 (Γ + 1 and Γ + 2 in the figures), related quentities inside the thin layer are assigned with the superscript (•) -, and the ones outsides attached with bulk medium is given by (•) + .

Interface model for thin layers

Sound propagation represented by Transfer Matrix Method (TMM)

The transfer matrix method is a frequently used method for wave propagating in laterally infinite stratified structures, where the solution of governing equations are represented by a matrix formulation. This method establishes a relationship of the state vectors 𝐬 between two sides of the layers by a so-called transfer matrix [𝐓] as:

𝐬(𝐱) = [𝐓]𝐬(𝐱 ′ ) (8)
where coefficients in [𝐓] are obtained under the assumption of plane wave excitation, and are functions of the incident angle. Consequently, the accuracy of TMM depends on the incident angle, which will be addressed in detail in section 4. 𝐱 and 𝐱 ′ are coordinates close to the forward and the backward faces of the layers, and the state vector 𝐬 includes variables to represent the acoustic field:

𝐬(𝐱) = {𝑢(𝐱), ũ(𝐱)} 𝑇 , (9)
indeed, 𝑢 and ũ correspond to the primal and dual variables of modelled media, respectively. In two-dimensional problems, the number of variables in the state vector is equal to 2m (m being the number of waves, that depends on the nature of media). A detailed description of the state vectors for porous materials follows in section 3.1.2.

It is worth noticing that the state vector contains simultaneously the primal and dual variables on the two faces of a layer, connected by a matrix relationship: this is exactly what we are looking for when constructing an interface condition associated with thin layers. In addition, the relationship given in eq. ( 8) satisfies the governing equations, no simplification is made in this method unlike with the asymptotic model proposed by Bövik [START_REF] Bövik | On the modeing of thin interface layers in elastic and acoustic scattering problems[END_REF] or the pressure drop model [START_REF] Pierce | Acoustics: An Introduction to Its Physical Principles and Applications[END_REF] . Furthermore, TMM is also favourable to treat multi-layer problems, as the global transfer matrix for a multi-layer system can be obtained directly by multiplying the matrices associated with each layer when they are in the same nature. This benefit will be illustrated through a numerical example in section 6.

TMM for porous thin layers

The porous thin layers in the problems stated above are considered as stratified plane layers, since their thickness is much smaller than the other dimensions. Thus, the TMM is used to establish our interface model for thin films. Concerning the TMM for porous materials, two different state vectors and transfer matrices are specified here for our two governing equations.

The state vector 𝐬 𝑓 for materials modelled as fluid has two components, since only one wave propagates in the media. As shown in Fig. 2(a,c), pressure 𝑝 𝑓 and normal (𝑥 component in our geometric representation) total fluid displacement/velocity 𝑢 𝑡 𝑛 of fluid at the interface are the primal and dual variables:

𝐬 𝑓 (𝐱) = {𝑝 𝑓 (𝐱), 𝑢 𝑡 𝑥(𝑛) (𝐱)} 𝑇 , (10)
Thus, we would obtain a 2 × 2 transfer matrix [𝐓 𝐟 ] for layers modelled as a fluid. For PEMs governed by Biot's equations, there exists three waves. Six quantities are thus required to represent the solution in the layers, as shown in Fig. 2(b,d). The state vector for the considered PEMs 𝐬 𝑃 (𝐱) reads:

𝐬 𝑃 (𝐱) = [ σ𝑦𝑥 (𝐱), 𝑢 𝑠 𝑥 (𝐱), 𝑢 𝑡 𝑥(𝑛) (𝐱), σ𝑥𝑥 (𝐱), 𝑝 𝑓 (𝐱), 𝑢 𝑠 𝑦 (𝐱)] 𝑇 , ( 11
)
where 𝑢 𝑠 𝑥 and 𝑢 𝑠 𝑦 are the two components of the displacement of the elastic frame which are primal variables, σ𝑦𝑥 and σ𝑥𝑥 correspond to the normal and shear in-vacuo stresses of the solid phase denoting the dual variables in PEMs. Together with fluid pressure and normal total displacement, three primal and three dual variables are retrieved. With this state vector, the transfer matrix for the PEMs layer is of size 6 × 6 and is denoted as

[𝐓 𝐏 ].
It is seen that the dimension of transfer matrices is decided by the nature of thee media, while the coefficients and their values in transfer matrices [𝐓 𝐟 ] and [𝐓 𝐏 ] depend on the specific material model. In this paper, the JCA model [START_REF] Johnson | Theory of Dynamic Permeability and Tortuosity in Fluid-Saturated Porous Media[END_REF][START_REF] Champoux | Dynamic Tortuosity and Bulk Modulus in Air-saturated Porous Media[END_REF] , one variant of the Limp model, see [START_REF] Panneton | Comments on the Limp Frame Equivalent Fluid Model for Porous Media[END_REF] and a Biot matrix representation [START_REF] Martinez | Derivation of the State Matrix for Dynamic Analysis of Linear Homogeneous Media[END_REF] are considered. We specify and discuss the coefficients in these models in section 4. In fact, the concrete coefficients have no effect on the following derivation of interface conditions, but the variables in the state vectors do.

As expected, the variables for fluid media and PEMs in their state vectors (10) and (11) correspond exactly to the dual and primal variables in the variational formulations eqs. ( 6) and (7). This coherence facilitates the deduction and enforcement of the interface conditions to be presented in the following.

Up to this moment, we established the relationship between Γ - 1 and Γ - 2 in the thin porous layers. It needs to be combined with the interface conditions attached to the bulk subdomain to obtain the relationship between faces Γ + 1 and Γ + 2 , which is related to the coupling configuration.

Generalized interface conditions

This section presents four considered interface conditions generated from the four coupling configurations shown in fig. 2. To make all related expressions slightly simpler, a compact notation is used for variables at the interface, for example:

𝑝 𝑓 (Γ - 1 ) ∶ = 𝑝 𝑓 - 1 , 𝑝 𝑓 (Γ + 2 ) ∶= 𝑝 𝑓 + 2 σ𝑥𝑥 (Γ - 1 ) ∶ = σ- 𝑥𝑥1 , σ𝑦𝑥 (Γ + 2 ) ∶= σ+ 𝑦𝑥2 𝑢 𝑡 𝑛 (Γ - 1 ) ∶ = 𝑢 𝑡 𝑛 - 1 , 𝑢 𝑠 𝑥 (Γ + 2 ) ∶= 𝑢 𝑠 𝑥 + 2

Fluid-Fluid coupling with fluid layer (f-f-f)

When the film and bulk porous materials are treated as equivalent fluids as Fig. 2(a), pressure and normal displacement along the interface are the only primal and dual variables. With the conservation of pressure and total fluid displacement across two faces, we obtain:

{ 𝑝 𝑓 - 𝑖 = 𝑝 𝑓 + 𝑖 𝑢 𝑡 𝑛 - 𝑖 = 𝑢 𝑡 𝑛 + 𝑖 , 𝑖 = 1, 2 (12a) (12b)
Thus, the relationship between Γ + 1 and Γ + 2 is written as:

{ 𝑝 𝑓 + 1 𝑢 𝑡 𝑛 + 1 } = [ 𝐓 f-f-f ] 2×2 { 𝑝 𝑓 + 2 𝑢 𝑡 𝑛 + 2 } , ( 13
)
where coefficients in 𝑇 f-f-f 𝑖𝑗 are the same as in the aforementioned transfer matrix for the fluid modelled layer [𝐓 𝐟 ]. The relation exposed by eq. ( 13) and following deduced relations eqs. (15), (17) and (20) are the ones we call generalized interface conditions. Both pressure and total displacement (function of the gradient of the pressure) are potentially discontinuous through the interfaces. The discontinuity (jump) for these variables (primal and dual one) are controlled by every coefficient in the associated transfer matrices.

Fluid-Fluid coupling with PEMs layer (f-P-f)

Compared to the previous case, PEMs controlled by Biot's equations involves 6 variables instead of two due to the solid motion, as illustrated in Fig. 2(b). The coupling of PEMs with fluid media results in a free surface on the elastic frame sides, a free surface condition is thus added on both sides inside the film layer as:

σ- 𝑦𝑥𝑖 = 0 and σ- 𝑥𝑥𝑖 = 0, 𝑖 = 1, 2 (14)
These conditions do not change the structure of the generalized conditions between Γ + 1 with Γ + 2 since we still have two variables (𝑝 𝑓 and 𝑢 𝑡 𝑛 ) in the state vector of the two bulk media:

{ 𝑝 𝑓 + 1 𝑢 𝑡 𝑛 + 1 } = [ 𝐓 f-P-f ] 2×2 { 𝑝 𝑓 + 2 𝑢 𝑡 𝑛 + 2 } , ( 15
)
where the coefficients in the above matrix 𝐓 f-P-f are function of the coefficients in the transfer matrix for Biot's layer [𝐓 𝐏 ].

Even though the relationship between the two sides of the interface remains similar to the previous case (same state vector with a 2 × 2 matrix), the coefficients in the matrix are much more complicated as they are obtained by introducing the free surface conditions in [𝐓 𝐏 ] (the procedure for obtaining matrix 𝐓 f-P-f is given in appendix B.1). Note that it is more common to formulate the transfer matrix condition for fluid-fluid coupling using pressure and total velocity, where total velocity converts to displacement as 𝑣 𝑡 𝑛 = j𝜔𝑢 𝑡 𝑛 under the harmonic convention chosen in this paper. The interface conditions of "f-f-f" and "f-P-f" can be classified into one type, since only Helmholtz equations are involved for two bulk media associated by a 2 × 2 matrix at the interface. However, when porous bulk media Ω 2 is considered as full PEMs governed by Biot's equations, the number of variables is not compatible between the two sides of the film (2 against 6 when coupling with an acoustic fluid). As a result, non-square matrix relationships for the interface will be observed.

Fluid-PEMs coupling with fluid layer (f-f-P)

This configuration is described in fig. 2(c). Contrary to the previous case, the free surface condition on the stresses are only applied on the side of the PEM bulk subdomain Γ + 2 :

σ+ 𝑦𝑥2 = 0, and σ+ 𝑥𝑥2 = 0 (16)
As this condition is on the PEMs' side, it is required to be enforced as an interface condition in the equation of solid phase. Then, the generalized conditions in this case are cast as:

{ 𝑝 𝑓 + 1 𝑢 𝑡 𝑛 + 1 } = [ 𝐓 f-f-P ] 2×4 ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑝 𝑓 + 2 𝑢 𝑡 𝑛 + 2 σ+ 𝑥𝑥2 = 0 σ+ 𝑦𝑥2 = 0 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ , ( 17
)
where the coefficients associated with 𝑝 𝑓 + 2 and 𝑢 𝑡 𝑛 + 2 in the matrix above are the same as in eq. ( 13), which is because these two cases both employ the equivalent fluid model for the embedded thin film.

Fluid-PEMs couling with PEM layer (f-P-P)

Replacing the fluid model in the previous case by a Biot's medium for the film, we obtain an interface between two PEMs as shown in fig. 2(d). The continuity of all the variables are retained at the interface Γ 2 as:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ σ- 𝑦𝑥2 = σ+ 𝑦𝑥2 and σ- 𝑦𝑥2 = σ+ 𝑥𝑥2 , 𝑢 𝑠 𝑥 - 2 = 𝑢 𝑠 𝑥 + 2 and 𝑢 𝑠 𝑦 - 2 = 𝑢 𝑠 𝑦 + 2 , 𝑢 𝑡 𝑛 - 2 = 𝑢 𝑡 𝑛 + 2 , 𝑝 𝑓 - 2 = 𝑝 𝑓 + 2 , ( 18a
) (18b) (18c) (18d)
and a free surface condition on Γ - 1 :

σ- 𝑦𝑥1 = 0 and σ- 𝑥𝑥1 = 0 (19)
Consequently, the generalized conditions that will be prescribed in the variational formulation are:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ σ- 𝑦𝑥1 = 0 σ- 𝑥𝑥1 = 0 𝑝 𝑓 + 1 𝑢 𝑡 𝑛 + 1 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = [ 𝐓 f-P-P ] 4×6 ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ σ+ 𝑦𝑥2 𝑢 𝑠 𝑥 + 2 𝑢 𝑡 𝑛 + 2 σ+ 𝑥𝑥2 𝑝 𝑓 + 2 𝑢 𝑠 𝑦 + 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ , ( 20
)
where [𝐓 f-P-P ] is a 4 × 6 matrix, the two lines related to the displacement

𝑢 𝑠 𝑥 - 1 and 𝑢 𝑠 𝑦 - 1
in the original transfer matrix [𝐓 𝐏 ] are eliminated due to the fact that the fluid media does not have a solid phase displacement.

For these four coupling types, three different forms are derived for interface conditions. They will be imposed weakly in the variational formulations of bulk subdomains.

Generalized interface operators in variational formulations

Now, we impose the generalized interface conditions obtained previously in the coupled variational formulations eq. ( 6) and eq. ( 7). Three interface conditions types are enforced: one interface conditions for two fluid coupling and two for fluid-PEMs coupling. As only the variables on surfaces Γ + 𝑖 are retained, the superscript (•) + is removed to simplify following expressions.

Fluid-fluid coupling

In addition to the two formulations eq. ( 6) for media 1 and 2, coupling interface integrals on Γ 1 and Γ 2 are prescribed:

𝓘 f -f = -∫ Γ 1 q𝑓 1 𝑢 𝑡 𝑛1 d𝑆 + ∫ Γ 2 q𝑓 2 𝑢 𝑡 𝑛2 d𝑆, (21)
where 𝑞 𝑓 𝑖 are the test function of fluid pressure in each subdomain. The generalized interface conditions derived from TMM relations can not be used directly, admittance relationship with matrix [𝐀] needs first to be obtained. Relationships expressed with transfer matrices in eq. ( 13) and eq. ( 15) need to be rewritten with a generalized impedance matrix [𝐀 𝐟 𝐟 ] as:

{ 𝑢 𝑡 𝑛1 𝑢 𝑡 𝑛2 } = [ 𝐴 f f 11 𝐴 f f 12 𝐴 f f 21 𝐴 f f 22 ] { 𝑝 𝑓 1 𝑝 𝑓 2 } , ( 22
)
where the coefficients in matrix [𝐀 𝐟 𝐟 ] are expressed by the coefficients in

[𝐓] ([𝐓 f-f-f ] and [𝐓 f-P-f ]): 𝐴 f f 11 = 𝑇 22 𝑇 12 , 𝐴 f f 12 = 𝑇 12 𝑇 21 -𝑇 22 𝑇 11 𝑇 12 , 𝐴 f f 21 = 1 𝑇 12 , 𝐴 f f 22 = - 𝑇 11 𝑇 12 (23)
Substituting eq. ( 22) into eq. ( 21) and combining the interface Γ 1 and Γ 2 as Γ, the final interface operator is obtained in a compact form as following:

𝓘 f -f = 𝑛 ∑ 𝑖=1 (-1) 𝑖 𝐴 f f 𝑖𝑗 ∫ Γ q𝑓 𝑖 𝑝 𝑓 𝑗 d𝑆, 𝑛 = 2 (24)
Remind that configurations "f-f-f" and "f-P-f" lead to the same interface operator but only with different coefficients 𝐴 f f 𝑖𝑗 .

Fluid-PEMs coupling

Interface terms for fluid-PEMs coupling are obtained by adding eq. ( 6) for media 1 and eq. ( 7) for media 2 together as:

𝓘 f -P = -∫ Γ 1 q𝑓 1 𝑢 𝑡 𝑛1 d𝑆 + ∫ Γ 2 q𝑓 2 𝑢 𝑡 𝑛2 d𝑆 + ∫ Γ 2 v𝑠 𝑖 𝒏 𝑗 σ𝑠 𝑖𝑗 (𝐮 𝑠 ) d𝑆 (25)
Compared to fluid-fluid coupling, a stress surface term appears due to the dynamic equation of the elastic solid phase.

For configuration "f-f-P", prescribing interface conditions eq. ( 17) directly conduces the stress terms to vanish. Then, the remaining conditions are the same as fluid-fluid coupling eq. ( 22), resulting in the following form as eq. (24).

𝓘 f -f -P = 𝑛 ∑ 𝑖=1 (-1) 𝑖 𝐴 f f 𝑖𝑗 ∫ Γ q𝑓 𝑖 𝑝 𝑓 𝑗 d𝑆, 𝑛 = 2, (26)
it can be seen that we obtain the same interface terms as in the previous case. This is because only pressure fields are involved in the interface conditions when the film is modelled by an equivalent fluid.

To prescribe conditions eq. ( 20) of configuration "f-P-P", as previously, we rewrite it with an admittance matrix [𝐀 𝐏 ] to express dual variables by means of primal variables as following:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑢 𝑡 𝑛1 𝑢 𝑡 𝑛2 σ𝑥𝑥2 σ𝑦𝑥2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = [𝐀 𝐏 ] 4×4 ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑝 𝑓 1 𝑝 𝑓 2 𝑢 𝑠 𝑥2 𝑢 𝑠 𝑦 2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ , ( 27
)
where [𝐀 𝐏 ] is not easy to be obtained from [𝐓 𝐟 -𝐏-𝐏 ]. The procedure for deriving this admittance matrix from [𝐓 𝐟 -𝐏-𝐏 ] is detailed in appendix B.2. Once the matrix [𝐀 𝐏 ] is determined, the interface terms in eq. ( 25) are rewritten by substituting eq. ( 27) as the following compact form:

𝓘 f -P-P = ∫ Γ [-q𝑓 1 , q𝑓 2 , v𝑠 𝑥 , v𝑠 𝑦 ][𝐀 𝐏 ][𝑝 𝑓 1 , 𝑝 𝑓 2 , 𝑣 𝑠 𝑥 , 𝑣 𝑠 𝑦 ] 𝑇 d𝑆, (28)
where the negative sign in front of q𝑓 1 still stems from the choice of outward normal at the interface as with previous couplings. It is worth emphasizing that all these interface operators ( 24), ( 26) and ( 28) introduce strong discontinuities for both primal and dual variables, which is non-trivial as an implementation point of view, for conventional FEM. The discretization for such particular problems will be presented in section 5.

Symmetry of variational systems

The variational formulation with generalized interface operators yields a linear system after the finite element discretization, to be presented in section 5. The structure of the resulting systems is worth being discussed, as it has an impact on the choice of linear solver that needs to be considered for the following. The symmetry could be identified directly through the variational formulations stated in this section in continuous space.

Conventionally, the bulk parts render the symmetry to the system, while the proposed generalized interface operators may bring some non-symmetry. From variational interface operator (24), it can be recognized that the fluid-fluid coupling system is symmetric if [𝐀 𝐟𝐟 ] 𝟐×𝟐 satisfies:

𝐴 f f 𝑖𝑗 = (-1)𝐴 f f 𝑗𝑖 𝑖, 𝑗 ∈ 1, 2, (29a)
which results in an anti-symmetric matrix for fluid-fluid coupling. With the conversion from transfer matrix [𝐓] to impedance matrix [𝐀] for this coupling as illustrated in eq. ( 23), the anti-symmetry of

[𝐀 𝐟 𝐟 ] (𝐴 f f 12 = -𝐴 f f 21 ) gives: 𝑇 12 𝑇 21 -𝑇 22 𝑇 11 𝑇 12 = - 1 𝑇 12 , ( 30
)
that leads to the coefficients in [𝐓] satisfying :

𝑇 22 𝑇 11 -𝑇 12 𝑇 21 = 1, 𝑇 12 ≠ 0 (31)
which means that the resulting linear system will be symmetric when the determinant of the corresponding transfer matrix [𝐓] is equal to one. The same conclusion applies for operator (26) of the configuration "f-f-P". When it comes to configuration "f-P-P", symmetry requires that coefficients in the matrix [𝐀 𝐏 ] 𝟒×𝟒 for interface operator (28) fulfil:

{ 𝐴 𝑃 𝑖𝑗 = (-1)𝐴 𝑃 𝑗𝑖 when 𝑖 = 1 𝐴 𝑃 𝑖𝑗 = 𝐴 𝑃 𝑗𝑖 , otherwise , (32a) (32b)
which results in an anti-symmetry for the terms associated with 𝑝 𝑓 1 . In this case, the link between transfer matrix in porous film [𝐓 𝐏 ] and admittance matrix [𝐀 𝐏 ] is not direct. Their relationship concerns the reconstruction and inverse of matrices as illustrated in appendix B.2, depending on the specific problem. We cannot give an analytical condition for symmetry of [𝐓 𝐏 ] as eq. ( 31). However, [𝐀 𝐏 ] can be evaluated numerically a posteriori to verify if the corresponding system is symmetric, in order to choose a suitable linear solver.

DISCUSSION ON POROUS THIN LAYERS

This section gives more specific details associated with thin porous materials in time-harmonic problems. We introduce three commonly used transfer matrices with concrete coefficients for porous media and discuss their influence on the proposed formulations. Furthermore, we compare our model to the simplified pressure drop model. At last, to illustrate the aforementioned effects from the assumption of plane wave in TMM, the proposed model is evaluated by a set of sensitivity tests. In this paper, two representative woven and non-woven films already investigated in a previous work [START_REF] Gaborit | A Simplified Model for Thin Acoustic Screens[END_REF] are both considered. The properties of the two films that will be used there-under are listed in table 1.

Transfer matrix method for porous films 4.1.1 JCA model

When porous material is considered as fluid, the classical transfer matrix for fluid layers [START_REF] Allard | Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials[END_REF] can be used directly but with complexvalued JCA properties [START_REF] Johnson | Theory of Dynamic Permeability and Tortuosity in Fluid-Saturated Porous Media[END_REF][START_REF] Champoux | Dynamic Tortuosity and Bulk Modulus in Air-saturated Porous Media[END_REF] . The matrix associated with the state vector of pressure 𝑝 𝑓 and total velocity 𝑣 𝑡 for a layer of thickness 𝑑 is written as:

[𝐓 𝐟 ] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ cos ( 𝑘 eq 𝑑 ) j 𝜔 ρeq 𝑘 eq sin ( 𝑘 eq𝑥 𝑑 ) j 𝑘 eq𝑥 𝜔 ρeq sin ( 𝑘 eq𝑥 𝑑 ) cos ( 𝑘 eq𝑥 𝑑 ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 33
)
where 𝑘 eq𝑥 is the equivalent wave number of component 𝑥 (normal component at the interface). As seen that, the thickness of media 𝑑 is considered as well in the matrix. It is easy to verify that the determinant of this transfer matrix |𝐓 𝐟 | = 1, which satisfies the condition for a symmetric linear system prescribed previously eq. (31).

Limp model

When a porous film is attached to a PEM that is controlled by Biot's equations eq. ( 7), the JCA parameters in the equivalent fluid model are not adequate. Therefore, a modified JCA model considering the inertia of the elastic frame is introduced, so-call Limp model [START_REF] Panneton | Comments on the Limp Frame Equivalent Fluid Model for Porous Media[END_REF] where the density of porous materials in Limp model is defined as:

ρLimp = ρ ρeq ρ + ρeq γ2 , ( 34
)
where γ is the coupling parameter which appears in Biot's equations eq. ( 7). The Limp model does not change the structure and characteristics of the transfer matrix (33) except that the ρeq terms are all replaced by ρLimp . Thus, a symmetric linear system is recovered as well for this model.

Biot model

Biot's model considers the interaction between fluid and elastic solid phases in the material, which leads to a larger transfer matrix than the two fluid models presented before. The transfer matrix for Biot's model [𝐓 𝐏 ] associated with the state vector (11) in this paper is extracted from Dazel et al. [START_REF] Dazel | A Stable Method to Model the Acoustic Response of Multilayered Structures[END_REF] which gives a complete description for stratified PEMs layers in two dimensions:

[𝐓 𝐏 ] = exp(-𝑑𝛂), (35)
with state matrix 𝛂 specified as:

𝛂 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 j𝑘 𝑦 Â P j𝑘 𝑦 γ - Â2 -p2 p 𝑘 2 𝑦 -ρ𝜔 2 0 0 0 1 P 0 j𝑘 𝑦 Â P 0 0 0 0 -1 Keq + 𝑘 2 𝑦 ρeq 𝜔 2 -j𝑘 𝑦 γ j𝑘 𝑦 -ρs 𝜔 2 -ρeq γ𝜔 2 0 0 0 0 ρeq γ𝜔 2 ρeq 𝜔 2 0 0 0 1∕ N j𝑘 𝑦 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 36
)
where 𝑘 𝑦 represents a vertical wave number in the layer under the dependence of 𝑒 -j𝑘 𝑦 𝑦 . Other materials' parameters present in the matrix are associated with JCA models and Biot's theory (the detailed presentation of these parameters can be found in appendix A or in the reference [START_REF] Allard | Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials[END_REF] ). Compared to the previous fluid models, the transfer matrix for Biot's layers is much more complicated, even an exponential of matrix 𝛂 is involved. As mentioned before, it is difficult to give an a priori analysis for such matrices and the associated admittance matrices in the interface model. The corresponding impedance matrix [𝐀 𝐏 ] derived from eq. ( 35) has been verified numerically as not fulfilling the symmetry requirement eq. ( 32) (the transformation procedure from [𝐀 𝐏 ] to [𝐓 𝐏 ] is detailed in appendix B.2). As a result, solvers handling non-symmetric systems must be used in this case. Notice that there exists other Biot's transfer matrices such as a simplified representation [START_REF] Gaborit | A Simplified Model for Thin Acoustic Screens[END_REF] , the symmetry of the resulting system needs to be verified a posteriori for each model.

Comparison with the simplified pressure drop model

The simplified pressure drop model is frequently used to model thin layers in acoustic problems. When the fluid passes across a sample of material, a pressure difference occurs due to the flow resistance. This difference is proportional to the mean normal flow velocity and material resistivity as:

𝑝 𝑓 1 -𝑝 𝑓 2 = 𝜎𝑑 v, ( 37
)
where 𝜎, 𝑑 and v are the flow resistivity of the materials, material thickness and mean flow velocity, respectively. As layers considered here are very thin, the total velocity is assumed to be conserved across layers as:

𝑣 𝑡 𝑛1 = 𝑣 𝑡 𝑛2 = v𝑛 ( 38
)
Combining eq. ( 37) and eq. ( 38), the simplified pressure drop model is rewritten in a matrix form as:

{ 𝑝 𝑓 1 𝑣 𝑡 𝑛1 } = [ 1 𝜎𝑑 0 1 ] { 𝑝 𝑓 2 𝑣 𝑡 𝑛2 } ( 39
)
As it is expressed in the same form as general interface conditions with a transfer matrix, the simplified relationship can also be imposed weakly using the proposed framework. The variational interface operator corresponding to eq. ( 39) in fluid-fluid coupling eq. ( 21) is written as:

𝓘 = j 𝜎𝑑𝜔 ∫ Γ q𝑓 𝑝 𝑓 d𝑆, (40)
where • corresponds to the jump operator, that is 𝑝 𝑓 = 𝑝 𝑓 1 -𝑝 𝑓 2 and q𝑓 = 𝑞 𝑓 1 -𝑞 𝑓 2 . This interface operator is similar to the spring-like interface in composite, as investigated in [START_REF] Zhu | Three-Dimensional Numerical Modelling by XFEM of Spring-Layer Imperfect Curved Interfaces with Applications to Linearly Elastic Composite Materials[END_REF] . In fact, the relationship exposed in eq. ( 39) is only suitable for the low frequency limit (𝜔 ←→ 0). Its flow resistance is frequency independent, which is in contrast with the aforementioned models for porous films. To identify the difference between each thin layer model, we calculate the evolution of the coefficient that corresponds to the pressure drop (𝑇 𝑓 12 in matrices ( 33) and ( 34) for two equivalent fluid models, 𝑇 f-P-f 12 in eq. ( 15) for Biot's model and 𝜎𝑑 in eq. ( 39)) at the interface with an increasing frequency for the two films presented above. The variable 𝜃 = 0 °(Normal incidence) and thickness 𝑑 = 1 mm are chosen in these related matrices.

From fig. 3, the coefficients from the simplified pressure drop model stay constant with respect to frequency. The imaginary part is zero for both types of film, which corresponds to the fact that no dissipation is considered in this model. In contrast, the other models lead to a change in the pressure drop coefficient when increasing frequency, especially the imaginary part. At low frequencies, the imaginary part of the JCA model tends to zero. When taking the vibration of the elastic frame into account (that is Limp and Biot's models), an obvious gap in the imaginary part is noticeable compared with the simplified model, even in the low-frequency range. This behaviour confirms that the simplified pressure drop model is only valid in the limit of low frequency for porous materials. On the other hand, the Limp and Biot models exhibit a similar behaviour regarding the pressure drop, even in the medium-high frequency range where a superposition of the curves is observed. This is because with the large density of the films considered here, as seen in table 1, the displacement of the elastic frame considered in Biot's model has minor effects on the poro-fluid behaviour.

To conclude, the proposed general framework presented in section 2 can be readily applied to the pressure drop model, yielding the interface operator eq. ( 40). On the other hand, we illustrate that this simplified interface model is not accurate enough to model thin porous films in a wide range of frequency. A more comprehensive model with the generalized interface operator proposed in this paper is necessary to be used to account for thin resistive porous films. Furthermore, if we are interested in the pressure field of the solution only, which is the case in most acoustic analysis, the Limp model is preferable to Biot's model. Indeed, the Limp model is able to provide a similar and sufficiently accurate description of Biot's equations to model thin films in TMM. Additionally, it results in a symmetric linear system, which will accelerate the computational procedure.

Sensitivity of transfer matrix model for thin films

We demonstrated the necessarity of using complete transfer matrices as eq. ( 33)eq. ( 34) and eq. ( 35) rather than the simplified pressure-drop model to represent porous materials. However, as aforementioned, these transfer matrices are formulated under the assumption of plane wave propagation. This assumption is difficult to be ensured in the general case, especially when the problem is excited by a point source where waves propagate in all directions or if the geometry of the interface is complex. This is due to the fact that a mismatch in the incident angle 𝜃 will result in an incorrect wave number and effective propagation To address this aspect in our specific case where TMM represents relatively thin porous layers, the solution of a classical transfer matrix without any simplification to be examined. Here, we carry out a sensitivity study with respect to the incidence angle 𝜃.

The test case is performed in a one-dimensional setting as illustrated in fig. 4 where a horizontal incident plane wave propagates in an infinite plane of air (with wave number 𝑘 𝑎 ) in which a layer modelled by TMM is embedded in the middle (𝑥 = 0). Eqivalent fluid models with JCA parameters eq. ( 33) are used for the transfer matrix of the middle layer. To quantify the sensitivity of the interface model with respect to incidence angle, we compare the TMM solution with varying angle to the exact solution with 𝜃 = 0. The varying angle changes the exact wave number 𝑘 and thickness 𝑑 of thin layers to 𝑘 ′ and 𝑑 ′ as:

𝑘 ′ = √ 𝑘 2 -𝑘 2 a sin(𝜃) 𝑑 ′ = 𝑑∕ cos(𝜃) (41a) (41b)
The solution for this problem is expressed with incidence 𝑝 𝑖 and transmission pressure 𝑝 𝑡 :

𝑝 𝑖 = exp ( -j𝑘 a 𝑥 ) + 𝑅 exp ( j𝑘 a 𝑥 ) 𝑥 ≤ 0 𝑝 𝑡 = 𝑇 exp ( -j𝑘 a 𝑥 ) 𝑥 > 0, (42a) (42b)
where the reflection and transmission coefficients 𝑅,𝑇 are obtained through solving the generalized interface conditions eq. ( 13) associated with the transfer matrix (33). The error against the reference solution is calculated numerically in 𝐿 2 norm as

𝜖 = ∫ 𝑥≤0 | | | (𝑅 -𝑅 𝜃=0 ) exp ( j𝑘 a 𝑥 ) | | | 2 d𝑥 + ∫ 𝑥>0 | | | (𝑇 -𝑇 𝜃=0 ) exp ( -j𝑘 a 𝑥 ) | | | 2 d𝑥 = | | | 𝑅 -𝑅 𝜃=0 | | | 2 + | | | 𝑇 -𝑇 𝜃=0 | | | 2 (43)
As can be seen, the error depends on 𝑅&𝑇 which are functions of incident angle 𝜃 in the transfer matrix model for middle thin layer.

Three TMMs for air gap, woven and non-woven films with different 𝑘𝑑 are assessed. 𝑘𝑑 characterizes the number of waves propagating through the layer.

The relative error as a function of 𝜃 through sweeping 𝑘𝑑 are plotted in fig. 5. The largest 𝑘𝑑 in the figure corresponds to a wave propagating in a layer of 2 mm under an excitation of 5000 Hz, which could represent an upper bound for the use of thin films in a practical situation.

When the TMM models an air gap (fig. 5(a)), the more waves propagate in the layers (larger 𝑘𝑑), the more sensitive to the angle mismatch 𝜃 of the solution (larger error in same angle mismatch). At 90 °, a drastic increase of the error is observed when 𝑘𝑑 is small and the maximal error occurs for the highest 𝑘𝑑. These behaviours are consistent with the fact that when the layer is extremely thin, the TMM has obvious influence only for a grazing incidence. And, when 𝑘𝑑 becomes larger, more error from the mismatch is accumulated in the TMM layer. Remind that an air gap is a real-valued acoustic media without energy loss.

Fig. 5(b) and Fig. 5(c) exhibit the same assessments for porous films, indicating that the sensitivity against angle mismatch is less dependent of 𝑘𝑑. As seen in table 1, the two porous films have large flow resistivity and solid densities, far more rigid than air. This kind of materials decouple the infinite plane and produce a non-negligible pressure difference at the interface even when 𝑘𝑑 is small. Moreover, energy damping (complex-valued properties) exists in the porous films. When the layer is resistant enough or when the pressure dissipation in the layer is large enough beyond specific values, the error will be barely influenced by the angle mismatch. This is why for a fixed 𝑘𝑑, the error for the more resistant woven film increases more slowly with respect to change of angle 𝜃 than the non-woven one. Nevertheless, for both considered films, the error stays mostly below 20% even until an angle mismatch of 70 °. The angle mismatch for a target error below for instance 1% has a margin of at most 20 °for both films, which is not a harsh condition to satisfy in general cases. It is admitted that here a quite simple problem is treated to analyse the sensitivity. It might be different in other cases with different materials, boundary conditions and source excitation, in which a little mismatch may lead to an unacceptable error.

To summarize this part, being conscious firstly that errors exit in our interface model produced from the angle mismatch 𝜃 in TMM. On the other hand, thin porous layers considered in this paper are less sensitive to this mismatch. In most situation of angle mismatch, less error is reproduced compared to the air layer. Even in the cases where this angle mismatch cannot be avoided, the proposed approach is still able to provide reliable results for modelling porous films. We will demonstrate this aspect in the last numerical example of section 6.

SPATIAL DISCRETIZATION

In this work, the eXtended Finite Element Method (X-FEM) [START_REF] Moës | High Order X-FEM and Levelsets for Complex Microstructures: Uncoupling Geometry and Approximation[END_REF] is adopted to discretize the computational domain Ω for two reasons: (a) Prescribing the proposed interface conditions which lead to strong discontinuities; (b) Meshing independently of the geometry of the interface, meaning that the elements do not need to coincide with interface Γ.

Interface definition and enrichment strategy

Within X-FEM discretization, a level-set function interpolated on the mesh is used to classify elements, for example in our two-media domain:

𝜙(𝐱) = ∑ 𝑖 𝜙 𝑖 𝑁 𝑖 (𝐱) ⎧ ⎪ ⎨ ⎪ ⎩ < 0 on Ω 1 = 0 at Γ > 0 on Ω 2 , ( 44
)
where 𝑁 𝑖 (𝐱) are the standard finite-element shape functions, 𝜙 𝑖 are the signed distances to the interface at the nodes of the mesh. The geometry of the interface is defined implicitly as the iso-zero of the level-set function. The subdomain of different media is classified through the sign of the level-set function. With the interface represented by the level-set, the whole mesh is divided into three types of elements: standard elements as in conventional FEM, enriched elements which contain the interface boundary and are fully enriched, and blending elements which are partially enriched as illustrated in fig. 6. Geometrical errors exist when curved boundaries are defined by a piece-wise linear level-set function. These errors will prevent solution to achieve the optimal rate of convergence in high order approximation. To tackle such numerical issue, the space-tree quadrature approach presented in Legrain et al. 46 is considered, allowing the use of coarse approximation meshes but still retaining optimal convergence. FIGURE 6 Typical X-FEM mesh with interface in triangular elements (red nodes: enriched nodes; green: enriched elements; yellow: blending elements; white: standard elements) With regard to the enrichment strategy for the approximation, a strong discontinuity of both primal (pressure) and associated dual variables occurs across the interface in all coupling configurations. Therefore, the pressure field needs to be enriched so that this behaviour is correctly modelled. The pressure is approximated in an appropriate finite dimensional space 𝑈 ℎ ⊂ 𝑈 (Ω 1 ∪ Ω 2 ) as:

𝑝 𝑓 ℎ (𝐱) | | |Ω = 𝑛 ∑ 𝑖 𝑁 𝑖 (𝐱)𝑝 𝑓 𝑖 + 𝑛 𝑒𝑛𝑟 ∑ 𝑗 𝑁 𝑗 (𝐱)𝐻(𝐱)𝑎 𝑗 , ( 45
)
where 𝐻(𝐱) represents the generalized Heaviside function which introduces the strong discontinuity in the approximation:

𝐻(𝐱) = { -1 𝜙(𝐱) < 0 1 𝜙(𝐱) ≥ 0 , ( 46
)
and 𝑁 𝑗 (𝐱) are the shape functions satisfying the partition of unity. 𝑎 𝑗 are the additional degrees of freedoms (dofs for the following) to represent the discontinuity in the elements. For the fluid-PEM coupling, we also have a displacement field defined only in the porous subdomain (for instance Ω 2 ) which is expressed as in standard finite element approximation:

𝐮 𝐬 ℎ (𝐱) | | |Ω 2 = 𝑛 ∑ 𝑖 𝑁 𝑖 (𝐱)𝐮 𝑠 𝑖 , ( 47
)
where 𝐮 𝐬 (𝐱) is a vectorial field having two or three components in 2∕3D problems. In this paper, all conventional shape functions 𝑁 𝑖 and 𝑁 𝑗 are based on high-order Bernstein polynomial basis in order to recover an optimal convergence for Biot's equations at high frequencies. Interested reader can refer to our previous work [START_REF] Wu | High-Order X-FEM for the Simulation of Sound Absorbing Poro-Elastic Materials with Coupling Interfaces[END_REF] for additional details on discretization strategies for domains involving PEMs within the X-FEM.

Discretization for different couplings

As three typologies of interface operator were introduced for coupling the variational formulations, three discrete systems will arise. Substituting the enriched fields eq. ( 45) into the previous variational formulations eq. ( 6) together with interface operators (21) results in the following discretized algebraic system for fluid-fluid coupling:

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐾 𝑏 𝑝1 + 𝐾 𝐼 𝑝1 𝐾 𝐼 𝑝12 𝐾 𝑏 𝑝1𝑎 + 𝐾 𝐼 𝑝1𝑎 𝐾 𝐼 𝑝21 𝐾 𝑏 𝑝2 + 𝐾 𝐼 𝑝2 𝐾 𝑏 𝑝2𝑎 + 𝐾 𝐼 𝑝2𝑎 𝐾 𝑏 𝑎𝑝1 + 𝐾 𝐼 𝑎𝑝1 𝐾 𝑏 𝑎𝑝2 + 𝐾 𝐼 𝑎𝑝2 𝐾 𝑏 𝑎𝑎 + 𝐾 𝐼 𝑎𝑎 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑝 𝑓 1 𝑝 𝑓 2 𝑎 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑓 1 𝑓 2 𝑓 𝑎 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ , ( 48
)
where 𝑝 𝑓 𝑖 are the pressure dofs for each subdomain and 𝑎 are the enriched dofs at the interface. 𝑓 𝑖 and 𝑓 𝑎 represent external forces on the bulk and enriched elements. Terms 𝐾 𝑏 and 𝐾 𝐼 are the generalized stiffness matrices from bulk and interface contribution, respectively.

Configuration "f-f-P" has the same interface operator as fluid-fluid coupling, and additional coupling contributions in Biot's bulk part are taken into account in "f-P-P" case. Therefore, only the case with embedded PEM film is presented here for fluid-PEM coupling. When eq. ( 45) and eq. ( 47) are inserted in variational formulations eq. ( 6) and eq. ( 7) with interface terms (27), the resulting discretized linear system is written as:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐾 𝑏 𝑝1 +
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝑝 𝑓 1 𝑝 𝑓 2 𝑢 𝑠 𝑥 𝑢 𝑠 𝑦 𝑎 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝑓 𝑝1 𝑓 𝑝2 𝑓 𝑢𝑥 𝑓 𝑢𝑦 𝑓 𝑎 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ , ( 49
)
where all the terms including 𝑢 𝑥 and 𝑢 𝑦 are originated from the displacement field in PEMs and 𝑓 𝑢𝑥 , 𝑓 𝑢𝑦 are the external forces applied on the elastic phase. Additional bulk contribution 𝐾 𝑏 𝑝 𝑖 𝑢 represent the coupling terms between 𝑝 𝑓 and 𝑢 𝑠 in Biot's mixed formulations. the detailed expressions of the generalized stiffness matrices 𝐊 presented in these discrete systems are provided in appendix C for completeness.

Reminded that every term defined in enrichment blocks in the system above is integrated using a modified Gauss quadrature, as proposed in [START_REF] Moës | A Computational Approach to Handle Complex Microstructure Geometries[END_REF] using partitioned sub-elements at the interface. In addition, the leading discrete systems being potentially non-symmetric ("f-P-P" with the proposed Biot's interface operator as discussed in section 4), a direct solver (MUMPS [START_REF] Amestoy | Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures[END_REF] ) is adopted in our work using a well-suited assembly graph and algorithm for each system to improve the solving procedure.

NUMERICAL EXAMPLES

In this section, the proposed approach is tested by means of several numerical examples to study its accuracy and advantages in different simulation scenarii. Our approach is compared against analytical or reference FEM solutions if exact one is unavailable.

To assess the accuracy of the proposed formulations with interface operators ( 24), ( 26) and ( 28), the error between numerical and reference solutions is measured in 𝐿 2 norm. Be interested in acoustic problems, the pressure field is of our prime interest. We define a mean global relative error for the pressure of fluid 𝑝 𝑓 in 𝐿 2 norm as:

𝜀 𝐿 2 = ( ∫ Ω | | | 𝑝 𝑓 XFE -𝑝 𝑓 ref | | | 2 dΩ ) 1∕2 ( ∫ Ω | | | 𝑝 𝑓 ref | | | 2 dΩ ) 1∕2 × 100%, ( 50
)
where Ω is the entire computation domain Ω 1 ∪ Ω 2 . 𝑝 𝑓 XFE denotes the results obtained with our approach and 𝑝 𝑓 ref represents a reference solution (analytical one and FEM). As a pressure drop will appear at the interface due to the generalized interface conditions, the local accuracy also needs to be quantified with help of the following error:

𝜀 𝑑 𝐿 2 = ( ∫ Γ | | | 𝑝 𝑓 XFE -𝑝 𝑓 ref | | | 2 dS ) 1∕2 ( ∫ Γ | | | 𝑝 𝑓 ref | | | 2 dS ) 1∕2 × 100%, (51)
where 𝑝 𝑓 represents the pressure gap between the two sides of the interface Γ. Note that the associated dual variable (total displacement 𝑢 𝑡 𝑛 ) is also discontinuous over the interface, but as mentioned before, we are more interested in the pressure field, only accuracy of pressure is evaluated

In addition, when analysing convergence, we require an indicator illustrating the resolution of the mesh. In vibro-acoustic problems, the number of dofs per wavelength 𝐷 𝜆 is commonly used:

𝐷 𝜆 = 2𝜋 𝑘 * √ 𝑁 dof area(Ω) (52)
where 𝑘 * is the largest absolute value of wave numbers in the media, since there exists three complex-valued wave numbers when PEMs are involved in the problem. 𝑁 dof denotes the total number of dofs and area(Ω) is the measure of the computational domain.

Plane wave propagation in infinite plane

The objective of this first test case is to verify the implementation of the variational formulations for the three interface operators and coupling types within X-FEM. Plane waves propagating in infinite planes are considered, where both planar and curved interfaces are accounted for. Therefore, this numerical example allows for testing different acoustic fields such as reflection, transmission and scattering. The analytical solutions associated with each test case are available. They can be obtained by solving the reflection 𝑅, transmission 𝑇 𝑖 (three transmission coefficients for PEMs subdomains) and scattering coefficients at the interfaces with generalized interface conditions derived before. The reader can refer to [START_REF] Allard | Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials[END_REF] for the detailed procedure on the obtention of the analytical solutions. In this experiment, a plastic foam is chosen as porous bulk media whose properties are given in appendix C: table A.1. Fig. 7 depicts the three proposed plane wave problems and boundary conditions settings accompanied by representative X-FEM meshes. Operator (24) with equivalent fluid matrix (33) and operator (28) with Biot's matrix given in section 4 are tested in fig. 7(a) and fig. 7(b), respectively, while the scattering problem fig. 7(c) takes operator (26) with Limp model presented in section 4.1.2 into account. The computational domain is a square of length 0.2 × 0.2 m 2 for the two first cases and 0.5 × 0.5 m 2 for the scattering problem. The interfaces for the thin film layers are defined by level-set functions (𝑥 = 0 for planar interface problems and radius 𝑅 = 0.1 m for the scattering problem). The thickness of the thin layer is fixed to 𝑑 = 10 -3 m in the transfer matrices in all cases. It is observed that the arbitrary interfaces are all embedded inside the elements instead of lying along the boundaries as it would be requested for classical FEM.

For the planar interface problem, the calculation is performed using unstructured meshes to illustrate the robustness of the proposed X-FEM formulations no matter how the elements are cut by the interfaces. A simple structured mesh is used for the curved interface. Exact Neumann boundary conditions are prescribed on the four outer boundaries: acoustic velocity (displacement) 𝑣 𝑒𝑥 ∕𝑢 𝑡 𝑒𝑥 for fluid media (green line) and total displacement and in-vacuo stress for poro-elastic media (yellow line). Fig. 8illustrates numerical solutions of pressure field 𝑝 𝑓 for plane waves with oblique incidence (45 °, 30 °) and a normal incidence propagating under an excitation of 3000 Hz and 2000 Hz. A gradual decrease in pressure amplitude and a change of incidence angle are first observed in the porous bulk domains. Moreover, apparent strong discontinuities are observed across the interfaces: these pressure drops are the effect of accounting for the thin film layer through the interface operators. As exhibited in the results, these discontinuities inside the elements are all well captured at the location of the implicit interfaces, indicating that the X-FEM implementation works properly. Thanks to high order polynomial approximation (𝑝 = 4 here), these examples give a quite smooth solution even with such a relatively coarse mesh.

Convergence

To evaluate the accuracy of the implemented formulation, the convergence of the solution needs to be analysed 50) and eq. ( 51)) as a function of the number of dofs per wavelength 𝐷 𝜆 . The ℎ-convergence depicted here are performed with linear, quadratic and cubic approximations (𝑝 = 1, 2, 3). In general, all the curves are close to the so-called "optimal" rate of convergence (ℎ -(𝑝+1) ) in each interpolation degree. The 𝑝-refinement (red curves) provides a more efficient convergence rate, in the sense that smaller 𝐷 𝜆 are needed to achieve the same level of accuracy. Even so, this efficiency is not noticeable in certain configurations, which will be explained below. The formulation implemented within X-FEM provides a similar solution accuracy (magnitude of error are in the same order) in terms of the global solution of pressure and pressure jump at interfaces.

The comparison between fluid-fluid coupling and fluid-PEM coupling can be illustrated from the four first figures (fig. 9(a,b,c,d)). As seen that larger 𝐷 𝜆 are needed to arrive to the "optimal" convergence stage using the same meshes (preasymptotic range of convergence are more obvious) in curves in figs. 9(c) and 9(d) than the ones in figs. 9(a) and 9(b) . This is due to the fact that the displacement field is involved in fluid-PEMs coupling, leading to a large disparity of scale between pressure and displacement field in 𝑢 𝑠 -𝑝 𝑓 mixed Biot's formulation, which prevents a global convergence as reported in [START_REF] Rigobert | Investigation of the Convergence of the Mixed Displacement-Pressure Formulation for Three-Dimensional Poroelastic Materials Using Hierarchical Elements[END_REF] . This is also why the efficiency of 𝑝-refinement with fluid-PEMs coupling is not obvious as red curves in figs. 9(c)-9(f), which was already observed in our previous work [START_REF] Wu | High-Order X-FEM for the Simulation of Sound Absorbing Poro-Elastic Materials with Coupling Interfaces[END_REF] .

The above convergence curves verify the implementation of the proposed formulation within the X-FEM. No matter what types of coupling and how the elements are cut by the interface, a desirable accuracy with a close to optimal convergence is obtained. The convergence has been assessed for other porous systems at different frequencies during this work. For the sake of brevity, only these representative results are provided here, but the same conclusions apply to the other cases.

Conditioning

The assessment of the conditioning of the system stiffness matrices is also important, as it is directly related to the stability and efficiency of the solution, especially when the system is solved through iterative solvers. Configurations "f-f-f" and "f-P-P"

FIGURE 10

Evolution of the condition number with respect to 𝐷 𝜆 for configuration "f-f-f" (top) and "f-P-P"(bottom): (a,c) comparison between perfect interface and generalized interface conditions with 𝑑 = 10 -3 m using interface conforming mesh; (b,d) Influence of values in interface operators by varying thickness of film 𝑑 using unstructured interface non-conforming mesh including symmetrical interface operator (24) and non-symmetrical one (28) are studied. Interpolation degree 𝑝 = 3 is used for all experiments. Fig. 10shows the variations of condition number with respect to 𝐷 𝜆 . We first examine how the interface operator affects the conditioning of the linear system. Therefore, a perfect interface condition without any discontinuity is first solved within standard FEM as a reference. In order to eliminate the influence of X-FEM on cut elements (as only the influence of the interface operator is characterized here), the implicit interface is positioned along the boundaries of elements in a structured mesh. With this mesh, the X-FEM performs like the standard FEM with double nodes along the interface. The parameters for the film used in the interface operators are the same as in the solutions figs. 8(a) and 8(b). Fig. 10(a) illustrates that the conditioning for a fluid-fluid coupling involving the proposed interface operator is close to the one with the perfect interface, yielding to a rate close to (ℎ 2 ). The similarity between these curves indicates that the proposed approach remains relevant, even if implemented within existing commercial FEM codes (the thin layers model does not degrade the conditioning of the original system for typical porous films). For "f-P-P" as shown in fig. 10(c), the interface operator does not either deteriorate the conditioning compared to a perfect interface, even though we have a non-symmetric linear system in this configuration. It is remarked that the conditioning for this configuration is far larger than the one for "f-f-f", which stems from the characteristics of the mixed Biot's formulation (large scale disparity between pressure and displacement). This fact also slows down the overall convergence speed, as shown previously in the convergence study. Even with a so large conditioning number, a converged solution could be still obtained with the used direct solver.

Then, the effect of [𝐀] on the conditioning within X-FEM is examined. The same setup and meshes are used as in figs. 8(a) and 8(b), but the components of the matrix [𝐀] are varied by changing the thickness 𝑑 of the layer.

It is first observed in figs. 10(b) and 10(d) that the magnitude of the condition number is at least five orders larger and increases more quickly than the one in figs. 10(a) and 10(c) under the same 𝐷 𝜆 . This worse conditioning stems from the elements cut by the interface, as an unstructured mesh is used. This behaviour can be improved by considering relevant pre-conditioners such as [START_REF] Axelsson | Superlinear Convergence Using Block Preconditioners for the Real System Formulation of Complex Helmholtz Equations[END_REF] . In practice, we did not encounter any instability issue when enforcing interface conditions under such "irregular" cut elements. Secondly, a noticeable difference is observed for the first points of configuration "f-f-f": the smaller the value of the coefficients in the matrix, the larger the condition number of the system. That is because, under coarser meshes, the condition number of the linear system is dominated by the interface terms associated with the value of coefficients 𝐴 𝑖𝑗 . As transformed from matrix [𝐓] to [𝐀], smaller thickness in [𝐓] leads to larger value in [𝐀]. When increasing the mesh resolution, contribution to the global conditioning from the bulk parts' integrals becomes more important than the interface terms. As a consequence, for the last considered mesh points, the condition numbers of different operators approach to each other. In regard to "f-P-P" in fig. 10(d), no apparent influence is observed with respect to thickness: all the curves are intertwined together. As in our formulations, the final conditioning of the FE system depends on the magnitude of the bulk and interface contributions. The similarity between these conditioning curves is because with air-PEMs coupling, the contribution to conditioning from generalized interface terms is much smaller than the one from bulk parts in Biot's equations (it already has a worse conditioning), even with a very small thickness. Conditionings shown in figs. 10(b) and 10(d) imply that for a converged solution (with sufficiently large 𝐷 𝜆 ), no evident effect to the global conditioning is produced from the different proposed interface operators within X-FEM.

From the studies in section 6.1.1 and 6.1.2, we summarize that for the problems of interest, the variational formulations with the proposed interface operators implemented within the X-FEM are able to give accurate solutions with a limited impact on the conditioning of discrete systems.

Validation for the impedance tube

The second example is the simulation of an impedance tube, which is frequently used to characterize acoustic systems. Here, a sound absorption system composed of a 2 mm sandwich film layer coated on a 20 cm thick polyurethane (properties available in table A.1) is introduced in a rigid, rectangular air cavity tube. The tube has a dimension of 1 × 0.2 m subjected to a unit velocity source at the left end, as shown in fig. 11. This example aims at validating the proposed model for thin layers systems in a typical application situation through comparison with conventional finite elements where all the material layers are discretized.

Without loss of generality, two geometries of sound absorption systems are considered: an oblique surface at 45 °as well as a curved one with a radius of 0.27 m. For the standard FE models, the meshes for these two geometries have to be refined near the thin layer, as shown in fig. 12(a) and fig. 12(b), in order to avoid distorted elements with large aspect ratios and to capture wave propagation in all directions. Additionally, for curved boundaries, such fine meshes are required to describe the curvature of the contour properly, weakening the geometrical error. Note that for certain implementations of 𝑝-FEM as in [START_REF] Bériot | Efficient Implementation of High-Order Finite Elements for Helmholtz Problems[END_REF] , coarser and anisotropic meshes could be used to discretize the planar thin layers while maintaining a high accuracy. Here a more conventional mesh is considered for the standard FEM model in order to remain suitable even for low interpolation orders.

In contrast, when the problem is solved within the proposed approach, a unique simple background mesh is sufficient for the different geometries, see fig. 12(c), as the thin layers are condensed into interfaces that are implicitly defined by a level-set function. Besides, the number of dofs is reduced when using such simple non-conforming mesh compared to the ones used in standard FEM.

Furthermore, the sandwich of acoustic films used here is composed of two woven films on the two sides and one non-woven film in the centre with thickness of 𝑑 1 = 𝑑 3 = 0.6 mm and 𝑑 2 = 0.8 mm respectively. This kind of multi-layer system is straightforward to be modelled under our approach thanks to the TMM technique through changing the global transfer matrix. The global transfer matrix for three layers that have the same physical nature is directly written as:

[𝐓] sand = [𝐓] 1 × [𝐓] 2 × [𝐓] 3 (53)
where [𝐓] 𝑖 are the transfer matrices for each layer. If stratified layers are not of the same nature, the global matrix is obtained by multiplying interface matrices containing the continuity equations at the interfaces (more detail can be found in 32, section 11 ). It is noticed that in the case of an oblique 45 °surface, the thickness 𝑑 𝑖 and 𝑘 eq𝑥 in each matrix are readily corrected to maintain the angle match between source excitation and orientation of film surface. When standard FEM solves such multiple layers problems, at least three dofs per node for one element in two-dimensional problems will be added in the global system to represent PEMs layer. It introduces far more unknowns to be solved compared to the proposed method. Therefore, treating multiple layers is another strength of our approach compared to standard FE models. Fig. 13 presents the pressure field for the two sound absorbers depicted in figs. 12(a) and 12(b). We observe that the pressure fields computed by the (XFEM+TMM) approach are in excellent agreement with the ones computed by full FEM discretization. A significant drop and an obvious pressure "discontinuity" is captured across the thin films in both numerical models. Moreover, it is worth pointing out that 𝑝-refinement is used to ensure proper convergence of the solution for both models. It is an efficient approach in acoustics due to its favourable behaviour against pollution, and it also allows for keeping a fixed mesh during the convergence procedure. However, it leads to a highly large linear system in the case of the full discretized FEM model involving thin layers. Note also that 𝑝-adaptative methods such as the work presented in [START_REF] Bériot | Anisotropic Adaptivity of the P-FEM for Time-Harmonic Acoustic Wave Propagation[END_REF] were developed for acoustics and could be used to improve the efficiency when using conventional meshes. Nevertheless, our approach is quite direct: one fixed background mesh can be used for any geometry and interpolation degree without any adjustment. To further quantify the accuracy of the proposed approach against standard FEM, the absorption coefficient 𝛼 of the porous system is estimated. We evaluate this coefficient numerically in the models by calculating the local average pressure and its gradient at the interface for each frequency. Here, an ideal normal incidence (planar interface perpendicular to the incident direction) is considered.

The absorption coefficient with respect to the frequency is plotted in fig. 14 for the considered sound absorption system. The solution calculated by a semi-analytical package "pymls" (based on work of [START_REF] Dazel | A Stable Method to Model the Acoustic Response of Multilayered Structures[END_REF] ) serves as a reference. This absorption exhibits multiple resonances from 1 Hz to 5000 Hz. Seventeen points are mostly chosen at the resonance frequencies, to estimate coefficient 𝛼 for the full FE model and the proposed approach. The converged solution for each frequency is obtained by increasing the interpolation order while keeping the same mesh. As illustrated, the two numerical models are in good agreement and stay close to the reference in such a wide frequency range. This numerical example reveals twofold benefits of the proposed model: (a) Flexibility with respect to the geometry, as the mesh is independent of the surface location. This benefit provides a convenience to test multiple geometries (configurations) of the sound absorption system without re-meshing; (b)The proposed method ensures high accuracy in a wide range of frequencies.

Car cavity with multiple porous seats

As a last example, let us consider a more realistic case in which the interior sound pressure of a full-size car is predicted by the proposed method. The dimensions of the vehicle are approximately 2.67 m × 1.1 m. Two driver seats made of XFM foam (properties available in table A.1) with a 1 mm thin layer of non-woven porous film are presented in the cavity. XFM foam is governed by Biot's mixed equations, and the film is controlled by the Limp model. An arbitrary normal velocity (𝜕𝑝∕𝜕𝑛 = 1) is prescribed on the front wind-shield with a 1000 Hz excitation frequency. As the acoustic field is too complex to correct the incidence angle 𝜃 in the TMM (keeping 𝜃 = 0), this example allows us to account for the error due to the angle mismatch as well.

In the standard FEM model, the geometry of the two seats and films are represented explicitly: elements need to conform to the interfaces between the different media and the film is required to be discretized. To this aim, neighbour elements need to be refined in order to reduce the geometrical error and capture its acoustic behaviour as shown in fig. 15(a). Within our approach, the geometry and mesh are shown in fig. 15(b) where the contour of the seats are defined implicitly by level-set functions and the films are reduced as interfaces with Limp transfer matrices (34). It is reminded that the mesh in the proposed approach is independent of the seats and attached films. Such coarse mesh is in stark contrast with the standard FE model, and results in fewer dofs in the solution.

The problem solved using standard FEM and the proposed approach are compared. Again, a 𝑝-refinement is performed for both methods to obtain converged sound pressure fields, as fig. 16(a) and fig. 16(b). These two solutions exhibit a perfect visual agreement, where the pressure damping in the seats and a noticeable pressure drop close to the thin films are observed. The difference between these two solutions is evaluated in fig. 16(c). Most of the distinction is below 1dB(≤ 1%), which stems partly from the angle mismatch exits everywhere in the film model and a re-interpolation of the pressure fields on two non-identical meshes. For engineering problems, this difference is completely acceptable compared with the full FE model, especially when the computational cost is considered.

For further analysis, a clipping is done along the line depicted in figs. 16(a) and 16(b). The corresponding results are illustrated in fig. 16(d) in solid blue line for the FEM model and dashed red line for the proposed approach. Herein, the pressure drop induced by the thin film is revealed more clearly at the boundaries of the seats, and an obvious sound dissipation is observed in the first seat. The two solution curves are nearly overlapped along the clipping line.

To illustrate the influence of the film and porous seats on the sound field, solution without film and solution without seat/film (fig. 16(d) green and orange solid line) are provided for comparison with the previous configurations. As expected, the solution without film is continuous across the boundaries of the seats. Besides, the pressure in the air cavity is slightly lower than the one with films. This is because the film was selected with a large flow resistivity and density, it can be considered as a rigid material so that more energy is reflected in the air. Because of the lack of dissipation from PEMs made seats, the empty car cavity (orange solid line) exhibits a larger pressure amplitude at the resonance zones. This comparison thus indicates as well that porous thick materials are able to stabilize and decrease the pressure level in the cavity. In order to have geometrically adapted meshes, two more re-meshing are needed when using the standard FEM model for the last two configurations. In contrast, the mesh depicted in fig. 15(b) can be used directly within the proposed approach for any change of interior geometry. This is quite beneficial for simulation engineers to carry out geometric or parametric investigations on complex geometries.

CONCLUSION

In this work, we established a computational framework to solve time harmonic problems involving thin acoustic layers. Largescale problems including such layers are highly time-consuming in terms of pre-processing and solving when treated by standard finite element methods (especially when the boundaries of the layers are complex). It is even more complicated when coupling problems involving various (fluid-solid interaction) fields.

To avoid the discretization of the thin layers, we replace its volume by a surface along with appropriate interface conditions based on the transfer matrix method (TMM). We proposed a general variationally consistent formulation enforcing such conditions weakly, which is able to comprise different physics of media such as fluid, elastic solid and poro-elastic materials. We applied the proposed formulations to sound absorption systems composed of acoustic fluid and poro-elastic materials, leading to four coupling configurations and three interface operators. Then, we demonstrated the relevance of our thin layer model through a comparison with the simplified pressure drop model. We discussed as well the sensitivity of the TMM due to the incidence angle mismatch, exhibiting that the proposed interface model is able to maintain an acceptable accuracy for thin resistive film in most cases. As the interface conditions lead to strong discontinuities, and to tackle the geometry restriction, the eXtended Finite Element Method (X-FEM) is employed to discretize the computational domain.

The proposed generalized variational formulations are assessed by different plane wave benchmarks. The solution and convergence verify the implementation and illustrate the accuracy of the proposed formulations within the X-FEM. Optimal rates of convergence in 𝐿 2 norm ((ℎ -(𝑝+1) )) are obtained. The conditioning of the resulting algebraic system remains comparable to standard FE models. The impedance tube example demonstrates the convenience of the proposed model for problems involving multiple layers, and its flexibility in the case of variable geometries. This numerical example shows as well that the proposed model is in good agreement with the full standard FE model in an extensive range of frequency with much fewer degrees of freedoms. The last example provides an engineering scenario where multiple inclusions are coated by thin layers in a real size car. Again, significant benefits on pre-processing and size of the discrete system are shown compared to standard FE models. The reliability of the proposed method for predicting such complex acoustic fields is also demonstrated in the case where angle mismatch exists anywhere in the TMM interface model.

The proposed model is proven competitive and advantageous in comparison with the conventional FE model for problems containing thin layers. It enables to capture the solution accurately in such problems using simple, coarse and even uniform background meshes which are independent of the location and geometry of thin layers. To the authors' knowledge, shape or topological optimization for such sound absorption systems would be a relevant application of the proposed approach. It would also be interesting to apply the proposed general framework for other media, such as elastic solids or thin elastic plates, which have an important vibro-acoustic impact as well. Finally, remind that although it is proven efficient and accurate for modelling resistive thin porous layers, the current approach may lead to a non-negligible error for other materials in some complex wave problems due to the plane wave assumption in TMM, as discussed in section 4. In addition, the proposed variational formulation involving the Biot's matrix eq. ( 35) is non-symmetrical. The symmetrization of this formulation would be an interesting and open subject. More thoughtful effort could be given to address this aspect in further study. The effective density of an equivalent fluid ρeq is written as:

ACKNOWLEDGEMENTS

APPENDIX A USED POROUS MATERIALS WITH ASSOCIATED JCA AND BIOT'S COEFFICIENTS

ρeq = 𝜌 a 𝛼 𝜙 (1 + 𝜔 0 j𝜔 √ 1 + j𝜔 𝜔 ∞ ) (A.1)
Under the continuity condition at the interfaces Γ - 𝑖 and Γ + 𝑖 , a transfer matrix relationship based on Biot film layer can be obtained with arbitrary coefficients as following:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ σ- 𝑦𝑥1 = 0 𝑢 𝑠 𝑥 - 1 𝑢 𝑡 𝑛 + 1 σ- 𝑥𝑥1 = 0 𝑝 𝑓 + 1 𝑢 𝑠 𝑦 - 1 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑇 𝑃
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ σ- 𝑦𝑥2 = 0 𝑢 𝑠 𝑥 - 2 𝑢 𝑡 𝑛 + 2 σ- 𝑥𝑥2 = 0 𝑝 𝑓 + 2 𝑢 𝑠 𝑦 - 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ (B.1)
where 𝑇 𝑃 𝑖𝑗 are the coefficients for the PEM layer that could be the matrix

B.2 From transfer matrix to Admittance matrix for Biot's model

Let us consider a generalized interface condition for the case f-P-P as presented eq. ( 20), which can be re-expressed by a transfer matrix with arbitrary coefficients: To obtain the admittance matrix which expresses the relationship between dual variables and primal variables, these linear equations need to be rearranged. We put dual variables on the left-hand side and the primal variable in the right-hand sides, leading to the following equations:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ σ- 𝑦𝑥1 = 0 σ- 𝑥𝑥1 = 0 𝑝 𝑓 + 1 𝑢 𝑡 𝑛 + 1 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
[𝐁] 4×4 Thus, the final equation for the dual variables expressed by primal variables on the positive side of Γ 2 is written as following:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑢 𝑡 𝑛1 𝑢 𝑡 𝑛2 σ𝑥𝑥2 σ𝑦𝑥2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = [𝐁] -1 4×4 [𝐂] 4×4 ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑝 𝑓 1 𝑝 𝑓 2 𝑢 𝑠 𝑥2 𝑢 𝑠 𝑦 2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ (B.11)
where all the superscript (•) + are removed to be consistent to eq. ( 27). The generalized admittance matrix [𝐀 𝐏 ] in eq. ( 27) is obtained as:

[𝐀 𝐏 ] = [𝐁] -1 [𝐂] (B.12)
The methodology presented here to obtain the generalized admittance matrix is also suitable for other case, for instance the fluid-fluid coupling. Even though, derivation of admittance matrix for fluid-fluid coupling is quite direct.

C DISCRETIZED TERMS IN LINEAR SYSTEMS

We offer all related discretized terms in the linear system (48) and (49) where 𝐍 represents the vector of the shape function as defined in eq. ( 45) and eq. ( 47), 𝛁𝐍 is its gradient. Compared to the previous system, all the terms in (49) involving only pressure 𝑝 𝑖 remain the same, new terms including displacement of solid phase 𝑢 𝑖 are given here. The bulk contribution including displacement are written as: The symmetry of discrete linear system can be recognized with this discretized expressions. From the interface terms in eq. (C.2), the conclusion in eq. ( 29) is obtained. The same conditions for the system of "f-P-P" are obtained under the expressions provided above eq. (C.4).

𝐾 𝑏 𝑝𝑖𝑢 = 𝐾 𝑏 𝑢𝑖𝑝 = -γ ∫ Ω 𝑒

FIGURE 2

 2 FIGURE 2Four interface coupling configurations.

FIGURE 3

 3 FIGURE 3 Variation of coefficient related to pressure drop (𝑇 drop with normal incidence 𝜃 = 0 °and thickness of 𝑑 = 1 mm) for different film models as function of frequency: (a) Non-woven film, (b) woven film

FIGURE 4 FIGURE 5

 45 FIGURE 4 Configuration for testing the sensitivity of the transfer matrix

FIGURE 7

 7 FIGURE 7 Three test configurations with X-FEM meshes and boundary conditions: (a) coupling fluid-fluid; (b) coupling fluid-PEMs with embedded PEM film; (c) cylinder scattering problem with embedded Limp film

FIGURE 8

 8 FIGURE 8Pressure solutions (Pa) for the proposed approach: (a,d) configuration "f-f-f" (45 °, 3000 Hz), (b,e) configuration "f-P-P" (30 °, 3000 Hz), (c,f) configuration "f-f-P" in cylinder scattering problem (normal incident 2000 Hz)

FIGURE 9

 9 FIGURE 9 Global convergence (left) and local convergence (right) for:(a-b) fluid-fluid coupling; (b-d) fluid-PEMs coupling with an embedded PEM film; (e-f) fluid-PEMs coupling in scattering problem with an embedded Limp film.

FIGURE 11

 11 FIGURE 11Kundt's tube problem in two-dimensions.

FIGURE 12

 12 FIGURE 12 Meshes within standard FE model (a) for oblique planar films; (b) for curved boundary films; (c) Meshes within the proposed approach for the above two geometries.

FIGURE 13 FIGURE 14

 1314 FIGURE 13 Pressure solution (Pa) for standard FEM and the proposed approach (XFEM+TMM) for sandwich absorption systems under different geometries and an excitation of 2000 Hz

FIGURE 15

 15 FIGURE 15 Boundary condition and meshes for the car cavity problem: (a) standard FEM model; (b) proposed approach.

FIGURE 16

 16 FIGURE 16 Example of solution (pressure in dB) at 1000 Hz (a) within standard FEM; (b) with the proposed approach; (c) Difference between (a) and (b); (d) Clipping of the solutions along the dashed lines for different configurations.

TABLE 1

 1 Film parameters

	Parameters	Woven	Non-woven
	Porosity 𝜙 (-)	0.04	0.72
	static flow resistivity 𝜎 (N.m -4 )	775 × 10 3	87 × 10 3
	Tortuosity 𝛼 (-)	1.15	1.02
	thermal characteristic length Λ ′ (m) 230 × 10 -6 480 × 10 -6
	Viscous characteristic length Λ (m) 230 × 10 -6 480 × 10 -6
	frame density 𝜌 1 (kg.m -3 )	809	171
	Young's modulus 𝐸 (Pa)	260 × 10 6	50 × 10 3
	Poisson ratio 𝜈 (-)	0.3	0.3
	Loss factor 𝜂 𝑠 (-)	0.5	0.5

TABLE A . 1

 A1 Porous materials parameters for bulk partThis appendix gives a brief explanation of the coefficients appearing in the JCA model and Biot's equations. One can obtain the coefficients with the provided material's parameters in table 1 and table A.1. All the definition presented here can be found in other literature as well.

	Parameters	Plastic foam Polyurethane	XFM
	Porosity 𝜙 (-)	0.97	0.98	0.98
	static flow resistivity 𝜎 (N.m -4 )	57 × 10 3	3.75 × 10 3	13.5 × 10 3
	Tortuosity 𝛼 (-)	1.54	1.17	1.7
	thermal characteristic length Λ ′ (m) 73.8 × 10 -6	742 × 10 -6	160 × 10 -6
	Viscous characteristic length Λ (m) 24.6 × 10 -6	110 × 10 -6	80 × 10 -6
	frame density 𝜌 1 (kg.m -3 )	46	22.1	30
	Young's modulus 𝐸 (Pa)	214 × 10 3	70 × 10 3	200 × 10 3
	Poisson ratio 𝜈 (-)	0.3	0.39	0.35
	Loss factor 𝜂 𝑠 (-)	0.4	0.265	0.05

  35. Since displacements do not exist in the fluid bulk parts, two lines associated with the 𝑢 𝑠 𝑥 and 𝑢 𝑠 𝑦 at the interface Γ - 1 of PEM layer need first to be eliminated, then 𝑢 𝑠 𝑥 and 𝑢 𝑠 𝑦 at the interface Γ - 2 should be re-expressed by the 𝑝 𝑓 and 𝑢 𝑡 𝑛 . And, based on the relationship in eq. (B.1), we obtain four linear equations, respectively as: the two equations associated with 𝑢 𝑠 𝑥 and 𝑢 𝑠 𝑦 at the interface Γ - 2 eqs. (B.2a) and (B.2c) are re-written in matrix form as: Finally, only the variables 𝑝 𝑓 and 𝑢 𝑡 𝑛 at the faces Γ + 𝑖 are remained in above equation, we write it in a matrix form as:

						0 = 𝑇 𝑃 12 𝑢 𝑠 𝑥	-2 + 𝑇 𝑃 13 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 15 𝑝 𝑓 + 2 + 𝑇 𝑃 16 𝑢 𝑠 𝑦	-2	(B.2a)
					𝑢 𝑡 𝑛	+ 1 = 𝑇 𝑃 32 𝑢 𝑠 𝑥	-2 + 𝑇 𝑃 33 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 35 𝑝 𝑓 + 2 + 𝑇 𝑃 36 𝑢 𝑠 𝑦	-2	(B.2b)
						0 = 𝑇 𝑃 42 𝑢 𝑠 𝑥	-2 + 𝑇 𝑃 43 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 45 𝑝 𝑓 + 2 + 𝑇 𝑃 46 𝑢 𝑠 𝑦	-2	(B.2c)
					𝑝 𝑓 + 1 = 𝑇 𝑃 52 𝑢 𝑠 𝑥	-2 + 𝑇 𝑃 53 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 55 𝑝 𝑓 + 2 + 𝑇 𝑃 56 𝑢 𝑠 𝑦	-2	(B.2d)
				[ -𝑇 𝑃 12 -𝑇 𝑃 16 -𝑇 𝑃 42 -𝑇 𝑃 46 ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ ]	{ 𝑢 𝑠 𝑥 𝑢 𝑠 𝑦	-2 -2	}	=	⏟⏞⏞⏞ ⏟⏞⏞⏞ ⏟ [ 𝑇 𝑃 15 𝑇 𝑃 13 𝑇 𝑃 45 𝑇 𝑃 43 ]	{ 𝑝 𝑓 + 2 𝑛 2 𝑢 𝑡 +	}	(B.3)
						𝑇 1				𝑇 2
	Thus, we have
						{ 𝑢 𝑠 𝑥 𝑢 𝑠 𝑦	-2 -2	}	= [𝐓 𝟏 ] -1 [𝐓 𝟐 ] ⏟⏞⏞⏞⏟⏞⏞⏞⏟	{ 𝑝 𝑓 + 2 𝑛 2 𝑢 𝑡 +	}	(B.4)
										[𝐓 𝐮 ]
	the 𝑢 𝑠 𝑥 equations eqs. (B.2b) and (B.2d), we get: -2 and 𝑢 𝑠 𝑦 -2 are expressed by the 𝑝 𝑓 -2 and 𝑢 𝑡 𝑛	-2 with the matrices 𝑇 𝑢 . Then, replacing terms 𝑢 𝑠 𝑥	-2 and 𝑢 𝑠 𝑦	-2	in the left two
	𝑢 𝑡 𝑛	+ 1 = 𝑇 𝑃 32 (𝑇 𝑢 11 𝑝 𝑓 + 2 + 𝑇 𝑢 12 𝑢 𝑡 𝑛	+ 2 ) + 𝑇 𝑃 33 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 35 𝑝 𝑓 + 2 + 𝑇 𝑃 36 (𝑇 𝑢 21 𝑝 𝑓 + 2 + 𝑇 𝑢 22 𝑢 𝑡 𝑛	+ 2 )	(B.5a)
	𝑝 𝑓 + 1 = 𝑇 𝑃 52 (𝑇 𝑢 11 𝑝 𝑓 + 2 + 𝑇 𝑢 12 𝑢 𝑡 𝑛	+ 2 ) + 𝑇 𝑃 53 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 55 𝑝 𝑓 + 2 + 𝑇 𝑃 56 (𝑇 𝑢 21 𝑝 𝑓 + 2 + 𝑇 𝑢 22 𝑢 𝑡 𝑛	+ 2 )	(B.5b)
	{ 𝑝 𝑓 + 1 𝑢 𝑡 𝑛 + 1	}	=	[ 𝑇 𝑃 32 𝑇 𝑢 11 + 𝑇 𝑃 35 + 𝑇 𝑃 36 𝑇 𝑢 21 𝑇 𝑃 32 𝑇 𝑢 12 + 𝑇 𝑃 33 + 𝑇 𝑃 36 𝑇 𝑢 22 𝑇 𝑃 52 𝑇 𝑢 11 + 𝑇 𝑃 55 + 𝑇 𝑃 56 𝑇 𝑢 21 𝑇 𝑃 52 𝑇 𝑢 12 + 𝑇 𝑃 53 + 𝑇 𝑃 56 𝑇 𝑢 22	]	{ 𝑝 𝑓 + 2 𝑛 2 𝑢 𝑡 +	}	(B.6)
					⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟
										[𝐓 f-P-f ]

the matrix [𝐓 f-P-f ] is obtained, whose coefficients are all original from [𝐓 𝐏 ].

  section 3, this matrix is a subset of the original Biot's matrix [𝐓 𝐏 ] where two lines (index 2 and 6) related to the solid displacement in the fluid media are eliminated. We are able to reproduce four linear equations from this relationship as:

				𝑇 𝑃 11 𝑇 𝑃 12 𝑇 𝑃 13 𝑇 𝑃 14 𝑇 𝑃 15 𝑇 𝑃 16 𝑇 𝑃 21 𝑇 𝑃 22 𝑇 𝑃 23 𝑇 𝑃 24 𝑇 𝑃 25 𝑇 𝑃 16 𝑇 𝑃 31 𝑇 𝑃 32 𝑇 𝑃 33 𝑇 𝑃 34 𝑇 𝑃 35 𝑇 𝑃 36 𝑇 𝑃 41 𝑇 𝑃 42 𝑇 𝑃 43 𝑇 𝑃 44 𝑇 𝑃 45 𝑇 𝑃 46 𝑇 𝑃 51 𝑇 𝑃 52 𝑇 𝑃 53 𝑇 𝑃 54 𝑇 𝑃 55 𝑇 𝑃 56 𝑇 𝑃 61 𝑇 𝑃 62 𝑇 𝑃 63 𝑇 𝑃 64 𝑇 𝑃 65 𝑇 𝑃 66	⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦	⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩	σ+ 𝑦𝑥2 𝑢 𝑠 𝑥 + 2 + 𝑢 𝑡 𝑛 2 𝑦 2 𝑢 𝑠 + 2 𝑝 𝑓 + 𝑥𝑥2 σ+	⎫ ⎪ ⎪ ⎪ ⎭ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬	(B.7)
	as explained in 0 = 𝑇 𝑃 11	σ+ 𝑦𝑥2 + 𝑇 𝑃 12 𝑢 𝑠 𝑥	+ 2 + 𝑇 𝑃 13 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 14	σ+ 𝑥𝑥2 + 𝑇 𝑃 15 𝑝 𝑓 + 2 + 𝑇 𝑃 16 𝑢 𝑠 𝑦	+ 2	(B.8a)
		0 = 𝑇 𝑃 31	σ+ 𝑦𝑥2 + 𝑇 𝑃 32 𝑢 𝑠 𝑥	+ 2 + 𝑇 𝑃 33 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 34	σ+ 𝑥𝑥2 + 𝑇 𝑃 35 𝑝 𝑓 + 2 + 𝑇 𝑃 36 𝑢 𝑠 𝑦	+ 2	(B.8b)
	𝑝 𝑓 + 1 = 𝑇 𝑃 41	σ+ 𝑦𝑥2 + 𝑇 𝑃 42 𝑢 𝑠 𝑥	+ 2 + 𝑇 𝑃 43 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 44	σ+ 𝑥𝑥2 + 𝑇 𝑃 45 𝑝 𝑓 + 2 + 𝑇 𝑃 46 𝑢 𝑠 𝑦	+ 2	(B.8c)
	𝑢 𝑡 𝑛	+ 1 = 𝑇 𝑃 51	σ+ 𝑦𝑥2 + 𝑇 𝑃 52 𝑢 𝑠 𝑥	+ 2 + 𝑇 𝑃 53 𝑢 𝑡 𝑛	+ 2 + 𝑇 𝑃 54	σ+ 𝑥𝑥2 + 𝑇 𝑃 55 𝑝 𝑓 + 2 + 𝑇 𝑃 56 𝑢 𝑠 𝑦	+ 2	(B.8d)

  in this appendix. The bulk contribution 𝐾 𝑏 𝑝𝑖 and the enriched one 𝐾 𝑏 𝑝𝑗𝑎 are written as: 𝐾 𝐼 𝑝21 = 𝐴 21 ∫

	𝐾 𝑏 𝑝𝑖 =	1 𝜔 2 𝜌 𝑖 ∫ Ω 𝑒,𝑖	𝛁𝐍 𝑇 𝛁𝐍dΩ -	Ω 𝑒,𝑖 K ∫ 1	𝐍 𝑇 𝐍dΩ	(C.1a)
	𝐾 𝑏 𝑝𝑗𝑎 = 𝐾 𝑏 𝑎𝑝𝑗 = (-1) 𝑗-1 1 𝜔 2 𝜌 𝑖 ∫ Ω 𝑒,𝐼	𝛁𝐍 𝑇 𝛁𝐍dΩ -	Ω 𝑒,𝐼 K ∫ 1	𝐍 𝑇 𝐍dΩ,	(C.1b)
	and the other matrices arise from the interface terms as:
	𝐾 𝐼 𝑝𝑖 = (-1) 𝑖 𝐴 𝑖𝑖 ∫	𝐍 𝑇 𝐍d𝑆	(C.2a)
				Γ
	𝐾 𝐼 𝑝12 = -𝐴 12 ∫							𝐍 𝑇 𝐍d𝑆	(C.2b)
											Γ
	𝐾 𝐼 𝑝1𝑎 =	( -𝐴 11 + 𝐴 12	)	∫	𝐍 𝑇 𝐍d𝑆, 𝐾 𝐼 𝑎𝑝1 =	( -𝐴 11 -𝐴 21	)	∫	𝐍 𝑇 𝐍d𝑆	(C.2c)
									Γ		Γ
	𝐾 𝐼 𝑝2𝑎 =	( 𝐴 21 + 𝐴 22	)	∫	𝐍 𝑇 𝐍d𝑆, 𝐾 𝐼 𝑎𝑝1 =	( -𝐴 12 + 𝐴 22	)	∫	𝐍 𝑇 𝐍d𝑆	(C.2d)
								Γ			Γ
	(						)
	𝐾 𝐼 𝑎𝑎 =	∑ 𝑗	(-1) 𝑗 𝐴 𝑖𝑗		Γ ∫	𝐍 𝑇 𝐍d𝑆,	(C.2e)

Γ

𝐍 𝑇 𝐍d𝑆,

  ,𝑖 that has been explained in appendix A is the material's modulus for solid phase. The interface terms involving displacement are written as:

					𝛁𝐍 𝑇 𝐍dΩ					(C.3a)
		𝐾 𝑏 𝑢𝑖 =	P ∫		𝐍 𝑇 𝐍dΩ					(C.3b)
						Ω 𝑒,𝑖
	where P 𝐾 𝐼 𝑝𝑖𝑢𝑗 = -𝐾 𝐼 𝑝𝑗𝑢𝑖 = (-1) 𝑖 𝐴 𝑖𝑗 ∫						(C.4a)
	𝐾 𝐼 𝑢𝑥𝑎 =	( 𝐴 31 -𝐴 32	)	∫	𝐍 𝑇 𝐍d𝑆, 𝐾 𝐼 𝑎𝑢𝑥 =	( -𝐴 13 -𝐴 23	)	∫	𝐍 𝑇 𝐍d𝑆	(C.4b)
				Γ					Γ
	𝐾 𝐼 𝑢𝑦𝑎 =	( 𝐴 41 -𝐴 42	)	∫	𝐍 𝑇 𝐍d𝑆, 𝐾 𝐼 𝑎𝑢𝑦 =	( -𝐴 14 -𝐴 24	)	∫	𝐍 𝑇 𝐍d𝑆	(C.4c)
				Γ					Γ

Ω 𝑒,𝑖 𝛁𝐍 𝑇 𝛁𝐍dΩ -𝜔 2 ρ ∫ Γ 𝐍 𝑇 𝐍d𝑆 𝑖 = 1

with 𝜔 0 and 𝜔 ∞ expressed as:

where 𝜌 a , 𝜂 a are the density of air and dynamic viscosity of the fluid. The coupling coefficient γ in Biot's equations eq. ( 7) and solid equivalent density are expressed as:

where ρ11 , ρ12 and ρ22 are the coefficients related to the geometry of the frame, the detailed relationship between 𝜌 a and density of solid phase 𝜌 1 are given:

The dynamic effective fluid compressibility Keq is expressed as:

where 𝛾 represents the ratio of specific heats of air, 𝑝 0 is the referred ambient pressure, 𝑃 𝑟 denotes the Prandtl number.

Concerning the solid phase in the porous materials, the constitutive relationship between stress in vacuum and strain in solid frame is defined classically as:

where strain 𝑒 𝑖𝑗 is related to the gradient of the solid displacement as:

and the structural Lamé coefficient 𝑁 and shear coefficient read:

and we define:

where 𝐸 and 𝜈 are the Young's modulus and Poisson's coefficient of the elastic solid in vacuum, 𝜂 𝑠 is the loss factor of the elastic solid phase with unit imaginary number j.

B METHODOLOGY AND PROCEDURE FOR OBTAINING SPECIFIC TRANSFER AND ADMITTANCE MATRICES

B.1 Obtention of matrix [𝐓 𝐟-𝐏-𝐟 ]

In this appendix, we would find a matrix of 2 × 2 that connects the pressure and total displacement for faces Γ + 1 and Γ + 2 in the coupling f-P-f.