Model-Assisted Guided-Wave-Based Approach for Disbond Detection and Size Estimation in Honeycomb Sandwich Composites

One of the axioms of structural health monitoring states that the severity of damage assessment can only be done in a learning mode under the supervision of an expert. Therefore, a numerical analysis was conducted to gain knowledge regarding the influence of the damage size on the propagation of elastic waves in a honeycomb sandwich composite panel. Core-skin debonding was considered as damage. For this purpose, a panel was modelled taking into account the real geometry of the honeycomb core using the time-domain spectral element method and two-dimensional elements. The presented model was compared with the homogenized model of the honeycomb core and validated in the experimental investigation. The result of the parametric study is a function of the influence of damage on the amplitude and energy of propagating waves.


Introduction
Honeycomb Sandwich Composites (HSCs) are a type of multi-layered structure that are composed of the mid-core with the geometry of honeycomb sandwiched between thin skins. They are widely used in the aerospace, marine and automotive industries due to the high strength-to-weight ratio, high energy absorption capability and effective acoustic insulation. However, these complex structures are exposed to various types of damage that are not found in metal alloy materials, e.g., hidden disbonds between the skin and the core, delamination of the skin plates, or the core impact damage. They can occur either during a manufacturing process, storage or in-service life. Therefore, advanced methods are required for on-line damage detection.
The Guided Waves propagation method is a high-potential approach in SHM for damage detection in HSCs [1][2][3][4][5]. GW are mechanical waves that propagate in a bounded elastic medium, e.g., bars, beams, rods, plates and shells. An excitation and sensing of the GW can be realised by the lightweight and inexpensive piezoelectric transducers (PZT) [6]. The compact PZT can be surface-bonded to the inspected structure or even embedded between the composite plies so that the measurements can be conducted in situ.
Among numerous GW-based techniques developed for damage detection and localisation, the most popular are pitch-catch [7,8], pulse-echo [9,10], phase array [11,12] and time-reversal mirror [13,14]. For damage identification, some of them require a baseline to be determined. Due to the costs and time-consumption, experimental investigation is an inefficient approach to obtain references. Lonkar et al. presented new model-assisted diagnostics for SHM [15]. The numerical model was used to determine the exact velocity of the wave propagating in the stiffened panel, which is essential for accurate damage identification. A combination of 3D scanning laser vibrometry measurements and the numerical model to reconstruct or update baseline signals for damage detection with guided waves was proposed by Aryan et al. [16].

The Spectral Element Method
The general concept of the SEM is based on the idea of the FEM. The similarity of both methods lies in the fact that the modelled domain is divided into non-overlapping finite elements, and external forces and arbitrary boundary conditions are imposed in the particular nodes. The main difference between those methods is a choice of the shape function N = N(ξ), which is interpolated by a Lagrange polynomial that passes through the element nodes. The nodes are localized on the endpoint of an interval, ξ ∈ [−1, 1], and the roots of the first derivative of Legendre polynomial P of degree p − 1: The approximation of an integral over the elements is achieved according to Gauss-Lobatto-Legendre (GLL) rule at points coinciding with the element nodes, and the weights w = w(ξ) calculated as: This approach guarantees a diagonal mass matrix. The shape functions and the weights for 2D or 3D elements are obtained by the Kronecker product of vectors of individual axes, denoted by ⊗ as follows: (3)

2D Spectral Modelling
According to the first-order shear deformation theory [37,38], the displacement field is expressed as: where u e 0 , v e 0 and w e 0 are nodal displacements, ϕ e x , ϕ e y are the rotations of the normal to the mid-plane with respect to the axes x and y, respectively.
The nodal bending strain-displacement relations are given in the form: The nodal shear strain-displacement relations are given in the form:

3D Model of the PZT Transducer
The displacement vector of the PZT transducer is composed of three translational displacements and is defined as: where u e , v e and w e are displacements of the element nodes in ξ, η and ζ direction. The nodal strain-displacement relations are given as [39]: The electromechanical coupling is governed by the linear constitutive equation of piezoelectric material according to [6,40], and this is defined as: where σ and S are the stress and strain components, respectively, c E is the stiffness coefficient matrix measured at zero electric field, e is the piezoelectric coupling tensor, S is the electric permittivity, and E and D are the electric field and electric displacement measured at zero strain. The superscript T denotes a transpose matrix. The electric field is defined as: where φ e is a nodal voltage of the transducer.

Displacements Coupling at the Substructures Interface
The present model of the sandwich panel consists of 2D and 3D elements. Moreover, there are non-matching grids between two adjacent substructures. These involve connecting them by imposing the compatibility of the displacements at the interface, see Figure 1. This type of connection is implemented through the interface elements based on Lagrange multipliers, which are interpreted as forces responsible for determining the appropriate displacements of nodes. The coupling can be expressed as: where s i1 and s i2 are substructures connected by the interface Γ i . For the whole structure, the Equation (12) can be written in the matrix form: where G is the coupling matrix, which contains the equations to interpolate the substructures displacements at the interfaces, and d is a global displacement field for nS number of substructures, composed as: The general formulation of the matrix G is presented in Algorithm A1 from Appendix B. The main task of the algorithm is to calculate shape functions for each adjacent substructures at the points X p = (x k p , y k p ), which are projections of the interface nodes onto these substructures.
The shape function can be calculated after finding an owner element and local coordinates of the points. The owner element is a spectral element in the domain of the substructure s ij , which contains the interface node, for example, interface node k Γ = 36 (see Figure 1a) is located in the element e I 3D and e I I I 2D for the substructures s 11 and s 12 , respectively. This can be found in two ways: using MATLAB's built-in function inpolygon or more the efficient procedure proposed by Silva et al. [41], which was used in the current implementation.
The transformation from global to local coordinates was realised by the iterative method presented in the work of Li et al. [42]. The computational effectiveness of Algorithm A1 can be easily improved if certain precautions are taken. First, the mesh of the interface has to be based on the mesh from one of the substructures s i1 , s i2 , which may be referred to as a slave. Thus, the shape function takes only the values of one and zeros. Moreover, the code can be implemented in vectorized form rather than using for-loops.

Elementary Governing Equations of Motion
The classical equations of motion Md + Cḋ + Kd = F known from FEM are complemented by piezoelectric and interface coupling. Thus, the governing equations are defined as: where M dd , C dd and K dd are the structural mass, damping and stiffness matrices, respectively; K φd = K dφ T are piezoelectric coupling matrices; K φφ is the dielectric permittivity matrix, d is the vector of unknown nodal displacements, φ is the electric potential vector, F is the nodal external force vector, Q is the nodal charge vector, λ is the Lagrange multiplier vector, and G is the interface coupling matrix (˙) = ∂ ∂t . The formulae of the matrices are provided in Appendix A. The coupling is realised by imposing the traction forces as represented by a vector of Lagrange multipliers.

Parallel Implementation of the Internal Force Vector Calculation
The presented HSC model occupies much more operating memory than the homogenized one; thus, in order to achieve the solution in a reasonable time, the computation is performed using a multicore graphics processing unit (GPU). The most time-consuming operation in the Equation (15) is calculation of the internal force vector as: It should be noted that the stiffness matrix K dd occupies a large amount of memory. Instead of allocating matrix K dd , Kudela proposed a parallelized computation of the internal force vector [32].
The calculation is performed in three steps. First, the strain vectors are determined by multiplying the vectors of local node displacements and a sparse matrix containing local shape function derivatives. Then, the local internal force vector is obtained by multiplying the strain vectors by an appropriate material coefficient and a matrix of local shape function derivatives. Finally, the transformation of the local internal forces into the global forces is performed.

Transformation of the Core Elements
All core elements are rotated relative to both skins, and thus it is necessary to transform the degrees of freedom from the local coordinate system of the core to the global coordinate system. For this purpose, an additional sixth DOF is incorporated, i.e., rotation with respect to the z-axis: First, the displacement vector is transformed from the global to local coordinate system by the direction cosines as follows: where V e 1 ,V e 2 and V e 3 are direction cosines of the core element. Then, internal forces are calculated according to guideline from Section 2.6 and transformed to a global coordinated system: Additionally, a part of the mass matrix accounted for rotary inertia has to be transformed, and, in contrast to the internal forces vector, this has to be done only once in pre-processing as follows: As the matrix becomes non-diagonal after transformation, some approximation is necessary. Off-diagonal terms in the matrix given in Equation (20) are neglected following the analysis performed in [43].

A Solution of the Equation of Motion
Assuming b and f represent order lists of the electrode nodes and free nodes of the PZT, respectively, the electrical potential vector is rewritten: Then, Equation (16) is expressed as: where the notation K(r, c) uses vectors r and c to extract rows and columns from the matrix K, respectively, and (:) means all rows or columns of K. The electrical potential of the free nodes can be extracted from Equation (22): Substituting Equations (21) and (23) into Equation (15), the equation of motion can be rearranged into the form: where . The unknown displacement vector d t is found using a central difference algorithm [39]. Thus, Equation (24) is rewritten as: where ∆t is the time increment. Imposing the constrain Equation (13), the vector of Lagrange multipliers λ t can be extracted from Equation (25): where L ± = 1 ∆t 2 M dd ± 1 2∆t C dd .

Experimental Validation
The presented model was validated with results from two experimental studies. The first one was performed for determination of the full wavefield of the propagating waves by the scanning laser Doppler vibrometer (SLDV, Polytec PSV-400). The second study was performed for wave acquisition by the PZT sensor. The schematic of the experimental setup is shown in Figure 2. The sample of interest was a not-regular hexagonal aluminium honeycomb bonded to one CFRP plate using the epoxy adhesive (Loctite EA3479B) as shown in Figure 3a).The subject of the parametric study was the effect of the disbond size on the propagating GW.  After a reference measurement was made on an intact sample, several measurements were taken for the subsequent damage introduced on the same specimen. The circular area of the core was detached from the adhesive at the centre of the plate using a sharp hooked tool. For this purpose, the bottom skin was omitted so that damage could be introduced. The damage size was controlled by its diameter Φ D = [10, 30, 50, 70, 90, 110, 130] mm.
The generation and reception of elastic waves were achieved with a pair of PZT transducers mounted on the skin top surface with the cyanoacrylate glue. The coordinates of the actuator were (x 1 , y 1 ) = (−100, 0) mm, and for the sensor, (x 2 , y 2 ) = (100, 0) mm. The dimensions of the sample components were as follows: • CFRP skin: L × W = 500 × 500 mm, H = 1.5 mm. The N c = 5 cycle Hann windowed signal at carrier frequencies f c = [75, 100, 125, 150] kHz was generated using an arbitrary waveform generator (National Instruments, PXI 5413). The signal was amplified 40 times and supplied to the piezo actuator (Noliac, NCE51). Each measurement was conducted in the room temperature and averaged 20 times in order to improve the signal to noise ratio.

Simulation Parameters
All structures used to create the sample were modelled in the simulation with the following elements: 2D for the core, epoxy adhesive and cyanoacrylate glue and 3D for the CFRP plate and PZT transducers. During the creation of the mesh, special attention was taken to reduce the number of non-zero values in the matrix G. While the inversion of the matrix GL −1 + G T is necessary to calculate the vector of Lagrange multipliers in Equation (26) and L + is a diagonal matrix, the sparsity of the matrix G has a significant effect on the computation cost.
One spectral element was intended for each wall of the honeycomb core, while the meshes of the skin plates and the adhesive layers were divided by three rhombus elements per area under the core cell. In this way, the interface nodes coincide with the nodes lying on the hexagon edges (thick line on Figure 4b). The mesh of the cyanoacrylate glue was generated using external software GMSH [44] (see Figure 4c) and joined to the plate by non-matching interface elements. The PZT mesh coincides with the glue mesh. The convergence of the solution requires time increment to be less than a critical value, above which the displacements go to infinity. The critical value of time increment depends on the mesh size and the wave mode velocity. In the present model, convergence was achieved for 3 × 10 −9 s. Additionally, the following number of nodes in the elements were used: the core 6 × 5, epoxy adhesive and cyanoacrylate glue 6 × 6, the plate 6 × 6 × 4 and PZT transducers 6 × 6 × 3.
As the cells in the damaged area become distorted during core separation Figure 5a, the damage was modelled by removing the core elements in the disbond area as shown in Figure 5b. The material properties used in the simulations are gathered in Appendix C.

Homogenized Model
For this paper, comparative studies were conducted between the current model and the homogenized one. In the simplified model, the values of the material constants of the panel core were calculated according to the method presented by Malek and Gibson [45]. The effective mechanical properties for an aluminium core are gathered in Table A2, while the properties for other structures, i.e., the skin, the epoxy adhesive, the cyanoacrylate glue, and the sensors remained unchanged. The core element has 6 × 6 × 4 nodes, and the mesh coincides with the plate mesh. The models of the other structures remain unchanged.

The Severity of Damage Estimation
The severity of damage was estimated based on the function determined with the numerical simulation. A simple flowchart given in Figure 6 represents a process for the sample assessment. When the structure model is developed, several computer simulations for various damage sizes must be conducted to determine the MADIF.
The MADIF indicates the damage size according to measured damage index I normalized by the value obtained for the pristine sample I re f . In the paper, two types of damage index I are considered: the energy I eng and the maximum value of the half-width of the first package arrived in the sensor I amp , and these are defined as: where T is a period of the signal. Ψ g (t, Φ D ) is for the damaged case scenario, whereas Ψ g (t, 0) is for the pristine sample and it is realized in the same way by windowing the full-length signals of the sensor Ψ(t) with a flattened Gaussian window g(t) as follows: where t 0 is the centre and w g = 0.5N c / f c is a half-width of the window. Windowing the signals ensures obtaining the signals without any reflections from the boundaries. The determination of Ψ g is pictured in Figure 7a.  In the time domain, an equivalent numerical signal to the signal registered by the PZT acquisition instrument is calculated as an average value of the electrical potential of the electrode surface where n = 1 and n = 2 correspond to the homogenized and presented model, respectively. The MADIF is achieved by approximating the inverse of the computed damage index that best matches the experimental one. Finally, the damage size Φ D is obtained from the MADIF curve for measuring the normalized value of I/I re f as it is presented in Figure 7b.

Comparison of the Models
The snapshots for the pristine and the damaged sample are shown in Figures 8 and 9, respectively. One can observe the wave reflections in the core cells for experimental measurements and the present model. Additional, the front of the incident wave is distorted for the measurements of 125 and 150 kHz. The wavefront distortion in the present model is observed in the full range of frequency. Such effects are not noticeable in the simplified model because the wave propagates smoothly through the structure. The wavefront improvement in the experiment and the present model is noticeable in the undamaged region and marked by the red curves in Figure 9. This is the effect of a lack of reflection with the core.

Model-Assisted Damage Identification Function
I eng and I amp were determined for experimental measurements and numerical calculations in the function of damage size and the carrier frequencies. The indices with the curve fitted by the polynomial interpolation of order three are shown in Figures 10 and 11. It can be seen that both indices increase with the damage size in all cases. This is the effect of the leaky GW phenomenon [26].
Waves propagating through the plate lose energy in contact with the core. While GW propagates in the damaged area, such an effect does not occur, and thus the signal amplitude arriving at the sensor is higher. The present model is in good agreement with the experimental results for the tested frequencies f c = [75, 100, 125, 150] kHz. The homogenized model is in good agreement with the experiment for the lowest frequency, and the differences increase for higher frequencies. This issue may be related to the fact that the a shorter wavelength wave can lose more energy due to reflections at the edge of the damage of a homogenized model compared with a wave reflecting off the core cells.
To qualify the index as the MADIF, it must meet the condition of matching the numerical results with the experiment. The matching condition is that the value of the mean absolute error (MEA) must be less than an assumed threshold. The MEA is defined as: where the superscripts n and e correspond to the numerical models and experimental measurements, respectively; and d is the number of damage cases. The threshold is the smallest to the most extensive damage ratio expressed as a percentage, i.e., threshold = Φ min D /Φ max D × 100.  It can be seen from Figure 12 that the following indices satisfy the MAE condition: honeycomb index I eng for frequencies f c = [100, 125] and I amp for f c = [75, 100, 125, 150] kHz, but none from the homogenized model satisfy the MADIF selection criterion. In the case under consideration, the best fitting index turns out to be honeycomb I amp in 125 kHz; therefore, it was chosen for the MADIF, which is shown in Figure 13.

Conclusions
This paper presents preliminary research on the possibility of using a model-assisted approach to identify the severity of damage in a composite structure using GW propagation. For this purpose, the HSC model was implemented with the actual geometry of the honeycomb core. In contrast to full-structure homogenization, which is the most common HSC model found in the literature, the interactions of the propagating wave with the core cell walls were visible in the current model. The MADIF determined by the present model was in better agreement with the experimental measurements compared with the homogenized one.
In future works, our model will be used for parametric investigation to determine the MADIF in varied environmental conditions. The extended model will be usable in developing SHM scenarios.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to due to privacy restrictions of the ongoing research.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The formulae of matrices for 3D elements are: where c is the stiffness tensor, ρ is the mass density, and V e is the element volume.
In the case of the 2D elements, the matrices are defined as: where h = h t + h b is the element thickness, while h t(b) is the distance between mid-plane and top (bottom) surface of the element, and Ω e is the element area: The dielectric conductivity matrix K e φφ and piezoelectric coupling matrix K e uφ are defined: