Wave Finite Element Formulation of the Acoustic Transmission Through Complex Infinite Plates

Summary A ﬁ nite element-based derivation of the transmission loss (TL) of anisotropic layered inﬁnite plates is presented in this paper. T he wave-ﬁnite element method (WFE) i s u sed to represent the plate with a ﬁ nite element model of a s ingle unit cell. The incident acoustic ﬁeld is a k nown plane wave, a nd the reﬂected and transmitted pressures are supposed to be plane waves w ith unknown amplitudes and phases. The periodicity conditions on the unit cell allow t o ﬁ nd a s imple matrix equation linking the amplitudes of the transmitted and reﬂected ﬁelds as a f unction of the incident one. This approach is validated for several cases against classical analytical models for thin plates and sandwich constructions, where the results agree perfectly for a r easonable mesh size. The method is then used to study the e ﬀ ect of stacking order in a l aminated composite plate. The main interest of the method is the use of ﬁnite elements, which enables a r elative easy modelling since most packages readily include di ﬀ erent formulations, compared to analytical models, where di ﬀ erent formulations have to be implemented for every kind of material.


Introduction
Sound transmission modelling is av ery important topic in the industry,y et it is often still difficult to calculate for complexg eometries or materials. Several methods have been proposed for the computation of transmission loss through infinite plate-likes tructures. Analytic models based on the constitutive equations for infinite thin plates [1], sandwich constructions [2,3,4], and even double plates with periodical stiffeners [5] have been proposed in the past. The Transfer Matrix Method (TMM) [ 6,7] is ag eneral framework used to compute acoustic transmission features of infinite plane isotropic layerings comprising elastic solids, poroelastic materials and air gaps. It relies on an analytical formulation of the problem in each layer,thus making it necessary to derive newequations for every kind of material that can be encountered (isotropic or anisotropic, elastic or poroelastic).
Finite Element (FE) modelling has the practical advantage overanalytical models that most constitutive laws and couplings between materials are already implemented in commercial codes. Howeverm eshing al arge and simple structure such as ap late can lead to high computational times. The Wave-Finite Element (WFE)m ethod is ar e-sponse to this issue in the case of large periodic structures. It couples finite elements and the periodic structure theory (also called Bloch'stheory), that has been successfully used in the past for the analysis of free wave propagation in mono-and bi-dimensional waveguides [8]. The WFE wasalso used to compute vibroacoustic indicators such as group velocities and modal densities for use in Statistical Energy Analysis (SEA) [9,10]. Extensions of this method towards the calculation of forced responses of waveguides have also been proposed [11], without considering fluidstructure coupling. An application of 1D WFE formulation wasused by Serra et al. [12] for finite multilayers comprising poroelastic materials, butthe transfer approach necessitated meshing the whole surface of the plate. An assessment of wave methods can be found in [13], comparing 3 methods, namely TMM, the WFE+SEA procedure of [9] and aRayleigh-Ritz procedure.
Since the periodic structure theory deals primarily with infinite structures, the WFE can be used to compute the transmission loss of infinite plane structures excited by a superposition of plane wavesd irectly,w ithout having to resort to SEA. The derivation of this application is proposed in this paper,m aking use of the WFE framework for aplate with homogeneous properties in its plane. This enables to model an arbitrary layering with the finite element model of as ingle unit cell, which is just 1e lement wide in each of the plane'sdimensions, as in [14] and [15]. This paper is structured as follows: the WFE model under forced plane wave is presented in Section 2, with the derivation of the formulation of the TL and absorption coefficient. The method is validated on several cases in Section 3: an isotropic thin plate, an orthotropic plate, asandwich construction with as oft viscoelastic core and as tiff sandwich plate with honeycomb core. Adiscussion on the element size closes this section. Finally,Section 4presents an application case in which the proposed method is used to study the effect of the order of the layers in acomposite plate.

Forced response to an incident plane wave
We consider amultilayered plate with one surface lying in the plane z = 0. Each layer is considered homogeneous, while this is of course not the case for the assembly.A s ah omogeneous structure is periodic for anyp eriod, the unit cell is ap arallelepiped. It is modelled with standard (8 nodes)b rick elements. The FE discretisation has only asingle element in each of the x and y directions, and as manyasneeded in the z direction.

Equilibrium formulation
Let D denote the dynamic stiffness matrix of the unit cell. This matrix is computed from the finite element stiffness and mass matrices K, M with the relationship Viscoelastic damping can be taken into account through a complexpart of the stiffness matrix K. The equilibrium equation of the unit cell then writes where f represents the forces due to the interaction with the neighbouring cells, and e the forces imposed on the plate by the acoustic pressure fields on each side. As can be seen in Figure 1, the displacement vector u can be partitioned as We nowsearch for solutions as plane waveswith the same wave characteristics as the projection of the incident wave on the plate, hence where λ x = exp(ik x d x )and λ y = exp(ik y d y ). The reduced variable u 1 can be used with the transformation where I is the identity matrix. The imposed pressure is condensed to forces exerted on the nodes, and can be partitioned in the same way, The internal forces exerted on each vertical segment are derivedf rom the stress field applied to the normal to the twofaces related to that segment. The linearity of the stress tensor with respect to displacement field u and the previous periodicity relations lead to 2) can then be premultiplied byT and rewritten as Let'sd enote A = 1 4T DT.W efi nally obtain the reduced equation, linking the forces and the displacements on one single segment instead of the initial four.

Acoustic coupling
We are concerned with an oblique incident plane wave impinging the plate on one side, with aknown amplitude p I . The interaction between this incident wave and the plate create ar eflected wave on the same side, and at ransmitted wave on the other one, whose respective amplitudes p R and p T we want to calculate. The sound field on the incident side is then the superposition of the incident and reflected waves, where the common factor exp(i(ωt − k x x − k y y)has been omitted for legibility.T he wavenumber components satisfy the relationship k 2 x + k 2 y + (k z,I ) 2 = k 2 = (ω/c I ) 2 and c I is the speed of sound in the fluid on the incident side. In the same way, the sound field on the transmission side contains only the transmitted wave propagating in the same direction as the incident wave, with the relationship k 2 where c T is the speed of sound in the fluid on the transmission side. Potential phase difference between the pressure fields are accounted for through the fact that the amplitudes p R and p T may be complex. The wavenumbers k x and k y are conserved across the plate, so only the k z component may vary with the nature of the fluid. We will consider in the following that the fluid is the same on both sides, thus ρ I = rho T = ρ 0 and c I = c T = c 0 .
In this case, the incidence angle θ is the angle between the wave vector and the normal of the plate, while the azimuthal angle φ givesi ts orientation in the plane. We therefore have k x = k sin θ cos φ, k y = k sin θ sin φ and k z = k cos θ.
The load imposed on the plate can be written from these twop ressure fields lumped on the nodes of the finite element model. As the pressure force is exerted along the normal to the structure, the only non-zero terms in the vector e 1 of forces imposed on segment 1will be those relative to the DOFs in z-direction on the incident side, e I ,and on the transmission side e T .For an elastic material, these quantities are scalar.Let u O and e O denote the displacement and force vectors on all other degrees of freedom in the segment. The force exerted on the nodes of segment 1isthen, as shown in Figure 1 where S is the free surface of the element, which is identical on both sides. All other vectors and matrices can be written following the same decomposition, allowing to rewrite (8) as The second line of (12) leads to which allows to condense this equation into The fluid-structure interaction is characterized by the continuity of normal particle velocity at the interface. This writes for the incident side, and for the transmission side: Introducing acoustic admittance Y 0 = cos(θ)/(iωρ 0 c 0 ), we obtain These expressions for u I and u T can be re-injected into (14),leading to twoscalar equations linking the unknowns p R and p T and the incident pressure p I , Solving this equation givesthe acoustic transparency τ = |p T /p I | 2 and the absorption coefficient α = 1 −|p R /p I | 2 . The transmission loss for aplane wave is then finally TL = −10 log 10 τ. (18) This leads also to the diffuse field transmission loss by integrating overa ll possible incidences between minimum and maximum bounds θ ∈ [θ min ,θ max ]a nd directions φ ∈ [0, 2π]. It is recommended to avoid the full range [0, π 2 ]for θ because numerical errors occur close to grazing incidences θ = π/2. However, θ min = 0c auses no issue. The diffuse field transparencyisthen and the diffuse field transmission loss

Validation results
In order to assess the validity of the method, acomparison with various analytical models is performed in this section. Adiscussion on the size of the elements to be used is proposed in Section 3.5.

Isotropic plate
The first case is an isotropic aluminium plate, under two different excitations: a45 o plane wave and adiffuse field. The material parameters are the Yo ung modulus E = 70 GPa, the Poisson ratio ν = 0.33, the density ρ = 2700 kg.m −3 .T he loss factor is η = 0.5%, which leads to ac omplexY oung'sm odulusẼ = E(1 + iη). The unit cell consists of 10 elements through the thickness, with an overall thickness of 2mmand lengths d x = d y = 0.1mm.
Twok inds of analytical models are used for comparison. The first one is based on the Kirchhoff-Love plate assumption, while the other is the Transfer Matrix Method (TMM)d escribed in [7]. All three model agree very well with each other,w ith ad i ff erence of less than 1dBo ver most of the frequencyrange, except in atinyregion around the coincidence frequency, which is slightly underestimated in the thin plate model. Because the TMM and the current model agree in this region, this difference must be due to the fact that the Kirchhoff-Love assumptions neglect shear,while the other twouse a3Delasticity theory.
The diffuse field integration has been performed over the range [0, 0.999π/2] to avoid numerical issues with grazing incidence. The range wasd iscretised with 1000 points. The three models are found in this case to agree up to 1dBo vert he whole frequencyr ange, the small difference at coincidence disappearing because of the averaging over1000 incidence angles.

Orthotropic plate
The diffuse field transmission loss of an orthotropic plate is compared against an analytical model [16]. The material parameters are giveni nT able I, where E stands for Yo ung'sm odulus, ν for Poisson'sr atios and G for shear modulus. The TL results are presented in Figure 4for frequencies ranging between 100 Hz and 20 kHz. The mesh of the unit cell has 10 elements in the z direction and lengths d x = d y = 0.1mm. Here again, the analytical and current models agree well with each other,w ith an error     smaller than 1dBoverthe whole range. There are twocritical frequencies due to the very different Yo ung'sm oduli in the x and y directions. An expression for these frequencies can classically be derivedfrom the analytical model, and is givenbythe formula where ν yx = (E y /E x )ν xy .T he other critical frequency f crit,y is obtained by replacing E x by E y in (21).I nthis case, the frequencies obtained by these formulas agree well with the dips observed in the TL curve.

Stiff sandwich plate
The 6 th order constitutive equation for the bending of a sandwich beam wasgiven by Mead and Markus [17] and later extended to aplate by Narayanan and Shanbhag [3].
The main assumptions are that the predominant deformation of the core is due to shear,while that of the skins is due to pure bending. This implies that the Yo ung modulus of the core is high compared to the out-of-plane shear modulus. The model has fivem ain parameters, namely skin bending stiffness D t ,overall bending stiffness B,damping η,surface mass m and shear parameter g.The constitutive equation reads where w = v/iω is the normal displacement of the plate.
In the considered frame where af orced wave is imposed on the plate with aw avenumber k = ω c sin θ,t he spatial derivative operator ∇ can be replaced by −ik.T his leads to the following expression of the impedance: Fors andwiches made of isotropic materials and identical skins, the skin bending stiffness is D t = E s h 3 s /[6(1 − ν 2 s )] the overall bending stiffness is B = E s h 2 c h s (1 + h s /h c ) 2 /2 and the shear parameter is g = 2G c /(E s h s h c ), where G c = G c,xz = G c,yz is the core transverse shear modulus, E s and ν s = ν s,xy are the in-plane Yo ung'sm odulus and Poisson ratio of the skins, and h s (resp. h c )isthe thickness of one skin (resp. the core). In-plane core shear modulus and all Yo ung'smoduli of the core are neglected in the analytical model.
The twom odels are compared between 1kHz and 20 kHz on as andwich plate with isotropic skins and orthotropic core impinged by a45 o plane wave.The parameters of the sandwich construction are giveni nT able II, and the TL results shown in Figure 6. With amesh using 3e lements for each skin and 45 for the core, and again d x = d y = 0.1mm, the twom odels differ also by less than 1dBi nt he whole range. Since the analytical model takes shear into account, the agreement is also excellent at coincidence.

Sandwich plate with soft core
The same calculation wasm ade with at hick sandwich plate made of ar ather soft isotropic core and thick skins. The material properties are described in Table III. The transmission loss curvei ss hown on Figure 6. Again this correlates well with the analytical TMM, with amaximum error of 1dBoverthe considered range. Three dips in TL     can be observed between 1kHz and 6kHz, corresponding to coincidence frequencies. The first twodrops around 1730 Hz and 5470 Hz are due to coincidences with awave having as ymmetric motion of the skin and compression  in the core, while the third one around 5620 Hz is due to coincidence with the bending wave.The dispersion curves for these twowaves are presented in Figure 7.
The displacement field can be calculated back from (13).T heya re displayed at several frequencies on Figure 8. The first one is related to the symmetric motion of the skins (Figure 8a), while the second is related to a global bending motion of the sandwich, or antisymmetric motion of the skins (Figure 8b). In both cases, it can be observed that ashearing motion of the skins is involved, due to their rather high thickness. In the antisymmetric case, the cross-section remain straight and vertical, indicating that only shear is at stake, with practically no bending motion.
The field represented on Figure 8c is at afrequencybetween the twocoincidences, showing adecoupling of the twoskins. On the latter example, it may be noted that the amplitude of the z-displacement is higher on the transmission side than on the incident side, yet the amplitude of the transmitted wave is of course smaller than that of the incident wave.This is due to the fact that the pressure field on the incident side is result of interference between incident and reflected waves, while the transmitted field consists of asingle plane wave.

Discussion on element size
In the validation cases presented in this section, the mesh parameters are chosen so as to ensure convergence with the TMM model, which is based on 3D isotropic elastic theory and therefore exact in this case. Ac oarser mesh leads to errors in the coincidence frequencyregion. Of all the wavesconsidered in acoustic transmission, the acoustic wave in the fluid has the shortest wavelength, at the highest frequency( 20 kHz)c onsidered in all these cases. Since it is imposed on the structure at grazing incidence, it is relevant to compare the mesh size to it. The element sizes in x and y directions must therefore satisfy about 100 elements per wavelength, which is much stricter than the usual 10 elements per wavelength used for finite elements. Since there is only one element in the in-plane directions, this has no impact on the number of degrees of freedom in the mesh. However, the results are also found to be quite sensitive to the aspect ratio of the elements, which means that the element size in the z-direction must be about the same as that in the plane. Forthick layers, asmall d x may thus lead to ahigh number of elements through the thickness. On the other hand, it also sets amaximum size for d x and d y ,since there must be at least one element per layer. Above 10 kHz, this latter limit is usually less restrictive than the "100 element per wavelength" criterion for plates with thickness in the order of 1cm or less.

Influence of stacking on the transmission loss
The current approach can be used to study the influence of the stacking sequence of amultilayer plate on the TL. We consider a1.6 mm thick composite plate made of 24 layers of the same transverse isotropic material described in Table IV in order to study the effect of the stacking sequence on the TL. Three different stacking sequences have been studied, which are described in Table V. All of these have the same proportion of layers with fibres oriented towards each direction, at -45, 0, 45 or 90 degrees of the global coordinate system, orientation 0being aligned with the x axis. The layers are supposed to be perfectly bonded to each other,sothe FE model ensures continuity of displacements at the interfaces. The three configuration are aclassical symmetric quasi-isotropic stacking ("initial" configuration in Table V),arandom permutation of it ("random" case), and agrouping of the same orientation in ascending order ("sorted" case). The TL wasc omputed for these three sequences with the method described above,u sing 2e lements per layer and amesh size of 0.1 mm in the x and y directions. The results are presented in Figure 9.
Concerning the lowfrequencybehaviour,itcan be seen that the stacking sequence has practically no influence for frequencies below2kHz, and little up to 4kHz. This is due to the fact that belowthis frequency, the plate behaves according to the mass law, and the actual stiffness has no effect on TL.
In high frequencies (above 8kHz), the random and initial stacking sequences have the same asymptotic behaviour.T his is due to the dominant effect of damping on the high-frequencyT L, above the highest critical frequency. These effects have been studied in [18]. The sorted case exhibits the same asymptotic behaviour,which is only visible for much higher frequencies, above 20 kHz and thus not represented on Figure 9 Between these, all configurations present acoincidence region, characterized by alower TL. Because of the presence of different layer orientations, the coincidence zone of all plates is strictly contained between the twoextreme critical frequencies of ap late where all layers are identically oriented. These bounds are givenb y( 21), which yields for this material f crit,x = 3770 Hz and f crit,y = 19 kHz. The random and initial configurations exhibit only one dip, indicating that theybehave closely to an isotropic plate in this respect. This is due to agood repartition of the layer orientations throughout the thickness. On the other hand, the sorted stacking exhibits aw ide coincidence re- gion with twocritical frequencies, abehaviour that is similar to that of an orthotropic plate.
The coincidence phenomenon is illustrated in Figure 10, which shows the free propagating wavenumbers as afunction of direction φ at afrequencyof7kHz, which is in the coincidence region for all 3cases. These results were obtained with aWFE calculation using the same mesh. The acoustic trace wavenumber k a = 2πf/c 0 is represented for ag razing incidence, and does not depend on φ.T he random stacking k-space has ashape close to acircle, and lies mostly below k a .The initial stacking has amore distorted k-space, partly above and partly belowt he acoustic wavenumber,b ut the deviation to the circle remains small, hence the narrowc oincidence region. Finally,t he sorted case has av ery stretched k-space, close to an ellipse with major axis oriented around −π/6radians, which also crosses the k a circle. This behaviour is related to the strongly inhomogeneous nature of this stacking case: for aw avew ith direction φ = 0, the stiffest layer (oriented along the fibers)i so nt he outer part of the plate, while a wave with φ = π/4" sees" the stiffest layer close to the central plane.

Conclusion
An ew method for computing the transmission loss of infinite plates made of arbitrarily layered materials has been presented in this paper.T he originality of the approach is that it couples the periodic structure theory to finite elements as in the WFE method, while taking into account the coupling with an acoustic fluid analytically. The method is validated against classical analytic Transfer Matrix Method computation for isotropic, orthotropic and sandwich constructions. The main interest of the method is that it natively takes into account complexbehaviour of the structure, such as symmetric motion of the skins with respect to the neutral plane. Besides, the use of finite elements makes it easy to retrieve the displacement field inside the structure. The computation time is much higher than that of analytical methods, butremains low, typically up to about 10ms per frequencyfor aplane wave.Apossible use of this method is the validation of analytical model for multilayers.
An application case is proposed to study the effect of the stacking sequence of alayered composite plate on the diffuse field TL. It is shown that the stacking sequence has an effect on the width of the coincidence range, and little influence outside it, the lowfrequencies being governed by mass lawand the higher by damping effects. The most extreme case of asorted stacking shows awide coincidence zone limited by twocritical frequencies, while the random permutation exhibits ad ip in an arrower zone around the lowest critical frequency( 7kHz in this case). In applicat ions where noise control is important at frequencies corresponding to the coincidence range, stacking sequence should therefore be considered for acoustic criteria too.