A Real-Time and Wireless Structural Health Monitoring Scheme for Aerospace Structures Using Fibre Bragg Grating Principle
INTRODUCTION
The population of artificial satellites and near-Earth objects in space has massively increased since the last decade. With an advancement of space industries all over the world, accumulations of expired and malfunctioning satellites along with natural space debris are presenting a rising concern for current and prospective aerospace structures. In the past, critical damages were reported multiple times due to impacts between rogue and victim aerospace structures leading to severe threats of catastrophic accidents arising from structural failures [1.2]. In this book chapter, a wireless remote monitoring scheme was presented to detect, identify, and locate such unwanted and sudden impacts in real-time using fibre Bragg grating sensors. The developed scheme was first analytically modelled, then implemented, and finally validated through in-lab measurements. For the sake of the discussion, the overall work was divided into two major sections. The first section presented identifying the location and estimation of any impact on a two-dimensional surface using fibre Bragg grating sensors on-board, while the second section presented a wireless monitoring scheme for transmitting the sensor data in realtime to a remote location. It was envisioned that such critical mission information will be extremely useful in (1) making real-time decisions on flight management information system (for example, continuing flight or returning to the base) and (2) understanding the effects of impacts on aerospace structures. Now, location estimations and characterizations of impact loads (IL) for various types of structures are already well established in current literatures [3]. On the contrary, simultaneous identifications of not only the IL locations but also analysing the impact load time history (ILTH) in realtime for space applications have not been properly discussed yet. In this research, a non-stochastic inverse analysis was used to address the uncertainty in detecting, identifying and locating ILs on a simply supported two-dimensional plate to obtain and simplify ILTH using a sinusoidal function and Whitney and Pagano's solution [3]. To do so, an effective search method was utilized to reduce the objective function along with a novel and simplified least square method in detecting the maximum limits of the ILTH and IL. Finally, a real-time monitoring system was proposed through optical sensor interrogators and wireless transmitters.
To further investigate the efficiency of the proposed algorithm, variable distances between sensors could be tested. Also, a study could be done on the minimum distance between sensors required to ensure the algorithm would work properly. Simultaneous multi-point and distributed impacts could also be taken into consideration; however, such studies were beyond the scope of the current discussion.
HAND ILTH ANALYSIS
Existing Techniques
Passive sensing is a well-established method for health monitoring of structural bodies experiencing mechanical stress over the time in civil and aerospace engineering. Analysis of IL and its identification constitute the primary step to correctly attribute the health of a structure in operation by locating potential damage. On the other hand, ILTH provides the information regarding the sequence of values as a time-varying quantity which is measured at a set of fixed intervals. Now, different techniques of inverse analysis exist for IL and ILTH assessment [4,5]. Inverse analysis can be executed using two methods; either by reducing the error between predicted data and measured data acquired using model-based forward solving of a mathematical model or by using a neural network that requires extensive training of any given dataset. The model-based approach is well accepted among researchers due to its convenience of implementation and simplicity of calculation compared to a neural network-based approach. Another advantage of the model-based method is that ILTH and IL can be obtained either separately by transfer matrix or simultaneously by back analysis [4]. Now, IL and ILTH can be determined from local sensor measurements by using different methods of inverse analysis. However, inverse (deconvolution) analysis is challenging to implement as the objective function at the minimum point does not always converge to zero as the process is noisy and measured responses result in instability. This challenge can be eradicated by the regularization of the measured responses. To do this, additional conditions are typically imposed to optimize the objective function by implicitly filtering the noisy measured data. For non-linear structure, Kalman filter and recursive least square estimator are used. In the case of linear structural systems Tikhonov regularization, truncation, and singular value decomposition can be used. In general, one common drawback of such classical regularization methods is the requirement of determining appropriate regularization parameters. Additionally, the presence of noise further degrades the overall performance of these algorithms. Unlike the deterministic methods, the stochastic analysis considers the uncertainties associated with structural properties, like measurement error and elastic modulus, this method also accounts for the noise in the measured data. Inverse analysis and stochastic analysis can work together, where the stochastic analysis will discover an affected area and inverse analysis will then work to detect the area of the impact location and the IL magnitude. The drawback for such statistical distribution as prior information on the data needs to be known. An alternative to stochastic analysis method is interval analysis where to rectify this disadvantage, interval analysis can be introduced instead of stochastic analysis. Interval analysis takes into consideration the minimum and maximum boundaries of the system response. However, it does not provide the insight into the distribution of ILTH and IL. Also, interval analysis alone or combined with the regularization method it is a costly process for real-time detection.
In this research, the IL and ILTH were concurrently obtained, and the uncertainty was addressed by the inverse analysis which is non-stochastic. To minimize the objective function, particle swarm optimization (PSO) technique was adopted as an efficient search method. A novel method was created to detect the maximum limits of the IL and ILTH which is based on the least square method. The ILTH was simplified by a sinusoidal function.
Three-Layer Identification Process
A three-layer identification process for IL and ILTH analysis is introduced in this section. It is an effective and simple scheme. These layers can be integrated to effectively detect IL and ILTH in an intelligent structure, which contains a network of at least three sensors. These layers are designed in such a way that, it will be able to differentiate IL and ILTH of a structure. Now' for this design to work, the network must have at least three sensors. At first, strain time histories are collected from all sensors and only three least times of arrivals related to the nearest sensors to the IL are selected for further analysis. Time of arrival (TOA or ToA) is the travel time of a (surface) wave from the point of impact to the sensor. IL can then be evaluated by using the triangulation method [4]. Here, IL is decided by minimizing the error between measured and the predicted TOAs. The output from the first layer helps in defining the sampling space. One significant challenge in this approach is finding TOAs under a noisy response. In addition, measurement of the speed of the surface wave is a complex process. For an anisotropic plate, that also depends on the wave frequency and propagation direction.
The next layer of the identification process contains detecting IL and reconstructing ILTH by reducing the predicted and measured strain time histories. Now the noisy measured data from the sensor locations were filtered using the moving average method. Those values are later applied to determine IL and ILTH. By subtracting the filtered strain time history from the noisy one, the error is measured. Finally, the identification process approaches the third layer based on the maximum observed error in the measured data and establishes the extreme ILs and ILTHs.
Theories
Estimation of Location
At first, TOAs were analysed by differential TOA (Figure 9.1) as described in [6] because the effect of noise must be reduced:
where ti is the TOA at the /th sensor and /, is the TOA at the firstsensor location. As three sensors are employed, i can be of values of two and three (for second and third sensors, respectively) when the reference was the first sensor and so on. The objective function was then defined as the following:

FIGURE 9.1 Time of arrival (TOA) based impact identification., where
Atjm: measured value of differential TOA.
ДГ1С: predicted value of differential TOA,
x, : ^-coordinate of the /'th sensor,
y, : у-coordinate of the /'th sensor,
£: .v-coordinatc of the IL,
//: у-coordinate of the IL,
V: wave speed and
S|, S2 and S,: three sensors located closest to the point of impact (f, q).
In Equation (9.2), the plate was considered to be isotropic. For an anisotropic plate, a direction-dependent wave speed must be incorporated. In a plate structure, the bending wave’s speed V can be calculated as follows:

where
E = elastic modulus, h = plate’s thickness, p = unit of weight, w = wave’s frequency and v = Poisson’s ratio.
Fast Fourier transform (FFT) on the sensor measurements can help measure wave frequency [4]. After calculating TOAs and wave speed at the sensor locations, the steepest descent method was used to minimize the objective function (Equation (9.2)).
Deterministic Identification of Impact Location and Load Characteristics
This layer consists of a plate dynamics model. This model is used for forward solving method. Layer 2 tries to make a prediction of the strain time history. Now, two classical methods can be used here to understand the plate problem: Kirchhoff-Love’s solution and Mindlin's solution. Kirchhoff-Love’s formulation only considers the in-plane shear. In the case of Mindlin’s solution, it also considers the out-of-plane shear. In this research, a different method known as Whitney and Pagano's approach [4] that also considers anisotropy was utilized for a simply supported rectangular plate. Now, to derive an expression for strain along the л-direction (ex), Soares [4] as w'ell as Carvalho and Dobyns [4] methodology can be adapted. The relationship between modal amplitudes, the time-de- pendent function of bending deflection would be used for the expression[4], which is

where z is the respective coordinate to mid-plane, and *F( is the bending curvature respective to the л-axis.
Next, strain respective to the л-axis can be derived using the following equation [4]:

where
i and j = mode numbers, a = length of plate, b = width of plate, q = л-coordinate of IL,
4 = v-coordinate of IL and o)f/ = modal frequency of plate.
Kc and K:/ are next defined as follows:
where и is the length and v is the width of the location (impact area) and KA is a transverse shear stiffness bending factor. Although small values of и and v were chosen to symbolize pointwise IL, they can be assumed as unknown quantities. Now, F(t) that symbolizes the ILTH can then be reduced using a half-cycle sinusoidal function as in (9.8):
where
F(l = IL amplitude, со,, = unloading frequency, ы, = loading frequency, t0 = impact start time, r, = time of impact and t2 = impact end time.
It should be noted here that if the measurement starts while impact occurs, t0 is then considered to be zero. In practice, this is different since there is always a time lag as the strain recording is initiated ahead of the impact occurrence for measurement. Hence, whenever the first signal was received, the time was then set to zero. The time lag was equalled to TOA of the closest sensor. Then, r,, t2, and f, can be defined as

Equation (9.5) can be solved analytically by replacing Equation (9.8) as an integral part of the equation, which leads to
where

the objective function can then be expressed as follows:
where eic(t) is the predicted strain and eim(t) is the measured strains at the f'th sensor, and the number of time steps is N. The above function is expressed using the forward solving model. The variables of the objective function are a 6-dimensional space in terms of X, expressed as
The layer 1 output was used to shorten the range of t and. The objective function was computed iteratively to converge the sample vectors and in computing the destination vector. The following equation represents how the particle speed and position updated after each iteration:
where
Xk = particle position after Ath iteration,
XM = particle position after the (A+ l)lh iteration.
Vk = the particle’s speed after the Ath iteration,
Vk+I = the particle’s speed after the (A+ 1 )lh iteration,
P, = the particle’s past best position,
P2 = the finest position experienced by other particles in the herd,
#), B2 = two factors relates to the personal influence and social influence of particles,
A = momentum factor, r, and r2 = stochastic variables
C, D = the weight factors for updating the particle’s position, considered as equal to one, and
Each term of Equation (9.17) represents a different aspect of the computation like the first term points out the momentum of motion of the particles in the search space. The second term indicates the speed of the particles based on their personal experience, and the third term updates, relying on the values of A, B, and B2. The long-term trend of the particles can show convergent (to the optimal position), harmonic oscillatory, and zigzagging behaviour [4]. To make the algorithm more efficient, these parameters can be tuned.
Verification
At Layers 1 and 2, finite-element modelling (FEM) was adopted to verify the inverse model and forward solving mathematical model. Thus, a comparison was being done between the actual IL and ILTH with detected IL and reconstructed ILTH.
Verifying the Forward Solving Model
To verify the proposed mathematical model, finite-element model (FEM) was used for a virtual experiment using ABAQUS. An aluminium plate of size 88 x 88 x 0.25 cm-’ was modelled. A sinusoidal IL was applied. This aluminium plate was given an IL at the centre at time zero. The IL had an amplitude of 1000N, with an unloading frequency of 300 Hz and a loading frequency of 100 Hz. The simulation was done using a quadrilateral shell element with reduced integration for small strains (S4RS). The FEM consisted of 17,161 nodes and 16900 shell elements. Only lOOOpS vibration time was used, and vibrating duration was divided into 100-time intervals. It is worth to mention that due to 1000pS vibration time, a significant damping effect may not be witnessed.
Figure 9.2 shows a comparison between the suggested mathematical model and FEM result at the impact location. The results obtained from both methods showcased a suitable match.
![Comparison of the mathematical model and FEM’s strain time histories is shown at impact location at the center of the plate [3]](/htm/img/39/2134/466.png)
FIGURE 9.2 Comparison of the mathematical model and FEM’s strain time histories is shown at impact location at the center of the plate [3].
TABLE 9.1
Plate Structure's Mechanical and Geometrical Properties [3]
Length of plate |
88 cm |
Width of plate |
88 cm |
Thickness of plate |
0.25 cm |
Elastic modulus |
69GPa |
Poisson’s ratio |
0.34 |
Unit weight |
2700kg/nv! |
Stiffness properties |
D„ = D,: = 101.59Pa: D12 = 34.54Pa; D66 = 33.52Pa; A,, = 195 MPa; A+, = A55 = 64.4MPa |
Verifying Process
Table 9.1 shows the exact plate with the properties. The exact IL that was used for IL and ILTH identification through Layers 1 and 2 was then used upon the plate. For the strain time histories and impact location, three different points were selected at nine different areas that were collected for IL and ILTH identification (Figure 9.3), and in Table 9.2, coordinates of impact and sensor locations were reported.
At the first layer, TOA needs to be determined at different sensor locations. The threshold method was used at each sensor location. It was expected that the level of threshold strain would be two times the maximum amount of noise. To find the dominant frequency content. FFT analysis was applied to the time histories from the three closest sensor locations. To calculate the wave speed, corresponding parameters such as the frequency of plate properties were provided to Equation (9.3). After
![Impact locations and Sensor configuration set for identification (represents sensor location and impact location. Please refer to Table 9.3 for coordinate information) [3]](/htm/img/39/2134/467.png)
FIGURE 9.3 Impact locations and Sensor configuration set for identification (represents sensor location and impact location. Please refer to Table 9.3 for coordinate information) [3].
TABLE 9.2
Impact Locations and Sensor Locations of X- and У-Coordinates Shown in Figure 9.2 [3]
Points |
X (cm) |
Y (cm) |
IL 1 |
25 |
50 |
IL 2 |
64 |
62 |
IL 3 |
53 |
30 |
SENS 1 |
22 |
77 |
SENS 2 |
44 |
77 |
SENS 3 |
66 |
77 |
SENS 4 |
22 |
55 |
SENS 5 |
44 |
55 |
SENS 6 |
66 |
55 |
SENS 7 |
22 |
33 |
SENS 8 |
44 |
33 |
SENS 9 |
66 |
33 |
SENS 10 |
22 |
11 |
SENS 11 |
44 |
11 |
SENS 12 |
66 |
11 |
having the parameters, Equation (9.2) was reduced to evaluate the IL. Table 9.2 shows the results obtained of IL estimation for the three ILs. As an example, from the TOA. it can be inferred that SENS 4, SENS 5 and SENS 7 were the nearest sensors to IL l (Figure 9.3), and X- and У-coordinates of IL were estimated by 20% and 9% error (origin was at the bottom left corner of the plate), respectively. In all those cases, it was observed that there was a relatively slight anomaly in identifying the exact impact location. This information helped to adjust the analysis of the second layer.
Verification
In this step, IL was refined and ILTH was calculated. From the previous layer, analysis of X- and У-coordinates of impact locations were obtained and they were utilized to define the sampling space. As the number of particles increases, computation costs also increase because of the large number of particles needs a lot of calculation of the objective function, so it helped to maintain a minimum number of particles with lower cost. Surrounding sensors formed the boundary and helped to refine the sampling space. As can be seen, when the TOA in Table 9.3 show that IL 1 was contained within a boundary formed by SENS 4, SENS 5, SENS 7, and SENS 8, it implied that the upper bound for £ was 0.40. After estimating £ = 0.31 m, the closest vertical gridline had X = 0.22 m. Hence, £ = 0.22 was considered as the lower bound of the sampling space. Hence, £ = 0.31 was the mean of £ = 0.22 - 0.40. Analogous with previous, q = 0.37 and 0.55 were used as the sampling space.
TABLE 9.3
Layer 1 IL Detection Results Summary [3]
TOA (mS)
Impact locations (IL) |
||||||
Sensor |
IL 1 |
IL 2 |
IL 3 |
|||
SENS 1 |
0.34 |
0.48 |
0.59 |
|||
SENS 2 |
0.40 |
0.31 |
0.51 |
|||
SENS 3 |
0.53 |
0.29 |
0.49 |
|||
SENS 4 |
0.13 |
0.48 |
0.40 |
|||
SENS 5 |
0.25 |
0.29 |
0.31 |
|||
SENS 6 |
0.42 |
0.30 |
0.34 |
|||
SENS 7 |
0.25 |
0.57 |
0.34 |
|||
SENS 8 |
0.33 |
0.44 |
0.13 |
|||
SENS 9 |
0.48 |
0.59 |
0.25 |
|||
SENS 10 |
0.60 |
0.67 |
0.41 |
|||
SENS 11 |
0.45 |
0.60 |
0.41 |
|||
SENS 12 |
0.58 |
0.54 |
0.29 |
|||
Dominant wave frequency | ||||||
Dominant frequency* (Hz) |
2600 |
3900 |
4160 |
|||
Impact location estimation | ||||||
Impact locations coordinates (ILC) |
IL 1 |
IL 2 |
IL 3 |
|||
Predicated |
Target |
Predicated |
Target |
Predicated |
Target |
|
zeta (cm) |
31 |
25 |
56 |
64 |
52 |
53 |
eta (cm) |
46 |
50 |
65 |
62 |
25 |
30 |
*For Table 9.3. note that the dominant frequency is reported as the average of the frequencies at each sensor location due to impact at the indicated point [3].
After parameter tuning. A, Bl, and B2 matrices were as follows:

Figure 9.4 illustrates the fitness convergence of the three analyses corresponding to the three impact locations. For each particle at each iteration fitness is calculated, fitness was determined. Convergence occurred for the analyses of IL 1 and IL 3 at 70 iterations and for IL 2 at about 80 iterations.
The identified impact locations for IL 1, IL 2, and IL 3 and the corresponding reconstructed load time histories are, respectively, demonstrated in Figure 9.5. First,
![At three impact locations, fitness value vs. number of iterations for particle swarm optimization (PSO) is shown here (please refer to Figure 9.2 for impact locations) [3]](/htm/img/39/2134/469.png)
FIGURE 9.4 At three impact locations, fitness value vs. number of iterations for particle swarm optimization (PSO) is shown here (please refer to Figure 9.2 for impact locations) [3].
an approximated location with an error of l.9% - 24% was obtained through the triangulation method (Table 9.3), and then in the second layer of analysis, it was adjusted with an error ranging between 0.02% and 2.7%. ILTH showed an error for the amplitude of about 5.3%, 3.1% and 1.7%, for unloading and loading frequencies, respectively, and 3.4% for t0 at the IL2 impact location, so it matched well with the actual time history.
Summary
In this discussion, a numerical scheme based on the inverse analysis technique for a simply supported plate is proposed for IL and ILTH identification. The model proposed here offers an economic and computationally efficient approach. It includes a simplified interval analysis and heuristic algorithm for optimization so that noise effect can be considered. Even though the suggested approach is based on a mathematical forward solving model for a simply supported plate, it can also be used for other types of structures, but a suitable mathematical model might have to be created for such structures.
For verification, a 88 x 88 x 0.25 cm3 aluminium plate was used to implement the numerical scheme after being subjected to an IL, which also established the accuracy and efficiency of the proposed IL and ILTH identification model.