1. Introduction
Deformation sensing in real-time is an essential technology for providing feedback to the actuation and control systems of smart structures, especially of next-generation aircrafts. Moreover, when the detailed state of structural deformation is known, other necessary response quantities, such as stress and failure criteria, can also be assessed [
1,
2,
3]. A research of deformation sensing focuses on the utilization of in situ strain measurements, captured from a network of strain sensors, to estimate the structural deformation, which is commonly referred to as shape sensing. Specifically, fiber Bragg grating (FBG) sensors have been extensively studied for shape sensing due to the corresponding lightness, accuracy and ease of embedding [
4,
5,
6,
7].
Nevertheless, the structural deformation sensing from in situ strain data always represents an inverse problem. The inverse problem is frequently ill-posed, since it does not satisfy conditions of existence, uniqueness and stability [
8,
9]. Although many types of inverse problems and their applications have been proposed, few researchers have dealt with the deformation sensing of the wing shape. The current methods used to sense the wing shape can be divided into two categories [
10]: one category focuses on the deformation sensing of the skin comprised of plate/shell structures to reflect the deformation of the wing shape and the other category emphasizes on the deformation sensing of a wing frame made of beam structures.
The effectiveness of the modal transformation scheme was analyzed and verified by Bogert et al. [
11], in which the in-situ surface strain is used to sense the bending plate deformation. Despite the advantage of this method, accurate mode shapes and extensive eigenvalue analysis are required; for eigenvalue analysis, the elastic and inertial material properties need to be described in detail for application in high-fidelity finite element methods. On the basis of the least-squares and first-order plate theory, which includes the membrane and the bending and transverse shear deformations, a mathematical framework was developed by Tessler and Spangler to conduct the deformation sensing of plate and shell—called the inverse Finite Element Method (iFEM) [
8,
12]. The advantage of iFEM is that the structural deformations are sensed only from strain measurements, without the finite element model and a priori knowledge, such as damping, loading and material properties [
13,
14]. Through the iFEM framework and the Refined Zigzag Theory (RZT), Cerracchio and Gherlone et al. explored an effective way to sense the shape and stress of multilayered composites and sandwich structures [
15,
16,
17]. Nevertheless, the kinematic variables in the above iFEM framework are derived from the three-node inverse-shell (i3-RZT) element or three-node inverse-plate (iMIN3) element, without including the hierarchical drilling rotation degrees-of-freedom (DOF). In order to extend the practical application of iFEM for the shape-sensing analysis of large-scale structures (such as the ship), a four-node quadrilateral inverse-shell element named iQS4 was formulated by Kefal et al. that utilized the kinematic assumptions of the first-order and transverse-shear deformation theory, to add the hierarchical drilling rotation DOF into the kinematic variables [
18,
19]. In a series of works [
14,
15,
16,
17,
18,
19,
20,
21,
22,
23], the capability of iFEM framework is verified through simulation or experimentation for 3D deformation sensing of plate/shell structures of aerospace vehicles and ships. However, for the aforementioned iFEM schemes, the strain sensors need to be attached on the top and bottom of the plate/shell structures. When the sensors are placed on the exterior surface of the structure, the protection mechanism are performed on the FBG or resistance strain sensors, which affects the aerodynamic performance of the wing. To avoid this problem, the strain sensors can be placed on the wing frame covered by the wing skin. If so, the wing frame deformation can be calculated to reflect the wing shape change.
In order to conduct the deformation sensing of beam structures, Ko et al. developed a load-independent scheme to approximate the beam curvature, named Ko’s Displacement Theory [
24]. The proposed scheme employs Euler-Bernoulli beam theory to account for the beam deformation, in which the discrete surface strain measurements are integrated by using piece-wise continuous polynomials to sense the beam deformation. The experimental tests show that Ko’s Displacement Theory is sufficiently accurate for the deflection sensing of aircraft wings [
25,
26]. Regardless of its advantage, the scheme is only suitable for one-dimensional deformation of a beam structure, requiring a high number of strain sensors. According to Timoshenko beam theory, Gherlone et al. analyzed the displacement field of a constant cross-section beam, consequently constructing the relationship between the displacement and the surface strain data through iFEM [
27,
28]. Furthermore, Gherlone et al. conducted a comparative research on the aforementioned three shape sensing approaches: the iFEM, the modal transformation method and the Ko’s Displacement Theory [
10]. It is found that the iFEM method is slightly more accurate and attractive than the other two methods for shape sensing. The experimental verifications demonstrate that the accuracy of iFEM for the beam deformation sensing is severely affected by the locations where the strain sensors are placed but the optimal placement criterion of the sensors is not mentioned [
29]. Meanwhile, the stability and existence of iFEM are also severely affected by the locations where the strain sensors are placed. For instance, all sensors are placed parallel to the centroidal axes of the beam element.
The existing optimal placements of the strain sensors are based on the mode shape analysis of the structure, such as the modal kinetic energy (MKE) method, the effective independence (EI), the modal assurance criteria (MinMAC) and the drive point residue (DPR) [
30,
31,
32,
33]. Based on a mode approach, Geng et al. proposed an optimization scheme and improved genetic algorithms (GA) for optimal FBG sensor placement; consequently, the best reconstruction effects are obtained for the reconstruction of flexible plate structures [
34]. Yang et al. developed a robust optimal sensor placement for uncertain structures [
35]. In their schemes, the fitness function based on the Fisher information matrix (FIM) is derived and extended to the interval parameters. Genetic algorithm is employed to obtain the optimal solution. Being different to the traditional sensors, FBG sensors possess different detecting probe length at different angles. In view of the above difference, Yi et al. developed a detective model for FBG sensor which fit its feature based on probability model and the optimal result of the FBG sensor placement is solved by Particle Swarm Optimization (PSO) algorithm [
36]. In Reference [
37], Soman et al. proposed an optimal placement model based on Modal Identification (MI) and Accurate Mode Shape Expansion (AMSE) to place the strain sensors and accelerometers simultaneously. The results indicate that the use of a multi-type sensor system can improve the quality of Structure Health Monitoring (SHM). Given the fact that the iFEM reconstructs the structure deformation without the mode shape, all previous optimal placement schemes of the strain sensors are not suitable for iFEM.
When the strain sensors are placed on the ill-suited locations, the coefficient matrix of the equations may be ill-posed since the essence of deformation sensing with iFEM focuses on the solution of a linear system of equations. In [
38,
39,
40], the optimal placement of sensors is regarded as the parameter identification of linear system of equations, using the condition number as the optimal object to construct the optimal model of sensors. When the condition number in a matrix is higher than 10
16, the solution algorithm can return results with no accuracy at all; thus, such a matrix is numerically singular and linear systems with this matrix are not solved [
39]. It is considered that a matrix is morbid when the condition number of the matrix is higher than 10
3. Nevertheless, in our research, it is found that the condition number of the coefficient matrix in iFEM tends to fall into [10
3, 10
16] when the coefficient matrix is nonsingular; consequently, the bigger condition number in this range does not mean that the coefficient matrix is morbid. Therefore, the condition number is not suitable for the optimal sensor placement for iFEM.
This paper focuses on the construction of an optimal sensor placement model, based on eigenvalue analysis of the relationship between the sensor placement and the stability of sensing of the beam deformation with iFEM. This is conducted with an aim to maintain the sensing stability of iFEM. The contents of this paper are: firstly, the beam deformation sensing process through iFEM based on Timoshenko beam theory was reviewed. Following, the instability of sensing is discussed; consequently, the optimal placement model of strain sensors based on eigenvalue analysis is constructed. The two optimal results are obtained through the model solution with PSO algorithm. Moreover, a high-fidelity finite element model of the beam is constructed with ANSYS, which generate the discrete strain data to replace the experimental strain measurements captured from strain sensors and which produce the deformation data to replace the experimental deformation. This is executed to verify the robustness of iFEM under the optimal placement schemes. Finally, to examine the iFEM sensing capability for the two optimal placement schemes, deformation sensing test is conducted on an aluminum alloy wing frame model.
2. Inverse Finite Element Method Algorithm for Beam Deformation
In the iFEM algorithm, the displacement fields of an isotropic, straight and circular cross section beam element can be described from Timoshenko beam theory [
27,
28,
29]:
where,
are the displacements along the
axes;
,
and
denote the displacements along the center axis (
);
are the rotations around the three coordinate axes (see
Figure 1). The six kinematic variables are grouped in vector form as
.
According to the small-strain hypothesis, the linear strains are given by Equation (1):
where, the section strains
are theoretically defined as:
The iFEM senses the beam deformation shape by minimizing a weighted least-squares function
, containing the section strains vector
computed from in situ surface strain data and the theoretical section strains vector
defined through Equation (3):
In the finite element framework, the kinematic variables vector
u(x) is obtained through interpolation from certain shape functions
and the nodal degrees of freedom
as following:
With regard to different loading cases, the orders of shape functions
differ. For the end-node loads, the order of
is
; whereas for the uniformly distributed loads, the order is
[
28]. Thus, the total least-square function is a sum of
N individual element contributions
:
Taking into account the effects of axial stretching, bending, twisting and transverse shearing, the element contributions
are given by the dot product as:
with
where,
are identically set as 1;
,
,
and
are the cross-section area, the second moments of the area according to the y- and z-axis and the polar moment of the area of the beam element, respectively.
is the length of the beam element;
and
n are the axial coordinate of the locations where the section strains are evaluated and the number of locations, respectively. For the end-node loads,
and for the uniformly distributed loads,
[
27,
28,
29,
41].
The invoking of Equations (3) and (5) results in the theoretical section strains as following:
where, the matrix
contains the derivatives of the shape functions
.
Substituting Equation (9) into Equation (8), Equation (4) yields the following quadratic form:
where,
is a constant vector;
and
are defined as follows:
Finally, the minimization of function
in Equation (10) in terms of
yields the sensing equation for the beam element deformation:
The two key steps in iFEM are: (1) the selection of suitable shape functions for the beam deformation sensing with iFEM; (2) the calculation of section strains from the measured surface strain data. For step (1), literatures [
27,
28,
29,
41,
42] provide the detailed derivation process and results. For step (2), the section strains are computed through the following equation [
27,
28]:
where,
is the in-situ section strains vector at location
, along the
x-axis.
R is the external radius of the beam element.
denotes the measured surface strain at the location
, which is signed by the cylindrical coordinate system (see
Figure 2).
indicates the transformation relationship between the surface strain measurements and the section strains.
3. Construction of Optimal Placement of Sensors
When the section strains vector
is calculated from the surface strain measurements only with Equation (13), the strain measurements
(
i = 1, 2…6) and the section strains vector
must be distributed in one cross-section. Consequently,
strain measurements are required to calculate the
section strain vectors. In the case of end-node loads, 12 strain measurements are required for the calculation of
Section 2 strain vectors; in the case of uniformly distributed loads, 18 strain measurements are required for the calculation of
Section 3 strain vectors. As aforementioned in the literatures [
27,
28,
29,
41], the number of strain measurements is reduced through the constitutive equations of Equation (14) (refer to
Figure 3) and the equilibrium equations of Equation (15).
where, the section forces (
) and moments (
) are related to the section strains
.
is the axial rigidity, where
is the area of the cross-section of the beam element.
,
are the shear rigidities with
and
denoting the shear correction factors, where
G is the shear modulus. For the thin-walled section,
=
= 0.531; for the thick-walled section,
=
= 0.62; and for the solid section,
=
= 0.887 [
28].
is the torsional rigidity.
and
denote the bending rigidities.
When the forms of distribution loads
,
and
are known, the forms of the section strains
can be estimated. For the end-node loads, the section stains
,
,
and
are constant and
are linear. For the uniformly distributed transverse loads, the section strains
are constant,
are linear and
are parabolic.
[
28,
41]. Through the aforementioned results, the section strains are calculated as:
For the end-node loads, the section strains are expressed as:
or in a matrix equation:
where,
is a constant parameters vector;
indicates the transfer matrix between
and the section strains vector
.
For the uniformly distributed transverse loads, the section strains are expressed as:
or in a matrix equation:
where,
is a constant parameters vector;
indicates the transfer matrix between
and the section strains vector
. Following, the relationship between the undetermined parameters vector
and any measured strain data
is expressed as:
When the parameters vector
is solved, the arbitrary section strains vector is calculated through Equations (16) and (18).
where,
is the location where the section strains vector is calculated;
is the number of sections stated in Equation (8);
and
is the minimum number of the strain sensors used to capture the surface strains. The value of
is different under different loading cases analyzed through Equations (16) and (18). For the end-node loads,
, whereas for the uniformly distributed loads,
[
28]. When the section strains vector
is obtained, the kinematic variables
can be directly calculated from Equation (12).
However, through in-depth research, it is found that the transformation matrices
and
are ill-conditioned or even singular, when the strain sensors are placed at the inappropriate locations along the beam surface, such as
in Equation (13) set to the same value
. Thus, the deformation reconstructed from Equation (12) does not necessarily satisfy the conditions of existence, uniqueness and stability [
8].
The constant parameter vectors
and
determination is regarded as the solution of a linear system of equations. The stability of Equation (20) solution depends on whether the product matrices
and
are well-conditioned or ill-conditioned. In actual measurements, the sensor placement is different from the pre-set placement. This means that one sensor is set at the node
ideally but the actual sensor location might be at the node
and the strain measurement contains the noise. Therefore, the strain input errors are added to Equation (20):
where,
indicate the strain input errors, which couple the errors of sensor placement with the noise of the strain measurement system.
and
are the ideal strain input at the node
and the actual strain input, respectively.
,
is the errors vector, resulting from the above errors. Following, the relative error between
and
is estimated as:
where,
indicates the matrix norm and
signifies the condition number of matrix
. To a certain extent, the condition number can reflect the morbidity degree of a matrix. From Equation (23), it is observed that the relative error is high and the matrix is ill-posed, when the condition number is high. Kunsoo Huh and J.L. Stein proposed that the reason for the morbidity of a matrix is that a good distribution of eigenvalues for the matrix does not exist. For well distributed eigenvalues, large differences among the matrix eigenvalues exist [
43], such as the high value in
,
and
are the eigenvalues of the matrix). Therefore, the optimal placement model of sensors is constructed as:
where,
(
i = 1, 2,…,
m) are the eigenvalues of the matrix
.
represents the locations where the sensors are placed and
is the number of the sensors used to capture the strain data of the structure. For the end-node loads,
, whereas for the uniformly distributed loads,
[
28]. In view of application environment in engineering, it is difficult to set FBG strain sensors on the clamped node and on the free end node of the beam element. Moreover, the curve of the beam surface affects the strain measurement accuracy when a non-zero angle exists between the sensor and the generatrix of the beam surface [
29]. Thus, it is difficult to set two or more FBG sensors at one node along the surface of the beam element, as the fiber sensors are not stacked (see
Figure 4). Consequently,
is set, where only one sensor is placed at
and the other sensors are placed at
.
The optimal placement model of Equation (24) is a multi-parameter optimization problem, which is solved by using the PSO algorithm with good convergence speed. The optimal result is efficiently obtained in hyperspace.
In a population that includes
N particles, the status of the
i-th particle is described by the current location and the current velocity.
is the current location of the
i-th particle and
is the current velocity of the
i-th particle
The fitness function of PSO algorithm is defined as:
and the updated velocity and location are expressed as:
where,
are constants;
are uniformly distributed random numbers between 0 and 1;
and
arethe current location and the current velocity of the
i-th particle in the
k-th iteration, respectively;
is the individual best location of the
i-th particle at step
k;
is the local best location of the population at step
k.
and
are determined as:
where,
and
are defined as the local extremum and the global extremum of the population, respectively. The key steps of the optimal model Equation (24) solution with PSO algorithm are described as:
- Step 1.
PSO algorithm initializing. The size of the particle swarm is set to N = 50 and the maximum iteration is set to ; the initial location and the initial velocity are set to the random values included in the constraint of Equation (24). The corresponding fitness value of each particle, which is calculated from the fitness function Equation (25), is set as the initial individual extremum , (i = 1,2,…,N). Following, the initial local extremum and the global extremum of the population at the initial phase are selected for Equation (27). The threshold of the algorithm breaking is set to . When the global extremum is equal to , the algorithm breaks and the current particle location is the optimal placement of the sensors. Else Step 2 will follow.
- Step 2.
Velocity and location update of each particle with Equation (26). The fitness value of each particle , which is calculated from Equation (25), is compared to the initial individual extremum . If > , the individual extremum is replaced with ), or else, the individual extremum is invariant. Subsequently, the local extremum () and the global extremum of the population at step 1 are selected with Equation (27). When the global extremum is equal to , the algorithm breaks; or else, the iteration continues and k = k + 1.
- Step 3.
PSO algorithm termination. The iteration will not break until the iteration number = or the current global extremum is equal to .
The flow chat of above optimization procedure is shown in
Figure 5.