A RealTime and Wireless Structural Health Monitoring Scheme for Aerospace Structures Using Fibre Bragg Grating Principle
INTRODUCTION
The population of artificial satellites and nearEarth 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 realtime using fibre Bragg grating sensors. The developed scheme was first analytically modelled, then implemented, and finally validated through inlab 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 twodimensional surface using fibre Bragg grating sensors onboard, 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 realtime 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 nonstochastic inverse analysis was used to address the uncertainty in detecting, identifying and locating ILs on a simply supported twodimensional 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 realtime 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 multipoint 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 wellestablished 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 timevarying 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 modelbased forward solving of a mathematical model or by using a neural network that requires extensive training of any given dataset. The modelbased approach is well accepted among researchers due to its convenience of implementation and simplicity of calculation compared to a neural networkbased approach. Another advantage of the modelbased 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 nonlinear 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 realtime detection.
In this research, the IL and ILTH were concurrently obtained, and the uncertainty was addressed by the inverse analysis which is nonstochastic. 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.
ThreeLayer Identification Process
A threelayer 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 t_{i} 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
At_{jm}: measured value of differential TOA.
ДГ_{1С}: predicted value of differential TOA,
x, : ^coordinate of the /'th sensor,
y, : уcoordinate of the /'th sensor,
£: .vcoordinatc of the IL,
//: уcoordinate of the IL,
V: wave speed and
S, S_{2} 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 directiondependent 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: KirchhoffLove’s solution and Mindlin's solution. KirchhoffLove’s formulation only considers the inplane shear. In the case of Mindlin’s solution, it also considers the outofplane 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 (e_{x}), Soares [4] as w'ell as Carvalho and Dobyns [4] methodology can be adapted. The relationship between modal amplitudes, the timede pendent function of bending deflection would be used for the expression[4], which is
where z is the respective coordinate to midplane, 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 = vcoordinate of IL and o)_{f/} = modal frequency of plate.
K_{c} and K_{:/} are next defined as follows:
where и is the length and v is the width of the location (impact area) and K_{A} 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 halfcycle sinusoidal function as in (9.8):
where
F_{(l} = IL amplitude, со,, = unloading frequency, ы, = loading frequency, t_{0} = impact start time, r, = time of impact and t_{2} = impact end time.
It should be noted here that if the measurement starts while impact occurs, t_{0} 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,, t_{2}, 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 e_{ic}(t) is the predicted strain and e_{im}(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 6dimensional 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
X_{k} = particle position after Ath iteration,
X_{M} = particle position after the (A+ l)^{lh} iteration.
V_{k} = the particle’s speed after the Ath iteration,
V_{k+I} = the particle’s speed after the (A+ 1 )^{lh} iteration,
P, = the particle’s past best position,
P_{2} = the finest position experienced by other particles in the herd,
#), B_{2} = two factors relates to the personal influence and social influence of particles,
A = momentum factor, r, and r_{2} = 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 B_{2}. The longterm 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, finiteelement 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, finiteelement 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 100time 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.
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: D_{12} = 34.54Pa; D_{66} = 33.52Pa; A,, = 195 MPa; A+, = A_{55} = 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
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,
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 t_{0} 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 cm^{3} 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.