Experimental identification of flexural and shear complex moduli by inverting the Timoshenko beam problem

Abstract This paper addresses the problem of estimating the local viscoelastic parameters of sandwich beams. An original procedure involving an inverse vibratory method (Force Analysis Technique) and the Timoshenko beam theory is detailed and applied experimentally on a sample presenting a honeycomb core. The major philosophy relies in considering multi-layer beams as equivalent homogeneous structures. This simplified approach is thought to be more representative of the global dynamic behaviour, in addition the reduction of degrees of freedom is obviously an improvement for modelling on Finite Element software. When compared to other usual approaches, the method developed in this paper shows a very good agreement between the experimental sandwich beam and the homogeneous model, which highlights interesting insights for applying it to industrial structures. The local aspect, the robustness and the self-regularization properties are verified on a wide frequency range, making the procedure possibly efficient for characterization of structures on a production line, flaw detection and Structural Health Monitoring.


Introduction
In order to reduce energy consumptions, transportation industries are increasingly focusing on new structures to diminish the mass of their vehicles. Thus composite materials appear to be ideal since they present a high stiffness and a low weight. Nonetheless the complexity in estimating and representing their dynamic behaviour relegates their use to secondary structures (car's bonnet, plane's tail, inner parts) which are subject to weaker disturbances. Although standard procedures are reliable for identifying elastic properties of common materials and alloys, their application on composite structures can be struggling.
Firstly measurement approaches involving modal analysis theory are widely used in engineering and in the literature. A verification of the theoretical resonance frequencies allows to identify the Young's modulus and the structural losses by applying the general half-power bandwidth method [1]. One can also mention standardized tests for estimating the elastic properties of multi layer composite beams, especially the shear modulus of the core, using the Oberst beam method [2]. Nonetheless these modal procedures imply specific boundary conditions that can be tricky to achieve, such as the cantilever beam. Montalvao et al. [3] demonstrated an important variability over the parameter estimation depending on the clamping condition.
In addition a potential fibre breakage or core flattening can happen when applying such condition. To avoid these issues a hard piece of metal is usually inserted in the core. Eventually, the so-called "base beam" must present a stable stiffness versus frequency, which is not always the case for composite material, needing some changes in the measurement procedure as proposed by Liao et al. [4].
Secondly the finite element method approaches provide realistic simulations and can couple various dynamic phenomena. A comparison between the predictions of a theoretical model and measurement on a real structure may lead to the estimation of elastic parameters [5,6]. Unfortunately the results strongly depend on user's decisions (assumptions, mesh refinement used for the finite element model, etc). One can also underline the difficulties to characterize but also to model some structures. For instance a bonding between two layers can be harsh to determine and even more to define accurately. The same trend can be extended to boundary conditions, coupling, and more generally to the trustworthiness of a model. Hence the misdiagnosis on a given structure can come from the measurement step and/or from the modelling stage.
Third some analytical approaches have been studied to represent the dynamic behaviour of composite structures. Amongst the numerous theories developed, always representative for one specific kind of sandwich organization, the zig-zag theory [7] appears as the most reliable method since it allows a global kinematic description of the multilayer material. Using a matrix transfer function, the vector of degrees of freedom (DOFs) of the bottom layer is projected onto the one above, whose computed vector of DOFs is applied to the upper layer, and so on until reaching the top layer whose vector of DOFs will represent the dynamic behaviour of the entire structure as an equivalent homogeneous material. Guyader et al. applied this approach to estimate equivalent elastic properties of homogenized multi-layer composites [8,9] and proved the reliability of this homogenization approach. Once more, it is of importance to perfectly know the mechanical properties of each layer to apply this analytical method.
Instead of analysing each layer independently (first plausible source of bias) in order to model a realistic multi-layer structure in a finite element simulation or analytical formulation (second possible source of bias), an original experimental approach is suggested in this paper using an homogenization of the whole composite. Such consideration of the sandwich beams has also been studied by Rak et al. [10] with the Euler-Bernoulli theory using a wave approach, and similar work has been conducted on plates by Ichchou et al. [11]. The proposed procedure herein involves the Timoshenko beam theory combined with the inverse vibratory method called the Force Analysis Technique.
To the authors' knowledge the use of Timoshenko beam theory to represent sandwich beams has not been an extended subject of research so far. Labuschagne et al. [12] compared the Euler-Bernoulli and the Timoshenko models, Lim et al. [13] developed an improved Euler-Bernoulli formulation to take shear effects into account. Dixit [14] focused on the differences between these beam theories to identify defects from the analysis of modal frequency shifts, while Mei et al. [15] expanded analytical expressions to represent vibrations of cracked beams using a wave approach and the Timoshenko theory. Concerning the vibrations of sandwich beams, the important shear and the associated damping effects induced by the core layer are often highlighted [16,17]. Besides the lack of experimental data regarding the Poisson's ratio (directly linked to the shear modulus) for composite structures may influence the analysis of their stiffness properties [18].
A method allowing to simultaneously characterize the Young's modulus and the shear modulus using few assumptions and easily reproducible measurement step is a requirement commonly aroused by engineering and industries.
The second originality of the proposed formalism lies in the use of an inverse vibratory method, the Force Analysis Technique (FAT), initially developed by Pézerat and Guyader to identify vibratory sources [19,20,21]. Based on a local verification of the equation of motion, this method allows analyses of various structures such as beams, plates, shells [22], and more complex structures using the finite element method [23,24]. Its application is valid for a wide frequency range and does not depend on the boundary conditions because of its local feature. Any unknown or undesired phenomena which is out of the measurement area can be simply discarded. For theses reasons numerous developments arose from the Force Analysis Technique such as the suppression of the filtering stage using a smart measurement mesh [25,26], detection of defects [27,28,29], localization of acoustic sources [30], or identification of acoustic part in turbulent wall pressures [31]. One of these developments focused on the characterization of the stiffness and the loss factor on plates as suggested by Ablitzer et al. [32,26] and showed promising results as a novel method to analyse local properties of structures. The work presented in this paper consists in an expansion of this approach to multi-layer composite beams and aims at characterizing the Young's modulus, the shear modulus and their respective damping simultaneously on a wide frequency range by only measuring the transverse displacement of the structure.
The paper is structured as follows. First, the general principle of the method is presented. Afterwards some simulation cases are considered in order to prove the feasibility of the procedure. The different steps to de-noise the data and stabilize the inverse problem are given in details. Finally an experimental validation is exposed to confirm the potential of applicability of the method.

Theoretical approach
Let us consider the harmonic vibration of a beam. Out of the exciting source area, the equation of motion for an Euler-Bernoulli beam and a Timoshenko beam respectively take the form : andẼ where W is the transverse displacement,Ẽ the complex Young modulus (Ẽ = E(1 + jη E )), η E its associated damping coefficient,G the complex shear modulus (G = G(1 + jη G )), η G the shear loss factor, S the cross section of the beam, I the second moment of area, ρ the density and ω the angular frequency.
The identification of the complex moduli lies on the knowledge of W (x), as seen in Eq.
(1) and (2). In practice the transverse displacement W (x) is simply estimated by using an accelerometer or a scanning laser vibrometer but the spatial derivatives cannot be assessed with a straightforward measurement.
This issue is bypassed through a fourth order finite difference scheme leading to the approximated spatial derivatives δ 2x (≈ ∂ 2 W (x) ∂x 2 ) and δ 4x (≈ ∂ 4 W (x) ∂x 4 ) of which expressions are detailed in Appendix A. When considering a matrix representation the previous equations of motion take the form : for the Euler-Bernoulli model and for the Timoshenko one.
One may notice that the complex elastic parametersẼ andG are isolated in the last expressions, allowing to find a unique value of E, η E , G and η G by mean of a matrix division in the least square sense as : . . .
The unknown parameterG is still present in the right hand side of Eq. (6) in which ρ 2 Ĩ G ω 4 expresses the rotatory inertia. An iterative procedure is therefore applied as presented below. The first iteration consists in removing the rotatory inertia which is neglected in first approximation, leading to estimation ofẼ 1 and G 1 . The latter is inserted in the second iteration in which the right hand side is complete, allowing a more precise estimation ofẼ 2 andG 2 , used for computation ofẼ 3 andG 3 , and so on until reaching a variation of G smaller than 0.1% between two iterations, which is generally achieved after 4 loops.

Proof-of-concept
This part will deal with simulations to verify the validity of the suggested method. First an analytical expression of the displacement field using the Timoshenko beam theory is developed. Then, the inverse problem is tested using exact data or noisy data. The latter implies a regularization procedure, but also specific post treatments for the Timoshenko problem that are explained in a last section.

Direct problem
Before applying the inverse approach detailed previously, it is necessary to define a convenient displacement field according to the Timoshenko beam theory. An analytical expression of the displacement field W is computed through a forced wave decomposition [33], for a simply supported beam of length L with F the amplitude of the source applied to the structure, k p the propagative wavenumber and k e the evanescent wavenumber defined as follows : The complex displacement field is computed on a thick aluminium beam with the parameters resumed in Table 1. The length to thickness ratio is about 33, which is sufficient to imply a difference between the Euler-Bernoulli and Timoshenko models because of the presence of shear effects, as it will be shown later.
The moduli of displacements are plotted in Fig. 1.

Parameters Values
Length

Inverse problem using exact displacements
The exact displacement and its discrete spatial derivatives estimated on the studied area x ∈ [0.1, 1] m, i.e out of the source region, are put into Eq. (5) and (6)  The aluminium beam is seen as a softer and more damped structure as frequency increases.

Inverse problem using noisy displacements
The displacement field computed in the last section is now blurred with some noise contribution as defined by Pereira et al. [34], using a Signal to Noise Ratio (SNR) of 40 dB : where W noisy is the noisy transverse displacement at a given frequency, W exact is the exact transverse displacement at a given frequency, α and γ are zero mean Gaussian random variables with unit variance, θ and φ are random variables uniformly distributed in [0, 2π], N s is the number of samples of the signal.

Regularization procedure
The presence of noise is a current drawback in inverse problems unavoidably implying a regularization step [35]. Most of them rely on a Singular Value Decomposition to properly inverse the matrix of data [24] or on probabilistic approaches [36]. Even if the approach suggested in this paper does not involve a matrix inversion, the extreme amplification of noise is still present in the terms δ 2x and δ 4x because a weak error in the displacement field implies important discrepancies in the estimation of the spatial derivatives. For these reasons a simple procedure using a windowing (Hanning window) and a low-pass filter (with cutoff wavenumber k c ) is applied on the noisy data as initially proposed by Pézerat et al. [19] for the Force Analysis Technique. More details about the regularization step can be found in Appendix B.

Simulation using noisy data with regularization
For this second simulation case the noisy displacement field is used to compute the spatial derivatives δ 2x and δ 4x . The data are treated through the windowing and low pass filtering procedures with a cutoff wavenumber k c = k p . Eq. (5) and (6) are then solved to identify the elastic parameters plotted in Fig. 3.
Apart from the presence of small variations on the damping, due to residual noise, the results obtained with the Euler-Bernoulli formulation ( Fig. 3(a) and (d)) remain the same than the ones computed from the exact case ( Fig. 2(a) and (d)). The filtering procedure on noisy data is efficient and allows the identification of the elastic parameters. On the contrary the Timoshenko model fails to identify the parameters despite the regularisation step (with exception for the Young's modulus ( Fig. 3(b)) at some antiresonance frequencies).
The Timoshenko equations allow to compute more viscoelastic parameters than the Euler-Bernoulli equations, and could explain the presence of a more important bias on each viscoelastic parameter. But in this case the behaviour of the estimated complex moduli is almost random as shown by Fig. 3. This instability of the Timoshenko equations, when using regularized data at a single frequency, is explained in the next section.

Phenomenon explanation
To make the instability easier to represent, let us consider the previous displacement field with only real variables, i.e. the damping parts are ignored (η E = η G = 0). When expressing the Eq. (4) as represented in Eq. (14) one may interpret the Timoshenko's equation of motion as the equation of a plane in the reference frame (W, δ 4x , δ 2x ), with A and B corresponding to the slopes depending on E and G.
To illustrate this, the exact data computed at f = 1530 Hz for W , δ 4x and δ 2x is plotted in Fig. 4(a).
The information gathered along the spatial positions on the beam forms a curved line, as shown by the side projections, which belongs to the equation of a plane as expected from Eq. (14). Nonetheless the regularization procedure vanishes this trend as presented on Fig. 4 The deterioration induced on the data by the regularizing step appears mostly on δ 2x . To bypass this difficulty it is necessary to use additional data belonging to the plane, which is achieved using a multifrequency analysis.

Multi-frequency analysis
A convenient solution to get stability back for the Timoshenko model is to estimate the parameters from The identification of the elastic parameters E, η E , η G and G is now carried out using the multi-frequency analysis and presented in Fig. 6 Fig. 3(a) and (d)) or a multi-frequency analysis ( Fig. 6(a) and (d)). Indeed as seen in Fig. 5 the rotation of W and   δ 4x for f ± ∆ f is equally distributed around the central frequency f . Hence the least square estimation of the complex Young's modulus is neutral regarding the side frequencies. Third the estimation of the elastic parameters is highly improved for the Timoshenko model. The Young's modulus E (Fig. 6(b)) is fitting the one defined in the direct problem, the estimated shear modulus G (Fig. 6(c)) and elastic damping η E (Fig.   6(e)) are spread around the expected values. The identification of the shear loss factor η G (Fig. 6(f)

Width criterion ∆ f
In order to stabilize the Timoshenko model without preventing frequency variations of the elastic parameters too much, a compromise has to be done on the width criterion ∆ f . This part only concerns the Timoshenko model since the Euler-Bernoulli one does not necessarily need the multi-frequency analysis.
The elastic parameters are now computed for every applicable multi-frequency analyses For instance at f = 1800 Hz the multi-frequency can be achieved with ∆ f ≤ 800 Hz otherwise unavailable frequency samples at f ≤ 1000 Hz would be needed. Similarly at f = 3400 Hz we must have ∆ f ≤ 600 Hz so the multi-frequency analysis is limited to the available frequency samples, i.e f ≤ 4000 Hz.
The identified elastic parameters are represented in Fig. 7. One can remark a rapid convergence with respect to ∆ f on the Young's modulus E (Fig. 7(a)) and a slower convergence for its loss factor η E (Fig.   7(c)). The shear modulus G ( Fig. 7(b)) is a lot more disturbed for small values of ∆ f , and the damping η G (Fig. 7(d)) is almost random, for the same reasons as explained before, and will not be used for the choice of the width criterion. The computed data for G presents an important drop for ∆ f < 100 Hz and an acceptable stability for ∆ f > 600 Hz. Nonetheless the smoothing effect on the elastic parameters towards frequency implied by such a width criterion and the important reduction of the frequency bandwidth are not satisfactory. For these reasons an objective function is needed. As explained previously the higher ∆ f the more the stability and the less the number of analysed frequencies. Since our goal is to reach coherent elastic parameters on the widest frequency band the objective function has to depend on the dispersion of the shear modulus G (the purpose to use the Timoshenko model) assessed using an Euclidean norm along frequency, at a given ∆ f , and on the number of output frequencies N f (∆ f ) (i.e the bandwidth of the estimated complex moduli after applying the multi-frequency analysis), as : An illustration of the computed objective function is available in Fig. 8. A transient zone can be seen and is followed by a slow increase that finally diverges. The best width criterion is chosen to be the one that minimizes the objective function apart from the transient zone. For this case the optimal width criterion is found to be ∆ f = 530 Hz. By using a multi-frequency analysis with such width criterion one can obtain the results that were presented in Fig. 6.

Optimized self-regularization
For a given frequency f a structural wavenumber k p is associated. In the Force Analysis Technique, a common way to adjust the low-pass filter is to choose a cutoff wavenumber proportional to the natural wavenumber of the structure. From past experience it has been noticed that for the Timoshenko model the optimal cutoff wavenumber k opt c must be tuned around the natural wavenumber such as k opt c = k p [37]. The application of a low-pass filter to the displacement field of an unknown structure can appear to be tricky. The cutoff wavenumber has to be tuned with respect to the natural wavenumber of the structure, which is itself unknown since it depends on the complex moduli to identify. To overcome the lack of a priori knowledge on the natural wavenumber the following procedure is applied.
Our goal is to find an optimal cutoff wavenumber k opt c such that k opt c = k p . To do so the noisy data is filtered with different candidate cutoff wavenumbers k c and the elastic parameters are identified as explained previously. Then for each candidate wavenumber the identified complex moduli are used to obtain an estimation of the natural wavenumber k p as given by Eq. (11). Finally the values of k p are compared to the ones of k c . Their coincidence means that the cutoff wavenumber is equivalent to the one estimated on the filtered data, indicating optimal regularization of k opt c = k c = k p . This case is considered as the best adjustment of filter and the associated elastic parameters are kept for this frequency f . An illustration of this procedure is available in Figs. 9 and 10 for the elastic parameters and the estimated wavenumber respectively, with the previous aluminium thick beam at f = 1990Hz with added noise. About 400 different candidate values of cutoff wavenumber between 12 and 48 rad/m are tested. It can be seen that the Euler-Bernoulli model does not really depend on the filter's severity ( Fig. 9(a) and (d)). It actually allows a higher filtering order depending on the noise level as already shown by Pézerat [38]. On the contrary the Timoshenko model suffers important variations especially for G and η G (Fig. 9(c) and (f)). But when the estimated complex moduli are used to compute k p the result is astonishingly constant, no matter the cutoff wavenumber, as shown in Fig. 10. This trend may be explained by the equilibrium between the elastic parameters : a deviation of E can be counterbalanced by a more important deviation of G. One can also notice that k p and k c are crossing each other at a location marked by a green dot. This specific coincidence actually depicts the optimal filtering coefficient as confirmed by the elastic parameter identification. Back to Fig. 9 it is clear that this situation is the only one for which all parameters are well-estimated with the Timoshenko model. Moreover this stable region is very tight and corroborates with the importance of an efficient filtering. The optimal cutoff wavenumber k opt c computed at each frequency is presented in Fig. 11 and compared to the propagative wavenumber k p defined in the direct problem. The high coherence between them confirms the robustness in the estimation of the optimal cutoff wavenumber k opt c .

Measurement set-up
The beam is mounted on a Bruël & Kjaer LDS V201 shaker excited by a periodic chirp in the frequency range [200, 4000] Hz at position x 0 = 0.048 m (Fig. 12). The displacement field is acquired on a mesh grid of 7 × 157 = 1099 points (spatial step ∆ x = ∆ y = 4.1 mm) using a PSV-500 scanning vibrometer and its respective acquisition software.

Identified elastic parameters
The displacement field W is kept only out of the vibratory sources (space domain for which x > 0.1 m) and a mean is performed over the y direction to discard the torsional modes. The spatial derivatives δ 2x and δ 4x are computed by mean of the finite difference scheme. The regularization procedure (windowing and filtering) is applied to these variables with various cutoff wavenumbers k c and are then put into Eq. (7) to (9).
The elastic parameters are computed for every combination of (f ,∆ f ,k c ). The optimized cutoff wavenumber k opt c (section 3.5) is obtained as explained previously and depicted in Fig. 13. It does not depend on the width criterion ∆ f and its linear behaviour, confirmed in Fig. 14, is typical of beams presenting an  Fig. 16. The Young's modulus estimated with the Euler-Bernoulli model decreases rapidly from 6 GPa to 1 GPa (Fig. 16(a)) while the loss factor η E grows up from 1% to 4% as frequency increases ( Fig. 16(d)). Regarding the Timoshenko model the Young's modulus varies around 7.5 GPa and finally decreases for f > 3000 Hz (Fig. 16(b)) but its respective loss factor η E cannot be estimated (Fig. 16(e)). It may be implied by a too small dissipation associated to bending, its estimation could be impossible because of measurement noise. On the other hand the shear modulus is found to be constant around 45 MPa in the whole frequency range (Fig. 16(c)) and its associated loss factor η G roughly grows up from 1% to 4% when frequency increases (Fig. 16(f)). One can also distinguish a singularity at f = 1390Hz for the shear modulus and the Young's modulus, due to a coupling between a torsional mode and a flexural mode, a case that is not taken into account in the considered beam theories ; the optimal cut-off wavenumber k opt c can't be well estimated (see Fig. 14) nor the elastic parameters.
The general trend of the elastic parameters leads to three major observations. Third the Young's modulus estimated with the Timoshenko problem tends to decrease for high frequencies although the shear modulus keeps a stable behaviour. This trend may be associated to potential Lamb waves that would imply an important dynamic in the core and a weak contribution of the glass fibre faces, accompanied with a flexural stiffness decrease.

Comparison with static tests
The elastics parameters estimated with the Timoshenko beam model are now confronted with classical methods implying a three points bending test.
The sandwich beam is placed on an INSTRON 8801 three points flexural bending test bench as the one presented in Fig. 17

Synthesis of displacement field with homogenized beam
To confirm the validity of the suggested method, the identified elastic parameters are used to compute the displacement field of a beam with free boundary conditions. The elementary stiffness and mass matrices are taken from the works of Zienkiewicz [42]. The frequency responses obtained from measurements and from the finite element model using the homogeneous elastic parameters identified from Fig. 16 are presented in Fig. 19. The Euler-Bernoulli model doesn't fit the experimental frequency response because of its simplicity. On the other hand one can remark a very satisfactory agreement between the measurement on the real structure and the reconstructed displacement of the equivalent homogeneous Timoshenko beam. The matching is almost perfect from 500 Hz to 2700 Hz, except around 1390 Hz because of singularities observed on the estimated elastic parameters.
Finally the divergence that appears slowly for f > 3000 Hz could be explained by a lack of convergence of the finite element model for high frequencies (even if the mesh refinement is at least 40 elements per wavelength) and/or by the presence of the accelerometer for the experimental case which could dissipate energy as frequency increases.
Still the general concept of modelling multi-layer structures as homogeneous one with equivalent elastic parameters appears as a smart alternative to static approaches, modal theory or FEM models. It reduces complexity in describing the dynamic behaviour of these materials, but also decreases the number of DOFs when implementing numerical models. Moreover the experimental set up is clearly simplified. The use of an inverse problem suppresses the consideration of boundary conditions and the homogenization approach removes the need of characterizing each layer separately, the sandwich material is considered as a whole.

Conclusion
The analysis of vibratory behaviour of composite structures is a struggling problem with classical approaches such as finite element methods or modal testing. Achieving specific boundary conditions or modelling interactions between the different layers may be tricky.
This study dealt with characterization of composite beams using an inverse vibratory method based on the local verification of the equation of motion applied to the Timoshenko beam. This novel approach con-  Since the method relies on a local identification, it does not depend on the boundary conditions, its application is simple and allows developments for structural health monitoring or implemented analysis on a production line. As perspective a development on plates could be efficient too. Finally an extension to flaw detection is of interest since the kind of defect may be identified through the cartographies of elastic parameters. For instance a fibre breakage could change only the Young's modulus, a delamination might mostly influence the shear modulus, while a Barely Visible Impact Damage would increase local dissipation and so the damping coefficients.

Appendix A. Finite differences approximation
The spatial derivatives involved in the equation of motion are estimated using a fourth order finite difference scheme and are defined as follows : (A.1) where ∆ x is the spatial step between the discrete measurement points on the structure, and W the displacement field. This approximation implies the loss of 3 points on each side of the beam, the computation of the inverse problem is then applied on a mesh of N x − 6 points. Secondly a low-pass filter, described by its spatial impulse response H (Eq. (B.2)), is applied to the windowed noisy data : H(x − where β x = f x λ c and f x defines the form factor of the impulse response (f x ∈ N). The higher f x the wider H(x). In this paper the regularizing procedure is always applied with f x = 1. The Fig. B.21 depicts the typical impulse response obtained for f x = 1.
Lastly the filtered data can be mathematically described by the following equation :