Nicolas Auquier
email: nicolas.auquier@entpe.fr

Kerem Ege

Emmanuel Gourdon

Equivalent dynamic model of multilayered structures with imperfect interfaces: application to a sandwich structured plate with sliding interfaces

Keywords: equivalent single layer, equivalent model, imperfect interface, sliding interface

This works aims to build an equivalent model of multilayered structures with imperfect interfaces for vibroacoustic modelling and characterization. To take into account imperfections, new interface conditions including constitutive equations, that describe the imperfections, are implemented. Once the displacement field with imperfect interfaces is obtained, the dispersion relation of the structure is derived from the equivalent model. The bending wavenumbers are then used to compute the equivalent flexural rigidity and the damping of the sandwich panel. In this paper, the methodology is applied to a sandwich panels with sliding interfaces. The equivalent parameters are computed and compared to the reference case, i.e. perfect interfaces model. The main impact of the imperfection implementation observed in the results is a shift towards the low frequency of the equivalent parameter curves.

Introduction

This work aims to develop an equivalent model able to characterize and optimize the behavior of multilayered structures with imperfect interfaces. The structures of interest here are sandwich panels, made up of three layers with a soft layer as core, usually thick, and hard skins from either side of the core. Sandwich panels are mainly used in the industry because of its high ratio stiffness over mass and its effectiveness to dissipate vibrations by using the high damping of the core. Hence, such structures yield many applications in several fields such as transport [START_REF] Rao | Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes[END_REF][START_REF] Subramanian | Optimization of Damping Treatments for Structure Borne Noise Reduction[END_REF] or civil engineering [START_REF] Caniato | Acoustic of lightweight timber buildings: A review[END_REF][START_REF] Jamal Arbaoui | Effect of Core Thickness and Intermediate Layers on Mechanical Properties of Polypropylene Honeycomb Multi-Layer Sandwich Structures[END_REF], mostly for customer comfort and energy savings in the audible spectrum.

This particular interest has found echoes in the scientific community, which has developed many models to predict the response of such materials. Equivalent plate models are increasingly developed in the scientific community these last years [START_REF] Carrera | An assessment of mixed and classical theories on global and local response of multilayered orthotropic plates[END_REF][START_REF] Carrera | Theories and Finite Elements for Multilayered, Anisotropic, Composite Plates and Shells[END_REF][START_REF] Ghinet | Modeling thick composite laminate and sandwich structures with linear viscoelastic damping[END_REF][START_REF] Dazel | A stable method to model the acoustic response of multilayered structures[END_REF][START_REF] Viverge | Model of highly contrasted plates versus experiments on laminated glass[END_REF], because of their several advantages for the calculation and the understanding of the modeling. For instance: (1) a minimum calculation time allows the testing of many configurations efficiently i.e. without using a finite element model mainly 3D with a significant mesh.

(2) A small amount of parameters which make the method easily implemented. (3) Physical meaning brought by the physical model. Carrera [START_REF] Carrera | An assessment of mixed and classical theories on global and local response of multilayered orthotropic plates[END_REF] has gathered the models into three categories depending on their implementation types. The first category is the Equivalent Single Layer (ESL) models, they describe the displacement of the whole structures through the displacement variables of a reference layer. The second category is the Layer-Wise (LW) models, they describe the behavior of the structure in each layer, which increase the calculation time but increase also the accuracy. For instance, LW description is more accurate than a corresponding ESL model, i.e. same hypothesis and order for z-expansion. The third category is for hybrid models, hereby providing the advantages of both previous categories. Once the model is chosen, the goal is similar and consist in developing a calculation methodology from the displacement field of the panel, to get the equivalent parameters of the multilayered system or its strains.

Guyader and Lesueur [START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF][START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part II: Transmission loss[END_REF] developed a hybrid model for n-layers structures based on the work of Sun [START_REF] Sun | Theories for dynamic response of laminated plates[END_REF]. The methodology is developed by assuming a displacement field composed of 3 types of displacement: transversal, shearing and bending displacement, which is related to the Reissner model [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF]. The latter is an extension of the Love-Kirchhoff plate theory [START_REF] Love | The small free vibrations and deformation of a thin elastic shell[END_REF] that takes into account shear strains throughout the thickness of the plate. Firstly, each layer has a specific displacement field as a LW model, but then interface conditions between the layers are implemented in order to propagate the kinematic field of the layer n into the layers. The propagation is applied until the kinematic field of the layer n is fully described with respect to the kinematic variables of the first layer, chosen as the reference layer in most cases. Once the kinematic field is fully propagated, the model is equivalent to an ESL model because the kinematic of all the layers are described thanks to the reference layer only. Thus, this kinematic field can be used to derive equivalent parameters of the multilayered structure, such as the equivalent flexural rigidity or the damping. The work of Guyader has been recently improved by Marchetti [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF][START_REF] Marchetti | On the condensation of thick symmetric multilayer panels including dilatational motion[END_REF], he extended the previous work to panel structures with anisotropic layers and it gathers many details about the implementation of the model.

However, the interface perfection between two consecutive layers in multilayered media is not guaranteed, which yields to uncertainties in the characterization of panels. Interface characterization for multilayered structures has found some interest in the scientific community. Schoenberg [START_REF] Schoenberg | Elastic wave behavior across linear slip interfaces[END_REF] showed the impact of a slipping interface over the propagation of an elastic wave through the media, especially at the interface between two layers, where the plane wave reaction (described by the reflection coefficient for instance), will depend on the interface properties. In this work, the effects of the sliding between two consecutive layers are modeled thanks to a zero-thickness interface, but interface can also refers to actual thickness inside a layer or between two layers where inhomogeneities are present as in [START_REF] Margetan | The Interaction of Ultrasound with Imperfect Interfaces: Experimental Studies of Model Structures[END_REF]. In the latter, the goal of using interface modeling is to simplify the model by considering a small thickness with specific parameters that will take into account the impact of the inhomogeneities. Larentyev [START_REF] Lavrentyev | Ultrasonic spectroscopy of imperfect contact interfaces between a layer and two solids[END_REF] experimentally highlighted pure interface effect due to interface quality. He analyzed the interaction of ultrasonic waves sent through two roughened aluminum layers in contact. Since the layers have the same properties, the interface properties are only function of the coupling quality. Throughout the previous works, it is shown that an interface stiffness is a good tool to describe interface properties, low interface stiffness is related to a low bonding quality and a very high interface stiffness is similar to a perfectly bonded system.

The goal of this work is to extend the latest improvement of the Guyader's model for multilayered plates, by taking into account imperfect interfaces, and more specifically, as a first step, sliding between the layers. The modeling of sliding interfaces is performed as in the work of Larentyev [START_REF] Lavrentyev | Ultrasonic spectroscopy of imperfect contact interfaces between a layer and two solids[END_REF], which consist of assuming interface stress resulting of sliding and interface stiffness. Thus, an equivalent model is derived with new assumptions for the interface conditions between the layers, and several material parameters are calculated to check the impact of such modifications. So far, the effects of imperfect interfaces in equivalent plate models observed in simulations are mainly about internal stresses of the multilayered plate [START_REF] Liu | Stress analysis of imperfection composite laminates with an interlaminar bonding theory[END_REF][START_REF] Cheng | Theory for multilayered anisotropic plates with weakened interfaces[END_REF][START_REF] Massàbo | Assessment and correction of theories for multilayered plates with imperfect interfaces[END_REF][START_REF] Massàbo | An efficient approach for multilayered beams and wide plates with imperfect interfaces and delaminations[END_REF][START_REF] Lu | Interlayer shear slip theory for cross-ply laminates with nonrigid interfaces[END_REF]. The interest is mainly in the evaluation of the maximum stress for each layer of a panel under a static load. It is shown that an equivalent model with perfect interfaces will underestimate the stresses applied to the skins, which could be a problem if the skins are too thin for example.

One original feature of this work is to implement imperfect interfaces into a dynamic model and to highlight its effect on frequency-dependent equivalent parameters, such as the equivalent flexural rigidity. Moreover, the application chosen, instead of being a two-layered or three-layered plate with the same thickness and the same material properties as the previous cited models, is a sandwich panels, which implies three layers with two different thickness and two different materials. The paper is organized as follows: 2 describes the model development in details. The 3 describes the methodology used to derive dynamic parameters in order to show the imperfection impacts and to compare them with the perfect interfaces case. The 4 shows the results of the methodology application to a sandwich panel; in 5 is discussed the model assumptions and the results.

Models

The aim of this section is to detail the models used in the current work. It starts with the Guyader-Marchetti's model and its assumptions, in order to see clearly the modification carried out in the proposed model and also to have a better understanding of its starting assumptions. The modification of the model is done to take into account imperfections. In this work, the imperfections assumed in the model are jumps of the transverse displacement, which implies uncontinuous displacement field between the layers. This modification is implemented similarly to other imperfect interface models [START_REF] Liu | Stress analysis of imperfection composite laminates with an interlaminar bonding theory[END_REF][START_REF] Cheng | Theory for multilayered anisotropic plates with weakened interfaces[END_REF][START_REF] Massàbo | Assessment and correction of theories for multilayered plates with imperfect interfaces[END_REF][START_REF] Massàbo | An efficient approach for multilayered beams and wide plates with imperfect interfaces and delaminations[END_REF][START_REF] Lu | Interlayer shear slip theory for cross-ply laminates with nonrigid interfaces[END_REF] but in dynamics, and is summarized by the new displacement field. Throughout this paper, the model presented is developed in the isotropic case for the sake of simplicity. But the model with imperfect interfaces can be developed for an anisotropic case as in Marchetti's work for the perfect interfaces case [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF][START_REF] Marchetti | On the condensation of thick symmetric multilayer panels including dilatational motion[END_REF].

𝑥 𝑦 𝑧 h 𝐿 𝑥 𝐿 𝑦 𝑛 = 1 𝑛 = 2 𝑛 = 3 𝑅 𝑛=2 𝑅 𝑛=3
Figure 1: Sandwich panels with perfect interfaces.

Guyader-Marchetti's model

The beginning of our work settled in the equivalent model for multilayered plate of Guyader [START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF] and especially in the recent developments of Marchetti [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF] who extended the model to anisotropic plates. The Guyader's model is categorized as a "Zig-Zag" model, the main characteristic of such models is that they take into account a kinematic analysis for each layer and then use perfect interface conditions between the layers as in 1 to obtain an Equivalent Single Layer (ESL) modeling.

The Guyader's model starts with the displacement field introduced by Sun and Whitney [START_REF] Sun | Theories for dynamic response of laminated plates[END_REF]:

     u n x = u n 0x (x, y, t) -z n ϕ n x (x, y, t), u n y = u n 0y (x, y, t) -z n ϕ n y (x, y, t), u n z = W (x, y, t), (1)
with, u n i the displacement of the layer n for the coordinate i, z n the thickness coordinate with the origin located at the midplane of the nth layer and ϕ n i the rotation in the (Oxz) or (Oyz) plane.

However, the displacement field is written differently in [START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF] by adding a new kinematic variable to the model. The Guyader's displacement field is assumed to be described thanks to 3 displacement terms, namely: the bending, the shearing and the membrane displacement, respectively represented by W (x, y, t), ϕ(x, y, t) and ψ(x, y, t), which corresponds to a Reissner's [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] or Mindlin's [START_REF] Mindlin | Influence of rotary inertia and shear on flexural motions of isotropic elastic plates American Society of Mechanical Engineers[END_REF] plate.

     u n x = ψ n x (x, y, t) -(z -R n )( ∂W (x,y,t) ∂x + ϕ n x (x, y, t)), u n y = ψ n y (x, y, t) -(z -R n )( ∂W (x,y,t) ∂x + ϕ n y (x, y, t)), u n z = W (x, y, t), (2)
with, ψ x (ψ y ) the translation over the Ox n (Oy n axis), which simply describes the displacement yielded by the extensional waves in the layer n, R n is the z mid-plane coordinate of the layer n with respect to the global axis. ϕ n x (ϕ n y ) is the rotation around the Oy (Ox) axis, which is depending on shear of the structure. Firstly, all the layers have the same theoretical displacement field formulation but their index. Thus, to compute the displacement field of the whole plate, one has to propagate the displacement field into the layers thanks to interface conditions between the layers. For the Guyader model, it is done throughout continuity conditions of the kinematic field and the stresses. In the Guyader-Marchetti's case, the interface conditions are continuity conditions for the transverse shear stresses and the displacement field.

σ n αz = σ n-1 αz , U n = U n-1 , (3)
with, α = {x, y} and

U n = [u n x , u n y , u n z ]
. The interface conditions leads to a linear system of equations that can be rewritten by means of a transfer matrix of the interface n [T n ]. This linear system of equations expresses the kinematic field of the layer n (U n ) with respect to the kinematic field of the layer n -1 (U n-1 ), which yields

(U n ) = [T n ].(U n-1
). After all, [T n ] can be seen as an operator, its application to the kinematic field of the layer (U n ) yields the kinematic field of the next layer (U n+1 ). The operator has to be applied as much as there are layers, so one can derive the expression of the kinematic field of the layer n with respect the layer n = 1 and by doing so, derive an ESL model because the displacement field of the whole panel is described thanks to a single layer. If the reader is interested in the details of the derivation he can refers to the work of Marchetti [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF].

Once the displacement field is fully propagated, the displacement variables of the first layer are the only unknowns of the model:

     u n x = ψ 1 x (x, y, t) + F ω ∂W ∂x + F n xx ϕ 1 x + F n xy ϕ 1 y , u n y = ψ 1 y (x, y, t) + F ω ∂W ∂y + F n yx ϕ 1 x + F n yy ϕ 1 y , u n z = W (x, y, t), (4)
with,

F ω = R 1 -z, F n ij = α n ij (R n -2) + γ n ij .
α n ij and γ n ij are material parameters (mainly Young's modulus ratios) yielded by the transfer matrix [T n ]. The details are given in appendix of Loredo [START_REF] Loredo | A multilayer anisotropic plate model with warping functions for the study of vibrations reformulated from Woodcock's work[END_REF] or by Marchetti [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF].

Thus, the Eq. ( 4) are the same as Eq. ( 2) but propagated thanks to the interface conditions. Fundamentally, the next step is to obtain the equation of motions and then compute equivalent parameters. However, the next subsection will describe the other model used in this paper. It introduces the interface conditions that take into account bonding imperfections, which has an impact on the propagated displacement field seen in Eq. (2).

Imperfect interface implementation

Previously, the interface conditions are assumed to respect perfect continuity of the displacement field and the transverse and vertical stresses, describing hereby perfect interfaces. For example, Massàbo's model [START_REF] Massàbo | Assessment and correction of theories for multilayered plates with imperfect interfaces[END_REF] also starts with a Zig-Zag model, it takes into account transverse shear as previously but the mathematical formulation is slightly different and used only in a static application. A displacement field for each layer is also stated, but displacement variables come from mathematical assumptions sorted by smoothness types, defined by the number of continuous derivatives it has over a domain. This yields the global displacement to be a composition of different function classes: a global displacement with a continuous first derivative function over z also known as a C 1 function; and perturbations with a piecewise continuous function over z namely a C 0 function. This model is equivalent to the previous model because the displacement field takes into account the same phenomenon such as shearing, bending and membrane effects.

The imperfect interface is implemented by assuming discontinuous displacement field at the interfaces 2. Then, a constitutive equation is formulated for the interfaces in order to compute the sliding term in the new displacement field. A spring approach is used to represent imperfect interface effects, the springs take into account the coupling quality between two consecutive layers. If the equivalent stiffness approaches infinity the interface has a perfect interface behavior. However, if the equivalent stiffness approaches 0 the consecutive layers are totally debonded. The spring stiffnesses are piecewise linear functions, which makes such formulation quite adaptive if needed. On one hand, these functions are called "jumps" and are written ûx n or ûy n for interfacial slip in the x or the y direction. On the other hand, they are called "openings" and are written ûz n for a discontinuity in the z direction. These discontinuities can be computed by the difference between the displacement field of two consecutive layers:

Û i n = U i n -U i n-1 . (5)
It is noteworthy to specify that the interface indexation depends on the indexation of the layers. In this work, the interface index is chosen according to the upper layer, which implies that for n = 1 there is not an interface since it indicates the first layer. Mathematically, it is represented by writing simply Û i 1 = 0. Moreover, openings are not considered in the model, they are assumed to be negligible ûz n ≈ 0. Thus, this new term in the displacement field is implemented into the Guyader-Marchetti's model and is represented in 2. To do so, the interface conditions are changed to discontinuous interface conditions for the displacements but still continuous for the stresses since contact is still occurring. The new interface conditions taken into account are modeled as:

σ n αz = σ n-1 αz , U n = U n-1 + Û n i . (6)
Once the interface conditions are modified, the derivation of the equations is done similarly as Guyader-Marchetti's model. Fortunately, the propagated displacement field can be written almost the same way as in Eq. ( 4) by gathering the terms with respect to the kinematic variables:

     u n x = ψ 1 x (x, y, t) + F ω ∂W ∂x + F n xx ϕ 1 x + F n xy ϕ 1 y + n i=1 ûi x , u n y = ψ 1 y (x, y, t) + F ω ∂W ∂y + F n yx ϕ 1 x + F n yy ϕ 1 y + n i=1 ûi y , u n z = W (x, y, t). (7)
However, the so-called jumps are not directly quantified and need another equation in order to be physically valuable. Thus, an equation has to be used in order to write the jumps with respect to known parameters. So the following constitutive equations are used:

σn z=z n = K n . Û n , (8)
And in its inverse form:

Û n z=z n = B n . σn , (9)
with, "." the scalar product, σn = [σ n xz , σn yz ] the vector of the interface stresses between the layers n and n -1,

K n = [K n xx K n xy , K n yx K n yy ] the matrix of interface stiffness, B n = [B n xx B n xy , B n yx B n yy ]
the matrix of interface compliance. The stiffness matrix K n or the compliance matrix B n are complex matrices in order to take into account interface dissipation. Nevertheless, the dissipation can be implemented thanks to a dissipative term such as σn xz → 0, so ûn x → 0, which corresponds to a perfectly bonded condition such as the Guyader's case where no imperfections are assumed.

z=z n = K n . Û n +t n with t n = [t x ,
Moreover, since the interface stiffness (or compliance) is written as a piecewise linear function, several kind of interfaces can be described. Such functions can, by their formulation, illustrate many kinds of interface behaviors, for instance a classic linear function of û, a non-linear behavior or even non-elastic behavior. Fig. 3 shows two different cases: an elastic piecewise linear case; and a non-elastic case. Thus, the jumps ( Û n ) can be written thanks to Eq. ( 9), by means of the material parameters and the stresses. To link the kinematic variables to the interface variables, the Hooke's law has to be used. Then, the displacement field can be rewritten by using only the kinematic variables and some material parameters.

In the case of ûn

x for the most general case:

ûn x = B n xz (Φ 1 x (-α n yx Q n 45 -α n xx Q n 55 ) + Φ 1 y (-α n yy Q n 45 -α n xy Q n 55 ) -t n x ), (10)
with, Q ij stiffness values following Voigt's notation. After writing an equivalent equation for ûn y , the propagated displacement field Eq. ( 7) can be written as:

     u n x = ψ 1 x (x, y, t) + F ω ∂W ∂x + F n xxg ϕ 1 x + F n xyg ϕ 1 y - n i=1 B i xz t i x , u n y = ψ 1 y (x, y, t) + F ω ∂W ∂y + F n yxg ϕ 1 x + F n yyg ϕ 1 y - n i=1 B i yz t i y , u n z = W (x, y, t). (11)
Where,

F xx became F xxg = F xx -B n xz (α n yx Q n 45 + α xx Q n 55 ), F yy became F yyg = F yy -B n yz (α n yy Q n 44 + α n xy Q n 45 ), F xy became F xyg = F xy -B n xz (α n yy Q n 45 + α n xy Q n 55 ) and F yx became F yxg = F yx -B n yz (α n yx Q n 44 + α n xx Q n 45 ).
It is noteworthy to say that in this case, no openings are taken into account, which leaves u n z unmodified. If the layers of the multilayered plate are assumed to be isotropic, it yields F xy = F yx = 0 and

F xx = F yy .
Once the displacement field is obtained, most of the previous applications [START_REF] Liu | Stress analysis of imperfection composite laminates with an interlaminar bonding theory[END_REF][START_REF] Massàbo | Assessment and correction of theories for multilayered plates with imperfect interfaces[END_REF] show the impact of the imperfections on the displacement field or on the stresses under a static load. However, the presented model is a dynamic model for infinite plate, but the displacement can still be visualized as shown in 3.1 and particularly thanks to 9, where the impact of the imperfection can clearly be checked after computation: a discontinuity of the displacement at the interfaces. In this work, the impact on equivalent dynamic parameters is shown, such as the flexural rigidity and the equivalent damping. A new displacement field Eq. ( 7) is used, but the computation methodology follows the same steps as Guyader or Woodcock's work [START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF][START_REF] Woodcock | Free vibration of advanced anisotropic multilayered composites with arbitrary boundary conditions[END_REF]. Similarly to the previous model, the coefficients can be propagated into the equations without repeating the derivation from scratch. But, this has to be done carefully by considering every step and the changes it yields on the equations. For instance, a spatial derivative of the displacement field could impact the new coefficients if they are functions of space.

The model has been fully described in this section. In the light of the details given throughout the model description some limitations can be drawn:

• Dilatation motion is not yet considered in the presented model. The first dilatation modes are hereby the high frequency limit of the model.

• The vertical displacement is assumed to be linear with respect to z. This yields the shear stresses to be constant through the thickness. Marchetti has shown that the effect of this assumption has low impact on the dynamics of multilayered plate [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF] for the audible spectrum.

• The interface parameter, namely the interface stiffness K ij or the interface compliance B ij is assumed to be constant throughout this work for the sake of simplicity. But, the parameter can be implemented frequency dependent or space dependent if needed.

• Openings ûn z are not currently considered here. They are assumed to be negligible.

Methodology

This section details the methodology applied in this work to visualize the impact of the imperfections on the equivalent parameters, this gives insight into the panel behavior and particularly into the imperfection impact. It is done by considering the multilayered plate as a single layer plate for each frequency and compute its equivalent parameters as in [START_REF] Ege | Assessment of the apparent bending stiffness and damping of multilayer plates; modelling and experiment[END_REF], this will yield frequency-dependent parameters for the bending motion. In order to give evidence about the correctness of the method, the dispersion curves are compared with the Spectral Finite Element Method (SFEM) used by Shorter [START_REF] Shorter | Wave propagation and damping in linear viscoelastic laminates[END_REF].

The methodology detailed in this section can be summarized thanks to the workflows in 4. Firstly, the displacement field u n is used to compute the kinetic energy and the strain energy in order to use the least action principle. This yields a differential equation system that can be written thanks to matrices. A particular solution, which describe a propagative wave in the (Ox) and (Oy) directions, is assumed and leads to a system of equation with respect to the wavenumber k. Finally, the determinant of the equation yields the wavenumbers with respect to the frequency. Then, the bending wavenumbers are retrieved and equivalent dynamic parameters are computed by using a Love-Kirchhoff model. The main idea is to compute the bending wavenumber taken into account by the displacement field assumed Eq. [START_REF] Dazel | A stable method to model the acoustic response of multilayered structures[END_REF]. Then, compute the frequency dependent flexural rigidity by doing a thin-plate equivalence thanks to the Love-Kirchhoff model. All these parameters will finally be used in an application in the next section to illustrate the impact of the imperfections on the equivalent parameters and so the behavior of an imperfect panel. All the examples shown in the methodology section are obtained from the parameters of a three-layers multilayered plate. Their values are gathered in 3.

Displacement field ! "

Differential!equations

# $ % $& % + ' . ( = 0

Least action principle

Particular solution ( = ( ) * , -/ 12&3 det ' 4 5 % # = 0 Wavenumbers 6 7

Equation!system

K45 % # = ) (a)
Bending Wavenumber ! Flexural rigidity

" ! = #$% & ! '

Love

Wavenumber calculation

Once the displacement field is obtained Eq. ( 11), the next step is to write the equation of motions. To do so, the Lagrangian of the system is computed L = T -V , which needs the computation of the kinetic energy T with Eq. ( 12) and of the potential energy V with Eq. ( 13):

T = + Lx 2 -Lx 2 + Ly 2 - Ly 2 1 2 n zn+ hn 2 zn-hn 2 ρ n ∂u n x ∂t 2 + ∂u n y ∂t 2 + ∂u n z ∂t 2 dz dydx. (12)
V = + Lx 2 -Lx 2 + Ly 2 - Ly 2 1 2 n zn+ hn 2 zn-hn 2 σ n .(ϵ n ) T dz dydx. ( 13
)
To compute these two energy terms, one needs only to use the kinematic field previously shown in Eq. ( 7).

The energy terms are developed with respect to the kinematic variables. Then, the principle of least action is applied for 3 variables x, y and t, which yields 5 equations, one for each kinematic variables [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF]. The equations to solve is:

∂L ∂α - ∂ ∂t ∂L ∂α ,t - ∂ ∂x ∂L ∂α ,x - ∂ ∂y ∂L ∂α ,y + ∂ 2 ∂x∂t ∂L ∂α ,xt + ∂ 2 ∂y∂t ∂L ∂α ,yt + ∂ 2 ∂x∂y ∂L ∂α ,xy + ∂ 2 ∂t 2 ∂L ∂α ,tt + ∂ 2 ∂x 2 ∂L ∂α ,xx + ∂ 2 ∂y 2 ∂L ∂α ,yy = 0. ( 14
)
with, α = {W, ϕ x , ϕ y , ψ x , ψ y }, α ,i = ∂α ∂i and α ,ij = ∂ 2 α ∂i∂j . The set of 5 equations thus obtained can be written thanks to a matrix equation Eq. ( 15):

M ∂ 2 ∂t 2 + K W = 0. (15)
A particular solution W = W 0 e jωt e -j(kxx+kyy) , with k x = k cos(θ), k y = k sin(θ) and with θ the direction of propagation, is then injected to remove the time differentials and write the system of equations as a mass and spring system as in Eq. ( 16):

K -ω 2 M W 0 = 0. ( 16
)
The matrices K and M are written respectively with the parameters λ i and δ i , also called Woodcock's parameters. These parameters have been expressed relatively to other plate models in [START_REF] Loredo | A multilayer anisotropic plate model with warping functions for the study of vibrations reformulated from Woodcock's work[END_REF]. They are functions of the material parameters such as thicknesses of the layers or Young's modulus of the layers. Since the λ i describe strain energy they are intuitively function of the Young's modulus of the layers and of their thicknesses. The δ i describe inertia effects and are functions of the density of each layer and also of Young's modulus ratio of the consecutive layers. These parameters can be expressed with respect to the coefficients of the kinematic variables in the displacement field Eq. ( 7) such as F ω or F xxg . Moreover, the matrices K and M are computed with the updated coefficients that take into account the imperfections. However, only some parameters are impacted by the imperfection implementation, these modifications are written as δ g i and λ g i and are implemented by adding them to the right coefficient such as δi = δ i + δ g i . All the parameters impacted for an isotropic case are written below: [START_REF] Shorter | Wave propagation and damping in linear viscoelastic laminates[END_REF].

λ g 2 = n i=1   Q i 11 h i   -2γ i i j=2 α j Q j 55 B j + i j=2 α j Q j 55 B j 2     λ g 4 = -2 n i=1   Q i 11 h i β i i j=2 α j Q j 55 B j   ν η
λ g 6 = -2 n i=1   Q i 11 h i i j=2 α j Q j 55 B j   δ g 2 = n i=1   ρ i h i   -2γ i i j=2 α j Q j 55 B j + i j=2 α j Q j 55 B j 2     δ g 4 = -2 n i=1   ρ i h i β i i j=2 α j Q j 55 B j   δ g 6 = -2 n i=1   ρ i h i i j=2 α j Q j 55 B j  
To have details about the classic Woodcock's coefficient in a similar case as the present work, the reader is referred to [START_REF] Loredo | A multilayer anisotropic plate model with warping functions for the study of vibrations reformulated from Woodcock's work[END_REF].

In order to solve Eq. ( 16) according to the bending wavenumber k = k f , one needs to solve the determinant of Eq. ( 16) det(K -ω 2 M ) = 0. The equation can be solved for ω or for k, in this work the choice has been to solve it according to k, which finally yields an 8th order equation. The equation can be simply solved by using a symbolic math solver software.

The solutions yield all the wavenumbers accounted by the modeling as in 5, where they are compared with SFEM results by Shorter [START_REF] Shorter | Wave propagation and damping in linear viscoelastic laminates[END_REF]. In order to make the comparison between the two results, the properties of the multilayered plate were chosen according to the case studied by Shorter. The characteristics of the multilayered are gathered in 1. All the wavenumbers obtained in a two dimensions kinematic field are plotted in this figure, so 5 curves are displayed. The highest wavenumbers correspond to the bending wavenumbers. The two curves linear with the frequency are the extensional waves in the (Ox) and (Oy) directions. Then, at higher frequencies shear appears in the structure in both direction of the (Oxy) plane. The dilatation wavenumbers are not appearing in the presented model since the vertical displacement W (x, y, t) is assumed constant through the thickness Eq. ( 1)

For the present work, only the bending wavenumbers are of interest. Its extraction is usually simple, because of the bending wave nature, their wavenumbers at a fixed frequency is mainly higher to the other wavenumber types. So they can be retrieve by taking the maximum wavenumber for a frequency. Otherwise, an algorithm can be used such as Newton-Raphson to get only one solution at the Eq. ( 16). The first proposition works in most of the cases but is not always a reliable technique because the bending wavenumber curve can be sometimes lower to the membrane wavenumbers for low frequencies. Once the bending wavenumbers have been retrieved, two parameters can be derived: the equivalent flexural rigidity D and the equivalent damping ratio η.

Equivalent flexural rigidity

The next parameter detailed in this subsection is the equivalent flexural rigidity. The goal is to explore the imperfection impacts thanks to another parameter. The equivalent flexural rigidity derivation will be first detailed, then compared to different cases, i.e. with or without imperfect interfaces.

Firstly, to compute the flexural rigidity one needs the bending wavenumber. This parameter is retrieved by solving the determinant of Eq. ( 16) thanks to a Newton-Raphson algorithm. Since we are interested only in the bending of the panel, the Love-Kirchhoff theory can be applied for each frequency to model this effect. This implies that the multilayered structure studied is, for a single frequency, seen as an equivalent thin plate. The methodology is applied for each frequency in order to have frequency-dependent equivalent parameters of the multilayered structure.

The governing equation for vertical displacement of the Love-Kirchhoff [START_REF] Love | The small free vibrations and deformation of a thin elastic shell[END_REF] model in isotropic conditions is:

D f ∂ 4 W ∂x 4 = -ρh ∂ 2 W ∂t 2 , (17)
with, D f the equivalent complex flexural rigidity of the isotropic plate, h = h i and ρ = hiρi hi . It is noteworthy to say that the isotropic assumption is done for the sake of simplicity, the model can be developed for anisotropic layers if needed in the same way as Marchetti's work [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF]. A particular solution can be used to solve the equation, such as: W = W 0 e j(k f x-ωt) . So the Eq. ( 17) yields Eq. ( 18), which is a straightforward equation that computes the flexural rigidity by using the bending wavenumber and the frequencies corresponding:

D f = ρhω 2 k 4 f . ( 18
)
This equation is particularly valid for beam-like structures since only one axis of propagation is considered.

The results for a sandwich panels are shown in 6. These are typical results for this type of structure and they can be divided into 3 zones: firstly, in low frequencies the multi-layered structure is behaving as an equivalent plate with the same flexural rigidity as the skins, but weighted by the thickness of the layers; in higher frequencies a transition is occurring because shear is appearing in the multi-layered structure, this makes the core, which is in most application a low rigidity material, involved in the behavior of the panel energy; finally, in the high frequencies, the layers are unbounded and so the flexural rigidity is mainly governed by the skins. The analysis of the first and last zones are more detailed in the next subsection.

Asymptotic analysis

The system of ( 16) can be used also to derive another equation of the flexural rigidity Eq. ( 19) as in [START_REF] Guyader | Viscoelastic properties of single layer plate material equivalent to a multi-layer composites plate[END_REF][START_REF] Arasan | A simple equivalent plate model for three-layer sandwich panels with shearing core[END_REF]. The computation is rather straightforward so only the result will be discussed here. By rearranging the system of equations the following result can be found out:

A 4 D 3/4 f + A 3 D f -A 1 A 4 D 1/2 f -A 1 A 3 + A 2 = 0, (19)
With

A 1 = λ 1 - λ 2 5 λ3 , A 2 = ω √ M (λ 4 -λ5λ6 λ3 ) 2 , A 3 = ω √ M (λ 2 - λ 2 6 λ3
) and A 4 = λ 37 , where the λ i are the Woodcock's coefficient [START_REF] Loredo | A multilayer anisotropic plate model with warping functions for the study of vibrations reformulated from Woodcock's work[END_REF], and M the total mass of the panel.

Accordingly, the asymptotic values of the flexural rigidity can be found out by computing the limits of the equation with respect to ω. If ω → 0: lim ω→0 D = D low , A 2 → 0 and A 3 → ∞ so

D low = A 1 ≈ D 1 (8 + 12 h 2 h 1 + 6 h 2 2 h 2 1
).

To obtain this equation approximation, the skins are assumed to be identical, which implies that the skins have the same flexural rigidity D 1 = D 3 , this is the case in most of the sandwich panels. Thus, the rigidity at the low frequencies is governed by the flexural rigidity of the skins and the ratio of thickness of the core over the thickness of the skins. So, a study on the effect of the thickness of the core shows that a thicker core yields a higher flexural rigidity at low frequencies.

On the other hand, if ω → ∞:

lim ω→∞ D = D high , A 2 ≫ (A 1 , A 4 ) and A 3 ≫ (A 1 , A 4 ) so D high = A 1 - A 2 A 3 ≈ D 1 + D 3 . (21)
Thus, the high frequency limit is simply the sum of the flexural rigidity of the skins. And so, in the case of the sandwich panels studied here, D high ≈ 2D 1 .

In this work, the asymptotic analysis has been re-implemented with imperfect interfaces in order to check what is the impact of the imperfections on the model asymptotic behavior and so increase our understanding on the result obtained for the flexural rigidity. When the propagation of the imperfect interface parameters is checked, one can notice that only the even λ i are affected by the modification. So only A 2 and A 3 are different from the previous results. In other words, only the high frequency limit D high is impacted by the imperfections. A 1 is not changed by the imperfection conversely to A 2 and A 3 . They can be written with respect to this modification by isolating the terms of the imperfection such as:

Ã2 = ω √ M (λ 4 + λ 4g ) -λ5 λ3 (λ 6 + λ 6g ) 2 and Ã3 = ω √ M (λ 2 + λ 2g ) - (λ6+λ6 g ) 2 λ3
. Finally they can be expressed as below:

Ã2 = A 2 + A 2g + 2 A 2 A 2g 2 , (22)
Ã3 = A 3 + A 3g , (23)
with,

A 2g ≈ h 1 (h 1 + h 2 )Q 1 11 Q 1 55 B 2 xz 1B 1 xz 2 and A 3g ≈ h 1 Q 1 11 (Q 1 55 ) 2 ( h2 Q 2 55 (B 1 xz + B 2 xz ) + (B 1 xz +B 2 xz ) 2 2
). The ω (M ) terms have been omitted since they vanish in the computation of D high thanks to the ratio between A 2 and A 3 . The imperfection has an impact on the asymptotic value only for important imperfections that finally would stand for fully debonded cases. Nevertheless, the approximation D high ≈ D 1 + D 3 is still a good for small damages as it will be shown in the application section.

Damping factor modeling

Another way to explore the impact of the imperfect interfaces in the model is to derive the equivalent damping factor. The damping implementation in the model has been done through a complex Young's modulus: E = E 0 (1+iη) of the layers, which yields the wavenumber to be complex, and so all the equivalent parameters. A classic way to estimate the equivalent damping ratio of the bending is by using the flexural rigidity as in Eq. ( 24):

η eq = Im(D f ) Re(D f ) = - Im(k 4 ) Re(k 4 ) . ( 24
)
However, as pointed out by Marchetti [START_REF] Marchetti | On the structural dynamics of laminated composite plates and sandwich structures; a new perspective on damping identification[END_REF], this equation based on Love-Kirchhoff's model overestimate the damping of the system because this model does not take into account the shear of the structure. At low frequencies, the behavior of a multilayered structure is depending mainly on pure bending. At higher frequencies, the behavior of the structure is changing because shear is appearing. However, in Love-Kirchhoff model the shear is not taken into account. This yields the ratio c g /c ϕ = 2 for any frequency, with c g = dω/dk and c ϕ = ω/k. In this case, the behavior modeled corresponds to a pure bending situation and overestimate the damping value when shear is occurring. In case of pure shearing of the structure the ratio c g /c ϕ = 1. Thus, the ratio of group velocity over phase velocity must vary between 2 and 1 depending on the frequency.

Furthermore, the damping ratio can be written through different formulations. The equivalent damping, as seen previously in the current subsection, will overestimate the damping at some frequencies because of the assumptions it implies. There is the spatial damping, which is obtained by computing the damping resulting of the propagation into space by using the particular solution W = W 0 e jk f x . The energy damping, obtained from a spatial formulation or a time formulation, which should give the same results as the first equivalent damping under the same conditions because both are computed by the ratio of the imaginary part and the real part of the flexural rigidity. All these damping ratio formulation are gathered into the 2.

The result for a sandwich panels is shown in 7. The low values for the low frequency limit and the high frequency limit indicates the damping of the skins, which are the only layers involved in the panel in these frequency intervals. A maximum damping is observed which corresponds to the optimum shearing zone, i.e. all the layers are involved in the behavior of the panel and particularly the high-dissipating core thanks to the shear of the panel. According to the damping chosen for the core, the equivalent damping shows its limits here by overestimating the maximum damping, which cannot be higher than the damping of the core.

Damping type

Equation

Equivalent Damping 3.

η eq = - Im(k 4 ) Re(k 4 ) Spatial Damping γ = Im(k) Re(k) Energy Damping η = -2 c g c ϕ Im(k) Re(k)

Application on sandwich panels

In this section is shown the results generated from the current model. The results are obtained for a three layered plate with the same properties for all simulations in order to check similarities and behavior modifications. The characteristics of the panel are shown in 3. In the previous papers [START_REF] Liu | Stress analysis of imperfection composite laminates with an interlaminar bonding theory[END_REF][START_REF] Cheng | Theory for multilayered anisotropic plates with weakened interfaces[END_REF][START_REF] Massàbo | Assessment and correction of theories for multilayered plates with imperfect interfaces[END_REF], the application is mainly about a two-layer or three-layer plate with the same thicknesses and the same material properties. In the application section, the system studied is a common type of sandwich panel that one can find in the industry, namely thin skins made out of aluminum with a soft and dissipative core to represent a polymer. This type of structure has been also studied and used as an example by the scientific community [START_REF] Marchetti | On the structural dynamics of laminated composite plates and sandwich structures; a new perspective on damping identification[END_REF][START_REF] Shorter | Wave propagation and damping in linear viscoelastic laminates[END_REF][START_REF] Mejdi | Wave spectral finite element model for the prediction of sound transmission loss and damping of sandwich panels[END_REF][START_REF] Ege | Vibrational behavior of multi-layer plates in broad-band frequency range: comparisons between experimental and theoretical estimations InterNoise 2015[END_REF][START_REF] Boutin | Generalized plate model for highly contrasted laminates[END_REF][START_REF] Arasan | On the accuracy limits of plate theories for vibro-acoustic predictions[END_REF].

In the results of the application, the perfect interface condition is always used as a reference. shown corresponds to the methodology previously detailed in 3. It first begins with the modification on the wavenumbers, to continue with the equivalent parameters, namely the equivalent flexural rigidity and the equivalent damping. Finally a comparison of the imperfection interfaces with a perfect interfaces model but with a softer core is shown. The parameter of defect has been chosen high enough to be physically valuable and fit illustrations, which corresponds to B s = 5e-10 m.Pa -1 .

Wavenumber modifications

As seen in 3.1, the wavenumbers considered in the model are bending, shearing and extensional wavenumbers. At higher frequencies other types of strain can appear such as breathing modes but are not taken into account in the presented model because the vertical displacement is assumed to be constant through the thickness, which implies that the vertical strain ϵ z = 0, which is a good assumption up to high frequencies.

In 8 the wavenumbers of the structure are plotted for different type of imperfect panels. The wavenumbers in the (Oy) direction are omitted, since they do not add information already brought by the ones in the (Ox) direction, in order to enhance the readability of the figures. The impact of one imperfect interface whether it is the first or the second interface is the same, which has been made to check the symmetry validity of the model. This could have been predicted from the impact of the imperfection in the displacement u x in 9. Indeed, whether the interface 1 or the interface 2 is imperfect, the displacement u x is the same for the skins. Since the behavior of a sandwich panel is mainly governed by the skins and the thickness of the layers, the equivalent parameters must be the same. Besides, this impact is clearly seen as a shift of the transition frequency towards the low frequencies, particularly visible on the bending mode thanks to its two slope changes. It is also particularly visible on the shear mode frequencies, they appear at lower frequencies for the imperfect interface 1 or 2 in comparison with the reference. No changes can be noticed for the extensional modes, they stay unchanged for any value of sliding interfaces. This is intuitive since the interface stresses are transverse shear stresses and so have effects only on the rotation kinematic variables in this model. This can be verified by deriving the transverse strain, which yields Eq. ( 25):

ϵ n xz = -α n xx Φ 1 x -α n xy Φ 1 y (25a) ϵ n yz = -α n yx Φ 1 x -α n yy Φ 1 y (25b)
Finally if both interfaces are imperfects, the effect is even higher since the total amount of imperfections is higher too. These results give evidence that such imperfections facilitate the shear between the layers, which is seen in 8 but also in the following subsection. An example of the resulting displacement field with imperfect interfaces is displayed in 9, it shows the discontinuities introduced in the displacement field depending on which interface is impacted by the imperfection. The displacement displayed is obtained from the computed wavenumbers. It shows the displacement in the x direction u n

x for each layer n in case of free vibrations. Thus, the figure provides an insight of the imperfection impact on the displacement field. An interesting result is that whether the first or the second interface is imperfect, the strains of the skins are the same. One can see that the imperfections change the strains evolution in the core since they are more and more independent.

Equivalent flexural rigidity and damping shift

As seen in 3.2, the equivalent flexural rigidity for a three layered plate is changing with respect to the frequency domain according 3 regions namely: the low frequency domain; the transition domain; and the high frequency domain. This behavior is seen in 10 for two kinds of simulation: a perfect interfaces panel; and a panel with imperfect interfaces. Thus, the imperfections have a similar effect as seen previously, i.e. the whole curve is shifted towards the low frequencies, and for this kind of structure does not show any additional impact than seen on the bending wavenumbers. Since the transition frequency is depending on the shear of the panel, the same conclusion can be drawn from this figure.

The damping used in the current results will only be the energy damping, which as seen in 3.4, is corrected by the ratio between the group velocity and phase velocity of the waves in order to take into account the shear of the structure. As we could expect from the previous results, the frequency corresponding to the maximum damping is shifted towards the low frequencies, because the location of this frequency is depending on the shear of the structure. Moreover, it is noteworthy to say that the damping of the panel is only shifted towards the low frequencies. Even if the defect is particularly high, no damping is added to the structure nor removed, only the transition between the medium frequency behavior and the high frequency behavior happens more easily. This indicates that the imperfections does not change the whole behavior of the panel but makes the shearing easier to happen. This is intuitive since the imperfections taken into account are sliding in the transverse directions. However, some damping can be added to the system if needed thanks to the imperfect parameter B i , by increasing its imaginary part.

Comparison with perfect interfaces with different parameter variables

As depicted by the previous results, the main impact of this new modeling is that the imperfections shift the transition frequency towards the low frequencies. Even though this behavior is intuitive, it can be already modeled simply with a perfect interface equivalent model by using a softer core as shown in 11.

The advantages of keeping the previous model and to model imperfection by reducing the Young's modulus of the core are quite obvious since its results are already used and have been experienced, and its modeling is simpler because it does not require any more computations to describe such effects. Nevertheless, to model imperfection with perfect interfaces by reducing the Young's modulus of the core is not physically valuable and can lead to errors on material parameters of the core. Thus, using the model with imperfect interfaces gives more insight into the physics of multilayered plate. As in 11, a rather important sliding defect is equivalent to a core approximately 8 times softer, which is too high to be acceptable as experimental uncertainties. Conversely, small imperfections that yield a slight shift of the equivalent parameters are hardly detected experimentally. They can be mistaken for a softer core, mainly if the variation of the rigidity is of the same order than the uncertainties.

Moreover, the application in the present work has been chosen for its simplicity as a first step for the model in order to see if it contains some particularities. Further applications are considered, with new assumptions in the model such as 3rd order vertical displacement or non-linear interface constitutive equation.

Comparison with Ross, Kerwin and Ungar (RKU) model

In order to validate the asymptotic behaviors of the dynamic model, a comparison is done with a model developed by Ross, Kerwin and Ungar model [START_REF] Ross | Damping of plate flexural vibrations by means of viscoelastic laminae[END_REF], denoted here by RKU. The RKU model is a simplified model of the dynamics of a constrained layer system. It is used to model perfectly bonded layers dynamics and so, can be compared to the presented model in the case B = 0 m.Pa -1 . The comparison is made in Fig. 11 and shows good agreement between the RKU model and the presented model. For the fully debonded layers case, no multilayered plate dynamics are involved. Thus, the flexural rigidity of the system must be the sum of the flexural rigidity of the three layers D = D 1 + D 2 + D 3 , which implies the flexural rigidity to not be a frequency dependent parameter. If the interface stiffness tends towards infinity B → ∞, the model shows good agreement with this result. Furthermore, the equivalent dynamic model with imperfect interfaces is able to model the transition between these two asymptotic behaviors.

Application to various sandwich panels

In this section are gathered the results for various sandwich panels with different mechanical parameters and different geometry. The cases were chosen according to the work presented here, but also to the work of Ghinet and Atalla [START_REF] Ghinet | Modeling thick composite laminate and sandwich structures with linear viscoelastic damping[END_REF]. All the results are summarized in 4.

One can find the input of each case needed for the model to yields results. For each case, two results are always computed: The first result is obtained with perfect interfaces B = 0 m.Pa -1 ; and the second with B = 5e-10 m.Pa -1 similarly as the first application. Therefore, it is possible to compare the effect of the interface parameter depending on the sandwich panel characteristics. The effect is seen thanks to ∆f = f T,B=0 -f T,B=5e-10 . The frequency values correspond to the frequency where the maximum damping is obtained, i.e. where the maximum shearing is obtained. So, the difference of the frequency values ∆f represents the frequency shift of the dynamical parameters. Some conclusions can be drawn in the light of the results shown in 4. The impact of the interface compliance B is different depending on the sandwich panel properties. Higher is the stiffness of the skins E 1 = E 3 in comparison to the core, lower will be the impact of the interface parameter B. This can be checked by comparing the case #1 with the case #3, where the latter has stiffer skins comparing to the core. Similarly, thicker is the core h 2 in comparison to the thickness of the skins h 1 = h 3 , higher will be the impact of the interface parameter B. The flexural rigidity asymptotic values D low and D high are unchanged since the interface parameter impact only the coupling of the layers and not the rigidity of the layers.

Discussion

Since the 3D-displacement of the plate is assumed to be linear with respect to the vertical axis and the layers have the same properties, the vertical stress is constant through the thickness. It has been shown that this assumption has meaningless impact on the equivalent parameters for the bending motion in case of perfect interfaces as in Appendix of [START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF]. But, since the interface stresses are strongly impacted by the imperfections, especially for highly imperfect conditions such as fully debonded layers, 3rd order vertical displacement will be considered and its influence will be checked on the results.

Moreover, it has been shown in [START_REF] Marchetti | On the condensation of thick symmetric multilayer panels including dilatational motion[END_REF] that dilatational motion represents the high frequency limit of the present model but also have strong impact on the behavior of a multilayered structure even in low frequency range (f ∈ [200, 500] Hz). This is a current limitation of the presented model since no dilatation is assumed here. Thus, take into account dilatational motion would increase the high frequency limit of the model but would also be interesting to model openings ûn z in the panel in order to implement more imperfection types. Furthermore, the constitutive law, and especially the interface parameters B ij or K ij are assumed to be constant in the present work. This was an unavoidable step for the model development since it provides a preview on the impact of the imperfection implementation. But this limits the present model to describe only sliding imperfect interfaces linearly with frequency. However, in order to fully take advantage of the model, other types of interface parameters are considered to be explored, such as frequency varying interface parameter or non-linear constitutive law described through piece-wise linear interface parameter.

Conclusion

This paper proposes a method to model imperfect interfaces, and in particular sliding interfaces for a better characterization and optimization of the behavior of multilayered structures. This implementation is done into an equivalent layer model, the Guyader model, a Zig-Zag model, recently improved by Marchetti for anisotropic plates. The advantage of such models is to yield reliable results efficiently by describing the behavior of all the layers of the panel one by one, but thanks to a single reference layer. For example, this allows to show the global displacement through the layers as well as equivalent parameters to model the panel as if it was an unique layer with specific properties. It is noteworthy to say that the presented method is focused on sandwich panels but can be used as well for other types of multilayered structures.

The results of the method show the impact of imperfect interfaces when added to the model, that is mainly the shift towards the low frequencies of the typical behavior of the frequency-dependent parameters for a sandwich panels. Since the transition between the first and the second behavior of a sandwich panels is depending on the shear between the layers, it is intuitive that the imperfection, which here is an interface sliding, will facilitate such behaviors. This effect is illustrated throughout the application done for this paper with several parameters: firstly, the set of wavenumbers shows this shift towards the low frequencies where the most obvious shift is the one of the shear which appears in lower frequencies. Then, the flexural rigidity shows this same effect by the shift of the transition. Finally, the damping, which implies the imaginary part into the process, shows the same effect, that is to say, the maximum damping is also shifted.

The model could be used to characterize the coupling conditions in a multilayered panel. It could be challenging to apply the method experimentally on several sandwich panels and analyze them in the scope of the presented model. It would be interesting to apply the method for Structure Health Monitoring (SHM). Adhesive defects that impact the dynamical behavior of the multilayered structure could be detected and characterized thanks to the shift of the dynamical parameters. Furthermore, the aging of systems could be analyzed in the light of the model, by checking the dynamical characteristics evolution.

Figure 2 :

 2 Figure 2: Sandwich panel with imperfect interfaces modelled by interface jumps and interface stiffness

  t y ] the vector of interface losses. Such equations describe the bonding behavior between consecutive layers with 2 limit cases. Let us write σn xz = K n xz ûn x or in its inverse form ûn x = B n xz σn xz , which corresponds to an isotropic plate. If K n xz → 0, so σn xz → 0, which describes a fully debonded case because no energy communication is carried out between the layers. If B n

Figure 3 :

 3 Figure 3: Illustration of two different cases for the constitutive equation of imperfect interfaces. (a) purely elastic interface with piecewise linear interface parameter, (b) linear non-proportional interface with piecewise linear interface parameter.

Figure 4 :

 4 Figure 4: (a) Wavenumbers k i computation workflow with the displacement U Eq (7) as input, (b) Flexural rigidity D f and damping η workflow with the bending wavenumbers k f as input.

Figure 5 :

 5 Figure 5: Full set of wavenumbers in 3D obtained by solving the system of equation (16), comparison between the presented model and SFEM [28]. (a) Bending wavenumbers, (b)-(c) Extensional wavenumbers, (d)-(e) Shear wavenumbers, and (f) Dilatation wavenumbers (SFEM only).

Figure 6 :

 6 Figure 6: Flexural rigidity of a sandwich panel as a function of frequency. The characteristics of the panel considered are gathered in Table 1

Figure 7 :

 7 Figure 7: Damping of a sandwich panel as a function of frequency. Two models used: the Equivalent Damping and the Energy Damping, where the first overestimate the maximum damping. The characteristics of the panel considered are gathered in Table3.

Figure 8 :

 8 Figure 8: Full set of wavenumbers in 2D, with perfect interfaces and one or two imperfect interfaces.

Figure 9 :

 9 Figure 9: Displacement ux through the thickness z of the panel for different imperfection namely: Only perfect interfaces, first interface imperfect and second perfect, first interface perfect and second interface imperfect, all interface imperfects.

Figure 10 :

 10 Figure 10: (a) Equivalent Flexural Rigidity D and (b) Equivalent damping η with perfect interfaces, one imperfect interface and only imperfect interfaces.

Figure 11 :

 11 Figure 11: Equivalent flexural rigidity comparison between an imperfect interfaces equivalent model and a perfect interfaces equivalent model but with a softer core.

Case # 1 :Case # 3 :Case # 4 :Table 4 :

 1344 E 1 = E 3 = 71e9Pa, E 2 = 30e6Pa, h 1 = h 3 = 1mm, h 2 = 5mm, ρ 1 = ρ 3 = 2700kg.m -3 , ρ 2 = 48kg.m -3 , ν 1 = ν 3 = 0.33, ν 2 = 0.2. f T,B=0544 Hzh 1 /h 2 0.2 D low 1448 N.m f T,B=5e-10 42 Hz E 1 /E 2 2.3e3 D high 13 N.m ∆f 502 Hz Case #2: E 1 = E 3 = 69e9Pa, E 2 = 2.1e6Pa, h 1 = h 3 = 1.52mm, h 2 = 0.127mm, ρ 1 = ρ 3 = 2800kg.m -3 , ρ 2 = 940kg.m -3 , ν 1 = ν 3 = 0.33, ν 2 = 0.2. E 1 = E 3 = 125e9Pa, E 2 = 5e6Pa, h 1 = h 3 = 1mm, h 2 = 3mm, ρ 1 = ρ 3 = 1600kg.m -3 , ρ 2 = 33kg.m -3 , ν 1 = ν 3 = 0.33, ν 2 = 0.2. E 1 = E 3 = 210e9Pa, E 2 = 15.8e6Pa, h 1 = h 3 = 0.453mm, h 2 = 0.035mm, ρ 1 = ρ 3 = 7850kg.m -3 , ρ 2 = 1000kg.m -3 , ν 1 = ν 3 = 0.33, ν 2 = 0.2. f T,B=01233Hz h 1 /h 2 13 D low 16 N.m f T,B=5e-10 0.03 Hz E 1 /E 2 1.3e5 D high 3.7 N.m ∆f 1233 Hz Multilayered plate mechanical parameters With, f T,B=0the frequency of maximum shearing, which correspond to the maximum damping value, for the perfectly bonded case B = 0 m.Pa -1 . f T,B=5e-10 the frequency of maximum shearing for the imperfect case with B = 5e-10m.Pa -1 . ∆f = f T,B=0 -f T,B=5e-10 the frequency shift of the dynamical parameters. D low the low frequency asymptotic value of the flexural rigidity. D high the high frequency asymptotic value of the flexural rigidity.

Table 1 :

 1 Multilayered plate mechanical parameters from Shorter's SFEM

	Skins	0.6	2700	71 000	0.33 0.5%
	Core	15	48	30	0.2	10%

Table 2 :

 2 Different equivalent damping parameters, the equivalent damping does not consider shear effects and is based on energy formulation, the spatial damping is based on spatial formulation and the energy damping is similar to the equivalent damping but takes into account shear effects. The definition are gathered from the work of Marchetti in[START_REF] Marchetti | On the structural dynamics of laminated composite plates and sandwich structures; a new perspective on damping identification[END_REF][START_REF] Marchetti | Modélisation et caractérisation large bande de plaques multicouches anisotropes[END_REF]

		0.2
			Energy damping (corrected)
			Equivalent damping
		0.15
	Damping	0.1
		0.05
		10 0 0	10 1	10 2	10 3	10 4	10 5	10 6
				Frequency [Hz]

Table 3 :

 3 Multilayered plate mechanical parameters

	The results

Acknowledgements

This work was supported by the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" operated by the French National Research Agency (ANR). We also express our gratitude to Fabien Marchetti (Research Engineer at Matelys) for fruitful discussions.