A likelihood-free Bayesian approach for characterisation of multiple delaminations in laminated composite beams

Abstract

This paper presents a probabilistic model-based optimisation framework for localising and characterising multiple delaminations in laminated composite beams using ultrasonic guided waves (GW). A likelihood-free Bayesian method, approximate Bayesian computation by subset simulation (ABC-SS), is implemented to determine the number of delaminations and identify the unknown damage characteristics and their associated uncertainties. The ABC algorithm provides a practical way to approximate the posterior distributions of uncertain damage parameters and select the most plausible damage model for determining the number of delaminations through direct comparison of the experimentally measured and numerically simulated GW signals without assuming any likelihood functions. To overcome the expensive computational cost of traditional finite element simulations, a higher-order laminated model (HOLM) is employed to model the GW propagation behaviour in the delaminated composite beams with satisfactory accuracy and acceptable computational efficiency. Benefiting from the accurate simulation, the dataset comparison utilises the original time-domain GW signals, thereby preventing the loss of any key information from the signals for damage identification. A comprehensive series of numerical case studies are used to demonstrate the accuracy, robustness and feasibility of HOLM simulation and the proposed multiple-delamination identification framework. The practicability and accuracy of the proposed ABC framework are further validated using two experimental datasets.

Introduction

Laminated composite beams

Laminated composite structures generally consist of an assembly of physical layers. These composite materials or structures have wide applications in many engineering fields because of their advantageous characteristics, such as high-strength, lightweight and anti-corrosion. One of the most common categories of defects in laminated materials is delamination.^1,2 This would cause a considerable reduction in structural stiffness and may lead to a high structural failure risk. A robust structural health monitoring (SHM) technique for identifying and characterising damages in their early state is essential to ensure structural serviceability and reliability, especially for high-value asset applications, such as aerospace and energy engineering. Depending on the diagnosed location, type and size, engineers can propose a specific maintenance plan and make recommendations for operation. ³ Delamination is considered a type of internal defect that cannot be easily detected using visual inspection. Non-destructive evaluation (NDE) techniques are, therefore, usually utilised in SHM for damage detection and assessment.^4–6 Compared with conventional NDE techniques, ultrasonic guided wave (GW) based methods are broadly used for damage detection in composite structures with an increasing trend due to their high sensitivity to local and incipient flaws, the capability of a relatively long propagation distance as well as the ability to inspect inaccessible areas.^7–11

Model-based GW approach

To determine the damage location and, in particular, to identify the damage geometry, the model-based GW approach has attracted much attention nowadays due to its baseline-free nature and capability to deal with more complicated damage situations. ¹² Specifically, an analytical or numerical model is used to simulate the propagation behaviours of incident waves and their reflection or transmission from the damages, boundaries and interfaces. The damage properties are considered unknown variables and identified by minimising the differences between experimentally measured and numerically simulated GW responses. ¹³

The ability of the numerical model to predict GW interaction with damages hence plays a vital role in model-based approach. GW propagation could be simulated accurately using the traditional finite element (FE) method with the use of an adequate small element size. Nevertheless, the computational cost is impractical, especially for composite structures, which usually need to be explicitly modelled using numerous elements per layer in order to obtain precise GW damage features.^7,9 Different modelling methods have been developed to overcome the limitations of the conventional FE method.^14,15 The frequency-domain spectral finite element (SFE) method is computationally efficient compared to the FE method, which utilises spectral analysis techniques to examine structural behaviours at discrete frequencies by solving the partial differential governing equations in the frequency domain.^16,17 Time-domain SFE is an alternative that employs Gauss-Lobatto-Legendre (GLL) nodes and high-order approximation polynomials to reduce element number in FE simulation while providing accurate prediction. Runge’s phenomenon can also be prevented by using GLL nodes.^18,19 Other advanced techniques, such as semi-analytical FE and hybrid/wave and finite element (WFE), involve the combination of analytical solutions from wave propagation principles and numerical solutions from FE method to obtain the GW scattering features of damage instead of the entire signal, such as reflection and transmission coefficients. Modelling complexity and cost can, therefore, be reduced while the simulation is relatively accurate.^20–22 These GW scattering features from damages can then be treated as numerical responses for damage quantification. Additionally, with the rapid development of machine learning and neural networks, surrogate models to predict GW propagation with damage features have also been widely investigated nowadays.^23–25 To summarise, while simulation accuracy is of paramount importance, a robust model-based approach must also take into consideration computational effort to remain feasible.

Bayesian GW model-based framework

In addition to the constraint of numerical simulation, the associated uncertainties tend to compromise the efficacy of damage identification. These uncertainties come from various aspects, such as measurement error due to noise from environment and data acquisition systems, modelling error due to model simplification and unspecified material/geometry properties, etc. To accurately define the damage properties and also quantify the involved uncertainties, the Bayesian inference approach was initially utilised in vibration-based SHM and was then extended to GW model-based SHM with a rising tendency.^26,27 Bayesian inference considers prior knowledge about the damage and current observations in the model updating procedure. It provides a useful platform to perform online SHM on structures.

Some Bayesian frameworks for GW damage identification can be found in the literature. ²⁸ Ng proposed a model-based Bayesian probabilistic framework with a particle swarm optimisation algorithm for identifying single crack in an isotropic beam structure based on different GW scattering features obtained from different data pre-processing techniques. ²⁹ Yan et al. ³⁰ developed a fast Bayesian scheme with the transitional Markov Chain Monte Carlo (TMCMC) algorithm for damage identification using an analytical probabilistic surrogated model with the WFE to simulate the scattering coefficients. Apart from parameter estimation, the Bayesian approach also provides an essential feature to determine the most plausible model. He and Ng introduced the Bayesian model assessment with TMCMC and subset simulation (SS) to define the unknown number of cracks in beam-like structures.^31,32 Different model classes, that is, models with different number of cracks, were assessed by comparing model evidence. Hibert transform of GW signal was used to enhance the sensitivity to identify the corresponding damage geometries. Recently, Cantero-Chinchilla proposed a multi-level Bayesian framework with Asymptotic independence Markov sampling algorithm for damage localisation in plate-like structures. ³³ The model class assessment was utilised at the first level to determine the most probable time-frequency (TF) transform technique from raw GW data. Damage localisation with uncertainty quantification was then processed based on the most plausible TF model.

Most studies focused on Bayesian GW model-based approach, despite their effectiveness and rigour, rely on extracting features from GW signals rather than using the entire signal information. Extracting these key features tends to reduce the computational burden as the signal is compressed so that some unnecessary ‘noise’ error can be ‘filtered out’. ²⁹ This works well when the GW damage reflections are simple, for example, in a single-crack situation. Conversely, in complicated damage situations, for example, multiple-damage situations, the signal can be complex as multiple scattered waves are induced by defects, and mode conversion also occurs. The damage reflections are always overlapped and difficult to be separated. Data compression or feature extraction may result in the loss of useful damage information from signals or inevitably introduce additional errors, leading to inaccurate or even incorrect damage identification. ³⁴ Thus, there is still a research gap in developing an accurate and computationally efficient model-based approach using the entire time-domain GW signal to identify damage characteristics.

Approximate Bayesian computation

In Bayesian approach, prior distributions are generally determined based on engineers’ prior knowledge or experience of the damage. Conversely, Bayesian inference also requires defining a likelihood function to measure the level of agreement between predicted and measured responses. ²⁷ The selection of the likelihood function can vary depending on the assumption of uncertainty about the underlying unknown data-generation distribution. In the same damage identification scenario, different likelihood assumptions may result in different posterior distributions of uncertain parameters, that is, different uncertainty quantification results. Additionally, for model simulation involving high-dimensional parameters or complex nonlinearity, the analytical formulation of the likelihood function might be elusive.^35–37 As a result, there is no consistent approach that can be applied to all circumstances.

The approximate Bayesian computation (ABC) method, which can bypass likelihood function evaluation, has recently attracted significant attention in the development of a more general and practical probabilistic framework for damage identification. ABC approximates the posterior distributions from the priors by evaluating the similarity directly between the observed and simulated responses.^35–37 In this way, ABC still obeys the Bayes theorem principle, while the difficulties of the likelihood function, as mentioned above, are avoided cleverly. ABC can be employed with any forward-simulated model. Some ABC applications can be found in vibration-based SHM for model parameter estimation and dynamic model selection,^38,39 while there is limited research in GW-based SHM for damage identification.^25,40 Numerous novel ABC algorithms have also been proposed in the last decade to enhance approximation accuracy and reduce computational costs. Some common examples include ABC with Markov Chain Monte Carlo (MCMC), ABC with sequential Monte Carlo (SMC), ³⁸ and approximate Bayesian computation by subset simulation (ABC-SS). ⁴¹

To summarise, achieving simulation accuracy is vital to precisely capture GW reflections and scattering. A robust model-based approach must also consider computational effort to remain feasible. It is still worthwhile to develop an advanced FE model that can accurately and efficiently simulate GW propagation within delaminated composites. Bayesian approaches are capable of quantifying unavoidable uncertainties. However, most GW-based Bayesian frameworks were developed depending on different likelihood assumptions, which might limit their applicability. Recently, ABC algorithms have been considered as a more general and practical alternative for model class selection and parameter estimation while still adhering to Bayes’ theorem.

In this paper, a likelihood-free Bayesian framework is developed to identify multiple delaminations in laminated composite beams. An advanced version, ABC-SS, is employed to enhance inference efficiency. In ABC-SS, the SS technique for rare-event simulation is combined with the ABC principle to automatically create a nested decreasing sequence of regions corresponding to an increasingly closer approximation of the observed data, that is, the posterior approximation. ABC-SS algorithm can also evaluate model evidence in a simple by-product manner, which enhances model class selection effectiveness. A higher-order laminated model (HOLM) is also proposed for modelling time-domain GW signals in composite beams. Higher-order deformation kinematics in both thickness and length directions are employed to guarantee simulation accuracy for damage identification. A sub-laminate scheme is used to provide a convenient, efficient and physical way for simulating multiple delaminations while maintaining computation efficiency. SFE and GLL nodes are also used in the model to further reduce computation costs.

The paper is arranged as follows: Firstly, the HOLM formulation and delamination modelling mechanism is provided in section ‘Higher-order laminated model’. ABC-SS and the likelihood-free framework for multiple-delamination characterisation are presented in section ‘Bayesian likelihood-free approach for multiple-delamination characterisation’. Numerical case studies under noise conditions and different numbers of delaminations are used to comprehensively demonstrate the robustness and accuracy of the HOLM and the probabilistic framework in section ‘Numerical case studies’. After that, in section ‘Experimental study’, experimental case studies are employed to further validate the accuracy and practicability of the proposed framework. Conclusions are given in section ‘Conclusion’.

Higher-order laminated model

Formulation of the laminated model GW propagation simulation

This section presents the proposed HOLM for modelling GW propagation in a multi-layered composite beam. The entire beam is simulated by representing it as a group of HOLM-based sub-laminates, with each sub-laminate containing several physical layers (see Figure 1(a) for an example).

Figure 1.

(a) An eight-physical-layer composite beam simulated by three HOLM sub-laminates. (b) $\mathrm{x}$ – $z$ plane orientation in a HOLM sub-laminate with four physical layers (using HOLM2 as an example).

To accurately capture GW propagation behaviour, higher-order through-thickness deformation kinematics is integrated to simulate each HOLM sub-laminate.^42,43 As shown in Figure 1(b), the deformation for a sub-laminate is assumed to have a two-dimensional (2D) plane stress condition in the $x-z$ plane. The in-plane ( $\mathrm{u}$ ) and transverse ( $v$ ) displacements of the HOLM are expressed as:

$$u=u_0-z\theta +z^2\alpha +z^3\beta$$

(1)

$$v=v_0+z\gamma +z^2\delta$$

(2)

where $u_0$ and $v_0$ are in-plane and transverse displacements at the mid-plane, respectively. $z$ is the coordinate along the thickness direction. $\theta$ is the sectional plane rotation normal to $x$ -axis. $\alpha ,\beta ,\gamma$ and $\delta$ are non-physical weighting terms that are used to provide higher-order variations.

These weighting factors can be expressed by incorporating both $u$ and $v$ displacements at the top ( $u_T$ , $v_T$ ) surface and bottom ( $u_B$ , $v_B$ ) surface. By substituting $z=\frac{h}{2}$ and $z=\frac{-h}{2}$ into Equations (1) and (2) and some formulation deduction of ( $u_T+u_B$ ), ( $u_T-u_B$ ) and ( $v_T+v_B$ ), ( $v_T-v_B$ ), these terms can therefore be expressed:

$$\alpha =\frac{2(u_T+u_B-2u_0)}{h^2}$$

(3)

$$\beta =\frac{4(u_T-u_B+h\theta )}{h^3}$$

(4)

$$\gamma =\frac{v_T-v_B}{h}$$

(5)

$$\delta =\frac{2(v_T+v_B-2v_0)}{h^2}$$

(6)

Both $u$ and $v$ displacement can then be simplified in the following forms:

$$u=A_u{u}_0+B_u\theta +C_u{u}_T+D_u{u}_B$$

(7)

$$v=A_v{v}_0+B_v{v}_T+C_v{v}_B$$

(8)

where $A_u=\frac{h^2-4z^2}{h^2}$ , $B_u=\frac{4z^3-h^{2}z}{h^2}$ , $C_u=\frac{2h{z}^2+4z^3}{h^3}$ , $D_u=\frac{2h{z}^2-4z^3}{h^3}$ , $A_v=\frac{h^2-4z^2}{h^2}$ , $B_v=\frac{\mathrm{hz}+2z^2}{h^2}$ and $C_v=\frac{2z^2-\mathrm{hz}}{h^2}$ . $h$ is the thickness of the HOLM sub-laminate.

GW propagation simulation is a dynamic problem that can be treated as solving a time-domain ordinary differential equation, also known as the equation of motion. The equation can be derived using the Lagrange equation as below.

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial L}{\partial \stackrel{·}{\mathrm{\Delta }}}\right)-\frac{\partial L}{\partial \mathrm{\Delta }}+\frac{\partial R}{\partial \stackrel{·}{\mathrm{\Delta }}}=0$$

(9)

where t represents the time. $L=T-V$ is the Lagrangian that consists of the kinetic energy $T$ and potential energy $V$ . $\mathrm{\Delta }$ and $\stackrel{·}{\mathrm{\Delta }}$ is the displacement and velocity vector. $R$ is the dissipation used to account for damping.

The main factors of $T,V$ and $R$ can be expressed as:

$$T=\frac{1}{2}\int \int ({\rho }_x{\stackrel{·}{u}}^2+{\rho }_z{\stackrel{·}{v}}^2)\mathrm{d}x\mathrm{d}z$$

(10)

$$V=\frac{1}{2}\int \int ({\mathrm{\epsilon }}^{\mathrm{T}}\mathbf{\sigma })\mathrm{d}x\mathrm{d}z-\int (u{q}_x+v{q}_z)\mathrm{d}x$$

(11)

$$R=\frac{1}{2}\int \int ({\mu }_x{\stackrel{·}{u}}^2+{\mu }_z{\stackrel{·}{v}}^2)\mathrm{d}x\mathrm{d}z$$

(12)

where $\rho$ and $\mu$ represent the mass density and the damping coefficient, respectively. $\mathrm{\epsilon }$ and $\mathbf{\sigma }$ are the strain and stress vectors. $\stackrel{·}{u}$ and $\stackrel{·}{v}$ represent the time derivatives of displacements, also known as velocity. $q_x$ and $q_z$ are the distributed in-plane load and transverse load that have no variations along the $z$ -axis.

The displacements $u$ and $v$ in Equations (7) and (8) can be expressed by a $z$ -axis dependent sectional displacement matrix ${\mathbf{H}}_d$ , and a one-dimensional (1D) displacement field $\mathbf{f}$ .

$$\left[\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right]={\mathbf{H}}_d\mathbf{f}$$

(13)

where

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{H}}_d}& =\left[\begin{array}{ccccccc}{A}_u\hfill & B_u\hfill & C_u\hfill & D_u\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & A_v\hfill & B_v\hfill & C_v\hfill \end{array}\right]\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}^{\mathrm{T}}}& =\left[\begin{array}{ccccccc}{u}_0\hfill & \theta \hfill & u_T\hfill & u_B\hfill & v_0\hfill & v_T\hfill & v_B\hfill \end{array}\right]\hfill \end{array}$$

SFEs with higher-order approximations of the $x$ -direction displacement components of $\mathbf{f}$ are employed to simulate GW propagation in the composite beam with small-scale defects. GLL nodes are employed to improve interpolation accuracy. All seven displacement factors in $\mathbf{f}$ are interpolated using the same number of GLL nodes $n$ per element; the $k\mathrm{t}\mathrm{h}$ displacement factor $f^k$ can be represented as:

$$f^k=\sum _{i=1}^n{N}_i{\Delta }_i^k\mathrm{for}k\in [1,7]$$

(14)

where $N_i$ is the 1D interpolation shape function. ${\Delta }_i^k$ is the nodal displacement at the $k\mathrm{t}\mathrm{h}$ factor.

The displacement vector at any nodes of the $x$ - $z$ plane within a HOLM sub-laminate can be further denoted by the substitution of Equation (14) into Equation (13).

$$\left[\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right]={\mathbf{H}}_d{\mathbf{B}}_d\mathrm{\Delta }$$

(15)

where ${\mathbf{B}}_d$ consists of 1D interpolation shape functions $N_i$ .

Following a similar manner, velocities $\stackrel{·}{u}$ and $\stackrel{·}{v}$ can be expressed using nodal velocity vector.

$$\left[\begin{array}{c}\stackrel{·}{u}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{v}}\end{array}\right]={\mathbf{H}}_v{\mathbf{B}}_v\stackrel{·}{\mathrm{\Delta }}$$

(16)

where ${\mathbf{H}}_v={\mathbf{H}}_d$ is the velocity sectional matrix and ${\mathbf{B}}_v={\mathbf{B}}_d$ .

The strain vector $\mathrm{\epsilon }$ in 2D plane can be expressed as a function of the derivatives of $u$ and $v$ . By utilising primary displacements in Equations (7) and (8), $\mathrm{\epsilon }$ can be denoted in terms of a sectional matrix $\mathbf{H}$ that depends on the z-axis, and a 1D strain vector $\overline{\mathrm{\epsilon }}$ that depends on the x-axis.

$$\mathbf{\epsilon }=\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_z\\ {\gamma }_{xz}\end{array}\right]=\left[\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial z}\\ \frac{\partial u}{\partial z}+\frac{\partial v}{\partial x}\end{array}\right]=\mathbf{H}\overline{\mathbf{\epsilon }}$$

(17)

where

$$\begin{array}{c}\mathbf{H}=\hfill \\ \left[\begin{array}{cccccccccccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & A_u\hfill & B_u\hfill & C_u\hfill & D_u\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{\partial A_v}{\partial z}\hfill & \frac{\partial B_v}{\partial z}\hfill & \frac{\partial C_v}{\partial z}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\partial A_u}{\partial z}}& \frac{\partial B_u}{\partial z}\hfill & \frac{\partial C_u}{\partial z}\hfill & \frac{\partial D_u}{\partial z}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & A_v\hfill & B_v\hfill & C_v\hfill \end{array}\right]\hfill \\ {\overline{\mathrm{\epsilon }}}^{\mathrm{T}}=\left[\begin{array}{cccccccccccccc}{u}_0\hfill & \theta \hfill & u_T\hfill & u_B\hfill & \frac{\partial u_0}{\partial x}\hfill & \frac{\partial \theta }{\partial x}\hfill & \frac{\partial u_T}{\partial x}\hfill & \frac{\partial u_B}{\partial x}\hfill & v_0\hfill & v_T\hfill & v_B\hfill & \frac{\partial v_0}{\partial x}\hfill & \frac{\partial v_T}{\partial x}\hfill & \frac{\partial v_B}{\partial x}\hfill \end{array}\right]\hfill \end{array}$$

By substituting all displacement vectors $f^k$ from Equation (14), $\overline{\mathrm{\epsilon }}$ can be expressed as a function of the nodal displacement vector $\mathrm{\Delta }$ :

$$\overline{\mathrm{\epsilon }}=\mathbf{B}\mathbf{\Delta }$$

(18)

where $\mathbf{B}$ is the strain-displacement matrix consisting of $N_i$ and the corresponding derivatives.

Since the HOLM assumes a plane stress condition in the $x-z$ plane that ${\sigma }_y={\tau }_{\mathrm{xy}}={\tau }_{\mathrm{yz}}=0$ , the dimension of the stress-strain relationship can be reduced from 6 (i.e., three-dimensional (3D) problem) to 3 (i.e., 2D problem). ⁴² The 2D strain-stress relationship of a composite ply can be expressed in terms of the constitutive matrix $\mathbf{D}$ and stresses $\mathrm{\epsilon }$ , as shown in the following form:

$$\mathbf{\sigma }=\mathbf{D}\mathrm{\epsilon }$$

(19)

By substituting Equations (15)–(19) into Equations (10)–(12), the formulation of the $T$ , $V$ and $R$ can be further derived and simplified.

$$T=\frac{1}{2}{\stackrel{·}{\mathrm{\Delta }}}^{\mathrm{T}}(\int {{\mathbf{B}}_v}^{\mathrm{T}}{\mathbf{D}}_m{\mathbf{B}}_v\mathrm{d}x)\stackrel{·}{\mathrm{\Delta }}=\frac{1}{2}{\stackrel{·}{\mathrm{\Delta }}}^{\mathrm{T}}\mathbf{M}\stackrel{·}{\mathbf{\Delta }}$$

(20)

$$V=\frac{1}{2}{\mathrm{\Delta }}^{\mathrm{T}}\int {\mathbf{B}}^{\mathrm{T}}{\mathbf{D}}_K\mathbf{B}\Delta \mathrm{d}x-{\mathrm{\Delta }}^{\mathrm{T}}\int {{\mathbf{B}}_d}^{\mathrm{T}}\mathbf{q}\mathrm{d}x=\frac{1}{2}{\mathrm{\Delta }}^{\mathrm{T}}\mathbf{K}\mathbf{\Delta }-{\mathrm{\Delta }}^{\mathrm{T}}\mathbf{F}$$

(21)

$$R=\frac{1}{2}{\stackrel{·}{\mathrm{\Delta }}}^{\mathrm{T}}(\int {{\mathbf{B}}_v}^{\mathrm{T}}{\mathbf{D}}_c{\mathbf{B}}_v\mathrm{d}x)\stackrel{·}{\mathrm{\Delta }}=\frac{1}{2}{\stackrel{·}{\mathrm{\Delta }}}^{\mathrm{T}}\mathbf{C}\stackrel{·}{\mathbf{\Delta }}$$

(22)

where

$$\begin{array}{c}{\mathbf{D}}_m=\int {{\mathbf{H}}_v}^{\mathrm{T}}\left[\begin{array}{cc}{\rho }_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\rho }_z\hfill \end{array}\right]{\mathbf{H}}_v\mathrm{d}z\hfill \\ {\mathbf{D}}_K=\int {\mathbf{H}}^{\mathrm{T}}\mathbf{DH}\mathrm{d}z\hfill \\ {\mathbf{D}}_c=\int {{\mathbf{H}}_v}^{\mathrm{T}}\left[\begin{array}{cc}{\mu }_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mu }_z\hfill \end{array}\right]{\mathbf{H}}_v\mathrm{d}z\hfill \\ \mathbf{q}={{\mathbf{H}}_d}^{\mathrm{T}}\left[\begin{array}{cc}{q}_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& q_z\hfill \end{array}\right]{\mathbf{H}}_d\hfill \end{array}$$

$\mathbf{M}$ , $\mathbf{C}$ and $\mathbf{K}$ is the element mass matrix, damping matrix and stiffness matrix, respectively. $\mathbf{F}$ is element load vector. It should be noted that $\mathbf{C}$ can be also expressed using the concept of Rayleigh damping to model the damping effect. In this case, $\mathbf{C}=a\mathbf{M}+b\mathbf{K}$ , where $a$ and $b$ are the damping coefficients. ⁴⁴

The equation of motion for an element of a HOLM sub-laminate can be obtained by deriving the Lagrange equation in Equation (9) after substituting Equations (20)–(22).

$$\mathbf{K}\mathrm{\Delta }+\mathbf{C}\stackrel{·}{\mathrm{\Delta }}+\mathbf{M}\stackrel{··}{\mathrm{\Delta }}=\mathbf{F}$$

(23)

where $\stackrel{··}{\mathrm{\Delta }}$ are the unknown acceleration vector. The same form of the equation of motion can be applied to the entire composite beam by assembling all element matrices of the equation from all SFEs and HOLM sub-laminates. The global stiffness, damping and mass matrices, as well as the load vector, can then be constructed.

Delamination modelling scheme

The proposed HOLM considers displacements at both top and bottom external surfaces, providing a simple but physical mechanism to simulate delamination and model multi-layer composite beam (as shown in Figure 2). Specifically, the connection of each sub-laminate HOLM is achieved easily by connecting the same degrees-of-freedom at both upper/lower planes, leading to just one set of nodes in the contact region. Delamination can then be easily inserted into the interface by taking double sets of nodes between two planes. It should be noted that the delamination region can be simulated by using only one or multiple spectral elements, mainly depending on the length of the defect as well as the convergence of the model. Moreover, the proposed composite beam model improves the simulation accuracy in both contact and delamination regions by adopting higher-order variation in the displacement vector and considering multiple HOLM sub-laminates in the thickness direction. The computational efficiency is also enhanced by adopting time-domain SFE and GLL nodes to reduce the number of elements of the entire model, as well as accommodating multiple composite layers (with no delamination) within each sub-laminate.

Figure 2.

(a) A diagram example of a composite beam with three HOLM sub-laminates and three delamination simulations. (b) The mechanism of contact and delamination simulation at the interface between two adjacent HOLM sub-laminates (incorporating with SFE and GLL nodes).

In this paper, the proposed HOLM is used to simulate GW propagation in composite beam with multiple delaminations. The point load vector in the form of a GW incident pulse is applied to the excitation location of the model. The GW time-domain signals are calculated by solving the equation of motion for the entire beam model using an explicit time integration method, the central difference method. ³²

Bayesian likelihood-free approach for multiple delamination characterisation

The procedure of detecting multiple defects using the Bayesian approach typically involves two main stages: mode class selection to determine the number of defects and model updating to identify damage properties, that is, delamination location and severity.^31,45 Recall Bayes’ theorem, the posterior distribution for predicting the uncertain delamination parameters $\theta$ , conditional on the measured GW signal measurement $D_m$ and an input-output model $M$ (i.e., the delaminated beam model), can be denoted as:

$$p(\theta |D_m,M)=\frac{p(D_m|\theta ,M)p(\theta |M)}{P(D_m|M)}\propto p(D_m|\theta ,M)p(\theta |M)$$

(24)

where $p(D_m|\theta ,M)$ is the likelihood function. $p(\theta |M)$ is the prior distribution that includes the initial knowledge about the unknowns. $P(D_m|M)$ represents the marginal likelihood or evidence.

At the model class level, given a set of model classes $\mathbf{M}=(M_i,i=1\dots ,n)$ (i.e., beam models with different numbers of delamination), Bayesian model selection is achieved by ranking the considered models based on their posterior plausibility. ⁴⁶

$$P(M_i|D_m,\mathbf{M})=\frac{P(D_m|M_i)P(M_i|\mathbf{M})}{\sum _{j=1}^{n}P(D_m|M_j)P(M_j|\mathbf{M})}$$

(25)

where $n$ is the total number of considered classes and $P(M_i|\mathbf{M})$ is the prior probability of $M_i$ . Assuming no initial knowledge, each model class is usually considered equally prior plausibility that $P(M_j|\mathbf{M})=1/n$ and $\sum _{j=1}^{n}P(M_j|\mathbf{M})=1$ . Hence, $P(D_m|M_i)$ plays a significant role in defining the most plausible model.

Approximate Bayesian computation

The general idea of ABC is to bypass the likelihood evaluation by directly assessing the similarity between measurements $D_m$ and simulations $D_s$ from the model $M$ . ³⁶ By applying Bayes’ theorem, the approximate posterior from ABC is expressed as:

$$p(\theta ,D_s|D_m,\epsilon ,M)\propto p(D_s|\theta ,M)p(\theta |M){\mathbb{I}}_{{\rho }_{\epsilon }(D_m)}\left(D_s\right)$$

(26)

where the model parameter and model output are simulated from the distribution $p(\theta |M)$ and $p(D_s|\theta ,M)$ , respectively. $\epsilon$ is the tolerance level. ${\mathbb{I}}_{{\rho }_{\epsilon }(D_m)}\left(D_s\right)$ is the indicator function that equals to 1 if the distance function $\rho (D_s,D_m)\le \epsilon$ and 0 otherwise.

For some convenience, for example, high-dimensional datasets or computational burden, summary statistics $S(.)$ can be employed to reduce the dimensionality of comparison and also computational difficulty, in such that $\rho (S(D_s),S(D_m))\le \epsilon$ . ³⁵ Similar to the posterior calculation in Equation (26), the ABC model evidence can also be approximated by:

$$\begin{array}{c}\multicolumn{1}{c}{P(D_s|D_m,\epsilon ,M)=\int p(D_s|\theta ,M)p(\theta |M){\mathbb{I}}_{\rho \left(\epsilon \right)}\left(D_s\right)\mathrm{d}\theta \mathrm{d}{D}_s}\\ \approx \frac{1}{N}\sum _{i=1}^N{\mathbb{I}}_{{\rho }_{\epsilon }(D_m)}\left(D_{s(i)}\right)\hfill \end{array}$$

(27)

where $N$ is the number of samples required for approximation.

It should be noted that the approximation quality of Equations (26) and (27) mainly relies on the choice of distance function $\rho$ , the level of tolerance $\epsilon$ and the summary statistic $S(.)$ . If $S(.)$ is sufficient to include all damage information and the tolerance is sufficiently small, the approximation error of posterior and evidence can be eliminated. Nevertheless, assuming such a small tolerance in ABC approximation always leads to a rare-event problem. Many samples are required to achieve an acceptable sample size, staying within the approximate region that meets the indicator function requirement. In other words, the computational costs are prohibitive. In addition, the selection of $S(.)$ may only be able to guarantee sufficiency for a few model classes, which may threaten ABC performance for parameter estimation, especially model class selection. ⁴⁷ In recent years, many advanced ABC algorithms have been developed to circumvent these difficulties, such as ABC-SMC and ABC-SS. The reader can refer to the comprehensive literature. ³⁷

Approximate Bayesian computation by subset simulation

In this paper, ABC-SS will be employed for multiple-delamination identification. It provides a computationally effective way to estimate the evidence of each model class and the posterior of uncertain damage parameters by combining the SS technique with the ABC principle. ⁴¹ Specifically, ABC-SS converts the Bayesian inference to a reliability analysis for solving a rare-event problem. The basic procedure involves a sequence of nested decreasing data-approximating regions $D_{{\epsilon }_d}$ to define the final region $D_{\epsilon }$ . The intermediate regions in the SS can be defined as $D_{{\epsilon }_d}=[(\theta ,D_s):\rho (S(D_s),S(D_m))\le {\epsilon }_d]$ , which is similar to the indicator function in the ABC principle, as shown in Equation (26). The probability $P(D_{\epsilon })$ can be calculated as a product of the conditional probabilities from all intermediate regions.^48,49

$$P\left(D_{\epsilon }\right)=P\left(D_{{\epsilon }_1}\right)\underset{d=2}{\overset{J}{\Pi }}P(D_d|D_{{\epsilon }_{d-1}})\mathrm{for}d\in [1,J]$$

(28)

where $J$ is the total number of simulation levels. The intermediate tolerance is a set with a reducing sequence that ${\epsilon }_1>\dots >{\epsilon }_{d-1}>{\epsilon }_d>\dots >{\epsilon }_J$ .

At the first stage of ABC-SS ( $d=1$ ), N independent samples are generated from the priors. These samples are substituted into the distance function $\rho (D_s,D_m)$ for rearrangement in such that ${\rho }_{1,(1)}\le \dots \le {\rho }_{1,(N)}$ . The probability $P\left(D_{{\epsilon }_1}\right)$ is approximated based on the chosen ${\epsilon }_1$ .

$$P\left(D_{{\epsilon }_1}\right)=\frac{1}{N}\sum _{i=1}^N{\mathbb{I}}_{{\rho }_{{\epsilon }_1}\left(D_m\right)}\left(D_{s_{1,i}}\right)$$

(29)

At the simulation stage ( $d\ge 2$ ), $N$ candidate samples are drawn from the conditional PDF using the modified Metropolis algorithm (MMA), which was proposed in the original version of the SS. ⁵⁰ In a similar manner, the conditional probability at the $d\mathrm{t}\mathrm{h}$ stage is estimated based on the chosen ${\epsilon }_d$ :

$$P(D_{{\epsilon }_d}|D_{{\epsilon }_{d-1}})=\frac{1}{N}\sum _{i=1}^N{\mathbb{I}}_{{\rho }_{{\epsilon }_d}\left(D_m\right)}\left(D_{s_{d,i}}\right)$$

(30)

The tolerance level ${\epsilon }_d$ in each simulation stage is adaptively defined, such that the $P(D_{{\epsilon }_1})$ and $P(D_{{\epsilon }_d}|D_{{\epsilon }_{d-1}})$ is always equal to an assigned and fixed conditional probability $P_0$ . The value of ${\epsilon }_d$ can therefore be defined as the $100P_0$ percentile of distance function set ${\rho }_{d,(1)}\le \dots \le {\rho }_{d,(N)}$ .

$${\epsilon }_d=\frac{({\rho }_{d,\left(N{P}_0\right)}+{\rho }_{d,(N{P}_0+1)})}{2}$$

(31)

It is observed that sample size generated in $D_{{\epsilon }_{d-1}}$ that lies in $D_{{\epsilon }_d}$ is always known in advance. These $N{P}_0$ samples from $D_{{\epsilon }_{d-1}}$ are employed in the MMA as seeds to create $N{P}_0$ Markov chains with the length of $1/P_0$ and hence $N$ new samples are generated in $D_{{\epsilon }_d}$ . The probability can hence be estimated as a function of a fixed $P_0$ ⁴¹ :

$$P\left(D_{{\epsilon }_d}\right)=P\left(D_{{\epsilon }_1}\right)\underset{d=2}{\overset{J}{\Pi }}P(D_{{\epsilon }_n}|D_{{\epsilon }_{n-1}})=P_0^d$$

(32)

ABC-SS approximates the posterior probability of the uncertain parameters through an adaptive tolerance procedure until the intermediate tolerance ${\epsilon }_d$ becomes smaller than the desired final threshold ${\epsilon }_{\mathrm{final}}$ or the simulation stage $d$ reaches the maximum level $m_{\mathrm{final}}$ .

It is also observed that evidence estimation of the model class at different tolerance levels using ABC-SS is a simple by-product. For a specific $\epsilon$ and $d$ ( ${\epsilon }_{d-1}>\epsilon >{\epsilon }_d$ ), the ABC evidence of the model class $M$ can be expressed as a function of tolerance^51,52:

$$P(D_{\epsilon }|D_m,\epsilon ,M)=P_0^{d-1}P(D_{\epsilon }|D_{{\epsilon }_{d-1}})=P_0^{d-1}{P}_d$$

(33)

where $P_d$ is the sample size generated in $D_{\epsilon }$ that lies in $D_{{\epsilon }_{d-1}}$ . If the same tolerance level $\epsilon$ is applied to all considered models, the posterior plausibility of each model class can be determined by substituting Equation (33) into Equation (25).

$$\begin{array}{c}\multicolumn{1}{c}{P(M_i|D_m,\epsilon ,\mathbf{M})=\frac{P(D_{\epsilon }|D_m,\epsilon ,M_i)P(M_i|\mathbf{M})}{\sum _{j=1}^{n}P(D_{\epsilon }|D_m,\epsilon ,M_j)P(M_j|\mathbf{M})}}\\ =\frac{P_0^{d_i-1}{P}_d}{\sum _{j=1}^n{P}_0^{d_j-1}{P}_d}\hfill \end{array}$$

(34)

where $d_i$ is the number of stages to reach the specific $\epsilon$ in the model class level $M_i$ .

Implementation of the proposed framework

Some critical factors need to be specified appropriately to guarantee ABC-SS performance. Firstly, ABC model selection is limited due to some constraints, particularly the need for sufficient summary statistics to guarantee sufficiency in all considered model classes. This limitation also directly affects posterior approximations. Different signal-processing techniques to extract GW features from original time-domain signals can be considered summary statistics. Even though data extraction can improve the sensitivity of damage identification, some essential information may be missing, and this would result in inaccurate results. One simple solution to overcome this difficulty is to use the entire measured signals. However, this has a critical requirement for simulation accuracy. ⁵² In this paper, benefiting from the HOLM, the laminated composite beam simulation can be accurate and close to the actual responses by taking higher-order variation in displacement expressions as well as employing a more physical delamination simulation mechanism. Euclidean distance is chosen as the distance function for assessing the similarity between numerical and measured GW signals, which is commonly used for comparing two same-dimensional data vectors.

$$\rho (D_s,D_m)=||D_m-D_s|{|}_2=||D_m-M(\theta )|{|}_2$$

(35)

The accuracy and efficiency of the ABC-SS are directly affected by MMA, while the Markov chain speed to the convergence state in MMA mainly depends on the spread of proposal PDFs. Assuming Gaussian proposal PDFs, the variance ${\sigma }^2$ should be properly chosen so that Markov chain samples can have the smallest correlation and statistical efficiency can be maximised. Zuev et al. ⁵⁰ proposed a strategy to adaptively select the optimal ${\sigma }^2$ by monitoring acceptance rate to stay within a suggested range of $[30\%50\%]$ . In a recent study, Vakilzadeh et al. ⁵¹ developed a more advanced strategy to self-regulate the proposal ${\sigma }^2$ for each $n$ chain at each stage, called self-regulating ABC-SS algorithm. This is achieved by learning the statistical properties of candidate samples so that the mean acceptance probability at each stage is always close to the target rate. Additionally, the algorithm provides an additional stopping heuristic rule to stop the calculation if the average acceptance rate falls significantly. In this study, self-regulating ABC-SS is employed to maximise the performance of MMA. This also provides an efficient way to determine the final tolerance and the required simulation stage for evidence estimation in Equation (31). The detailed implementation procedure can be found in the literature.^51,52

The selection of conditional probability $P_0$ also has a significant impact on the algorithm and there is a trade-off between efficiency and accuracy. A smaller $P_0$ requires less simulation stages to reach a fixed tolerance, but it also lowers the posterior distribution quality. Oppositely, a higher $P_0$ enhances the approximation accuracy but demands more simulation stages and computational effort. According to the sensitivity analysis, ⁵⁰ the recommended $P_0$ for optimising the algorithm is within the range of $[0.10.3]$ . Also, $1/P_0$ and $N{P}_0$ is always selected as integers for convenience.

In this paper, the ABC-SS is used to develop a likelihood-free Bayesian framework for multiple-delamination characterisation. A pseudocode implementation of the proposed framework is provided for clear clarification. The details of model class selection and damage parameter identification with evidence evaluation are shown in Figure 3 and Algorithm A, respectively. The HOLM described in section ‘Higher-order laminated model’ is employed to simulate the composite beam with different numbers of delaminations. Once the most plausible delamination number is determined, the corresponding delamination characteristics can be identified. The associated uncertainties can also be quantified by evaluating the posterior distributions.

Algorithm A.

Damage parameter identification using ABC-SS.

1. INPUT: $P_0,N,m_{\mathrm{final}}\&{\epsilon }_{\mathrm{final}}$ .

2. START: generate $[{\theta }_0^k,{D_s}_0^k{]}_{k=1}^N~p(D_s|\theta ,M_i)p(\theta |M_i)$ .

3. For $d=1:m_{\mathrm{final}}$ do

4. | Rearrange $[{\theta }_{d-1}^k,{D_s}_{d-1}^k{]}_{k=1}^N$ so that ${\rho }_{(1)}\le \dots \le {\rho }_{\left(N\right)}.$

5. | Calculate ${\epsilon }_d$ using Eq. (31).

6. | Select $N{P}_0$ seeds that lies in $D_{{\epsilon }_d}$ : $[{\theta }_d^n,{D_s}_d^n{]}_{n=1}^{N{P}_0}=[{\theta }_{d-1}^n,{D_s}_{d-1}^n{]}_{n=1}^{N{P}_0}~p(|D_{{\epsilon }_d})$ .

7. | Generate $1/P_0$ Markov Chain states lying in $D_{{\epsilon }_d}$ using MMA:

8. |  $[({\theta }_d^{n,1},{D_s}_d^{n,1}),\dots ,({\theta }_d^{n,1/P_0},{D_s}_d^{n,1/P_0}){]}_{n=1}^{N{P}_0}~p(|D_{{\epsilon }_d})$ .

9. | Renumber samples: $[({\theta }_d^{n,1},{D_s}_d^{n,1}),\dots ,({\theta }_d^{n,1/P_0},{D_s}_d^{n,1/P_0}){]}_{n=1}^{N{P}_0}=>[{\theta }_d^k,{D_s}_d^k{]}_{k=1}^N$ .

10. | If $({\epsilon }_d\le {\epsilon }_{\mathrm{final}})$ or $(d\le m_{\mathrm{final}})$ or (average acceptance rate drops significantly) then:

11. |  Algorithm END:evaluate the evidence using Eq. (33).

12. | end if.

13. End for.

Figure 3.

Flowchart of Bayesian model assessment for determining the number of delaminations.

Numerical case studies

Details of the case studies

Different numerical case studies were used to investigate the performance of the likelihood-free probability approach and the HOLM for GW-based multiple-delamination identification. A laminated composite beam with 300 mm of length and 2 mm of total thickness was used. The beam consists of eight equal-thickness composite layers with different fibre orientations: [ $0/90/0/90$ ]_s. The material properties of each composite layer are identical and summarised in Table 1. Three delamination properties, including lengthwise location measured from the beam left end, delamination length and layer-wise location measured from the top layer, are treated as unknown for identification. Details of numerical case studies are summarised in Table 2. Regarding layer-wise delamination location, it should be noted that due to the symmetry, the delamination at the first, second and third interface has the same effect at the seventh, sixth and fifth, respectively.

Table 1.

Material properties of each identical layer of the laminated composite beam in numerical study.

$E_1\left(\mathrm{GPa}\right)$	$E_1=E_2\left(\mathrm{GPa}\right)$	$G_{12}=G_{13}\left(\mathrm{GPa}\right)$	G 23 $\left(\mathrm{GPa}\right)$	${\upsilon }_{12}={\upsilon }_{13}$	${\upsilon }_{23}$	$\rho (\mathrm{kg}/{\mathrm{m}}^3$ )	Layer thickness (mm)
128.75	8.35	4.47	2.90	0.33	0.44	1517.00	0.25

Table 2.

Delamination details of damage cases in numerical study.

Case study	Number of delaminations	Delamination lengthwise location $L_i\left(\mathrm{mm}\right)$	Delamination length $D_i\left(\mathrm{mm}\right)$	Delamination layer-wise location $T_i$	Signal-to-noise ratio (%)
Case 1	1	$L_1=100.00$	$D_1=6.00$	$T_1=3/5$	2.5/5.0/7.5
Case 2	2	$L_1=100.00$	$D_1=10.00$	$T_1=4$	2.5
		$L_2=200.00$	$D_2=6.00$	$T_1=3/5$
	3	$L_1=60.00$	$D_1=4.00$	$T_1=2/6$	2.5
		$L_2=120.00$	$D_2=6.00$	$T_1=3/5$
		$L_3=220.00$	$D_3=10.00$	$T_3=4$

Two case studies were considered for different purposes. Case 1 aims at investigating the convergence of the proposed framework under different noise levels, while Case 2 aims at investigating the performance under multiple-delamination situations. The incident was selected as the fundamental asymmetrical mode (A₀) of GW because of its sensitivity to defects. ³¹ The excitation was applied at the top surface of the left end of the beam in the transverse ( $v$ ) direction, while measurement was also captured at the same location and direction. The excitation frequency in all cases was selected as a 60 kHz narrow-band five-cycle sinusoidal tone burst pulse modulated by a Hann window and with an amplitude of 1.0 N.

The HOLM was used for simulating numerical GW signals. Based on an initial convergence study, 15 mm-long SFE was used for beam modelling. ⁴² Eight GLL nodes per element and time step of $0.3\times {10}^{-7}$  s, were also used in the simulation. The total time duration of $4.3\times {10}^{-4}$  s was selected to ensure that all reflections from delaminations can be captured before the right-end boundary reflection. To guarantee simulation accuracy and computational efficiency, at least three HOLM sub-laminates were used to simulate the laminated beam in this study. The number of physical layers in each sub-laminate was automatically adjusted based on the delamination through-thickness locations. Adjacent Layers with no delamination were accommodated in an individual HOLM sub-laminate (see Figure 2(a) for an example).

A commercial Finite element modelling (FEM) software, ABAQUS (Version 6.14, DASSAULT SYSTEMS, Paris, French), was also employed to build laminated models for generating synthetic experimental GW signals. Assuming that the beam is 6 mm wide and 300 mm-long. 3D solid brick elements (C3D8R) were used to simulate composite beams. The excitation and measurement locations were same as the HOLM setting, and the load was acting across the entire width. The delamination region and intact region between layers were simulated using ‘untie’ and ‘tie’ connection constraints at the corresponding locations. The mesh size was 0.2 mm to ensure simulation convergency in all directions, and the time step was automatically defined by the dynamic explicit solver. Measurement error for GW experiment under the laboratory condition was also considered by adding different signal-to-noise levels to the synthetic signal. All captured time-domain signals from both HOLM and ABAQUS were normalised based on the maximum amplitude of signal for comparison.

In terms of prior distributions for each uncertain parameter, assuming early stage of damage, the prior distribution of location $L_i$ and length $D_i$ were independently and uniformly distributed within the range of $[10290\mathrm{mm}]$ and $[220\mathrm{mm}]$ . The prior distribution for location $T_i$ was an integer and has equal probability from 1 to 7. The number of samples $N$ and the conditional probability $P_0$ at each simulation stage was set as 500 and 0.2, respectively. The final tolerance threshold ${\epsilon }_{\mathrm{final}}$ was set as 1.00 × 10⁻³. The maximum simulation stage $m_{\mathrm{final}}$ was chosen as 35 in the numerical study based on some initial trial tests on the case study. The target acceptance rate in self-regulating ABC-SS algorithm for automatically calibrating the variance was set as 44.0%. ⁵¹

Case 1: single delamination with different measurement noise levels

This section investigates the ABC framework performance under various measurement conditions. The noise level from 2.5% to 7.5% is assumed to represent normal to adverse laboratory conditions. Figure 4 shows the time-domain GW signals generated from both the commercial FEM and the HOLM. In Case 1, only one delamination was considered. As can be observed, GW damage reflections can occur twice in the delamination region when the wave enters and exits the damaged area. The amplitudes of damage reflection tend to depend on defect size and also on the layer-wise location. Even though the length in Case 1 is only 6 mm, the reflection energy reflects that A₀ is sensitive to detect delamination in composite. According to Figure 4, there is little shifting error between the two signals due to modelling error (2D modelling from HOLM and 3D modelling from ABAQUS), but overall, there is good agreement between the signals. This simulation illustrates the accuracy of the proposed HOLM.

Figure 4.

Time-domain signals from FEM (in red with solid line) and HOLM (in blue with dot line) at measured location in Case 1.

The tolerance level is determined adaptively in the ABC-SS algorithm. For clear illustration, the distance function evaluation at each stage for M1 in Case 1a is presented in Figure 5(a). Samples were rearranged in ascending order based on the value of the distance function. The dynamic tolerance in each stage ${\epsilon }_d$ was automatically calculated using Equation (31). As the iteration increased, the tendency of distance function set from a ‘curve’ became a ‘flat shape’, which indicated that the sample movement was gradually approaching the global optimum. Additionally, the first $N{P}_0$ number of samples, that is, 100 in this case (see Figure 5(b)), was selected into the MMA as the seeds to generate the new sample set, such that the conditional probability for calculating the evidence is always a constant value.

Figure 5.

(a) Rearranged samples in ascending order based on distance function values at each stage at M1 in Case 1a. (b) Zoom-in details (red rectangle) of the 100 selected seeds at each stage in ABC-SS.

When either of the stopping criteria was reached, the algorithm was terminated, and the required simulation stage and stopping tolerance could be determined. The posterior probability was then calculated from Equation (34). As shown in Table 3, since the first model class has the highest probability, only two model classes, that is, single delamination and two delaminations, were involved in the model class assessment. It was observed that more stages were required to reach the optimum with the increasing noise effect. The stopping tolerance level was also increased due to more uncertainty involved in parameter estimation. Despite that, as highlighted in Table 3, the delamination number was accurately predicted.

Table 3.

Posterior probability with required simulation stage and the stopping tolerance level at different model classes in Cases 1a–1c.

Case 1	1a: 2.5% Noise level		1b: 5.0% Noise level		1c: 7.5% Noise level
Model class	$M1$	$M2$	$M1$	$M2$	$M1$	$M2$
Simulation levels $(d)$	9	21	12	18	14	24
Tolerance level $(\epsilon )$	1.52 × 10−3	2.33 × 10−3	5.22 × 10−3	6.00 × 10−3	1.14 × 10−2	1.51 × 10−2
Probability	1	0	1	0	1	0

As shown in Table 4, delamination properties in Cases 1a–1c were accurately identified. Since the GW reflection was affected by a higher noise level, a larger uncertainty level was quantified in the parameters, especially in the delamination length. The error percentage between the prediction and the exact value was increased from 3.17% to 4.17%, as well as the coefficient of variation (COV) value was increased from 0.02 to 0.29.

Table 4.

Identified results of delamination properties in Cases 1a–1c.

Case 1	Delamination location $\left(\mathrm{mm}\right)$	Delamination length $\left(\mathrm{mm}\right)$	Delamination though-thickness location (layer-wise)
	$L_1:\mathrm{mean}[\mathrm{error}\%](\mathrm{c}.\mathrm{o}.\mathrm{v}.\%)$	$D_1:\mathrm{mean}[\mathrm{error}\%](\mathrm{c}.\mathrm{o}.\mathrm{v}.\%)$	$T_1$
1a	$100.00[0.00](0.00)$	$6.19[3.17](0.02)$	$3$
1b	$99.99[0.01](0.00)$	$6.23[3.83](0.17)$	$3$
1c	$99.98[0.02](0.00)$	$6.25[4.17](0.29)$	$5$

Similar behaviours were observed in the scatterplots of the 500 samples (see Figure 6). These samples transitioned from the prior distribution and gradually to the posterior distribution. As the noise level increased from 2.5% to 7.5%, there was a relatively larger region for the global optimum, particularly for the length parameter. However, the overall quantified uncertainty in each case remained acceptable and within expectations. This was attributed to the large A₀ GW reflection energy from delamination, which reduced the impact of noise effects.

Figure 6.

2D scatterplot of the 500 samples (delamination lengthwise location $L_1$ vs length $D_1$ ) from the initial simulation stage to the final stage in (a) Case 1a, (b) Case 1b and (c) Case 1c.

Case 2: multiple delaminations with 2.5% of measurement noise level

This section further examines the proposed framework performance under multiple-delamination situations. The number of simulation stages and the final stopping tolerance levels in Cases 2a–2b are summarised in Table 5. Given the increased complexity of the GW signal due to additional damage reflections from multiple delaminations, it was expected that more simulation stages would be necessary than Case 1 to calibrate the unknowns and minimise discrepancies between simulations. More model classes were also involved in both cases until the highest posterior probability was determined. For clear illustration, Figure 7 displays the ABC model evidence calculated from Equation (33) in Case 2a against the dynamic tolerance level, along with the comparison of posterior probabilities calculated from Equation (34).

Table 5.

Posterior probability with required simulation stage and the stopping tolerance level at different model classes in Cases 2a–2b.

Case 2	2a: Two delaminations case			2b: Three delaminations case
Model class	$M1$	$M2$	$M3$	$M1$	$M2$	$M3$	$M4$
Simulation levels $(d)$	13	25	29	10	25	20	35 (Max reached)
Tolerance level $(\epsilon )$	1.31 × 10−1	9.44 × 10−3	1.41 × 10−2	1.25 × 10−1	3.08 × 10−2	1.76 × 10−2	—
Probability	0	1	0	0	0	1	—

Figure 7.

(a) ABC evidence estimation and (b) the posterior probability of all three considered model classes against the dynamic tolerance level in Case 2a from the initial to the final simulation stage.

Referring to the information-theoretic point of view for log evidence, the posterior probability computation considers a trade-off between two terms. The first term measures the posterior average data fit of the model class. The second term is defined as the relative entropy between the prior and posterior distributions, and it measures the extracted information about the model class complexity.^53,54 In short, an optimal model class tends to have an acceptable average data-fitting accuracy while extracting less information about its model parameters. It was found that posterior probability estimation using ABC-SS follows a similar concept to define the preferred model. Specifically, the tolerance level determined from the distance function provides a similar way to measure the goodness of fit. At the same time, the number of simulation stages reflects the overall penalty on the complexity of the model class, that is, more simulation stages are required in a more complex damage scenario.

Following the aforementioned rule, it was observed from Table 5 and Figure 7(a) that M1 in Case 2a has the smallest complexity penalty, as it requires the fewest simulation numbers to reach the final stage. The phenomenon makes sense as M1 only considers a single-delamination scenario. Nevertheless, M1 was unable to achieve good data-fitting accuracy, such that the algorithm is stopped at a relatively large tolerance level compared to other considered model classes. The posterior probability of M1 dropped and eventually became zero as the dynamic tolerance level decreased (see Figure 7(b)). In contrast, M2 was accurately selected as the optimal class due to achieving the best data-fitting accuracy (i.e., the smallest tolerance) while only extracting less information about the model class complexity than M3. In a similar manner, the number of delaminations in Case 2b was also accurately detected as M3, that is, three delaminations. It should be noted that M4 was excluded from model class assessment, as the algorithm was stopped at the maximum stage.

The uncertain parameters in both cases are summarised in Table 6. In terms of quantified uncertainty, the error percentage and COV values increased compared to Case 1. This is because the GW reflection signal was more complicated, and more reflections were generated from defects. This led to more noise effects or uncertainties involved in damage identification. Nevertheless, all delamination parameters were accurately identified with acceptable errors. The time-domain signals from FEM and HOLM (using the mean value of identified parameters in Table 6) are also presented in Figure 8 for comparison. Some reflection pulses, as observed were overlapped in an irregular wave shape, containing combined damage information from different delaminations. There is overall good agreement between the two simulations, indicating accurate multiple-delamination identification.

Table 6.

Identified results of delamination properties in Cases 2a–2b.

Case 2	Delamination location $\left(\mathrm{mm}\right)$	Delamination length $\left(\mathrm{mm}\right)$	Delamination through-thickness location (layer-wise)
	$L_i:\mathrm{mean}[\mathrm{error}\%](\mathrm{c}.\mathrm{o}.\mathrm{v}.\%)$	$D_i:\mathrm{mean}[\mathrm{error}\%](\mathrm{c}.\mathrm{o}.\mathrm{v}.\%)$	$T_i$
2a	$L_1=100.50[0.50](0.04)$	$D_1=10.64[6.40](0.40)$	$T_1=4$
	$L_2=201.14[0.57](0.07)$	$D_2=6.18[3.00](1.17)$	$T_2=5$
2b	$L_1=59.70[0.50](0.19)$	$D_1=3.98[0.50](0.03)$	$T_1=2$
	$L_2=119.98[0.03](0.04)$	$D_2=6.20[3.33](0.98)$	$T_3=5$
	$L_3=220.52[0.24](0.02)$	$D_3=10.32[3.18](0.88)$	$T_3=4$

Figure 8.

Time-domain signals from FEM (in red with solid line) and HOLM (in blue with dot line) at measured location in (a) Case 2a and (b) Case 2b.

Figure 9 presents the histograms for delamination lengthwise locations and lengths at the final stage in Case 2b for illustration. The posteriors were distributed close to the exact value in each parameter. It was also observed that the posteriors approximated by ABC algorithm showed different types of statistical distribution, reflecting that uncertainty error has different impact on different parameters. The associated uncertainties can be quantified by evaluating distribution shapes. The histogram range was consistent with the corresponding error percentage and COV value shown in the summary Table 6.

Figure 9.

Histograms for all delamination lengthwise locations $L$ and lengths $D$ at the final stage in Case 2b.

Experimental study

Experiment setup

Two fibre-reinforced polymer beams with 6 mm width were utilised in the experimental study to validate and investigate the practicability of the proposed likelihood-free framework. The Beams were made by eight HexPly M21/IM7 (Hexcel) carbon and epoxy physical layers with an assembling orientation of $[0/90/0/90]$ _s. Table 7 summarises the material properties. The delaminations were artificially generated during the manufacturing process by inserting a thin Teflon film between two laminates at the designed lengthwise location.

Table 7.

Material properties of each identical layer in experimental study.

$E_1\left(\mathrm{GPa}\right)$	$E_1=E_2\left(\mathrm{GPa}\right)$	$G_{12}=G_{13}\left(\mathrm{GPa}\right)$	G 23 $\left(\mathrm{GPa}\right)$	${\upsilon }_{12}={\upsilon }_{13}$	${\upsilon }_{23}$	$\rho (\mathrm{kg}/{\mathrm{m}}^3$ )	Layer thickness (mm)
160.00	8.50	4.20	2.70	0.35	0.53	1580.00	0.184

A similar experimental setup as He and Ng ³² was conducted in this study. Specifically, both ends of the beam were fixed using rigid clamps. The length of the beam between the two clamps is 300 mm. A plate piezoceramic transducer (PZT) with dimension $6\times 6\times 3$  mm³, was bonded to the left end surface using conductive glue to generate A₀ GW excitation. The selection of excitation location was motived by the goal of maximising the wave propagation distance before the incident wave reaches the beam ends and minimising the effect from boundary reflections. To enlarge the excitability, a metal backing mass with size $6\times 6\times 4$  mm³ was also bonded on the top of PZT. Because of the best signal-to-noise ratio in the experiment, the excitation frequency was selected as a 50 kHz narrow-band six-cycle sinusoidal tone burst pulse modulated by a Hann window.

The excitation was generated using a signal generator NI PXIE-5412 and then amplified by a high-voltage amplifier before applying to the PZT. The transverse displacement $v$ at the measured point was measured using a laser Doppler vibrometer (1D) Polytec PSV 400 with a controlling system OFV5000. The measured location was selected 10 mm away from the PZT due to its relatively flat surface, which could provide good A₀ observation quality during measurement. Five hundred times of signal averaging and band-pass filter were also applied to further reduce environmental noise effect in the GW signals. The same excitation and measured location were applied in the HOLM simulation. Figure 10 shows the detailed flowchart of experimental setup. Details of delamination in both composite beams can be found in Table 8 and Figure 11.

Figure 10.

Detailed experimental setup of GW excitation using transducer and measurement using 1D laser.

Table 8.

Details of delaminations in experimental study.

Case 3	No. delaminations	Delamination lengthwise location $L_i\left(\mathrm{mm}\right)$	Delamination length $D_i\left(\mathrm{mm}\right)$	Delamination layer-wise location $T_i\left(\mathrm{mm}\right)$
Beam 1	1	$L_1=100.00\pm 10.00$	$D_1=6.00\pm 1.00$	$T_1=3/5$
Beam 2	2	$L_1=100.00\pm 10.00$	$D_1=10.00\pm 1.00$	$T_1=4$
		$L_2=200.00\pm 10.00$	$D_2=6.00\pm 1.00$	$T_1=3/5$

Figure 11.

Details of (a) PZT installation and measurement location, (b) delamination locations in beam 1 and beam 2, (c) zoom-in detail of 6 mm delamination in beam 1 and (d) zoom-in details of 10 mm and 6 mm delamination in beam 2.

Results and discussions

The same HOLM and ABC-SS algorithm settings in numerical study were employed in both experimental studies. Table 9 summarises the number of simulation levels and the corresponding final tolerance level for both experimental cases. Different model classes were automatically considered based on the proposed flowchart in Figure 3 until the highest probability was obtained. M1 in Beam 1 and M2 in Beam 2 were selected as the most plausible model due to the best goodness of fit (i.e., the smallest stopping tolerance) and less model complexity penalty (i.e., less simulation stages). The number of defects in both cases was correctly identified as one delamination in Beam 1 and two delaminations in Beam 2.

Table 9.

Posterior probability with required simulation stage and the stopping tolerance level at different model classes in beams 1 and 2.

Case 3	Beam 1		Beam 2
Model class	$M1$	$M2$	$M1$	$M2$	$M3$
Simulation levels $(d)$	11	27	16	25	30
Tolerance level $(\epsilon )$	6.99 × 10−3	2.24 × 10−2	3.73 × 10−1	1.41 × 10−1	2.24 × 10−1
Probability	1	0	0	1	0

A summary of the identified delamination properties can be found in Table 10. The results are presented as sample mean at the final stage. The corresponding percentages of error and the sample COV are also included. Although the layer-wise location in Beam 2 differed by one layer from the actual location, most of the delamination properties were identified accurately. A comparison of the time-domain GW signal for Beam 1 between the measurement and the simulation (using the mean value of identified parameters in Table 10) is presented in Figure 12 for illustration. In addition to the excitation pulse, two GW damage reflections due to the wave entering and exiting the delamination zone were clearly captured in the measurement. As observed, both GW signals from simulation and measurement were highly matched, while the uncertainty level was also well-controlled during measurement. The quantified uncertainties, as shown in Table 10, were acceptable and within expectations. The accuracy of multiple-delamination identification in both experimental studies further validated and demonstrated the practicability of the proposed likelihood-free framework.

Table 10.

Identified results of delamination properties in beams 1 and 2.

Case 3	Delamination location $\left(\mathrm{mm}\right)$	Delamination length $\left(\mathrm{mm}\right)$	Delamination through-thickness location (layer-wise)
	$L_i:\mathrm{mean}[\mathrm{error}\%](\mathrm{c}.\mathrm{o}.\mathrm{v}.\%)$	$D_i:\mathrm{mean}[\mathrm{error}\%](\mathrm{c}.\mathrm{o}.\mathrm{v}.\%)$	$T_i$
Beam 1	$L_1=99.97[0.03](0.02)$	$D_1=6.37[6.01](0.33)$	$T_1=3$
Beam 2	$L_1=100.23[0.23](0.06)$	$D_1=10.69[6.94](0.73)$	$T_1=3$
	$L_2=198.65[0.68](0.03)$	$D_2=6.19[3.17](1.34)$	$T_2=5$

Figure 12.

Time-domain signals from experiment (in red with solid line) and HOLM (in blue with dot line) at measured location in Beam 1.

Figure 13 shows the 2D scatterplots for the delamination lengthwise locations and lengths in Beam 2 from the initial stage to the final stage. With the number of iterations increasing, 500 samples were moved from a wide initial range to a small search area and eventually arrived at the global optimum. The size of the optimum region mainly depends on the uncertainty effects of the unknowns from both measurement and simulation. It was observed that the final scattering areas in both delamination lengths were relatively larger than the ones in both locations, demonstrating that more uncertainties were involved in the length parameter estimation procedure. But overall, the uncertainties were well-controlled, which allowed the optimisation to identify delamination parameters precisely within acceptable error margins and refine global optimums.

Figure 13.

2D scatterplot of the 500 samples (delamination lengthwise locations vs lengths) (a) $L_1$ versus $D_1$ and (b) $L_2$ versus $D_2$ from the initial simulation stage to the final stage in Beam 2.

The posterior distributions at the final stage in Beams 1–2 are also shown in Figure 14 using histogram plots for better observation. The posterior samples were well-established around the actual value within acceptable error margins. By evaluating the statistical properties of each posterior distribution, the associated uncertainties were quantified, and also the error information provided reliability of the identified results. Figures 13 and 14 further demonstrated the accuracy and practicability of the proposed framework for multiple delaminations in composite in laboratory practical application.

Figure 14.

Histograms for all delamination lengthwise locations $L$ and lengths $D$ at the final stage in (a) Beam 1 and (b) Beam 2.

Conclusion

This paper proposed a GW model-based approach for characterising multiple delaminations in composite beams. The proposed framework incorporates a likelihood-free method, the ABC-SS algorithm, which provides a straightforward and practical alternative to Bayesian model selection and parameter inference without assuming any likelihood functions. Model evidence estimation of a model class is a simple by-product of ABC-SS algorithm. This provides a clear understanding of selecting the most plausible model using ABC approach. Instead of using any GW signal-processing techniques, damage identification in this study relies on the entire GW time-domain signal. A HOLM with a sub-laminate scheme is used to enhance the simulation efficiency and accuracy of GW propagation in composite beams. A cubic variation of in-plane displacement and a quadratic variation of transverse displacement expression over the beam thickness were employed. This allows delamination to be simulated by simply connecting or disconnecting the corresponding nodes at the top and bottom surfaces of the adjacent sub-laminates. The numerical studies have comprehensively demonstrated the procedures to identify multiple delaminations using ABC model selection and parameter estimation. The HOLM has also proven that precisely simulate GW propagation behaviours from delaminations compared to traditional FE methods, minimising the impact of modelling uncertainties on identification results. Two experimental studies have further verified the accuracy and practicability of proposed ABC framework in laboratory applications.

The current approach has great potential for identifying delaminations in any composite structures. However, different stacking materials may lead to different GW damage responses and more complex reflection signals. More associated uncertainties in real situations can also be involved compared to controlled laboratory conditions. The presented HOLM is capable of simulating any composite beam by applying relevant material properties and orientations. Further research is necessary to examine the proposed ABC framework in real applications. The current study only considers Euclidean distance as a similarity measure, while other distance functions could be explored to account for various aspects of captured wave signals. Further research is needed to examine the impact of ABC damage identification on the selection of different distance functions. The proposed framework can also be used with any numerical model to identify specific defect types of interest. In addition to defining the number of defects, ABC model selection can be used to classify defect types if sufficient damage information is available from measurement. Further research is worthwhile to enhance the applicability and practicability of the ABC framework in GW model-based damage identification.