Transducer Placement Option of Lamb Wave SHM System for Hotspot Damage Monitoring

: In this paper, we investigated transducer placement strategies for detecting cracks in primary aircraft structures using ultrasonic Structural Health Monitoring (SHM). The approach developed is for an expected damage location based on fracture mechanics, for example fatigue crack growth in a high stress location. To assess the performance of the developed approach, ﬁnite-element (FE) modelling of a damage-tolerant aluminum fuselage has been performed by introducing an artiﬁcial crack at a rivet hole into the structural FE model and assessing its inﬂuence on the Lamb wave propagation, compared to a baseline measurement simulation. The efﬁcient practical sensor position was determined from the largest change in area that is covered by reﬂected and missing wave scatter using an additive color model. Blob detection algorithms were employed to determine the boundaries of this area and to calculate the blob centroid. To demonstrate that the technique can be generalized, the results from different crack lengths and from tilted crack are also presented.


Introduction
Non-Destructive Testing (NDT) [1] has been implemented in many industries to ensure structural safety and reliability. Structural Health Monitoring (SHM) as a complementary solution with respect to the already existing NDT techniques has been a subject of interest in the past decade due to its potential economic benefit, particularly in structural maintenance [2,3].
Along with fiber-optic techniques such as Fiber-Bragg-Grating sensors [4,5], an ultrasonic guided waves-based solution such as those making use of Lamb waves is one of the promising techniques for SHM due to its long-range inspection capability up to several meters [6,7], which is suitable for monitoring large structures such as an aircraft fuselage. Also, in comparison to the FBG technique, Lamb wave SHM sensors typically require less complicated installations. This is particularly useful for monitoring aircraft which are already in-service, since aircraft operators generally tend to be reluctant to make larger modifications.
Lamb wave damage monitoring relies on the detection of interactions between the Lamb wave and the damage [8]. Various properties of propagating wave are modified by the crack and various properties of the signal such as amplitude, frequency, phase, etc. are assessed to calculate the damage index (DI). The most commonly used sensors for generating and capturing Lamb waves are piezoelectric transducers (PZT) [9]. Since PZTs are permanently attached to the surface of the structure, sensor placement becomes tremendously important because the quantification of DI relies heavily on the sensor positions. A poor transducer placement will result in weak or undetected capture of wave scatter which in turn decreases the damage detectability of the SHM system.
The literature describes different approaches for sensor placement options, such as prioritizing the sensor location based on a detectability limit [10], using a modal analysis parameter determination for damage localization assessment on a truss structure [11], and by using global search and greedy algorithms [12]. Li et al. [13] proposed a sensor optimization algorithm based on maximum energy consumption on sensor candidate location of a civil engineering structure.
Fendzi et al. [14] proposed a novel approach for sensor placement by using geometric dilution of precision (GDOP), which is based on Lamb wave ray tracing method for known damage locations. Haynes [15] proposed sensor placement by minimizing the Bayesian cost to select locally optimal sensor locations. However, if the damage occurs outside of the designated area, the system might fail to detect it. A similar approach by using Bayesian experimental design has been conducted by Flynn and Todd [16].
Furthermore, Sun et al. [17] performed discrete optimization by using the artificial bee colony algorithm to optimize their objective function which is based on modal assurance criterion (MAC). Similar approaches by using search metaheuristics were also delivered by Yi et al. [18], Zhao et al. [19], and Shan et al. [20]. In more recent study, Capellari et al. [21] proposed an optimal sensor placement by employing Polynomial Chaos Expansion and stochastic optimization to maximize the gain in Shannon information.
Thiene et al. [22] introduced the DI-free sensor placement optimization based on the fitness function which maximizes the coverage area of the sensor network. They calculated the coverage of each pixel in the geometry based on pitch-catch technique, so that every pixel that contributes to the probability that a damage in random location is being detected is counted. Venkat et al. [23] used a Finite Element (FE) simulation platform to build differential images between the undamaged and damaged structure. They plotted the summed-up energy captured by all sensors and determined the most optimal sensor location by the highest captured energy.
From all these studies, we concluded that most of the sensor placement strategies are for civil engineering structures and that one of the most common features used for the objective function is the MAC, which is typically applicable for low-frequency vibrations. However, we have not seen any application using MAC in Lamb wave SHM frequency, which is normally above 100 kHz. Only [14][15][16]22,23] are related to sensor placement for guided wave SHM.
Summarizing references [14][15][16]22,23], there are two streams of sensor placement research that are specialized for guided wave-based SHM: (1) The approach exhibiting known area of damage location; and (2) The random damage location approach. If a damage location can already be approximated by analytical methods from fracture mechanics, FE simulation and fatigue testing, then the optimum sensor location can be determined by maximizing the DI of the sensor network around the predicted damage location. However, for random location damage such as hail impact, a different objective is needed.
From this current state of the art, our work is performed to develop the work presented in two most recent articles [22,23]. For this paper, the objective is to propose an image processing for known damage location approach, known as hotspot SHM [22].
Due to the expected number of different damage locations occurring in an aircraft fuselage multiplied by the numerous possibilities for sensor placement, conducting an experiment or even a numerical simulation to track the wave scatter is not practical. To address this issue, we propose a design for a sensor placement in a hotspot SHM system based on an additive color model and blob detection algorithm that can detect residual wave scatter from simulation data.
After the introduction in Section 1, we have structured the paper as follows: In Section 2, we briefly review the Lamb wave propagation theory whose effect is employed in Section 3 for determining the SHM sensor positions by using image processing for predictable critical crack size, location, and orientation. The results and concluding remark are presented in Sections 4 and 5, respectively.

Lamb Wave Propagation
First described by English mathematician Horace Lamb in 1917 [24], the description of acoustic wave propagation in solid plates has been evolving ever since. A theoretical analysis of Lamb wave propagation in metallic materials, composite, and hybrid materials is described in [25][26][27]. The elastodynamic wave equation for an anisotropic inhomogeneous medium in a d-dimensional bounded domain Ω ⊂ R d (d = 2, 3) is given in by: where u(x,t) is the time and space dependent displacement, ρ the material density, C ijkl the material stiffness tensor, ε kl the strain tensor, and f (x,t) the source function, respectively. After applying boundary conditions of two parallel surfaces, there are two solutions to Equation (1) for wave propagation in homogeneous material of given density, described by: which are well-known as symmetrical (S-Mode) and anti-symmetrical (A-Mode) Lamb wave propagation modes, respectively. The definitions of p and q are given by: where ω is the frequency, c L is the longitudinal bulk wave velocity, c T is the transversal bulk wave velocity, and k is the wavenumber. The numerical solution of Equations (2) and (3) can be drawn as a dispersion curve, which describes the wave velocity as a function of frequency-thickness product [28]. As they introduce potentially overlapping signals, higher-order Lamb modes are generally undesirable, and the maximum cutoff frequency is normally determined as before the A1 mode appears.

Sensor Positioning Approach for Hotspot SHM
This section describes the developed methodology for sensor placement in hotspot SHM. It considers: (1) Crack growth and critical crack size in a damage-tolerant aircraft substructure; (2) Numerical simulation of Lamb wave propagation in a plate-like structure; and (3) Image processing that involves displacement subtraction of damaged from undamaged structures.

Crack Growth in Damage Tolerance Structure
To understand deterministic sensor positioning for a predictable crack location, firstly, it is important to understand the concept of damage-tolerant design. The definition of damage tolerance is "the ability of the structure to sustain design limit loads in the presence of damage caused by fatigue, corrosion and other sources until such damage is detected and repaired" [29]. The important key elements in damage tolerance design are: (1) the assumption of initial damage existence; (2) damage growth in the material due to structural loading; and (3) The critical damage size up to which the structure endures the loading before catastrophic failure. These three key elements are synchronous to the regions I, II, and III in a typical da/dN curve [30], which describes crack propagation rate as a function of stress intensity factor (SIF) range during fatigue cycling (∆K).
The assumption of initial crack existence falls in region I, where the crack growth is typically slow. This region can be covered by some advanced (but also large and expensive) NDT and material characterization techniques such as X-ray computer tomography [31] or scanning electron microscopy (SEM) [32]. After passing the threshold stress intensity factor range ∆K th , region II begins where the crack propagation rate is stable, and crack growth normally follows the Paris-Erdogan law [30,33], which has Paris-Erdogan constants C and m. In this region, the crack becomes much larger and sometimes can be seen with naked eye [34]. After transition between region II and III, the crack propagation rate becomes rapid and unstable until final failure.
Lamb waves can interact with a crack which has a size at least the half of its wavelength [35,36], and since the critical crack tolerant size (which can be up to several hundred millimeters [22]) is generally larger than the wavelength (which is typically less than 100 mm), it is safe to assume that Lamb waves can also interact with the critical crack. Fracture mechanics and fatigue analysis can predict the most probable crack location for certain loading conditions [37].
The SIF is generally larger in the area around the notch, thus one can expect that a crack will occur around a notch. For example, the changes in cabin pressure in an aircraft fuselage can be approximated by the internal pressure in a thin-walled cylinder where the hoop stress is two times larger than the axial stress [38]. Therefore, the crack orientation is expected to be orthogonal to the direction of the hoop stress. By knowing the most probable crack orientation and location and the critical crack size for a certain geometry, two FE simulations scenarios of wave propagation can be performed [39]: (1) Lamb wave propagation in an undamaged structure as a baseline and (2) Lamb wave propagation in a critically damaged structure.

Simulation of Lamb Wave Propagation
The FE formulation for Lamb wave propagation is based on Hamilton's principle [40] as described in: where Γ and Φ are the surficial and volumetric integral areas, ρ is the material density, u and ü are the particle displacement vectors in the material and their corresponding accelerations, respectively, ε is the strain tensor and C ijkl is the stiffness matrix. The external forces can be classified as surface load F S and volume load F V . To numerically solve Equation (6), the geometry involved is divided into discrete mesh elements over which the equation can be approximated [41]. For a 3D-problem, there are four element types: prism, pyramid, brick and tetrahedron [42]. To reliably model an ultrasonic signal, the recommended number of mesh elements is 4 elements per A0-or 8 per S0-wavelength [43]. Beside the spatial discretization, time incrementation of Equation (6) is needed as well. The minimum requirement to ensure numerical stability of time integration is the Courant-Friedrich-Lewy (CFL) condition [40,[43][44][45].

Additive Color Model
Mathematical representations of color can be summarized as additive and subtractive color models [46]. The additive color model is used for combining light of different color and the subtractive color model is used for combining pigments and dyes, which are typically represented in 4 basis colors: cyan, magenta, yellow, and black (CMYK).
Both additive and subtractive color models are typically represented in vector form. Both color models relate to the perception of the human visual system. In the additive color model, the background is black, and color is composed of a red, green, and blue (RGB) tuple as depicted in Figure 1. This relates to the three cone cells of the human retina (long, medium, and short) which have peak sensitivities to light of wavelengths 559 nm, 531 nm, and 419 nm. These correspond to red, green, and blue light, respectively [46]. Normally, red, green, and blue color may be written in 8  When working with digital imaging, the additive color model (i.e., RGB) is a natural choice because what humans see in a measurement device/computer screen is normally encoded in RGB values. Computer scientist and electrical engineers have selected to represent the digital vision (e.g., in computer monitor) as an RGB array because of the trichromatic cone cells in human eye which are sensitive to the frequency spectrum of red, green, and blue photons that are falling into human retina. Therefore, the origins of the additive color model are physical and biological rather than purely technical.
In the 3-bit model, shown in Figure 1, other colors such as cyan, magenta, yellow, and white are obtained by adding basic colors. For instance, cyan is the addition of blue and green, while white is the addition of red, green, and blue. Colors may also be subtracted from one another, as given in Table 1. Note that in the convention used negative values are set at zero, as negative light amplitude has no physical meaning. In this paper, material displacements are represented as colors in the additive color model and mathematical operations are performed to identify and visualize the optimum sensor positions.

Simulation
The method for choosing the simulation parameters has been described in [23,43]. The plate dimensions are 600 mm × 400 mm × 2 mm. For this work, the following parameters are used: ABAQUS explicit, aluminum properties (Young's modulus of 70 GPa, Poisson ratio of 0.33, density of 2700 kg/m 3 ), quadratic brick mesh (C3D20) with a global mesh size of 1 mm, single node out-of-plane excitation with windowed 5 sine-cycle of central frequency of 250 kHz with 1 N concentrated force, dynamic implicit step, no boundary conditions imposed, time increment of 0.1 µs (which means a sampling frequency of 10 MHz), total time period of 500 µs and a single nodal output precision. The specification of the computer where the simulation was run was: Intel Xeon E5-1620 3.5 GHz (Quad-core 8-Threads), 32 GB DDR3-RAM, and NVidia NVS310M Graphic card (GPU acceleration was not activated).

Data Extraction
To capture the wave propagation image, the screenshot program which is called Lightshot can be downloaded from prntscr.com [47]. In the program settings, the option "Keep the selected area position" must be selected so that Lightshot remembers the X-Y position of the captured screenshot for every time increment. As mentioned previously, this technique is less exhaustive and more memory efficient as every saved image has a size of only around 500 KB. For 20-time increments and 2 cases (uncracked and cracked plate), only around 20 MB of space is needed in comparison to ODB extraction which takes around 3.2 GB of space.
In a computer, the color in a single pixel is represented as an RGB array and this can be used to represent different Lamb wave displacement amplitudes. Figure 2 shows an example simulation of Lamb wave propagation 30 µs after excitation. The displacement magnitude U is defined as: where U x , U y , and U z are the displacements in the x, y, and z-directions, respectively. It is obvious that no displacement (U = 0 nm) is shown as a blue pixel, while a displacement of 2.5 nm is shown in green, and a displacement of 5 nm is shown in red. The values in-between such as 1.25 nm and 3.75 nm are shown in cyan and yellow, respectively. This colormap 'rainbow' is the default colormap in ABAQUS. Note that this colormap is slightly different from the basic 3-bit RGB colormap depicted in Figure 2 as it has more color transitions, i.e., there are smoother transitions between blue and cyan, cyan, and green, and so on. For ease, we neglected displacements larger than 5 nm (shown in grey color), which are mostly due to the excitation, because in this simulation this only appears in few locations. While this image processing procedure offers less displacement information as all displacement values are translated into an RGB array, we find this alternative procedure much faster, and more memory efficient rather than extracting data directly from the ABAQUS ODB binary file.

Differential Images
The approach is demonstrated by taking screenshots from the ABAQUS Viewer User Interface. Figure 3a shows Lamb wave propagation at t = 100 µs in an uncracked Al-7075 plate. The simulation images are actually similar to those which are depicted in Figure 4 of [16], except that the simulation platform, color map, excitation signal, and the geometries are slightly different.
The plate contains 3 rivet holes as depicted in Figure 3a,b shows Lamb wave propagation in the same plate but with a symmetric crack (from tip-to-tip, including the hole diameter of 10 mm) of 28 mm length in the middle of the plate (marked by a yellow rectangle). The images captured have size of 1210 × 807 pixels, so the resolution is 2 pixel/mm. Images are stored as matrices with a size of 1210 × 807 × 3, where each pixel has 3 arrays, each containing a normalized floating value between 0 and 1 for each of the RGB colors.
A similar pattern of Lamb wave propagation can be observed if the crack length differs by +/− 10%, as shown in Figure 3c. In this case, the crack length is 30 mm instead of 28 mm. However, if the crack is much larger, a notable change in the wave propagation pattern can be observed, as depicted in Figure 4d. In this case, the crack length is 60 mm. Similarly, the wave propagation at t = 125 µs for an uncracked plate, and for plates with crack lengths of 28 mm, 30 mm, and 60 mm can be observed in Figure 3e-h, respectively. By subtracting the image of the cracked plate (Figure 3h) from that of uncracked plate (Figure 3e), the reflected wave scatter image can be obtained, as shown in Figure 4a. Note the RGB values are subtracted rather than the displacement values. Pixels, for which there is no change in wave scatter are shown as black. In this figure, the reflected wave scatter is highlighted by the yellow rectangle, while the corruptly transmitted wave is obtained as well (in the red rectangle) but is not clearly visible. In an analogous way, we can subtract Figure 3e from Figure 3h, and in this case, the corruptly transmitted wave scatter image is highlighted more (red rectangle) as depicted in Figure 4b, while the reflected wave scatter can still be seen (yellow rectangle) but is less visible. For the further sections in this paper, only the results from the uncracked plate and the cracked plate of 60 mm crack will be shown for conciseness. To highlight both the reflected and corrupted wave scatter, the two images shown in Figure 4a,b can be joined to form a composite image.
Two common ways combining images are described here, the first one is called image addition. In this case, both image matrices are just mathematically added. The second one is called image fusion and uses the 'imfusion' function in MATLAB, where the two images are firstly converted into greyscale mode, given a false color, and then added mathematically. Furthermore, it is possible to invert the composite images to obtain the inverted images. This step is not necessary, but some people subjectively may find the colored wave scatters easier to track if the background is white. The inverted images from image addition and image fusion are depicted in Figure 4c,d, respectively. For the sensor placement procedure described in Section 4.4, the inverse fused image (Figure 4d) is used because we can still see the representation of reflected and missing wave scatter and that enables us to cross check the results with the simulation data.

Sensor Placement
The wave scatter which was caused by both reflections and corrupted transmissions due to the crack front are represented by false color (i.e., green and magenta pixels in Figure 4d). These pixels have a different color from the background (white). Regions of interest are determined by using the Matlab blob detection function to locate areas of adjacent green and magenta pixels. The larger the blob is, the larger the area of wave scatter. The sensor should be placed in the center of the largest blob, so that it will have a high probability of capturing a portion of wave scatter from cracks.
The blob detection algorithm is based on the Laplacian of Gaussian [48][49][50] with a kernel of 8-pixel connectivity. Given the Gaussian function G of an input image f (x,y) and feature scaling σ as given in Equation (8). The Laplacian operator ∇ 2 is given in Equation (9).
By applying the Laplacian operator in Equation (9) to the Gaussian function in Equation (8), one obtains the Laplacian of Gaussian, commonly known as LoG, as described in Equation (10). Hence, the blob centroidx,ŷ with the scaleσ is the simultaneously local extremum of the LoG in Equation (11).
For the blob detection, the kernel of 8-pixel connectivity is used because it is more suitable for a larger area since the diagonal neighbor is counted as well, see Figure 5a, while 4-pixel connectivity (Figure 5b) is typically used for line and corner detection. In Figure 5a,b, the meaning of −1 and 0 are pixels which are counted and are not counted as a neighbor of the center pixel, respectively. The blob detections from various time increments are depicted in Figure 6a-d. As mentioned before, we would like to place the sensor in the location with the largest change in signal over time, i.e., the largest blob. Since it generally contains more than a single pixel, the centroid can be calculated to determine the average pixel location that would receive wave scatter. In Figure 6a-d, the largest and the second-largest blob centroids are marked by red and green dots, respectively. Blob boundaries are marked by the yellow polylines.  Table 2. The detailed procedure in MATLAB to find to trace the blob is described in Algorithm 1. Loop over all traced region from 1 to n: Store the area information in matrix 'S' Sort from the largest to smallest blob Store the blob boundaries Calculate blob centroids Assign X-coordinate from the centroid Assign Y-coordinate from the centroid End the loop After a certain time, the wave pattern becomes more chaotic due to multiple reflections from the crack front, rivet holes, and plate edges so that more smaller centroids will be born that are not exactly aligned with the mid Y-axis anymore. The smaller centroids imply that the potential energy capture by PZT is getting smaller and this will be aggravated by Lamb wave attenuation. This is the reason we recommend 'early wave scatter capture' for hotspot SHM design. Typically, one can decide the best sensor position by considering the movement of the centroid per time increment, also known as ray tracing [51].
Furthermore, it is possible to fuse Figure 6a-d into a single image. We performed this operation and the result is depicted in Figure 7, where each pixel contains information about the normalized intensity between 0 and 1 which is then mapped into the rainbow color scale. From Figure 7, it can be subjectively judged that the best sensor position is between X = 44 cm and 57 cm and the second-best sensor position is between and X = 34 and 38 cm, while the vertical coordinate for both positions remains at Y = 20 cm. In order to demonstrate that the image processing algorithm also works for a different case of simulations, we slightly modified the case for the (a) the critical crack length to 30 mm, and (b) the orientation of the crack by 8 • . The whole procedure was repeated, and the results are depicted in Figures 8 and 9, respectively.  From Figure 8, it can be seen that the areas with higher pixel intensity (colored in red with value between 0.7 and 1) is smaller than those of Figure 7. This is absolutely normal, since the wave perturbation at the crack front due to a crack length of 30 mm is smaller than those of 60 mm.
Meanwhile from Figure 9, it can be seen that the crack orientation changes the direction of the reflected wave scatter and it is still conform with Snell's law. However, it does not change the orientation of areas where the wave scatter is not present. After cross validating with the original simulation data, it can be confirmed that angled crack only heavily influences the reflected wave scatter portion, but not the missing scatter portion (Figure 10a,b, marked in yellow). For current work, we have demonstrated a new efficient method based on an additive color model to find the best two locations for placing the sensor. However, the number of sensors that are allocated for every case can be adjusted according to the manufacturer and/or operator requirement to achieve required crack detectability. Furthermore, not only the number of sensors, but also the excitation frequency can be changed depending on the size of the critical crack that must be detected. A higher frequency means a shorter wavelength, and this would mean the wave will be able to interact with smaller critical crack. However, such a higher frequency Lamb wave will also be more quickly attenuated than the lower frequency Lamb wave. Therefore, in order to stabilize the SHM network performance, more sensors will be required. Nevertheless, when more sensors are employed, higher procurement costs are also expected due to more weight, more data processing capability, etc. We regard this as the classical trade-off between SHM investment cost and SHM network reliability and would like to pass this decision back to the aircraft manufacturer or operator according to their needs.

Conclusions
In this work, we have demonstrated a novel technique to design the sensor network topology for hotspot SHM by using differential images and blob detection algorithm. While our image processing technique does not allow a quantitative approach to observe the nodal displacement, i.e., displacement from every single FE node, we believe this technique offers a more holistic view (Figure 7) of where to place the PZT sensors on the structure to be monitored.
Also, with this technique we think that the sensor placement can be done more quickly without exhaustive data processing from simulation file for each surface while not sacrificing too much spatial resolution. In practice, even the extracted nodal data from the simulation must be interpolated, since in reality a PZT sensor would always occupy more than a single node (e.g., a typical PZT in our lab has a diameter of 1 cm, so would theoretically occupy about 78 FE nodes). Therefore, as a concluding remark, we hope this technique will help further research in sensor placement.