A coupled FEM-MFS method for the vibro-acoustic simulation of laminated poro-elastic shells

Summary A new simulation method for the vibro-acoustic simulation of poro-elastic shells is presented. The proposed methods can be used to investigate arbitrary curved layered panels, as well as their interaction with the surrounding air. We employ a high-order finite element method (FEM) for the discretization of the shell structure. We assume that the shell geometry is given parametrically or implicitly. For both cases the exact geometry is used in the simulation. In order to discretize the fluid surrounding the structure, a variational variant of the method of fundamental solutions (MFS) is developed. Thus, the meshing of the fluid domain can be avoided and in the case of unbounded domains the Sommerfeld radiation condition is fulfilled. In order to simulate coupled fluid-structure interaction problems, the FEM and the MFS are combined to a coupled method. The implementation of the uncoupled FEM for the shell and the uncoupled MFS is verified against numerical examples based on the method of manufactured solutions. For the verification of the coupled method an example with a known exact solution is considered. In order to show the potential of the method sound transmission from cavities to exterior half-spaces is simulated.


Summary
A new simulation method for the vibro-acoustic simulation of poro-elastic shells is presented. The proposed methods can be used to investigate arbitrary curved layered panels, as well as their interaction with the surrounding air. We employ a high-order finite element method (FEM) for the discretization of the shell structure. We assume that the shell geometry is given parametrically or implicitly. In order to reduce noise levels which are harmful for the human health, lightweight poro-elastic materials are often used in a wide range of applications such as building acoustics, automotive, and aircraft interior noise. Typically, porous materials used as dissipative components are introduced in multilayered structures. This is the motivation to develop a new simulation method, where a laminated shell structure made of elastic and poro-elastic materials is discretized by the finite element method (FEM) and the surrounding fluid is discretized by the method of fundamental solutions (MFS). In particular, we focus on situations as shown in Figure 1. We assume that the upper half-space has a sound hard ground and is divided into a bounded interior domain and an unbounded exterior domain by a shell structure. Our main focus is to develop a simulation method in order to calculate the sound transmission from the interior to the exterior domain.
The modeling of the vibro-acoustic behavior of systems including porous materials is far from being trivial. Due to the fact that in most applications only small vibration amplitudes occur, it is usual to use linear models. Therefore, air F I G U R E 1 Problem setting: Sound hard ground surface Ω + , interior fluid domain Ω int , exterior fluid domain Ω ext , shell structure Ω s [Colour figure can be viewed at wileyonlinelibrary .com] volumes are modeled with the acoustic fluid, whereas solid structures are modeled with the linearized theory of elasticity. Poro-elastic materials are described with the dynamic Biot theory. This theory was published in References 1,2 and has been adapted to acoustic applications. Present-day descriptions can be found in References 3 and 4.
Many structural parts can be classified as of shell-type. Their thickness is very small when compared with the other dimensions. In such a situation, it is reasonable to describe the geometry by a curved two-dimensional (2D) surface in space. A recent review of laminated plate and shell models is given in Reference 5. In principle, two types of theories can be distinguished. In the first type the number of parameters is independent of the number of layers. Such theories are termed equivalent single-layer theory. The second type referred to as layer-wise theories, where the number of parameters depends on the number of layers. Such a theory is formed by packages of single-layer shell models coupled at the layer interfaces. Concerning the literature on poro-elastic plate and shell theories, we mention the thin poro-elastic plate theory for the consolidation problem. 6 This theory has been extended to the dynamic problem in Reference 7. In Reference 8, the displacement kinematics are extended to allow for shear deformations. Following the idea of a three-dimensional (3D) resolution, a series expansion in thickness direction by means of monomials has been utilized for single-layer poro-elastic plates in Reference 9. This approach has been extended to layered panels in Reference 10, utilizing a layer-wise modeling. In this approach possible air layers can also be considered. 11 A rigorous derivation of poro-elastic plate and shell theories have been published only recently. 12,13 The FEM is the most popular numerical discretization method for the analysis of shell structures. Due to the geometry, various locking phenomena can occur. To overcome this issue, a huge number of different techniques have been developed. In order to resolve the locking issues, many finite elements are based on mixed variational formulations. Common techniques are reduced integration, assumed natural strains, enhanced assumed strains, and the discrete strain gap method (cf. [14][15][16]. Another possibility is the use a high-order ansatz functions to reduce the locking effects. 17,18 Typically, the geometry of the shell reference surface is approximated. However, in the present work, we make use of the exact geometry descriptions within the FEM. For the case of a parametrized surface we mention, [19][20][21] whereas in References 21-23 exact geometry methods for implicitly defined surfaces have been developed. The FE discretization of poro-elastic media can, e.g., be found in Reference 3. For sophisticated solution techniques of the coupled poro-elastic system we refer to 24 and 25. We assume that in the geometry definition of a problem only the shell geometry structure is modeled directly. The geometry of the fluid domains is described indirectly by their boundary, that is, by the boundary of the half-space and the shell structure. In such a situation it can be a rather difficult task to generate a volume mesh. To tackle this problem, one possibility is to use embedded/fictitious domain methods. These methods are based on the idea of defining an auxiliary domain which can be meshed easily. As recent contributions in this field, we would like to mention the finite cell method, 26,27 the Cartesian grid method, 28 and the CutFEM. 29 Recent developments in the context of constructive solid geometry modeling can be found in References 30 and 31. An interesting application of this concept to shell analysis is given in Reference 32. Another possibility for the analysis of domains described by their boundary is to resort to boundary related methods, where no volume mesh is required. The most developed method of this type is the boundary element method (BEM), see e.g. 33 . In the direct BEM, the unknown boundary data are discretized, whereas in the indirect BEM, an auxiliary density function defined on the boundary is sought. This density function is used to describe the solution field. The indirect BEM can be seen as a Trefftz-type method. In such methods, the solution is approximated by a linear combination of basis functions, which fulfill the underlying partial differential equation. Reviews on Trefftz-type methods along with their classification can be found in 34 and 35. A Trefftz-type method which was developed with the aim to solve vibro-acoustic problems is the Wave Based Method (WBM) developed in Reference 36. Since its first publication, this method has been constantly developed further. 37 Within the WBM, the solution is determined by means of a variational formulation. Despite that the BEM would be optimal with respect to the model of the surrounding air, here, a more simpler approach, the MFS is used. This Trefftz-type method (see, eg, 38) uses fundamental solutions of the underlying partial differential equation for the solution approximation. Typically, the solution is determined by collocation at the boundary. Due to the use of fundamental solutions, the Sommerfeld radiation condition is fulfilled exactly. Thus, unbounded domains are easily treatable. We mention three applications of the MFS in the field of poro-elasticity. Nennig et al 39 applied the MFS to scattering problems from poro-elastic bodies in two dimensions. Wen and Liu 40 derived the fundamental solution for a poro-elastic plate in Laplace domain and applied the MFS for solving boundary value problems. Augustin 41 presented density results and a MFS for quasistatic poro-elasticity.
In the present work, we couple the FEM and the MFS in order to benefit from both worlds. The coupling of different numerical methods is a well known approach for acoustic-structure interaction problems. The FEM is perfectly suited for models with complex geometries while boundary related methods enable to account for the radiation of waves in domains of infinite extent. The coupling of FEM-BEM is well-known in literature, see, e.g. 42,43, 44 among others for structural-acoustic coupling. The coupling of the FEM and the WBM has been proposed for different configurations. We mention the cases of structural (FEM) -acoustic fluid (WBM) coupling 45, acoustic fluid (WBM) -poroelastic domain (FEM) coupling 46 and the coupling of two different poroelastic domains 47. For two-dimensional soil-structure interaction problems a coupled FEM-MFS schema has been developed in 48. Thus, two elastic sub-domains are coupled. This was extended to a 2.5-dimensional model for the prediction of vibrations due to underground railway traffic in 49. Furthermore, we mention 39, where an acoustic fluid domain is coupled with a poro-elastic domain. Both are discretized using the MFS. To our best knowledge, the coupling of a poro-elastic shell with an acoustic fluid was not considered before.
In the present work, the main novelty is the development of a layer-wise poro-elastic shell theory. The individual layer displacement kinematics are based on a seven-parameter shell model, 19,50,51 whereas the scalar fluid pressure field is assumed to be quadratic through-the-thickness in each poro-elastic layer. The shell reference surface can be given either parametrically or implicitly. In both cases the exact geometry is incorporated in the FEM following the developments in Reference 21. The field approximation for the shell is done by means of arbitrary order hierarchical shape functions. Furthermore, a variational variant of the MFS in order to discretize the fluid domains is used. Schemata for the coupling of the FEM and the MFS at different interfaces are presented. While the coupling formulation at fluid-elastic solid interfaces are known in literature, 45 to the best of our knowledge, the coupling formulation of a poro-elastic shell with fluid is new. The implementation of the method is verified with high rigor. To that end, the method of manufactured solutions is adapted to the case of poro-elasticity for the first time. Finally, the proposed method is used to simulate the sound transmission from the inside of two cavities bounded by a poro-elastic shell structure to the infinite outer domain.

GOVERNING EQUATIONS
In this section, the governing equations of an acoustic fluid, an elastic solid and a poro-elastic solid are briefly introduced.
Harmonic time dependency of all fields is assumed where the angular frequency is denoted with .

Acoustic fluid
The propagation of sound waves in an fluid is governed by the Helmholtz equation 52 where p stands for the sound pressure fluctuation around the equilibrium state, k = √ K is the wave number, g the source term, and Δ the Laplace operator. The material parameters are the density and the bulk modulus K. A pure acoustic boundary value problem Equation (1) is supplemented by the boundary conditions On the Dirichlet boundary Γ D , the pressure has the prescribed value g D , whereas on the Neumann boundary Γ N the normal velocity n p(x, ) i has the prescribed value g N . In the case of an unbounded domain, the solution has to fulfill the Sommerfeld radiation condition first considered in Reference 53.

Elastic solid
The propagation of waves in elastic solids is governed by the following equations Here, u denotes the displacement vector, stands for the stress tensor, b is the volume load density, is the infinitesimal strain tensor, C denotes the fourth-order elastic tensor, and s is the solid density. Here, we consider linear isotropic materials where the elasticity tensor is given by where I is the second-rank identity tensor, and I is the symmetric part of the fourth-rank identity tensor. E denotes the Young's modulus and is the Poisson's ratio. The boundary conditions are given by On the Dirichlet boundary Γ D the displacement vector has the prescribed value u D . On the Neumann boundary Γ N the surface traction vector has the prescribed value t N i .

Poro-elastic solid
Following Biot's theory, wave propagation in a poro-elastic solid is governed by the equations Here, tot stands for the total stress tensor =(1− ) s + f is the bulk density, f is the fluid density, denotes the porosity, R is a poro-elastic parameter, and is the so-called effective stress coefficient. 54 To obtain (9) the relative fluid to solid displacement with the fluid displacement U has been eliminated. In (9) and (11) the abbreviation is used, where a is the apparent mass density and b is the viscous drag force. For typical poro-elastic materials used in acoustic applications, the bulk modulus of the elastic solid K s is very large as compared with the bulk modulus of the fluid K f and the bulk modulus of the skeleton K. Thus, we assume an incompressible solid skeleton material which results in Following, 55 the frequency-dependent viscous drag is given as Here, denotes the static flow resistivity and ∞ is the tortuosity. Λ denotes the viscous characteristic length, and f is the dynamic viscosity. In order to take thermal effects into account, Champoux and Allard 56 introduced the thermal characteristic length Λ ′ . Following this approach, the bulk modulus of the fluid becomes frequency dependent Here, p 0 is the ambient pressure of air, is the ratio of specific heats and Here, Pr denotes the Prandtl number with the specific heat capacity at constant pressure c p and the thermal conductivity . The boundary conditions for a poro-elastic solid are given by On the Dirichlet boundary for the solid Γ u D the displacement has a prescribed value u D , whereas on the Neumann boundary Γ u N the total surface traction vector is prescribed with t N . Furthermore, on the Dirichlet boundary for the fluid Γ p D the fluid pressure has a prescribed value p D , whereas on the Neumann boundary Γ p N the normal flux F I G U R E 2 Coupled continua: Acoustic fluid continuum Ω a , elastic solid continuum Ω e , poro-elastic solid continuum Ω p can be two independent nonoverlapping decompositions of the boundary. Thus, on a boundary point two conditions have to be fulfilled.

Coupling conditions
In this section, we outline the coupling conditions between the physical models. In Figure 2, the abstract setting of coupled continua is depicted. The domains Ω a , Ω e , Ω p refer to an acoustic fluid, to an elastic solid and to a poro-elastic solid, respectively. The boundary of Ω a is denoted by The coupling boundary between Ω e and Ω p is Γ ep . Before we consider the coupling of different models, we examine the coupling of two elastic continua. The coupling of two elastic continua (Ω 1 ,Ω 2 ) over the common interface Γ yields two conditions. The first condition which has to hold is the continuity of the displacements The second condition is the equilibrium of forces, which results in In the case of an acousticfluid -elasticsolid interface Γ ae also two coupling conditions are necessary. According to the inviscid assumption for the acoustic fluid, no shear stresses occur in the fluid. Hence, particles of the fluid can move tangential to the interface without resistance. Therefore, only the normal displacements are continuous The fluid displacement can be expressed through the pressure by u a = a 2 ∇p. Thus, the coupling condition gets a 2 ∇p a ⋅ n = u e ⋅ n on Γ ae .
The traction on the interface resulting from the scalar pressure field in the fluid is t a = −p a n. Thus, the equilibrium of forces yields t e = −p a n on Γ ae , where the normal vector is the outward normal vector of the elastic domain.
On an acousticfluid -poro−elasticsolid interface Γ ap , three coupling conditions have to be fulfilled. The continuity of normal displacements implies where u p =(1− )u+ U is the displacement of a "poro-elastic particle." Using the relative displacement between solid and fluid defined in (11), the coupling condition gets The second condition is given by the equilibrium of forces. The surface traction induced by the scalar pressure in the acoustic fluid has to be balanced with the total traction in the poro-elastic solid − p a n = t tot on Γ ap , where the normal vector is the outward normal vector of the poro-elastic domain. The third condition ensures the continuity of the pressure fields. Hence, The remaining combination is an elasticsolid -poro−elasticsolid interface Γ ep . The elastic solid and the poro-elastic solid are able to resist shear forces. Therefore, no relative motion between the two solid phases are allowed The elastic domain represents an impervious interface for the fluid in the poro-elastic solid. Therefore, the relative mass flux normal to the interface has to be zero As a third condition, the equilibrium of forces demands

PORO-ELASTIC SHELL STRUCTURE
In this section, we describe a layer-wise theory for laminated poro-elastic shells. We consider shell reference surfaces Ω which are given parametrically or implicitly. In the former case we have a parametrization ∶Ū ⊂ R 2 → Ω available with given parameter domainŪ. In the latter case, the reference surface is given as the zero-level set of a function ∶ R 3 → R inside a cuboid B, For this setting, the normal vector to the surface is given bỹ In the numerical method, we will make use of a piecewise parametrization of the exact surface over a space triangulation. Therefore, the implicit description of the reference surface is turned into a parametric one (cf. 21,23). Given the parametrization , we can define the two covariant base vectors G ∶= , which span the tangent plane to Ω. Here and in the following, Greek indices take the values 1 and 2 and Latin indices the values 1, 2, 3. With the base vectors we can define the unit normal vector

F I G U R E 3
Parametrization of the shell. The parameter space on the left is mapped to the physical space on the right. The reference surface is parametrized by , whereas the shell volume is parametrized by g [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 4 Layup of the shell
Having defined the normal vector the parametrization of the 3D shell volume Ω is given by with the one-dimensional (1D) thickness interval T. The geometric setting is shown in Figure 3. We consider a layered shell structure with a layup as depicted in Figure 4. The total number of layers is L. Each layer has a thickness t i , i=1… L and is classified as elastic or poro-elastic. The distance from the reference surface to the bottom of the shell is denoted by t 0 . Furthermore, we define the individual layer thickness intervals The relations between the local thickness coordinate i ∈[0,1] of layer i and the global thickness coordinate 3 The present layer-wise theory is based on the assumed seven-parameter displacement field for layer ∈{1,… ,L} of the form Here, V( ( 3 )) are through-the-thickness functions and u i , u i , u are seven local displacement parameters and e i are the base vectors of a Cartesian coordinate system. The through-the-thickness functions are chosen as (1) Thus, the parameters u i = u i ( 1 , 2 ) and u i = u i ( 1 , 2 ) describe the displacement of the bottom ( =0) and the top ( =1) of layer , respectively. This is enhanced with the parameter u i = u i ( 1 , 2 ), which accounts for a quadratic variation of the displacement in thickness direction, which vanishes at the bottom and the top surface of the layer. In all possible combinations of elastic and poro-elastic layers, the continuity of the displacement field is required. Thus, for two subsequent layers, we set Therefore, the total number of parameters describing the displacement is 3(L+1)+L and is independent of the stacking sequence. Following, 8 the fluid pressure field in poro-elastic layers is approximated by a quadratic expansion through-the-thickness. The pressure field for layer is assumed to be of the form are three local parameters. Between two poro-elastic layers, we require the continuity of the pressure. Thus, if layer and +1 both are poro-elastic, we set Therefore, the total number of parameters describing the fluid pressure is 3L p −n p , where L p is the number of poro-elastic layers and n p is the number of interfaces between two poro-elastic layers. Since the considered coupling of elastic and poro-elastic layers is natural, 3 it is sufficient to enforce the continuity of the displacement field and the pressure field, as it is done in (40) and (42).

FEM for the poro-elastic shell
The FEM developed for the layered poroelastic shell structure is based on the weak formulation of the respective governing Equations (4) and (9) and on the through-the-thickness discretizations (38) and (41).

Weak form of elastodynamics
The weak form of elastodynamics is obtained by multiplying (4) with a vector-valued test function u ∈ V u 0 = {u ∈ [H 1 (Ω)] 3 | u = 0 on Γ D }, integrating over the domain Ω and applying integration by parts we obtain where Thus, the weak form of the problem reads: (43) is fulfilled for all u ∈ V 0 .

Weak form of poro-elasticity
The derivation of the weak form of poro-elasticity follows the same arguments as for elastodynamics. We multiply (9a) with vector-valued test function u ∈ V u 0 , integration over the domain and apply integration by parts, Furthermore, we multiply (9b) with a scalar test function p ∈ V p 0 = {p ∈ H 1 (Ω)|p = 0 on Γ p D } and perform an integration over the domain. This yields with some rearrangements Now, we apply integration by parts to the first integral and using (19) in order to obtain Thus, we can formulate the weak form of the poroelastic boundary value problem: are fulfilled for all test functions (u, p) ∈ V u 0 × V p 0 . The bilinear and linear forms in (48) are defined in (44) and

Discretization
The discretization is based on the standard reference element technique, that is, on the reference element hierarchical shape functions of arbitrary order and quadrature are defined. 57 Thus, the shell reference surface is subdivided into n e nonoverlapping geometric elements e such that Ω = ∪ n e e=1 e . In the case of an parametrically defined reference surface the parametric domain is subdivided into quadrilaterals Q e and the geometric elements in real space are given by e = (Q e ), see Figure 5. The standard affine mapping from the reference element R to Q e is denoted by Φ e . For implicitly defined reference surfaces the parametrization is not explicitly available. Therefore, the exact geometry is parametrized over a piecewise flat triangulation with elements e in real space, see Figure 6. The elements in real space are obtained by e =a e ( e ). Here, a e is a mapping which is only implicitly defined, 21,23 that is, for the evaluation a 1D nonlinear root finding problem has to be solved. For further details on the evaluation of the integrals in (44) and (49) we refer to Reference 58.

MFS for a single fluid domain
In this section, we will introduce a variant of the MFS, which is a Trefftz-type method. This method is used to discretize the acoustic fluid surrounding the shell structure. In general, a Trefftz method consists of a discrete Trefftz space and a Trefftz variational formulation. Usually, the MFS refers to a particular combination of discrete space and variational formulation. The ansatz functions are fundamental solutions (Green's functions) with source points placed along a curve (in 2D) or on a surface (in 3D) surrounding the computational domain Ω. It is necessary that the source points lie in the complement of Ω, since the fundamental solutions are singular at the source points. The coefficients are determined by collocation. Here, we want to solve the boundary value problem of acoustics (2). The well known fundamental solution of the Helmholtz operator in the full 3D space is given by The point x is called field point, whereas y is called source point. We restrict our further considerations to half-space problems. This means that we assume that the upper half-space Ω + = {(x, y, z) ∈ R 3 |z > 0} is divided by the shell structure Ω s into two fluid domains Ω int and Ω ext , such that Ω + =Ω int ∪Ω ext ∪Ω s . We assume a sound hard surface Ω + at z=0. Thus, the normal fluid velocity vanishes, v a (x) ⋅ e 3 = 0, for x ∈ Ω + .
By making use of the well known half-space fundamental solution, the condition at sound hard surface Ω + is fulfilled exactly. The half-space fundamental solution is given by In the present work, we use fundamental solutions placed on a surface embracing the fluid domain to approximate the pressure field in the acoustic fluid, as it is done in the classical MFS. The unknown coefficients are determined by a variational formulation like in the WBM 37. To this end, the residuals are introduced. The solution in the fluid domain is approximated by a linear combination of fundamental solutions  weighted with coefficients c i In order to treat nonhomogeneous problems, that is, where sources are located inside the acoustic fluid domain, we extend (54) to where F p are particular solutions of the in-homogeneous problem. In the present work, we consider point sources only. Thus, F p is a linear combination of fundamental solutions where d s is the strength of an acoustic source at y s . The discrete fluid displacement is given by where (x, y) = ∇ x (x, y). Thus, In the proposed MFS either the interior acoustic problem in Ω int or the exterior acoustic problem Ω ext can be tackled. In both cases, we assume that the boundary Γ of the respective domain can be decomposed into Γ=Γ D ∪Γ N ∪ Ω + . Due to the ansatz (55), the Helmholtz equation in Ω + and the hard wall condition at Ω + are fulfilled exactly. Nevertheless, the other boundary conditions cannot be fulfilled exactly, yielding the residua (53). Therefore, the unknown coefficients c i in (55) are determined in a weighted residual sense. To this end, the boundary residua are weighted with the complex conjugate of the test functions and their normal derivative, Applying a Galerkin approach, the test functions are chosen to be fundamental solutions with the same source points as in (54). This leads to the linear system of equations where matrix and vector entries are given by

Coupled method
In this section, we present a coupling approach for the FEM and the MFS.

Acoustic-elastic FEM-MFS coupling
We consider the case of an elastic solid which is in contact with an acoustic fluid at the interface Γ ae . The presented formulation is similar to the one in 45. The variational formulation of the governing equations for an elastic solid is rewritten to Here, we have assumed the case Γ=Γ N ∪Γ D ∪Γ ae . The case of an additional interface between an elastic solid and a poro-elastic solid was commented in Section sec::poroShell. The two coupling conditions for the interface Γ ae stated in (23) and (24) are rewritten to t e = −p a n, on Γ ae , and 1 a 2 ∇p ⋅ n − u e ⋅ n = 0, onΓ ae .
The incorporation of the first coupling condition in (63) yields Following the variational MFS approach in Section 4.2, the second coupling condition is weighted with the complex conjugate of the test function and integrated over the coupling interface The discretization of (65) and (67) leads to a system of equations in the form of The matrices M u and K u are the mass and the stiffness matrix and f V +f N is the load vector, which arise from the FEM discretization. In order to state the entries of the other matrices, we use the global index i=i(i 1 ,i 2 ,i 3 ,l), which refers to the respective finite element function N i 1 (x) Here, N i 1 denotes a 2D finite element shape function. Thus, Using Δ(x, y) = 0 ∀ x ≠ y and applying integration by parts two times, and we obtain that the matrix H is symmetric. Therefore, the system matrix of (68) is symmetric. We remark that the integrals are transformed to the reference element for their numerical evaluation by means of a quadrature rule.

Acoustic-poro-elastic FEM-MFS coupling
We proceed with the coupling in the case of an acoustic fluid where the pressure field is approximated by an MFS ansatz and a poroelastic solid, which is discretized by the FEM. Our formulation is different to the one in 46, with the feature of yielding a symmetric system matrix. We rewrite the variational formulation of the governing equations for a poroelastic solid stated in (48) to , u). Here, we have assumed that Γ=Γ N ∪Γ D ∪Γ ap . We incorporate the coupling condition (27) and (28) in the first equation of (72), which gives for the integral over Γ ap The coupling condition (26) is incorporated in the second equation of (72). For the integral over Γ ap , we obtain Following the variational MFS approach in Section 4.2, the coupling condition (28) is weighted with the gradient of the complex conjugate of the test function and integrated over the coupling interface After discretisation we have the symmetric system of equations where the newly introduced matrices and vectors are

Verification
The aim of this section is to verify the implementation of the methods developed in Section 4. In order to ensure the reliability of a numerical simulation software, verification and validation (V&V) are unavoidable tasks. 59,60

Verification of the poro-elastic shell FEM
In this section, we apply code verification to the poro-elastic shell FEM developed in Section 4.1 based on order-of-accuracy tests and on the method of manufactured solutions (MMS). [61][62][63][64][65] The necessary prerequisite to apply an order-of-accuracy test to a numerical schema is the knowledge of a formal order of convergence and exact solutions. Thus, one has to know an estimate of the type where C is a constant and h is a characteristic discretization parameter. Here, we refer to a characteristic element size. Then q is called the formal order of convergence with respect to the norm ||⋅||. For two meshes with characteristic element sizes h 1 and h 2 , the experimental order of convergence (eoc) is defined as where is the numerical error corresponding to the discretization h i . The code is verified, if the eoc matches the formal order of convergence within the asymptotic range. For the FEM with arbitrary ansatz order p applied here, we expect q=p+1 for the error in the L 2 norm for smooth solutions.
In the case of the elastic solid problem (8), the procedure above is straightforwardly applicable. The implemented finite element code needs the components b i of the source term b=b i e i with respect to the global Cartesian frame. In the case of the poro-elastic solid problem (9), we include artificial source terms in the formulation, in order to apply the MMS. We modify the respective Equation (9) We remark that these source terms have no physical meaning. They are incorporated in the variational formulation and in the FEM easily, leading to additional entries in the load vectors.

Verification examples
We have checked the order of convergence for a number of examples, considering different material parameters, frequencies, geometries, displacement fields, and pressure fields. In all examples, the optimal asymptotic order of convergence could be observed.
In this section, we show the results of four verification examples in total. In the first three examples, we use a parametric description of the reference surface, whereas in the fourth example, an implicit description is used. In the first example, we prescribe the trivial displacement solution and a nontrivial pressure solution. In the second example, we make it the other way round.
In both examples, we use the reference surface given by  Table A3 and an angular frequency =20s −1 . For the construction of the manufactured solution, we take the structure of the shell model into account. Therefore, the solution is defined by specifying the parameters in (38) and (41). For the presented examples, these parameters are given in Table 1. The errors at multiple meshes with n e elements each are shown in the Figures 7 to 9. In all graphs, the respective absolute error defined by is plotted. Therein, u and p denote the numerical solution, whereas u M and p M denote the manufactured solutions. We observe the optimal convergence rate in all examples. Furthermore, we see from Figure 8 that a small discretization error TA B L E 1 Parameters for the manufactured solutions in the displacement field leads to large error in the pressure field. This is due to the conditioning of the physical problem and depends on the material parameters. Therefore, we conclude that an accurate discretization of the displacement field is necessary in order to obtain an accurate pressure approximation. Next, we discuss the results of the fourth verification example for the poro-elastic shell FEM in the case of an implicitly given reference surface. The considered spherical reference surface is given by =x 2 +y 2 +z 2 −1 and B=[0,2] 3 . The shell volume has the extension t in the thickness direction and is symmetric around the reference surface. We prescribe the solution as It is important to note that this solution cannot be exactly represented by the shell model. Therefore, a modeling error and a discretization error occurs. We use quintic shape functions for the discretization of a series of problems with decreasing thickness T. The results are depicted in Figure 10. We observe that the modeling error dominates in the case of thick shells. Therefore, a mesh refinement cannot reduce the overall error in this case. However, with decreasing thickness, the modeling error decreases and the discretization error dominates for the coarse meshes. In this regime, we observe the expected order of convergence of the FEM. In order to verify the layered shell theory we consider a fixed shell thickness T = 0.2 and use a varying number of layers n L through-the-thickness and different numbers of finite elements. The results of the convergence study are depicted in Figure 11 With increasing number of layers the individual layer thickness t is reduced, which reduces the modeling error. We observe optimal quadratic convergence with respect to the layer thickness.

Verification of the MFS implementation
This section deals with the verification of the MFS for uncoupled acoustic fluid problems described in Section 4.2. To this end, we present the results of two examples. In both cases, the considered fluid is air with the material parameters given in Table A1. In the first example, an interior problem is considered. The domain Ω int is bounded by the plane Ω + and a parametrically given surface Γ, x = cos( 1 ) sin ( 2 ), where 1 ∈[0,2 ] and 2 ∈ [0, 2 ], see Figure 12. The constructed solution is obtained by means of the fundamental solution. Therefore, we specify the point y 0 =[1m,0m,1m], which lies outside the problem domain. This defines the sought solution according to The boundary data is derived from this solution as This boundary data is the input for the numerical method. The source points for the approximation are obtained by placing points on Γ and moving each 0.3m in the direction normal to Γ. The problem setting and a source point configuration are depicted in Figure 12. We investigate the acoustic problem at a frequency f =[200]Hz. In order to study the convergence behavior we introduce the relative error We evaluate this error for a number of solutions obtained with different number of source points. Furthermore, we study the influence of the numerical integration. The computed errors are plotted in Figure 13 for different number of quadrature points n g used. It is evident that the integration has to be sufficiently accurate in order to obtain an accurate result. Therefore, one has to increase also the number of quadrature points when increasing the number of source points in order to obtain a monotonic convergence. Nevertheless, very accurate solutions are possible with only a few source points for the solution approximation.
In the second verification example, we consider an exterior problem. The unbounded problem domain is given as {(x,y,z)| (x,y,z)>0}∩Ω + , where (x, y, z) = ((x 2 + y 2 − 1) 2 + ( The geometry of the problem and a source point configuration is depicted in Figure 14. We use y 0 =[0.7m,0.7m,0.25m] for the construction of the solution and the boundary data by means of (88) and (89). Again, we solve the acoustic problem at f =[200]Hz. The convergence behavior of the error (90) is plotted in Figure 15. We observe that the MFS is able to reproduce the solution accurately in the case of this geometrically complex problem as well. As in the previous example, we see that a sufficient number of quadrature points is necessary in order to obtain stable results and a monotonic convergened method. Therefore, the number of source points and quadrature points has to be chosen with care, when applying the MFS.

Verification of the coupled method
In this section, we are concerned with the verification of the coupled MFS-FEM developed in Section 4.3. To this end, we consider radial symmetric problems, which allow for a closed form solution. Thus, we consider a spherical shell structure separating the upper half-space into an interior and an exterior domain. The exact solutions of the 3D problems are derived in Appendix B1. In our FEM model we approximate the spherical shell structure with a shell theory, hence, a modeling error is introduced. The reference surface for the shell structure is a hemisphere with unit radius. First, we consider the case of an elastic aluminum structure with a thickness of t=0.025m. Therefore, the ratio of curvature radius and thickness is 40, classifying the structure as thin shell. The used material parameters are given in Table A2. In this example, an implicit geometry description of the hemisphere is used. The displacement of the shell is discretized with quartic finite element shape functions. The interior and exterior fluid pressure is discretized by means of the MFS. The error according to (90) is evaluated for the exterior acoustic fluid and plotted in Figure 16. The error is evaluated for varying finite element meshes and for varying number of MFS source points. The finite element meshes are identified by their number of elements n e . Depending on the used discretization, the error is dominated either by the FEM or the MFS. For horizontal lines the error is dominated by the FEM discretization. This verifies the coupled method for the case of the acoustic fluid-elastic solid coupling.
Next, we consider a poro-elastic polyurethane shell structure separating the interior domain from the exterior domain. The material parameters are given in Table A3. The spherical shell is described parametrically and has a thickness of t=0.002m. Quartic finite element shape functions are used for the discretization of the displacement and pressure field of the shell structure. The error according to (90) is evaluated for the exterior acoustic fluid and plotted in Figure 17. The numerical integration for all matrix entries is done on the finite elements. In the case of n e =4, we can see the instability, which arises due to the insufficient numerical integration. However, we do not observe such instability in the case of the other meshes. In the present example, we see that due to the modeling error the minimal achievable error is around 10 −9 . Up to this error, we observe the convergence of the coupled method. Hence, this examples verifies the coupled method for the case of acoustic fluid − poro-elastic solid coupling.

Sound transmission through poro-elastic shells
In the previous section, the implemented numerical methods were verified. Here, the capabilities of the developed methods are shown on the basis of more complex examples. We investigate the sound transmission from the inside of two cavities bounded by a poro-elastic shell structure to the outside. In one example the reference surface of the shell is given parametrically, whereas in the other it is given implicitly.
with the mapping The parameter domain is ( 1 , 2 )∈[0,1]×[0,1]. In Figure 18, the geometry of the problem is visualized. The shell structure is composed of two layers. We assign the material parameters of aluminum (see Table A2) to the first layer, which has a thickness of 0.01m. The second layer is a poro-elastic polyurethane (see Table A3) layer with a thickness of 0.03m. This poro-elastic layer is in contact with the interior fluid.
For the sake of solution verification, we consider the uncoupled dynamic response of the shell structure as a result of a surface load t = −e z 10 3 N∕m 2 applied on the free surface of the aluminum layer. For the analysis at 250 frequencies in the range [[0]Hz, [500]Hz] we use sextic ansatz functions. The vertical displacements u z at the point ( 1 , 2 )=(0.5, 0.25) are plotted over the frequency in Figure 19 for three meshes (each mesh has n e elements). For the better visualization, we have used the logarithmic measure L p (u) = 10 log 10 ( u 2 4 ⋅ 10 −10 ) dB.
Due to these results, we will further use the 64-element mesh. Next, the performance of the uncoupled MFS is studied. To this end, we calculate the interior pressure field due to a source with unit strength at the point [0m, −1m, 0.5m]. The structure is assumed to be rigid in this case. The exact solution to this problem has to be real valued. In Figure 20, the imaginary part of the solutions at the evaluation point [0m, 1m, 0.5m] obtained for a varying number of approximation source points are plotted. Based on this result, we proceed with 372 source points for the interior, as well as for the exterior fluid.
With the discretization described above, we study the sound transmission from the interior to the exterior. Again, we consider a source at the point [0m,−1m,0.5m] with unit strength for the excitation of the system. In Figure 21, the sound pressure level determined at the interior evaluation point [0m,1m,0.5m] is plotted. We compare the cases uncoupled MFS (rigid structure), aluminum shell (only the aluminum layer), and poro-elastic shell (the aluminum layer with the polyurethane layer). The results for the cases uncoupled MFS and aluminum shell virtually agree. In the case of the uncoupled MFS, no dissipation occurs in the system and the solution is infinite at the eigenfrequencies. In the case of aluminum shell a small structural dissipation effect is present. This can be seen in Figure 21, where the eigenfrequencies are damped. However, we conclude that the compliance of the aluminum structure has only little influence on the interior sound pressure field. Nevertheless, the effects of the dissipation introduced by virtue of the poro-elastic polyurethane layer are clearly visible. In particular, the calculated sound pressure levels are significantly reduced for frequencies above f =[200]Hz. The same conclusion is valid for the results at the exterior evaluation point [0m,0m,1.5m]. The sound pressure level determined at the exterior evaluation point [0m,0m,1.5m] is plotted over the frequency in Figure 22. In the case of an rigid structure no transmission can occur. Therefore, this case is not considered in Figure 22. Again, the sound pressure level is significantly reduced in the case poro-elastic shell compared with the case aluminum shell.

Sound transmission through an implicitly given shell
In this section, we consider the sound transmission through a shell structure described by an implicitly given reference surface. The corresponding level-set function is (x, y, z) = (2x) 4 + (2y) 2 + (2z) 2 − 0.25.
The geometry of the problem is depicted in Figure 23. The shell structure is composed of an aluminum (see Table A2) and a poro-elastic polyurethane layer (see Table A3). The aluminum layer is t=0.002m thick, whereas the polyurethane layer has a thickness t=0.01m. The poro-elastic layer is in contact with the interior fluid. In this example, we have chosen the discretization according to similar considerations made in Section 5.2.1. We use octic shape functions constructed on a 64-element mesh for the discretization of the parameters in the shell model. For the discretization of the interior and exterior fluid pressure fields, we use 145 MFS source points each. The sound pressure level at the interior evaluation point [0m, 0m, 0.1m] due to a source at the point [−0.25m, −0.05m, 0.1m] with unit strength is plotted in Figure 24. Therein, the case of a rigid structure (uncoupled MFS), an aluminum shell and the full poro-elastic structure is considered. We conclude that the compliance of the aluminum structure has only little influence on the interior sound pressure field. However, the effects of the dissipation introduced by virtue of the poro-elastic layer are visible. The same conclusion is valid for the results at the exterior evaluation point [0m, 0m, 0.3m], see Figure 25.

CONCLUSIONS
We have presented a new simulation method for vibro-acoustic analysis. In particular, we have treated the simulation of laminated poro-elastic shell structures and their interaction with the surrounding fluid. In order to face this complicated situation, a layer-wise shell model has been developed. The through-the-thickness variation of the displacement field is described by a seven-parameter model, which is assumed in each layer. The pressure field occurring only in the poro-elastic layers is described with a quadratic expansion through-the-thickness. The numerical model consists of a FEM formulation for the shell and the MFS for the acoustic domains. Hence, the radiation condition is implicitly fulfilled by the MFS due to the usage of the acoustic fundamental solutions. The FEM for the shell is based on the exact geometry given by either a parametric or an implicit representation of the reference surface. In the case of a given parametrization, the evaluation of the quantities from differential geometry are done with respect to this parametrization and not with respect to a superfluous geometry discretization. In the case of an implicitly defined surface, the exact parametrization is constructed by means of the level-set function. Contrary to most MFS solutions, here, the strength of the sources for the field approximation are determined by a variational formulation. This has the advantage that the need of collocation points is circumvented.
The implemented numerical methods are verified against solutions obtained from the MMS. In this method, a distinct solution is chosen and the corresponding source terms and boundary conditions are derived from the chosen solution. For the developed shell FEM, this approach has been implemented in the curvilinear coordinates induced by the shell geometry. For the verification of the MFS, the prescribed solution is constructed by a fundamental solution. Thus, the FEM and the MFS could be verified with high rigor. The coupled method has been verified against radial symmetric problems. The reference solutions were obtained by solving these problems analytically. The calculation of the sound transmission through geometrically complicated layered shells structure has shown the suitability of the proposed approach for real world examples.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.

APPENDIX B. RADIAL SYMMETRIC SOLUTIONS
We derive radial symmetric solutions of two coupled problems. In both cases a spherical shell structure separates a bounded interior fluid domain from an unbounded exterior fluid domain. In particular, we consider the cases acousticfluid -elasticlayer -acousticfluid and acousticfluid -poro−elasticlayer -acousticfluid. We employ the spherical coordinates (r, , ). Thus,