© The Author(s), Published by EDP Sciences, 2024

1 Introduction

Multilayered structures are materials typically made by stacking different plates or beams with various mechanical and physical parameters. They aim to optimize performance by taking advantage of the specific characteristics of each layer. This results in materials that effectively absorb sound while remaining lightweight, which is of particular interest in building engineering. Lightweight structures are also extensively used in the automotive industry, but in this field, stiffness is a key parameter. Carrera [1] has conducted an extensive review of analytical modeling of multilayered structures. Equivalent dynamic models are employed to represent the dynamics of a multilayered structure using a similar equivalent structure, like a thin plate made of an unique layer. Such models are widely used to compute mechanical parameters more efficiently compared to the Finite Element Method (FEM), for instance, because they do not require a detailed mesh. Most of the time, these existing models assume perfect coupling with continuity between two layers. In the following, we will refer to these types of interfaces as “perfect interfaces”.

However, the actual coupling conditions within manufactured multilayered materials can result in some unexpected dynamical behaviors. Due to aging, damage occurring under operating conditions, or the manufacture itself, the coupling conditions may produce a different behavior than expected. Discontinuities (in displacement, stress, etc.) can arise at the interface between two layers. In the following, we will refer to these types of interfaces as “imperfect interfaces”. To model such effects, it is necessary to account for imperfect interfaces in the modeling of multilayered materials. Several studies on imperfect interfaces can be found in a few multilayer modeling works, such as [2–8]. These works are mainly based on the model proposed by Pagano in 1969 [9]. Unfortunately, as described by Carrera in [1], these models are not equivalent single-layer models. Thus, it is challenging to apply some equivalence to effectively describe the dynamical behavior of a multilayered material with imperfect interfaces. To address this, the authors have previously proposed an equivalent dynamic model that incorporates imperfect interfaces [10].

Experimentally, several studies have provided evidence of imperfect interfaces in multilayers through ultrasonic wave measurements, such as [11–14]. These studies suggest that damage to the interface, such as inclusions or delaminations, intuitively reduces the coupling quality between the consecutive layers. However, the impact of such damage on the overall dynamics of a multilayered material has not yet been demonstrated. Likewise, no experimental dynamical characterization of the impact of an imperfect interface has been conducted.

This work aims to experimentally characterize the dynamical behavior of imperfect interfaces within a multilayered structure. The characterization uses an equivalent dynamic model developed by the authors [10] to identify a theoretical parameter that represents the interface coupling conditions. This model is briefly described in Section 2. It is shown here that this characterization helps assess the impact of imperfect interfaces on the overall dynamics of a multilayered structure through a single parameter. This parameter, denoted as B, aggregates the contribution of all the factors involved in the interface mechanics according to an unknown law. The experimental methodology is fully presented in Section 3. First, the setup used to perform the vertical displacement measurement is detailed. Next, the data processing employed in this study is briefly outlined to emphasize the main features and their implementation. Additionally, the filtering methodology is explained in this section. Finally, the results are presented and discussed in Section 4. These results are analyzed in terms of space and frequency to illustrate the impact of imperfect interfaces and characterize the interface parameter B.

2 Model overview

This section is a brief summary of [10], the reader is referred to this paper for more details. The displacement field equation (1) used for the model is the 3D displacement of the layer n in a multilayered plate with N layers [15]. A representation of a sandwich panel with imperfect interfaces (represented by springs to indicate that there is a coupling law between layers) is shown in Figure 1. According to Mindlin’s or Reissner’s plate theory [16, 17], three types of displacements taking into account the physics are described by this model: the extension ψ(x, y, t), the bending W(x, y, t) and the shearing ϕ(x, y, t). In this work, a 1st order displacement is assumed, and the vertical displacement is constant through the thickness z. However, Marchetti has shown that a 3rd order displacement field with respect to z yields only a slight increase in the prediction accuracy of the dynamics of the multilayered structure [18]. Since no breathing is taken into account in the model, because of the constant vertical displacement W(x, y, t), the high frequency limit of model validity is defined by the first breathing mode. Their frequency depends highly on the structure geometry and mechanical properties.

$$\{\begin{array}{l}{u}_x^n={\psi }_x^n\left(x,y,t\right)-\left(z-z_n\right)\left(\frac{\partial W\left(x,y,t\right)}{\partial x}+{\varphi }_x^n\left(x,y,t\right)\right),\\ u_y^n={\psi }_y^n\left(x,y,t\right)-\left(z-z_n\right)\left(\frac{\partial W\left(x,y,t\right)}{\partial x}+{\varphi }_y^n\left(x,y,t\right)\right),\\ u_z^n=W\left(x,y,t\right).\end{array}$$(1)

Figure 1 Sandwich panel with imperfect interfaces (represented by springs to indicate that there is a coupling law between layers). According to this representation, only relative motion in x and y directions is accounted for.

with, ψ_x (ψ_y) the translation over the Oxⁿ (Oyⁿ axis), which simply describes the displacement yielded by the extensional waves in the layer n, z_n is the z mid-plane coordinate of the layer n with respect to the global axis. ${\varphi }_x^n$ (${\varphi }_y^n$) is the rotation around the Oy (Ox) axis, which depends on the shear of the structure. The model assumes laterally infinite dimensions.

The continuity of the transverse shear stresses σ_αz (with α either x or y) and the discontinuity of the displacement U, shown in equation (2) which applies at interfaces, are used to implement the imperfect interfaces into the model as in [3, 4, 19].

$$\{\begin{array}{l}{\sigma }_{\alpha z}^n={\sigma }_{\alpha z}^{n-1},\\ {\mathit{U}}^n={\mathit{U}}^{n-1}+{\widehat{\mathit{U}}}^n,\end{array}$$(2)

with ${\widehat{\mathit{U}}}^n$ the sliding displacement field of the layer n. Even if ${\widehat{\mathit{U}}}^n$ is a vector here, for the experimental characterisation section only the in-plane components will be used.

These interface condition equations make it possible to write the displacement field of the layer n − 1 with respect to the layer n and so on. Thus, the displacement field Uⁿ can be written as Uⁿ = TⁿU^n − 1, with Tⁿ the transfer matrix from the layer n to n − 1. Finally, the reference layer i = 1 can be written according to the layer i = n by using the transfer matrix as much as there are interfaces as in equation (3). Conversely to [15, 20], the sliding field $\widehat{\mathit{U}}$ is added during this homogenization step in equation (3) since here the interfaces are imperfect. Naturally, the sliding field $\widehat{\mathit{U}}$ of the interface i = 1 is null because there is no layer before the layer i = 1 and so there is no interface i = 1.

$${\mathit{U}}^n={\mathit{T}}^n{\mathit{T}}^{n-1}\dots {\mathit{T}}^2{\mathit{U}}^1+\sum _{i=1}^n\mathrm{}{\widehat{\mathit{U}}}^i.$$(3)

Equation (4) shows the propagated displacement field with imperfect interfaces.

$$\{\begin{array}{l}{u}_x^n={\psi }_x^1(x,y,t)+F_{\omega }\frac{\partial W}{\partial x}+F_{\mathrm{xx}}^n{\varphi }_x^1+F_{\mathrm{xy}}^n{\varphi }_y^1+\sum _{i=1}^n\mathrm{}{\widehat{u}}_x^i,\\ u_y^n={\psi }_y^1(x,y,t)+F_{\omega }\frac{\partial W}{\partial y}+F_{\mathrm{yx}}^n{\varphi }_x^1+F_{\mathrm{yy}}^n{\varphi }_y^1+\sum _{i=1}^n\mathrm{}{\widehat{u}}_y^i,\\ u_z^n=W\left(x,y,t\right),\end{array}$$(4)

with, F_ω = z₁ − z, $F_{\mathrm{ij}}^n={\alpha }_{\mathrm{ij}}^n(z_n-z)+{\gamma }_{\mathrm{ij}}^n$. ${\alpha }_{\mathrm{ij}}^n$ and ${\gamma }_{\mathrm{ij}}^n$ are material parameters (mainly Young’s modulus ratios of the consecutive layers) yielded by the transfer matrix [Tⁿ]. More details are given in Appendix of Loredo and Castel’s work [21] or by Marchetti [18]. However, in order to do further computation, the sliding $\widehat{\mathit{U}}$ is written with respect to the stresses at the interface of the multilayered plate thanks to the constitutive equation (5).

$${\widehat{\mathit{U}}}_{z=z^n}^n={\mathit{B}}^n{\widehat{\mathit{\sigma }}}^n.$$(5)

This yields the displacement field to be written exclusively by the kinematic variables of the first layer. It is also used to implement the bonding condition laws. Equation (5) expresses the in-plane displacement discontinuity at interfaces as a function of (some) stress components. This is a matrix equation, but later on, we will use only a single scalar value of interface parameter B for each layer since only beams will be considered for practical experiments. If we consider three layers with same interfaces (symmetrical structure with same skins) then a single value of B will be considered for the whole structure. If Bⁿ = 0 m.Pa⁻¹ the interface is perfect since the displacement obtained is the same as in the former models. If B → ∞, or the interface stiffness K = 0 Pa.m⁻¹, the interface is fully debonded. It should be underlined that the value of B is not so easy to interpret. In fact, other choices could be done: for example a dimensionless parameter may be used. However, this would introduce a parameter with no physical sense which would be more difficult to interpret and the question arises to know by which value of the stiffness (stiffness of the core? stiffness of the skin?) the compliance B could be divided. So the choice has been made here to let this value as a proportionality law between stress sigma and sliding displacement. The aim of this parameter is to remain a physical parameter that can be experimentally measured in further studies: work is in progress to quantify it experimentally. It should be also underlined that even if this value can vary between 0 and ∞ for which perfect coupling corresponds to 0 and perfect decoupling corresponds to ∞, a “small” value of B can already lead to a “quasi-perfect” decoupling. The final displacement field used in this work is equation (6).

$$\{\begin{array}{l}{u}_x^n={\psi }_x^1(x,y,t)+F_{\omega }\frac{\partial W}{\partial x}+F_{x{x}_g}^n{\varphi }_x^1+F_{x{y}_g}^n{\varphi }_y^1,\\ u_y^n={\psi }_y^1(x,y,t)+F_{\omega }\frac{\partial W}{\partial y}+F_{y{x}_g}^n{\varphi }_x^1+F_{y{y}_g}^n{\varphi }_y^1,\\ u_z^n=W\left(x,y,t\right),\end{array}$$(6)

where F_xx became $F_{x{x}_g}=F_{\mathrm{xx}}-B({\alpha }_{\mathrm{yx}}^n{Q}_{45}^n+{\alpha }_{\mathrm{xx}}^n{Q}_{55}^n)$, F_yy became $F_{y{y}_g}=F_{\mathrm{yy}}-B({\alpha }_{\mathrm{yy}}^n{Q}_{44}^n+{\alpha }_{\mathrm{xy}}^n{Q}_{45}^n)$, F_xy became $F_{x{y}_g}=F_{\mathrm{xy}}-B({\alpha }_{\mathrm{yy}}^n{Q}_{45}^n+{\alpha }_{\mathrm{xy}}^n{Q}_{55}^n)$ and F_yx became $F_{y{x}_g}=F_{\mathrm{yx}}-B({\alpha }_{\mathrm{yx}}^n{Q}_{44}^n+{\alpha }_{\mathrm{xx}}^n{Q}_{45}^n)$, with Q_ij stiffness values following Voigt’s notation (The details are given in [18]). So values Q_ij depend on the Young’s modulus. It is noteworthy to say that in this case, no disbonds are taken into account, which leaves $u_z^n$ 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. To derive dynamical parameters, the lagrangian is computed from the displacement field. This yields a system of equations that can be solved using the particular solution of a propagating wave. Once solved, the solutions obtained are the wavenumbers of the multilayered system (bending, shearing, and membrane wavenumbers). Finally, the bending wavenumbers k_f are used to compute the flexural rigidity D_f of the system thanks to Love-Kirchhoff’s [22] thin plate theory similarly to equation (7):

$$D_f=\frac{\rho h{\omega }^2}{k_f^4},$$(7)

where h = ∑h_i is the total thickness with h_i the thickness of each layer i, $\rho =\frac{\sum \mathrm{}{\rho }_i{h}_i}{\sum \mathrm{}{h}_i}$ the equivalent mass density with ρ_i the mass density of the layer i, and ω the angular frequency in rad.s⁻¹. It should be noticed that the flexural wavenumbers are the highest wavenumber for a given frequency.

Damping can be introduced by different manners. As underlined in [10] interface losses may be introduced in the constitutive equation (5) or directly in the Young’s modulus. In the present paper, it was chosen to introduce damping directly in the Young’s modulus E leading to complex wavenumbers. That is why, in the following a complex Young’s modulus E^* will be considered:

$$E^{\mathrm{*}}=E(1+\mathrm{i\eta })$$(8)

where η is the loss factor and E the modulus. More details are given in [10] to identify then the complex equivalent bending stiffness. It should also be underlined that for characterising materials used in the present study, the Poisson’s ratio ν of each layer will be given.

A classic way to estimate the equivalent damping ratio η_eq of the bending is to use the flexural rigidity:

$${\eta }_{\mathrm{eq}}=\frac{\mathrm{Im}\left(D_f\right)}{\mathrm{Re}\left(D_f\right)}=-\frac{\mathrm{Im}\left(k_f^4\right)}{\mathrm{Re}\left(k_f^4\right)}.$$(9)

However, as pointed out in [10], this equation based on Love-Kirchhoff’s model overestimates 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_f the group velocity and c_ϕ = ω/k_f the phase velocity. In this case, the modeled behavior 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. So the following energy damping (which will be more interesting) is used in the present study:

$$\eta =-2\frac{c_g}{c_{\varphi }}\frac{\mathrm{Im}\left(k_f\right)}{\mathrm{Re}\left(k_f\right)}.$$(10)

The impact of the imperfect parameter, namely the interface compliance B, is shown in Figure 2. First, the flexural rigidity has a low-frequency asymptotic behavior, which corresponds to bending strains only. Secondly, the shearing effect of the constrained core layer is increasing, which causes the flexural rigidity of the system to decrease. Finally, the bending stiffness decreases with respect to the frequency until the system is fully debonded. This tends to the high-frequency asymptotic bending stiffness and is simply the sum of the flexural rigidity of the layers D_high = D₁ + D₂ + D₃. The impact of the interface parameter B is to shift towards the low frequencies the dynamical parameters of the system. This is intuitive since the transition occurs when the system is shearing, so an imperfect coupling of the layers in the transverse plane allows the system to shear more easily.

Figure 2 Flexural rigidity with respect to the frequency for two different interface conditions. For the blue curve, the simulation has been done with B = 0 m.Pa−1, which is a perfect interface structure. The red dotted line has been done with the same material but with imperfect interfaces with B = 10e−9 m.Pa−1 = 1e−8 m.Pa−1. The material properties for this example are gathered in Table 1.

An equivalent analytical model of multilayers with imperfect interfaces has been introduced in this section. The model introduced here is used in Section 4 to characterize the dynamical behavior of experimental samples. The main goal will be to characterize the interface parameter B.

3 Experimental methodology

In this section, the experimental methodology employed in this work is detailed. It was used to characterize samples with imperfect interfaces. First, the data acquisition setup is presented, showing how a vibrometer and piezoelectric buzzers were used to perform measurements over a wide frequency bandwidth, from 1 kHz to 100 kHz. The samples fabricated and characterized in this work are presented in Table 2 and shown in different pictures below. Two sets of beams were manufactured at two different time periods: aluminum beams and glass beams. A quasi-similar epoxy adhesive was applied for both sets, assuming same interface properties for upper and lower skin-core interfaces. Mass was measured, and the variations in epoxy mass resulted in variations in total weight. There were small mass variations in the beams due to imperfections. The density and Young’s moduli (for glass and aluminum) were taken from the literature, and the values for the adhesive were inverted from measurements (for the perfect B = 0 case).

Then, the data processing method, the Corrected Force Analysis Technique (CFAT) is briefly presented for a beam-like structure. Finally, the filtering methodology applied in this work to increase the frequency bandwidth of the results is presented. The results obtained thanks to this methodology are shown further in Section 4.

3.1 Experimental validation on manufactured multilayered beams

In this study, we are interested in the variations of the Young’s modulus due to the interfaces and not on the overall value of the rigidity. Moreover, from an experimental point of view, it is complicated to create plates with complex interfaces that is why in the following we chose to perform experimental validation on beams. We do an approximation.

The experiment allows us to estimate wave numbers and to determine a “structure” parameter $s=\frac{\mathrm{EI}}{\rho S}$ (CFAT) linked to an equivalent Young’s modulus of the sandwich beam (I is the moment of inertia, S the cross-section area).

A comparison between experiment and theory of “equivalent Young’s moduli” is displayed on the “result” figures of next part. The theoretical E is recalculated via the analytic equivalent wavenumber estimated in the first part.

The estimated analytical wave number k_N for a plate gives us:

$$\frac{\mathrm{EI}}{\rho S}=\frac{{\omega }^2}{k_N^4}.$$(11)

where k_N = k_f (with k_f the wavenumber calculated in Sect. 2 in the model review). So:

$$E=\frac{D_{f}S}{Ih}.$$(12)

As pointed out in [20], consistent results between plate and beam in the case of perfect interfaces have been obtained. So the same approach is used here.

3.2 Setup

Laser Doppler Velocimetry (LDV) was the technique used in this work to measure the displacement of the freely suspended samples. For this purpose, a vibrometer, the PSV-400 made by Polytec¹, operating in the visible range, has been used. The aim of the setup was to perform vibration measurements in order to characterize the dynamic behavior of a planar structure, particularly a beam-like structure. Figure 3 illustrates the set of tools and their connections that make up the bench. The various elements of setup are:

• The sample being measured, which is chosen from a set of multilayered beams manufactured in the laboratory.

• The frame housing the measured sample.

• The excitation system, which in the present work is a piezoelectric buzzer. The diameter of the buzzer is 12 mm and its resonance frequency is 9 kHz.

• The laser vibrometer to do the optical measurement of the displacement field.

Figure 3 Picture of the setup consisting of: the measured sample, the vibrating device and the laser vibrometer. The picture is shown next to a diagram of the sample from the laser vibrometer point of view, with a representation of the scanning measurement mesh.

The samples are suspended as hanging beams to simulate free-free boundary conditions. Both ends of the samples are tied with a thin polypropylene rope to prevent continuous pendulum motion. The structure is excited using piezoelectric buzzers. Their low mass and small size make piezoelectric buzzers an interesting, non-invasive vibrating device that has minimal impact on the behavior of the structure under test. Any potential lack of signal due to their low vibrating amplitude is compensated by using two or more piezoelectric buzzers. This increases the amplitude level of the sample displacement field being measured.

For the results presented here, the input signal is a coherent white noise for both piezoelectric buzzers, with a frequency bandwidth generally ranging from 1 kHz to 100 kHz. The mesh size target was 2.8 mm and the diameter of the piezoelectric buzzers were ∅ = 12 mm. All these parameters are enclosed in Table 3.

The spatial derivatives of displacements are estimated along the length of the beam (long dimension of the beam). Several meshes were used during measurements (in particular to estimate the global coherence on the whole structure) that is why the LDV scans show three points along beam width in Figure 3. For this study the choice was made to keep only the central line to estimate these derivatives along the long dimension of the beam.

3.3 Corrected force analysis technique

Several inverse methodologies have been developed by the vibroacoustical community to identify the structural parameters (Young’s moduli, loss factors) of beams, plates or multilayered composites structures, in wide frequency band. One can cite here the Virtual-field Method (VFM) that solves the principle of virtual work using virtual field functions [23–26], or the recent High-Resolution Wavevector Analysis (HRWA) extracting wavenumbers using 2D-ESPRIT algorithm [27]. Some techniques based on function fitting also exist, using for example plane waves such as the Inhomogeneous Wave Correlation method [28]) or, with more accuracy, Hankel functions [29, 30].

The Force Analysis Technique (FAT) derives the structural parameter from an estimation of the spatial derivatives of the equation of motion [31–34]. The inverse method chosen for the present work is the CFAT [35–38] that corrects the bias error of the finite difference scheme approximation in the FAT technique. The methodology is similar to a wave number analysis method, such as spatial Fourier transform analysis, but adds the assumption that the experimental structure studied corresponds to a theoretical plane structure. This leads the results to be more accurate as far as the sample fits the assumption of the structure studied. CFAT is a local characterization method based on the application of the Fundamental Principle of Dynamics to a thin beam in bending. For characterization purposes, only the structural parameters are of interest. They are gathered by the single parameter s, defined and computed in equation (13), called the structural parameter.

$$s=\frac{\mathrm{EI}}{\rho S}={\omega }^2{\Delta }^4{\left[\arccos \left(1-\frac{{\Delta }^2}{2}\sqrt{\frac{{\delta }_{\Delta }^{4x}}{w}}\right)\right]}^{-4},$$(13)

with, E the apparent Young’s modulus of bending, I the moment of inertia, ρ the density, S the cross-section area, ω the pulsation in rad.s⁻¹, Δ the distance between two measurement points (mesh size), w(x) the vertical displacement as a function of x with the vertical axis parallel to the material thickness, and ${\delta }_{\Delta }^{4x}$ estimation of the 4th derivative (bilaplacian) of the displacement w(x) through the finite different scheme.

For beams, a 5-point finite-difference scheme defined in equation (14) is used to estimate the 4th derivative of the displacement w(x).

$${\delta }_{\Delta }^{4x}=\frac{w(x-2\Delta )-4w(x-\Delta )+6w\left(x\right)-4w(x+\Delta )+w(x+2\Delta )}{{\Delta }^4}.$$(14)

The optimal mesh size Δ is computed to have a spatial mesh of three points per wavelength [18, 36]. Therefore, Δ depends on the frequency processed and must be adaptive since the wavelength is a frequency-dependent parameter. This concept is illustrated in Figure 4.

Figure 4 Adapting mesh size Δ in the finite-difference scheme of equation (14) for (a) a long wavelength and (b) a short wavelength. The optimal mesh size is the one that allows three points per wavelength [36].

Figure 5 Map of the 2D-Fourier transform (a) raw and (b) filtered, of the vertical displacement U2D(k, f) with respect to the wavenumber and the frequency done by using equation (15). Figure (a) corresponds to the input data to be filtered and Figure (b) to the final result obtained after filtering according to the procedure explained in Section 3.4.

Figure 6 Filtering flowchart with border-padding. Border-padding is achieved using zero-padding and a succession of filterings. To achieve this, the measured displacement field is re-injected at each cycle without modifying the edges, in order to gradually obtain border-padding solely from the measured data.

So, to estimate the structural parameter s for a single frequency, only five points in the measured displacement field are needed. However, a complete scan of the structure reduces measurement noise and uncertainties. Indeed, the use of the least-squares method on all spatial results has been performed to estimate the frequency dependence of spatially-averaged values of the structural parameter s. For a beam, the spatial least-squares calculation is straightforward. Assuming A = w/C where $\mathit{C}={\delta }_{\Delta }^{4x}$, the least-squares estimate of A is given by $A=\frac{{\mathit{C}}^H\mathit{w}}{{\mathit{C}}^H\mathit{C}}$, where each element of the column vectors w and C represents the displacement and bilaplacian estimator at one point, and where C^H is the Hermitian transpose (conjugate transpose) of C [36, 39].

Conversely, it is also possible to estimate a spatial mapping of the structural parameter s by averaging the results over a small frequency bandwidth. This bandwidth has to be sufficiently small to prevent significant variation of the structural parameters with respect to the frequency. Otherwise, the value of the average could not be representative of the behavior of the structure studied.

3.4 Filtering methodology

The results obtained by directly applying the methodology discussed so far show some discrepancies for certain frequency bands, mainly in high-frequencies, as in Figure 7, around f = 40 kHz. The main source of this discrepancy is the contribution of the other deformation modes that can be seen thanks to the displacement in k-space in Figure 5 obtained through 2D-Fourier transform (Eq. (15)).

$$U_{2\mathrm{D}}(k,f)={\int }_{t=-\mathrm{\infty }}^{t=+\mathrm{\infty }}\mathrm{}{\int }_{x=-\mathrm{\infty }}^{x=+\mathrm{\infty }}\mathrm{}w(x,t)e^{-\mathrm{ikx}}{e}^{\mathrm{i\omega }t}\mathrm{d}x\mathrm{d}t.$$(15)

Figure 7 (a) Equivalent Young’s modulus (b) Equivalent loss factor, with respect to the frequency with or without the application of the filtering methodology. This methodology is presented thanks to the workflow in Figure 6.

Figure 8 Samples A (a), B (b), and C (c), made with aluminum skins, with varying epoxy width along the beams.

Measurement noise and vibrations of the piezoelectric buzzer can also lead to other contributions. Even if CFAT is local and does not depend on lateral dimensions or boundary conditions, practically, there are contributions of other modes and of the buzzer excitation (assumption is done by considering the measurement window far from the excitation). Filtering methodology may be applied.

Therefore filtering the displacement field in wavenumber space is done. However, windowing prevents the displacement field from being rebuilt in a way that is representative of the original signal, particularly at low frequencies. To counter this, border-padding has been implemented to extend the spatial domain of measurement beyond what has been measured. The idea of the procedure comes from [40] and similar operations are performed. This technique allows the signal to be windowed spatially without altering the measured signal, as the windowing will be applied only to the added borders. The border padding of the filtering loop is generated from experimental values only. Since one loop is not enough to guarantee a continuous displacement field, the filtering and border padding processing is run several times. The loop is run until the border padding is converged. The steps are repeated until the extrapolation has converged. This methodology follows the workflow shown in Figure 6.

Finally, the border-padding and the measured signal are filtered, retaining only the bending wavenumbers. This yields the Figure 5b in k-space. Only the part of the signal corresponding to the measured area (without border-padding) is used to apply CFAT data processing. Applying the whole methodology leads to Figure 7, for a zero-padding of 80%, a window width of 23 rad.m⁻¹ and a filter of order 5. The real part is well filtered, extending the measurement range. Besides, for the imaginary part, several artifacts are removed from the data, but errors remain when the measured level is too low. This is particularly the case for high-frequency results, where the noise-sensitive loss factor becomes negligible. The results of the application of the filtering methodology are shown in the Figure 7.

4 Results

This section deals with the dynamic characterization of various sandwich beams. The aim is therefore to provide experimental results that reflect the dynamic behavior of a multilayered structure with imperfect interfaces. The dynamic parameters are identified using the CFAT method presented in Section 3. The mechanical and geometrical parameters of the structure are all assumed to be known, except for the Young’s modulus of the core. The latter is estimated a priori by taking the specimen which has the best quality of coupling. We consider this specimen as perfectly glued (so with perfect interfaces): B = 0 m.Pa⁻¹. We applied the model by considering continuity at the interfaces and we determine the core Young’s modulus by inversion. Then, for other specimens which have imperfect interfaces, the interface parameter B is characterized by using the model presented in Section 2.

In this section, the inverse method presented is applied to have two types of visualization: frequency visualization and spatial visualization, respectively, with results in two different subsections. These two types of visualization can play a complementary role, which is a very interesting feature of CFAT as opposed to modal analysis, for instance, which cannot directly combine these two approaches.

4.1 Frequency domain identification

As a reminder, the main equation of CFAT is equation (13) used for the characterization carried out in this section. This equation leads to many solutions, one for each given frequency and space coordinate. In this subsection, namely the “Frequency domain resolution”, the space solutions are used all together to find the optimal solution thanks to the least squares method. Therefore, the solution obtained is frequency dependent.

Figure 9 shows the results for three samples A, B and C that are presented in Figure 8 and Table 2. Photographs are obtained after hand made uncoupling/delamination of the samples to better see the coupling and interfaces (only two layers and one interface remains on the photos). Quality of interfaces (imperfect interfaces) are performed by varying epoxy width along the beams. In Figure 9, we can clearly see spatial evolution of the quality of the bonding rates with an increase of the quality by going to the right side.

Figure 9 Characterization of three samples: A, B, and C (made with aluminium skins). (a) Equivalent Young’s modulus and (b) Damping ratio with respect to the frequency. Comparison with the model is done here, by varying the value of the interface parameter B from perfect (B = 0 m.Pa−1) to imperfect (B ≠ 0 m.Pa−1).

Figure 10 Sample G1, a glass beam with a varying quality of the coupling.

Sample C is assumed to be perfect. As described previously Young’s modulus of the core is first extracted from the multilayer model with perfect interfaces and measurements on sample C assuming perfect bonding, and then the coupling parameter B is fitted for samples A and B assuming the same Young’s modulus for the core as the sample C. Therefore, the dynamical behavior computed from the displacement measurement is as expected by the model presented in Section 2, i.e. a shift towards the low frequencies of the dynamical parameters. The values of B correspond to in between situations. The case B = 0 corresponding to complete coupling is an assumption. The model is computed from the characteristics of the samples, given in Table 2. It can be seen in Figure 9 that naturally, as the value of B increases, the samples with the more imperfect interfaces are better represented. However, the imaginary part shows less accurate results than the real part for these samples. The aim of the characterization was mainly to match the upward slope between 1 kHz and 10 kHz, because beyond these frequencies the loss factor is highly variable. The discrepancy between the model and the experimental loss factor can be due to a number of parameters (temperature, age, layer alignment, boundary conditions, etc.) which, when added together, can have a significant effect on the overall behavior. The loss factor is particularly impacted, contrary to the equivalent stiffness, because the imaginary part is much more sensitive than the real part, more sensitive to measurement noise, but also to measurement conditions (such as temperature).

It is interesting to note that very low maximum loss factor could be estimated in this work, i.e. around 2%, with good accuracy for characterizing a material. It is relatively rare to obtain a variation in dynamic loss factor as a function of frequency at such low maximum values. In the bibliography, for example, we find measurements of maximum loss factor that are between 7% and 10% at least: 10% for bilayer plates in [41] and 7% for sandwich plates in [20].

4.2 Spatial domain identification

Leclère et al. showed in [36] that CFAT has the particularity of being able to account for the mechanical properties of a structure in the space domain. Equation (13) is usually solved as a function of frequency and as a function of space. The apparent bending stiffness is then derived as a function of space for each frequency (E(x, ω) for a beam, E(x, y, ω) for a plate). In practice, to limit measurement noise, this stiffness will be averaged over at least a short frequency band, allowing correct spatial visualization of the dynamical parameters. The value therefore represents the average behavior around the chosen central frequency. An example is shown in Figure 11 for another sample, called sample G1. This sample is similar to the A, B, and C samples but it is made from glass and is shown in Figure 10. The parameters of interest for this sample are given in Table 2. The two visualizations are shown simultaneously to show the width of the frequency band used for spatial visualization.

Figure 11 Comparison between the frequency domain approach and the spatial domain approach. G1 is the sample processed here.

Spatial domain identification is performed for the three sample A, B, and C and results are shown in Figure 12.

Figure 12 Equivalent Young’s modulus with respect to the position estimated thanks to CFAT for the three samples shown in Figure 8. These figures have been obtained on a frequency bandwidth from f = 10 kHz to f = 11.5 kHz.

Damping is not displayed for this visualization because of its high sensitivity to noise. Indeed, the high dependence of the dynamical behavior of the structure with respect to frequency, the low frequency resolution for this type of application, and the greater sensitivity of the loss factor to noise mean that the imaginary part is too widely scattered.

The results show stiffness varying as a function of the position along the length of the various beams. This is intuitive since the structures studied have a bond quality that increases linearly with respect to the position of measurement. Pictures of the samples were indeed processed to highlight the varying epoxy width along the beams and check if it fits the apparent stiffness trends. An example of this process is given in Figure 13 for sample A. The original image was turned to levels of gray only. Then, thresholds were added in order to obtain the value 0 (totally black) where glue is seen, and the value 1 (totally white) where there is no glue. As the glue thickness is constant all along the beam, the ratio $\frac{\mathrm{glue} \mathrm{width}}{\mathrm{total} \mathrm{width}}$ directly represents the amount of glue responsible of the coupling quality. Figure 13c shows this proportion graphically and provides an indicator for comparing the equivalent Young’s moduli with the quantity of adhesive.

Figure 13 (a) Sample A. (b) Processed picture with grey levels. (c) Adhesive width ratio.

This comparison is made on all the three samples in Figure 14, in order to check the reliability of this analysis. The purpose of this figure is just to show a general trend and not to fit the curves or show that they follow the same position dependence. That is why, to better see this trend the axis of adhesive ratio and Young’s moduli were slightly adjusted (y-scales). The results are very interesting for the three samples, as the overall trend is well described, but some spatial areas show less correlated results (between 0.075 m and 0.175 m for sample C, for example). The error may mainly be due to image processing, which could be improved by using edge detection, for example, to apply thresholds (either 0 or 1) to a geometric area directly. The discrepancy may also come from the averaged frequency region. Indeed, these results were obtained over a frequency range from f = 10 kHz to f = 11.5 kHz, and as the behavior of the multilayered structure is strongly dependent on the frequency, it is possible to have a global behavior that is not representative of the real dynamics. However, the need to obtain results implies averaging over a relatively wide spectrum, depending on the frequency resolution. In order to reduce this averaging bias, one solution would be to increase the frequency resolution and/or to perform more averaging during the measurement.

Figure 14 Equivalent Young’s modulus comparison with respect to the position and the adhesive width for (a) sample A, (b) sample B, and (c) sample C.

In order to confirm the effect of adhesive width on coupling, bevel-shaped samples were manufactured in a reproducible way. As with all the other samples, rulers were used to manage the adhesive width, but these were laid diagonally to produce the desired shape. These samples are shown in Figure 15 and the parameters of interest for those samples are given in Table 2. Sample “Bevel 1” has a glue width that varies linearly from 2.8 cm to 0.2 cm, sample “Bevel 2” has a width that varies linearly from 2.7 cm to 0.9 cm. These samples were manufactured together under the same conditions, and measured identically. The samples in Figure 15 were manufactured using the protocol developed in [42]. The results are shown in Figure 15, and have been obtained by averaging the dynamical bending stiffness from 10 kHz to 11.5 kHz. The trends induced by the different bonding widths are clearly visible. “Bevel 1” and “Bevel 2” samples have more or less the same maximum width, but the minimum adhesive width is different by a factor of 5, which has direct effects on coupling.

Figure 15 (a) Sample “Bevel 1”; (b) Sample “Bevel 2”; (c) Equivalent Young’s modulus with respect to the position for the two beveled samples. The sample “Bevel 1” width varies from 2.8 cm to 0.2 cm, the sample “Bevel 2” width varies from 2.7 cm to 0.9 cm.

It is interesting to note that for a glue ratio value (around 60%), the coupling impact seems to be stable. Indeed, the equivalent Young’s moduli for both “Bevel 1” and “Bevel 2” samples reach the same value if this glue ratio value is reached. This would indicate a limited sensitivity of skin coupling to the width of glue applied: for a threshold value of around 60% glue width, coupling no longer evolves significantly. In other words, for this threshold value, coupling seems to have converged towards perfect coupling. Indeed, increasing the proportion of glue no longer improves the coupling of the skins, which fits the definition of perfect coupling.

To check the sensitivity of the coupling plateau obtained at 60% glue width in Figure 15, the same calculations were repeated on Sample “Bevel 1” but for different frequency ranges. Indeed, since low frequencies represent well-coupled behavior, this plateau should be wider for low frequencies if that corresponds to a perfectly bonded Young’s modulus. These results are shown in Figure 16. Glue width thus seems to be a less and less influential parameter as the frequency values decrease. For the f = 10:11.5 kHz frequency band, the equivalent Young’s modulus stabilizes at 35 cm. Nevertheless, for the frequency band f = 2.3:2.5 kHz, the equivalent Young’s modulus already stabilizes at 25 cm. Finally, this means that long wavelengths (low frequencies) are less sensitive to the effect of glue width than short wavelengths (high frequencies). This seems natural, since the smaller the wavelength, the more sensitive the wave will be to spatial variations in mechanical parameters such as material heterogeneity, or coupling in this case. Moreover, since partial decoupling of the layers is a function of both the parameter B and the frequency f of the bending wave, it is natural that the variation in Young’s modulus, an indicator of coupling in the present work, should be affected differently according to coupling and frequency. In other words, this means that the best-coupled zone of the material (at the lowest B), will decouple less rapidly as a function of frequency f. The converse implies that the less well-coupled zone of the material (at the highest B), will decouple faster as a function of frequency f. To sum up, the high equivalent Young’s modulus plateau corresponds to a perfectly bonded plateau since it decreases slower with respect to the frequency, and also it gets wider when the frequencies are lower.

Figure 16 Equivalent Young’s modulus with respect to the position for the sample “Bevel 1” averaged for three different frequency bandwidths.

4.3 Combined frequency and spatial identification

Inverse methods in the frequency domain or the spatial domain were demonstrated in the two previous subsections. However, it is also valuable to use these two approaches in a complementary way. By doing so, experimental bias, such as that caused by temperature variations or sample aging, can be significantly reduced by comparing a sample to itself by dividing it into different areas.

First, a spatial computation is performed over a reasonably wide frequency band. Next, different regions can be identified if the sample does not exhibit a spatially homogeneous stiffness profile. Finally, these areas are separated, and the computation for frequency visualization is performed only within each region. Figure 17 illustrates the results of this method’s application.

Figure 17 Young’s modulus rigidity with (a) Spatial domain visualization divided in three regions on the frequency bandwidth f = 4500:5500 Hz, (b) Frequency visualization for each region shown in (a). The sample “Bevel 2” has been used to obtain these results.

The results in the space domain are still shown, but together with its division into three equal areas. The use of these three parts for the frequency computation is shown in Figure 17b. The reference (B = 0 m.Pa⁻¹) is chosen from within this structure. The spatial variation of the interface parameter B is thus highlighted thanks to these results. In fact, freed from numerous experimental biases (different samples, different measurement conditions, for example), the only parameter that will impact the dynamics of this multilayer is its interface.

The same type of analysis is carried out on a similar three-layered beam, but with a thin foil of cellophane at one aluminum-epoxy interface. It is spatially centered along the length of the beam. A schematic top-view of this type of sample is shown in Figure 18.

Figure 18 Scheme of a three layered aluminium-epoxy-aluminium, with a thin cellophane layer introduced between a layer of aluminium and the epoxy core. The cellophane layer is as large as the beam, but only 20 cm long (a third of the beam).

Cellophane foil was chosen to decouple the aluminum from the epoxy glue, as it does not adhere to the epoxy glue and thus enables decoupling. Another advantage of this material is that its thickness is negligible compared to the beam layers. Results for this structure are shown in Figure 19. The sample shows three typical areas, in the following order: partially bonded, fully debonded, partially bonded. In the first and third area the coupling is not perfect but “partially” bonded with a value of B not equal to 0 (0 is for perfect coupling). The first and last thirds have a relatively high apparent bending stiffness. They have the same amount of bonding and shows hereby similar spatial and frequency results. The area in the middle, which is the zone with the cellophane material, is much softer; that is a consequence of the almost total decoupling of the layers. Indeed, the value shown in the spatial visualization for this area (around 15 GPa) is the value obtained over the entire spectrum for its frequency calculation. This value is the case where B → ∞ m.Pa⁻¹. This shows that decoupling has already been reached for all measured frequencies. The low-frequency results for the 1st third and the 3rd third diverge due to the lack of length used in the wavenumber estimation. Indeed, since the length used to apply CFAT is three times smaller due to the division of the computed surface, the size of the area used in the wavenumber estimation is not enough in low-frequency because the wavelengths to be estimated are too high for the well-computation of the finite-difference scheme.

Figure 19 Young’s modulus rigidity with (a) Spatial domain visualization divided in three regions on the frequency bandwidth f = 10000:11500 Hz, (b) Frequency visualization for each region shown in (a). The sample with cellophane has been used to obtain these results.

In the case of the multilayered sample presented here, the significant spatial discontinuity of the bonding accounts for very different local dynamics. This type of behavior cannot be modeled globally by the current model. In other words, it is not possible to obtain an equivalent Young’s modulus value as a function of frequency alone, as this would break the model’s assumption of spatial homogeneity. Originally, this was also the case for CFAT. Indeed, these strong discontinuities make frequency results uncertain if CFAT is used on the whole sample. The global dynamics oscillate between a fully debonded behavior and a partially bonded behavior with respect to frequency. However, CFAT’s spatial visualization solves these problems as shown here, making the analysis of such materials straightforward.

5 Conclusion

In this work, the characterization of the dynamics of imperfect interfaces has been demonstrated. The model used for this characterization was initially introduced in Section 2. It is an equivalent dynamic model that incorporates imperfect interfaces through a single parameter, the interface parameter B. This parameter represents the bonding conditions within the multilayered structure. The impact of partially debonded layers results in a softening of the apparent flexural rigidity of the entire structure under vibrations. The behavior of imperfect interfaces was experimentally demonstrated in this work and characterized using the interface parameter B presented in the model.

The experimental methodology employed in this work is detailed in Section 3. It enables the dynamic characterization of multilayered beam-like structures through displacement measurements using the LDV technique and the CFAT. The filtering methodology used to extend the results into the high-frequency range is also described in this section.

Finally, the results of the dynamic characterization of imperfect interface samples are presented in the Section 4. Three approaches are outlined: in the frequency domain, in the spatial domain, and a complementary approach combining both. This enabled the characterization of the interface parameter B with respect to both frequency and space. It also allowed for the mechanical characterization of highly discontinuous samples. The model is based on a global B parameter of interface coupling; once this value is adjusted, the spatial and frequency evolution of the measured and predicted Young’s modulus are in good agreement along the frequency spectrum. However, the results for the imaginary part could still be improved. If not through data processing, greater accuracy could be achieved by using samples with a higher loss factor. This could be accomplished by selecting another core material or by significantly increasing its thickness.

Acknowledgments

The authors thank the two anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions. The authors would also like to express their gratitude to Fabien Chevillotte, François-Xavier Bécot and Fabien Marchetti (engineers at Matelys Research Lab) for fruitful discussions and scientific advices on several experimental choices done in this work. We also thank Quentin Leclère (professor at LVA), Julien Chatard (technician at LVA) and Teddy Sacco (student at INSA, GM department) for their precious help during the set up of the experimental protocol.

Funding

This work has been done during the PhD of the first author at the LVA and LTDS financed 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).

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

Data are available on request from the authors.

References

All Tables

All Figures

	Figure 1 Sandwich panel with imperfect interfaces (represented by springs to indicate that there is a coupling law between layers). According to this representation, only relative motion in x and y directions is accounted for.
In the text

	Figure 2 Flexural rigidity with respect to the frequency for two different interface conditions. For the blue curve, the simulation has been done with B = 0 m.Pa−1, which is a perfect interface structure. The red dotted line has been done with the same material but with imperfect interfaces with B = 10e−9 m.Pa−1 = 1e−8 m.Pa−1. The material properties for this example are gathered in Table 1.
In the text

	Figure 3 Picture of the setup consisting of: the measured sample, the vibrating device and the laser vibrometer. The picture is shown next to a diagram of the sample from the laser vibrometer point of view, with a representation of the scanning measurement mesh.
In the text

	Figure 4 Adapting mesh size Δ in the finite-difference scheme of equation (14) for (a) a long wavelength and (b) a short wavelength. The optimal mesh size is the one that allows three points per wavelength [36].
In the text

	Figure 5 Map of the 2D-Fourier transform (a) raw and (b) filtered, of the vertical displacement U2D(k, f) with respect to the wavenumber and the frequency done by using equation (15). Figure (a) corresponds to the input data to be filtered and Figure (b) to the final result obtained after filtering according to the procedure explained in Section 3.4.
In the text

	Figure 6 Filtering flowchart with border-padding. Border-padding is achieved using zero-padding and a succession of filterings. To achieve this, the measured displacement field is re-injected at each cycle without modifying the edges, in order to gradually obtain border-padding solely from the measured data.
In the text

	Figure 7 (a) Equivalent Young’s modulus (b) Equivalent loss factor, with respect to the frequency with or without the application of the filtering methodology. This methodology is presented thanks to the workflow in Figure 6.
In the text

	Figure 8 Samples A (a), B (b), and C (c), made with aluminum skins, with varying epoxy width along the beams.
In the text

	Figure 9 Characterization of three samples: A, B, and C (made with aluminium skins). (a) Equivalent Young’s modulus and (b) Damping ratio with respect to the frequency. Comparison with the model is done here, by varying the value of the interface parameter B from perfect (B = 0 m.Pa−1) to imperfect (B ≠ 0 m.Pa−1).
In the text

	Figure 10 Sample G1, a glass beam with a varying quality of the coupling.
In the text

	Figure 11 Comparison between the frequency domain approach and the spatial domain approach. G1 is the sample processed here.
In the text

	Figure 12 Equivalent Young’s modulus with respect to the position estimated thanks to CFAT for the three samples shown in Figure 8. These figures have been obtained on a frequency bandwidth from f = 10 kHz to f = 11.5 kHz.
In the text

	Figure 13 (a) Sample A. (b) Processed picture with grey levels. (c) Adhesive width ratio.
In the text

	Figure 14 Equivalent Young’s modulus comparison with respect to the position and the adhesive width for (a) sample A, (b) sample B, and (c) sample C.
In the text

	Figure 15 (a) Sample “Bevel 1”; (b) Sample “Bevel 2”; (c) Equivalent Young’s modulus with respect to the position for the two beveled samples. The sample “Bevel 1” width varies from 2.8 cm to 0.2 cm, the sample “Bevel 2” width varies from 2.7 cm to 0.9 cm.
In the text

	Figure 16 Equivalent Young’s modulus with respect to the position for the sample “Bevel 1” averaged for three different frequency bandwidths.
In the text

	Figure 17 Young’s modulus rigidity with (a) Spatial domain visualization divided in three regions on the frequency bandwidth f = 4500:5500 Hz, (b) Frequency visualization for each region shown in (a). The sample “Bevel 2” has been used to obtain these results.
In the text

	Figure 18 Scheme of a three layered aluminium-epoxy-aluminium, with a thin cellophane layer introduced between a layer of aluminium and the epoxy core. The cellophane layer is as large as the beam, but only 20 cm long (a third of the beam).
In the text

	Figure 19 Young’s modulus rigidity with (a) Spatial domain visualization divided in three regions on the frequency bandwidth f = 10000:11500 Hz, (b) Frequency visualization for each region shown in (a). The sample with cellophane has been used to obtain these results.
In the text