Boundary element modelling of ultrasonic Lamb waves for structural health monitoring

In this paper, a novel boundary element plate formulation is proposed to model ultrasonic Lamb waves in both pristine and cracked plates for structural health monitoring (SHM) applications. Lamb waves are generated and sensed by piezoelectric discs. An equivalent pin-force model is newly proposed to represent the actuation effect of piezoelectric discs, which is more accurate than the classical pin-force model. The boundary element formulation is presented in the Laplace-transform domain based on plate theories, which allows three-dimensional analysis of Lamb wave behaviours, such as propagation and interaction with cracks, in thin-walled structures. A damage detection algorithm is used for crack localization alongside the BEM-simulated data. The BEM solutions show excellent agreement with both 3D finite element simulation and experimental results.


Introduction
Ultrasonic Lamb waves have been widely used for structural health monitoring (SHM) over recent decades. Recently, small-sized transducers made of piezoelectric Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. lead-zirconate-titanate (PZT), which are lightweight, inexpensive and could be bonded onto structural surfaces permanently, allow actuating and sensing Lamb waves conveniently [1]. The piezoelectric-mechanical coupled dynamic system can be modelled by the finite element method (FEM) [2], spectral element method [3] and local interaction simulation approach (LISA) [4]. A review of different numerical methods can be found in [5] .
The boundary element method (BEM) is a very promising alternative and has been demonstrated to be efficient and numerically stable for the health monitoring of threedimensional structures using piezoelectric transducers [6]. In contrast to domain-discretized methods, boundary element meshing is only limited to the boundary of a domain with the help of fundamental solutions and as a result significant reduction in the degrees of freedom is achieved. This advantage has also been shown in the modelling of ultrasonic Lamb waves [7], which is particularly attractive for high-frequency cases.
Furthermore, when various damage scenarios need to be simulated, only the boundary of damage, rather than all domain boundaries, need to be remeshed, resulting in less data preparation. Further discussions regarding the advantages and disadvantages of the BEM in modelling Lamb waves are given in reference [7]. Since the three-dimensional BEM used in reference [6] suffers some numerical problems encountered in modelling thin-walled structures, the BEM formulations for plate theories are preferred in terms of modelling Lamb waves, in which only plate edges are required to be meshed using 1D line elements. When it comes to the modelling of cracks, the well-known dual boundary element method (DBEM) [8] is more suitable and allows treating the crack surfaces as extra boundaries. The DBEM provides a general framework, which is applicable to both isotropic and anisotropic materials, and to both static and dynamic problems.
In order to incorporate the PZT transducers into the BEM model, a direct strategy is to model the transducers using a numerical method, such as a semi-analytical finite element method [6], and to solve the system in a fully coupled manner. However, when it comes to the boundary element plate formulations, this introduces additional plate-domain points to couple the transducers to the plate, which significantly deteriorate the efficiency of the boundary-element-only formulations. This drawback can be avoided by replacing the actuator by prescribed tractions obtained directly from the normal and shear stresses at the actuator-plate interface. The tractions are able to be calculated using analytical or numerical methods. When it comes to the analytical models, perhaps the simplest replacement is the pin-force model [9,10], which considers that the shear traction is concentrated at the edge of the actuator assuming the perfect bonding between the actuator and the plate. Taking the effect of the adhesive layer into account, Crawley et al [11] and Roy et al [12] gave the analytical expression for the distributed shear traction generated by rectangular and circular actuators, respectively. However, they neglected the dynamics of the system, and therefore their models are not suitable for high-frequency excitations. Huang et al [13] and Kapuria et al [14] improved the analytical model for rectangular actuators by including the inertia term in the analysis. After the improvement, some limitations still remain. First, these dynamic models cannot be applied to circular actuators. Second, the normal traction at the interface, which has been shown to be as important as the shear traction near the resonant frequencies of the actuator [15], was disregarded. Thus, in order to deal with general actuation scenarios and obtain more accurate traction distribution, numerical methods, such as the FEM, should be relied on. Such distributed tractions can be directly inserted into the boundary element plate formulations. However, this direct method is not numerically efficient. In order to improve the efficiency, an equivalent pin-force model, in which the distributed tractions are simplified into line tractions, is proposed in this study. As for the PZT sensor, one common model [10,12,[16][17][18], in which the sensor is treated as a capacitor, can be easily integrated into the proposed BEM model.
In this paper, Laplace-transform boundary element formulations based on Kane-Mindlin [7,19] and Mindlin [20,21] plate theories are used to model the S0 and A0 Lamb modes in an isotropic plate below the cut-off frequency of the A1 mode, respectively. The fundamental solutions in reference [20] are modified slightly to represent the dispersion characteristic of the A0 mode more accurately and are rewritten in a more compact form. Cracks are modelled using the DBEM based on the above-mentioned plate theories [22,23]. An equivalent pin-force model is newly proposed to represent the PZT actuator, and an existing sensor model is implemented numerically. The accuracy of the proposed BEM-based model is validated using 3D FEM and experimental results. In addition, this model provides data for a damage detection algorithm to localize the crack, and the accuracy of this data is checked against the experimental data. Finally, some parametric studies are carried out using this BEM model.

Boundary integral formulations
Throughout this paper, Greek indices (α, β, γ) take the values 1 to 2, while Latin indices (i, j, k) run from 1 to 3; repeated indices are summed (summation convention) except where otherwise indicated.
The mid-plane of a linearly elastic isotropic plate is located at the x 1 − x 2 plane, and the x 3 -axis is the normal to the midplane. The displacements can be approximated by taking the constant and linear terms from their Taylor series about x 3 = 0: in which the superscripts S and A denote the motions symmetric and antisymmetric with respect to the mid-plane of the plate, respectively; t denotes time. Based on this assumption, Mindlin and Kane [24,25] have given the governing equations for both symmetric and antisymmetric cases: ; the mass density, Young's modulus, Poisson's ratio and plate thickness are denoted by ρ, E, ν and h, respectively; κ S and κ A are adjustable shear factors; the generalized surface tractions are given by: where δ ij is the Kronecker delta and the summation convention is suspended with regard to i. The Laplace transform with respect to time can be expressed as: Applying equation (5) to the above governing equations, Equations (2) and (3) can be rewritten in the following form: On the basis of Betti's reciprocal work theorem, the displacements at the point X ′ in a plate domain Ω can be calculated as follows:ũ where X and x are the field points in the domain Ω and on the boundary Γ respectively;Ũ Y ij andT Y ij are dynamic fundamental solutions listed in appendix A; the superscript Y denotes S and A for the symmetric and antisymmetric motions respectively; the tractions are defined as: in which the summation convention only applies to the subscript β; n β denotes the components of a unit outward normal vector to the plate boundary.
When the domain point X ′ is moved to the boundary Γ, denoted by x ′ , the well-known displacement boundary integral equation (DBIE) can be derived using equation (8) and written as: where − stands for Cauchy principal value integral and the free term c ij is equal to δ ij /2 when x ′ is located on a smooth boundary. In order to model crack problems in a more efficient way, the dual bonudary element method (DBEM) is adopted. In the DBEM, the DBIE is used for the modelling of one crack surface, while the other crack surface is modelled using the traction boundary integral equation (TBIE). Utilizing the DBIE and the stress-displacement relationship, the TBIE on a smooth boundary is obtained: where = stands for Hadamard finite-part integral and the dynamic fundamental solutionsD Y iβk ,S Y iβk are given in appendix A.
In this paper, the plate edges are divided into continuous quadratic isoparametric line elements. When it comes to the crack surfaces, discontinuous quadratic line elements are adopted in order to ensure the existence of the Hadamard finite-part integrals in the TBIE. The shape functions and nodal locations of these elements can be found in reference [8]. After collocating at every nodal point and introducing boundary conditions, the discretized boundary integral equations can be written in matrix form as: in whichÃ Y is the coefficient matrix (the superscript Y denotes S and A for the symmetric and antisymmetric motions respectively);x Y is comprised of unknown nodal displacements and tractions;f Y is the known vector obtained by using the boundary conditions;F Y is due to the effect of the generalized surface tractions (domain loads). More details regarding the numerical implementation can be found in reference [26]. The timedomain solutions are computed using the numerical inversion of Laplace transforms, which is given in appendix B.

Wave dispersion characteristics
Since the proposed formulations are based on plate theories, which are approximations of the 3D theory of elasticity, the guided waves described in the Mindlin and Kane-Mindlin (K-M) theories do not represent all Lamb wave modes. Thus, the applicable range of the proposed formulations should be known before the numerical analysis is carried out. It is shown from equations (A16) and (A37) that three symmetric and three antisymmetric waves can be modelled in the K-M theory and Mindlin theory respectively. The dispersion relations of the above six waves are given explicitly in appendix C. In the K-M theory, the three waves are similar to the S0, S1 and SH0 modes, while A0, A1 and SH1 modes are represented approximately using the three waves in the Mindlin theory [27]. The exact dispersion relations of Lamb waves can be calculated numerically from the well-known Rayleigh-Lamb frequency equation [28]. In order to achieve the desired approximation to the exact dispersion relations, the shear factors κ S and κ A in the K-M and Mindlin theories should be adjusted.
In most SHM applications, only S0 and A0 modes are excited and the operating frequency is below the cut-off frequency of the A1 mode [29,30]. If this frequency range is required, it has been found that the satisfactory approximation can be achieved when κ S = 0.73 and κ A = 0.94 for the aluminium plate. Such shear factors were determined using the method given in reference [7]. The group-velocity dispersion curves of the Lamb wave modes (S0, S1, A0 and A1) on the basis of the plate theories and those obtained by solving the Rayleigh-Lamb equation are compared in figure 1. As can be seen on this figure, the K-M and Mindlin plate theories can be used to model the S0 and A0 modes with excellent accuracy at frequencies below the cut-off frequency of the A1 mode.

Actuator model
A new PZT actuator model, named an equivalent pin-force model, is developed here. First, the tractions applied to the host plate by the PZT actuator are obtained from the finite element simulation. Then, such distributed tractions are projected to the excited S0 and A0 modes using an integral-transform technique. Lastly, an equivalent pin-force, which leads to the same projection, is calculated.
The commercial finite element software ABAQUS is used to get the applied tractions in the frequency domain, and the time-domain results are obtained through the inverse Fourier transform. Since the PZT discs are polarized along the thickness direction and the material properties of the plate and adhesive layer are isotropic, the actuation problem can be solved using a 2D axisymmetric model (see figure 2). In the finite element analysis, only a small portion of the host plate is modelled, including an artificial damping zone at the end of the plate section to avoid wave reflections. The damping coefficients in this zone can be determined using the method given in reference [31].
The distributed tractions applied to the plate due to the actuator are equal to the normal and shear stresses at the adhesive-plate interface. In the equivalent pin-force model,  the distributed shear traction, τ r (r, ω), and normal traction, τ z (r, ω), are replaced by simple concentrated line tractions at r = R, τ r|EPF (R, ω)δ (r − R) and τ z|EPF (R, ω)δ (r − R), in which δ denotes the Dirac delta function. Thanks to this simplification, the corresponding 2D surface integral resulting from the actuator in the boundary integral equations is reduced to 1D line integral.
When the plate is subjected to the distributed tractions τ r and τ z , only S0 and A0 modes are generated at the frequency used in this study. According to Quaegebeur et al [15], the analytical solutions for the displacements of the S0 mode, u S r , where Y denotes S and A for the S0 mode and A0 mode, respectively; ω and k Y represent frequency and wavenumber respectively; i denotes the imaginary unit; H 1 is the Hankel function of the second kind of order one;τ Y r andτ Y z are the projections of the distributed tractions on the S0 and A0 modes, which are computed by the Hankel transforms of order one and zero (with first-and zeroth-order Bessel functions J 1 and the other terms involving the wavenumber are given by: in which It can be seen from equations (13)-(15) that the influence of the distributed tractions on the Lamb wave response is determined solely by the projectionsτ Y r andτ Y z . In other words, two different sets of tractions lead to the same response of a specific Lamb wave mode, provided that they have the same projections on this mode. Thus, if the projections of the equivalent line tractions are equal to those of the distributed tractions, these two sets of tractions are able to be interchangeable.
Assuming the equivalent line tractions τ Y r|EPF and τ Y z|EPF are applied at r = R, the projections on the S0 and A0 modes are given by : By letting the right-hand side of equation (16) equal to the right-hand side of equation (14), the frequency components of the unknown equivalent line tractions can be obtained as follows: If the frequency spectrum of the excitation voltage V(t) is denoted by V(ω), the time-domain equivalent line tractions can be computed by taking the inverse Fourier transform of After substituting the equivalent line tractions into equation (4), the loads due to the PZT actuator are able to be integrated into the boundary integral equations and the resulting line integral is estimated using the composite trapezoidal rule, which is the same as the rule for θ in equation (19). Finally, the actuation effect of the PZT actuator is represented by the domain-load termF Y in equation (12).
It should be noted that, unlike the conventional pin-force model (see reference [10]), the proposed equivalent line tractions are not necessarily applied at the circumference of the actuator disc, and furthermore the radius R to which the equivalent line tractions are applied is not unique. For different R, the amplitudes of such tractions, which are determined by equation (17), are different. Since these different equivalent line tractions lead to the same projections, they have the same actuation effect on the corresponding Lamb wave mode. When the numerical difficulty, caused by the zeros of the Bessel functions in the denominators of equation (17), is encountered, it can be remedied by choosing R which is smaller than the radius of the actual actuator disc.

Sensor model
One common PZT sensor model [10,12,[16][17][18]] is adopted in this paper. In such model, the output voltage of the sensor, which is directly proportional to the total in-plane strains in the sensing area, is estimated using the following equation: (18) in which E pzt and ν pzt denote the in-plane Young's modulus and Poisson's ratio of the PZT; h pzt and A pzt are sensor thickness and the area of the sensor bottom planar surface respectively; d 31 and ε 33 are the piezoelectric coefficient and out-ofplane dielectric permittivity respectively;ũ α,α represents the total in-plane strains at the sensor-plate interface.
Since the integration domain for the sensor disc is circular, the integral in equation (18) is computed numerically in polar coordinates by using the quadrature rule suggested in reference [32]. A composite trapezoidal rule is applied for θ, in which the interval [0,2π] is subdivided into N equal intervals, and the m-point Gauss quadrature rule is applied for r. Thus, the integral can be estimated as follows: where R s denotes the radius of the PZT sensor; w m and r m are the weights and nodes of Gaussian quadrature over the domain [0, R s ], respectively. Differentiating equation (1) with respect to x α in the Laplace domain, the total in-plane strains at the sensor-plate interface can be calculated as follows: After obtaining all nodal displacements and tractions by solving equation (12), the strain components for the symmetric and antisymmetric motions are computed directly using the following integral equations: in which Y denotes S for the symmetric motion and A for the antisymmetric motion;Ũ Y j andT Y j are fundamental solutions listed in appendix A.

Specimen description
A 2-mm-thick square plate is considered in this study. The plate is made of aluminium which has a Young's modulus of 69 GPa, a Poisson's ratio of 0.33 and a density of 2700 kg m −3 . There is a through-thickness crack located at the center of the plate. Four PZT disc transducers, 0.5 mm thick, are attached to the top surface of the plate. The plate dimensions and layout of PZT transducers are illustrated in figure 3. Figure 4 shows the specimen for experiments. A throughthickness slit of 0.4 mm width is machined with the help of electrical discharge machining (EDM) to represent the crack. The PZTs are made of PIC 255 produced by PI Ceramic GmbH. The material properties of PIC 255 are given in table 1. The PZT discs were bonded to the plate using thermoplastic adhesive films.

Actuation signal and frequency selection
A sinusoidal tone burst modulated by Hann window is chosen as the excitation signal because its limited bandwidth can alleviate the dispersion effect of Lamb waves. Thus, the actuation voltage V (t) has the following form : where f c and N w denote the central frequency and the number of cycles respectively; V m represents the amplitude and H is the Heaviside function.
In the SHM applications, it is desirable to generate a pure Lamb wave mode to simplify the signal analysis. For the disc actuators used in this study, both the S0 and A0 modes were excited. The A0 and S0 modes dominate the sensor response around 50 kHz and 300 kHz, respectively. Since the A0 mode is highly dispersive at 50 kHz, the actuator is excited by a Hann-windowed tone burst centred at 300 kHz with 5 cycles.

Experimental setup
The experimental setup shown in figure 5 is mainly composed of a LabVIEW-based control panel, a waveform generator card and a data acquisition card. An arbitrary wave generator National Instruments (NI) PXI-5412 was used to generate the actuation signal. The response signals of the PZT sensors were recorded with a NI PXI-5105 digitizer card at the sampling frequency of 60 MHz.

Sensor and plate responses
The proposed BEM model is used to generate sensor signals in the plate. The generated BEM signals are compared with those obtained from 3D finite element simulation as well as experimental measurements. It is worth mentioning that the thermoplastic adhesive film should be included in the numerical simulation for accurate modelling when the excitation frequency is around the local resonant frequency of the coupled system consisting of the PZT disc, adhesive film and host plate. Willberg et al [33] have found that the adhesive film is able to shift this resonance frequency and it has a great influence on the system response amplitude at the frequencies near the resonance. Hence, although the film does not change the group velocity of the Lamb-wave packet propagating in the host plate, it affects the shape of the wave packet.
In order to take the adhesive film into consideration in the numerical modelling, the material properties and thickness of the cured thermoplastic film are essential. However, it is difficult to accurately measure these parameters. Instead, they can be determined by calibrating the numerical model with experimental data. Kapuria et al [14] pointed out that the film stiffness, rather than its density, dominates the stress-transfer effect between the PZT disc and the host plate. Thus, in the calibration process, the stiffness is adjustable while the density remains the same. Furthermore, according to reference [15], the stiffness is not sensitive to the Poisson's ratio, and the thickness and Young's modulus of the film contribute the most to the stiffness. Since the thickness is easier to measure than the Young's modulus, the Young's modulus is recommended to be the sole parameter for the calibration. The value of the calibrated Young's modulus is able to be determined by comparing the shape of the direct-arrival wave packet in the simulated sensor signal with that in the experimental signal. The density and Poisson's ratio can be found in the data sheet provided by the manufacturers. The errors due to these two parameters, as well as the measurement error of the thickness, are not important because the calibrated Young's modulus can compensate for them. After the calibration, the final value of the Young's modulus for the numerical simulation does not necessarily represent the actual Young's modulus of the adhesive film. It is more appropriate to treat it as a generalised elastic modulus which reflects the adhesive effect on the system response, and this value could vary in different numerical models. In the 3D FEM simulation, both actuator adhesive and sensor adhesive are modelled, while only actuator adhesive is taken into account in the BEM model. Thus, the Young's modulus used in the BEM, which should include the contribution from the sensor adhesive, is different from that input to the 3D finite element simulation. In the following numerical simulations, the adhesive film has a density of 1200 kg m −3 , a Poisson's ratio of 0.44 and a thickness of 20 m. The Young's moduli are 1.2 Gpa and 2 Gpa for the BEM and 3D FEM models, respectively.
The equivalent pin-force model in the boundary element analysis requires the tractions which are applied to the host plate by the PZT actuator. These tractions are able to be obtained using a 2D axisymmetric analysis conducted with ABAQUS. As shown in figure 2, it is only necessary to model a small portion of the plate (l p = l d = 50 mm). In the damping zone, the mass proportional damping coefficient has the value α R (x r ) = 11 × 10 6 × [(x r − l p ) /l d ] 3 . The plate and the adhesive layer were discretized with 5000 and 200 elements (CAX8R), respectively, while 75 piezoelectric elements (CAX8RE) were used to model the actuator. The connection of these three components was achieved by tie constraints.
The convergence properties of the Laplace-transform boundary element method are shown in figure 6. The term 'Path 14' refers to that PZT1 (the first number '1') is the actuator and PZT4 (the second number '4') is the sensor. The scattered signal is obtained by subtracting the response of the sensor on the intact plate from that on the cracked plate. Satisfactory results can be obtained when each plate edge was divided into 50 boundary elements and the crack was discretized using 12 elements (6 elements along each crack surface). In addition, 40 Laplace terms are sufficient for the convergence. Thus, the BEM solutions based on these parameters are used in the following analysis. The proposed BEM model has low memory requirements. The amount of memory required to store the global coefficient matrix (double precision) in the discretized boundary integral equation (12) is 88 Mb and 105 Mb for the intactplate and crack-plate cases, respectively. Furthermore, if the sparsity of this matrix is utilized, such memory requirement can be halved. When it comes to large-scale problems, such as assembled plate structures, the parallel implementation of this BEM can be achieved in a simple way. Since both the creation and the solution of equation (12) for every Laplace term are independent, these computational tasks for different Laplace terms are able to be conducted simultaneously on multiple CPUs. Alternatively, a fast solver based on hierarchical matrices [34] can be used to reduce the computational costs for solving the large-scale problems.
The 3D finite element simulation is based on the ABAQUS implicit solver, which provides the piezoelectric elements. Tie constraints were used to connect different components. The PZT actuator was discretized with 192 piezoelectric elements (C3D20RE), and this mesh was also applied to the PZT sensor. The adhesive layer was modelled using 9,904 linear elements (C3D8R). The group velocities of the S0 and A0 modes at 300 kHz are 5268.5 m s −1 and 2974.8 m s −1 , respectively. In order to capture the behaviours of these short-wavelength waves accurately, the plate should be modelled using sufficiently fine FE mesh. The intact and cracked plate were discretized using 164,738 and 179,800 quadratic hexahedral elements of type C3D20R, respectively, and such mesh density ensures that there are 37 nodes per wavelength of the S0 mode and 21 nodes per wavelength of A0 mode at 300 kHz. The finite element meshes for the intact plate and cracked plate lead to 2,731,968 and 2,980,428 degrees of freedom respectively. By contrast, the degrees of freedom required for the boundary element discretization of the intact plate and cracked plate are 2,400 and 2,616 respectively. It is shown clearly that the proposed BEM model is able to reduce the total degrees of freedom significantly.
In SHM, sensor signals from the intact plate are often used as the baseline for damage detection. Thus, the peak amplitude of the first wave packet in these signals is adopted to normalize the sensor responses from both the intact and cracked plate. When PZT1 acts as an actuator, the normalized responses of the sensors (PZT2, PZT3 and PZT4) on the intact and cracked plate are presented in figure 7. It can be seen that the BEM results agree very well with both the experimental and FEM signals.
The waves sensed by the sensors on the intact plate are comprised of two parts, the direct-arrival waves excited by the actuator and the reflected waves from plate edges. When it comes to the cracked plate, in addition to these two components, the sensor responses include the scattered waves caused by the crack. In order to visualize the wave propagation and interaction with the crack, the displacements within the plate are calculated using equation (8). The snapshots of the displacement field (u 2 ) for the intact and cracked plate over two time frames are illustrated in figure 8. It is shown that the excited S0 mode reaches the plate edges which are close to the actuator first, followed by its interaction with the crack. The scattered waves by the crack can be observed comparing figures 8(c) and (d). Figure 9 compares the normalized scattered signals computed based on BEM, FEM and experimental data. When it comes to PZT3 and PZT4 signals, it can be seen that the BEM results are in good agreement with both FEM and experimental results, except for some noise in the experimental data (see figure 9(b)). This noise (the first wave packet) is mainly due to the reflected S0 mode wave from the right plate edge. When the measurements were taken on the intact plate and the cracked plate, it was difficult to ensure that the boundary conditions in these to cases were the same, which could cause such signal difference. Apart from the boundary-condition difference, the noise could also originate from the environmental change and electronic devices. When the crack-scattered signal is weaker than the noise level, it is prone to be contaminated by the noise and thus cannot be identified accurately. This happens to the sensing path 12 (see figure 9(c)), in which the amplitude of the scattered wave is only around half of the noise amplitude. Such poor-quality signal could make some crack-detection algorithms, such as a method of triangulation [35], ineffective, in which the first wave packet in the scattered signal is considered as the scattered wave directly from the crack. In order to improve the detection accuracy of these algorithms, higher sensing amplitude of the scattered wave is required and BEM simulations can be conducted to provide guidance for the layout of PZT sensor network.
The snapshots of the crack-scattered wave field (total strain u α,α on the plate surface x 3 = h/2) at t = 45 s and 60 s are illustrated in figure 10. It can be seen that the scattered wave does not reach the sensors at t = 45 s, and thus the first wave packets from the experimental signals in figures 9(b) and (c) do not represent the crack-scattered wave but are caused by the noise. This scattered wave is being sensed by the sensors at t = 60 s, which corresponds to the first wave packets in the BEM-simulated results shown in figure 9.

Crack detection
The proposed BEM model can be used for many SHM applications, such as the evaluation of damage detection algorithms. In this section, a reconstruction algorithm for probabilistic inspection of defects (RAPID), which was proposed to detect small defects by Zhao et al [36], is utilized alongside the proposed BEM model for crack localization. The aim of this section merely is to assess the application of the simulated signals (in particular the damage scattered signals) as input to damage detection algorithms, which can then be used for development purposes.
In the RAPID algorithm, the difference between the baseline sensor signal X and any monitored signal Y is described using a correlation coefficient ρ cor = CXY σXσY , in which C XY is the covariance of the baseline data set X and the monitored data set Y, and σ X and σ Y denote the standard deviations of X and Y. Such correlation coefficient is linked to the crack probability at any position x = (x 1 , x 2 ) in the plate domain by the following expression: in which the PZT network consists of N PZT transducers and an adjustable scaling parameter β scale > 1 is used to control the size of effective crack-probability area around every sensing path (β scale = 1.02 in this study); R nm is given by: representing the ratio of the sum of the distance between the point x and an actuator center x n and the distance between x and a sensor center x m , to the distance between the actuator and the sensor.
The crack distribution probability image is plotted using the probability function P (x 1 , x 2 ) and the pixel with maximum probability indicates the estimated crack location. Figure 11 shows the constructed images on the basis of simulated and experimental data. It is shown that the RAPID algorithm succeeds in estimating the crack location using the data from both BEM simulation and the experiment, and these two images are nearly the same. This indicates that this BEM model can be used to evaluate a damage detection algorithm and check its applicability before conducting experiments, which helps to improve the algorithm and optimize sensor network layouts.

Parametric study
The scattered field signal due to the crack can be fully described using the proposed BEM model. In order to quantify   the strength of the scattered signal ( ⌢ R scatt ) received by the sensors, the following definition is used: (25) in which ⌢ V peak scatt represents the peak value of the envelop of the scattered signal and ⌢ V peak direct denotes that peak value of the signal directly from the actuator. The envelope of the sensor signal can be calculated as follows: where H denotes the Hilbert transform and i is the imaginary unit.
First, one actuator is supposed to be positioned along a quarter circle (125 mm radius) which is centred at the crack center and it accommodates different incident angles (θ inc ) ranging from 0 • to 90 • with 15 • angular spacing. Seventytwo artificial sensors are evenly distributed along the the circle centred at the crack center with a radius of 100 mm. The angular distribution diagrams of the scattered-signal strength under different incident angles are illustrated in figure 12. Each diagram focuses on one incident angle and the angular distributions under other incident angles are shown with a transparent effect to facilitate comparison. Generally speaking, it can be seen that the scattered signals become stronger as the incident angle increases. In addition, the strength of the transmitted signals is higher than that of the reflected signals, which confirms the basic assumption in the RAPID algorithm [36].
Next, the distance between the sensor and the crack center decreases from 100 mm to 75 mm, and then to 50 mm.
Such decrease causes the variation of the angular distribution of the scattered-signal intensity, which is shown in figure 13. The scattered signals become stronger when the sensors are placed closer to the crack. This indicates that in order to minimize noise interference in experiments and real SHM applications, the covered area by the sensor network can be shrunk to increase the intensity of the scattered signals. When the noise strength is much higher than the scattering intensity, the 'scattered' signals calculated from the measurement data include more noise than the scattered-wave information and as a result the RAPID algorithm could fail to localize the defect. This indication is especially important for the damage detection algorithms which depend highly on the time-of-flight of the scattered wave, such as the method of triangulation, because the time-of-flight obtained from the noise-dominated 'scattered' signal could deviate quite a lot from the actual value.
Some comments regarding the limitations of the proposed BEM model should be made. First, it is difficult to consider the non-linearity due to the contact between crack surfaces. Time-domain BEMs based on time-stepping technique are more suitable to solve this non-linear problem. Secondly, since the PZT transducers are not incorporated into the boundary element formulation in a fully coupled manner, the electromechanical impedance responses of the system cannot be computed based on the proposed BEM model. This coupling can be achieved based on the PZT-transducer modelling technique used in reference [6]. Finally, this BEM model is only suitable to isotropic plates. In order to extend it to anisotropic plates, new fundamental solutions should be derived.

Conclusions
In this paper, the first boundary element formulation based on plate theories has been proposed for the modelling of Lambwave structural health monitoring of thin plates. In order to incorporate the piezoelectric actuators into the boundary element formulation efficiently and simply, an equivalent pinforce model, a new actuator model, has been proposed. Existing dynamic fundamental solutions used for such formulation have been rewritten in a more compact form, including some modifications to represent the dispersion behaviours of Lamb waves more accurately. Some parametric studies were conducted based on the newly-developed BEM model, which showed that in order to increase the intensity of the sensor signals due to the scattered waves by a crack, sensors should be placed closer to the crack. The BEM-simulated results were shown to be in excellent agreement with the 3D FEM solutions and experimental results. The proposed BEM formulation is very promising in modelling Lamb waves and wave-crack interaction because of the advantages of boundary-only discretization and low labor-cost requirement in mesh generation. This formulation is especially suitable for providing a large amount of data for a plate with many different crack parameters, such as length, orientation and location, because only the crack need to be remeshed and only coefficients associated with the crack in the global system coefficient matrix need to be updated.

Acknowledgment
This research was financially supported by China Scholarship Council.

Appendix A. Fundamental solutions in the Laplace domain
The dynamic fundamental solutions for the symmetric case can be written as follows: [(δ βγ r ,α + δ αγ r ,β ) r ,n + (n α r ,β + n β r ,α ) r ,γ ] + 4 where the normal derivative r ,n = r ,α n α ; the three Laplacedomain wave numbers can be expressed as: ] and the cut-off frequency ω cS = 2 √ 3κ S c 1 /h; and The fundamental solutions for calculating total in-plane strains are given by: The dynamic fundamental solutions for the antisymmetric case are given by: [(δ βγ r ,α + δ αγ r ,β ) r ,n + (n α r ,β + n β r ,α ) r ,γ ] (r ,α n β + r ,β n α ) } (A34) (δ βγ r ,n + r ,γ n β ) } (A35) where the three Laplace-domain wave numbers can be expressed as: in which Q A = √ (S 1 − S 2 ) 2 − 4S 1 S 2 (ω cA /p) 2 , S 1 = κ 2 A µh/D, S 2 = 12/h 2 and the cut-off frequency ω cA = 2 √ 3κ A c 2 /h; and The fundamental solutions for calculating total strains can be written as: In the above fundamental solutions, the functions ψ , χ , η and ξ, and their derivatives are given by: where z Y i = ζ Y i r; K 0 and K 1 are the modified zeroth-and firstorder Bessel functions of the second kind, respectively; the superscript Y denotes S and A for the symmetric and antisymmetric cases, respectively.