# The Application Research of Inverse Finite Element Method for Frame Deformation Estimation

Authors: Yong Zhao; Hong Bao; Xuechao Duan; Hongmei Fang
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

A frame deformation estimation algorithm is investigated for the purpose of real-time control and health monitoring of flexible lightweight aerospace structures. The inverse finite element method (iFEM) for beam deformation estimation was recently proposed by Gherlone and his collaborators. The methodology uses a least squares principle involving section strains of Timoshenko theory for stretching, torsion, bending, and transverse shearing. The proposed methodology is based on stain-displacement relations only, without invoking force equilibrium. Thus, the displacement fields can be reconstructed without the knowledge of structural mode shapes, material properties, and applied loading. In this paper, the number of the locations where the section strains are evaluated in the iFEM is discussed firstly, and the algorithm is subsequently investigated through a simple supplied beam and an experimental aluminum wing-like frame model in the loading case of end-node force. The estimation results from the iFEM are compared with reference displacements from optical measurement and computational analysis, and the accuracy of the algorithm estimation is quantified by the root-mean-square error and percentage difference error.


---

## Body

## 1. Introduction

Aircraft flexible wings with embedded conformal antennas and large frame structures that carry antennas require accurate real-time deformation estimation to provide feedback for their actuation and control systems [1–3]. Using the measured strain to reconstruct the shape of the structure is a key technology in the accurate real-time deformation estimation, which has been studied by many researchers.The computation of the displacement field of the deformed structure is commonly performed on the basis of strain data measured in real time by a network of strain gauges [4–6]. For example, fiber Bragg grating (FBG) sensors have been extensively researched for deformation estimation due to their lightness, accuracy, and easy embedment. The strategies that reconstruct the deformed shape displacement field of the structure with in situ strain data can be divided into two kinds. One of the two kinds trains the mapping relation between the displacement field and the measured strain by model learning algorithms, for example, neural network model and fuzzy network algorithm [7, 8]. When the training system has enough measured strains and displacement field data of the structure, a certain relation matrix can be determined and the stable relationship between the measured strain and the displacement field can be obtained. But the strategy requires a large number of training data, and the mapping relation is easy to fail when the actual loading is beyond the range of the training cases.The other establishes the mapping relation between the displacement field and in situ strain data without the model learning. In the literature [9–11], the global or piecewise continuous basis function methods were employed to fit the surface-measured strain into the structure strain field, and then, the structure deformation displacement was obtained from the strain-displacement relationship. These methods are easy to implement, but the reconstruction accuracy of deformation estimation depends on the appropriate selection of basis function and weight coefficients. Mode shapes have been used as basis function in [12, 13]. The deformation displacements are reconstructed from measured strains by using the modal transformation method. However, there exist the following disadvantages in this method. (1) The detailed material elasticity and inertial parameters are needed to precisely construct mode shapes. (2) The accuracy of deformation reconstruction is severely limited to the modeling precision of the structure, and it is quite difficult to precisely model the complex structure. Maincon [14] developed a finite element-based methodology involving an inverse interpolation formulation that employs the surface-measured strain to determine the loads and structural response of aerospace vehicles, while this algorithm needs an appropriate quality function to adapt the different loading cases and this function is constructed based on a mass of computer simulation and experimental statistics.Based on the Euler-Bernoulli beam equation, Jute et al. evaluated the deflection of beam by the integration of discretely measured strains directly [15]. The algorithm applied the classical beam equation and piecewise continuous polynomials to approximate beam curvature through integration. Derkevorkian et al. [16] compared the algorithm with the modal transformation method and demonstrated the additional benefits of the algorithm to achieve a robust method for monitoring ultralightweight flying wings or next-generation commercial airplanes. Though this one-dimensional scheme has displayed high accuracy in predicting deflection, it fails to estimate the element deformation under multidimensional complex loads.Tessler and Spangler [17] proposed the inverse finite methodology (iFEM), which can be used to reconstruct the displacement field of shear-deformable structures, not only beam structure but also plate and shell structures. The main idea is reconstructing a three-dimensional displacement field of beam structure from the surface-measured strains according to a least squares approach. Due to the fact that only the displacement-strain relationship is used, the deformation reconstruction can be accomplished by the methodology without the prior knowledge of loads, materials, and inertial and damping properties. To model arbitrary plate and shell structures [18, 19], Tessler and his partners developed the iFEM algorithm using the first-order shear deformation theory and a three-node inverse shell element. FBG sensors were applied to measure the surface strains on the slender beams, and then, the deformed displacement was reconstructed by using an iFEM shell model.The beam-deformed displacement and cross-sectional torsion were reconstructed by Gherlone et al. who employed the inverse finite element formulation to achieve high reconstruction accuracy of deformed displacement [20–21]. The authors first used the Timoshenko beam theory to model the beam kinematic accurately, then used the C0 or C1 inverse frame elements and least squares formulation to establish the relationship between the measured strain and the arbitrary node displacement field of the beam element with no prior knowledge of the finite element model and loads. The reconstruction equation is nonsingular when the boundary conditions are applied; that is, the status of one end node of the beam element must be known. As the displacement field of the beam element in iFEM is constructed based on an isotropic straight beam structure, whose cross section is invariable along the whole beam, the deformation of the tapered beam cannot be estimated. Meanwhile, in view of different load cases, the authors discussed different types of shape functions that are used to interpolate the kinematic variables in the beam element; then, the number of surface strain measurements and the deformation reconstruction equation are confirmed. Unfortunately, the minimum number of locations where the section strains are evaluated, which is critical for obtaining a correct solution to the algorithm, is not discussed clearly in different load cases. If the minimum number of the locations where the section strains are evaluated is inappropriately set, the solution will be nonunique; that is, the reconstructed displacement may be incorrect.The contribution of this paper is twofold. (i) The relationship between the loading case and the number of the section strain locations is discussed in detail, and the following conclusions are verified; on the one hand, the minimum number of locations should be 2 in the loading case of end-node forces; on the other hand, the minimum number should be 3 under the uniform distributed loading. (ii) In order to examine the minimum number which is 2, estimation studies are carried on a beam structure and a wing-like three-dimensional frame structure, in the loading case of end-node force. Six fiber optic strain sensors are placed at two nodes of each beam to capture in situ strain dates. To assess the estimation error of the iFEM algorithm, the RMS and difference percent parameters are employed.


## 2. iFEM Algorithms for Beam Estimation

A beam deflection estimation algorithm was developed by Gherlone and his collaborators [20]. In the algorithm, the Timoshenko beam theory is first used to analyze the expression of the displacement field in a straight beam element (Figure 1).
(1)uxx,y,z=ux+zθyx−yθzx,vyx,y,z=vx−zθxx,uzx,y,z=wx+yθxx,where ux,vy, and wz are the point displacements along the x, y, and z axes, respectively. ux,vx,andwx denote the displacements at y=z=0; θxx, θyx, and θzx are the rotations about the three coordinate axes. The six kinematic variables in the middle axes can be grouped in vector form as follows:
(2)u=u,v,w,θx,θy,θzT.Figure 1
Beam geometry and kinematic variables.The arbitrary section strainseu=e1,e2,e3,e4,e5,e6T can be obtained by (1).
(3)e1x=uxx,e2x=θy,xx,e3x=−θz,xx,e4x=wxx+θyx,e5x=vxx−θzx,e6x=θx,xx.The six kinematic variablesu can be interpolated by the right shape functions
(4)u=Nxue,where Nxandue denote the shape function and nodal degrees of freedom, respectively. Substituting (4) into (3) gives arbitrary section strains in terms of the nodal degrees of freedom as follows:
(5)eu=Bxue,where the matrix Bx=B1x,B2x,…,B6x contains the derivatives of the shape functions Nx. Once the section strains eu are obtained, the nodal displacement ue is determined; then, the kinematic variables u can be acquired by (4). However, the section strains eu are derived from the kinematic variables u theoretically, rather than the strain measurements. So, iFEM uses in situ section strains eε computed from the measured strains to replace eu when the least squares error function φu reaches the minimum.
(6)φu=eu−eε2,eε=e1ε,e2ε,e3ε,e4ε,e5ε,e6εT.In view of the effect of the axial stretching, bending, twisting, and transverse shearing, the improved least squares error functionalφeu is obtained by the dot product of the weighting coefficient vector W and the original vector φu=φke.
(7)W=wk=w10,w20IyeAe,w30IzeAe,w40,w50,w60IpeAe,φeu=W⋅φu,φke=Ln∑i=1nekxi−ekεxi2,k=1,2,…,6,where wk0k=1,2,…,6 denote dimensionless weighting coefficients whose initial values are identically set as 1; Ae, Iye and Ize, and Ipe are, respectively, the cross-sectional area, second moments of the area according to the y- and z-axes, and polar moment of the area of the beam element. L is the length of the beam element; xi0≤xi≤L and n are, respectively, the axial coordinate of the locations where the section strains are evaluated and the number of locations, that is, the axial coordinate of sections where section strains are distributed in, and the number of sections.For a straight beam member of constant circle cross section, the in situ section strains atx=xi, eεxi=e1εxi,e2εxi,…,e6εxiT, could be derived from the measured strains of the beam surface with appropriate strain-tensor transformations from the θ,x,r to xX,xY,xZ coordinates [22].
(8)ε2∗xi,θ,β=e1εxicβ2−vsβ2+e2εxicβ2−vsβ2sθR+e3εxicβ2−vsβ2cθR+e4εxicβsβcθ−e5εxicβsβsθ+e6εxicβsβR,with cβ≡cosβ,sβ≡sinβ,cθ≡cosθ,andsθ≡sinθ.In (8), v is the Poisson ratio, ε2∗xi,θ,β is the in situ strains that are obtained from strain sensors. R denotes the radius of the beam cross section (Figure 2).Figure 2
Location and coordinates of a strain gauge placed on the beam external surface.Substituting (5) into (7) results in the following quadratic form:
(9)φeu=12ueTkeue−ueTke+ce.Herein, ce is a constant, and ke and fe are indicated as follows:
(10)ke=∑k=16wkkke,kke=Ln∑i=1nBkTxiBkxi,fe=∑k=16wkfke,fke=Ln∑i=1nBkTxiekεxi.Finally, the relationship between the deformation and in situ section strains, shown in (11), can be confirmed when the minimization of functional φeu is performed.
(11)keue=fe.Once the appropriate shape functions and the problem-dependent displacement boundary conditions (e.g., setting the displacement of one end point to zero, which means that one of end nodes of the beam is fixed) are given,ue can be derived from a nonsingular system, and the vector fe depends on the measured strain values that change during deformation. Once the nodal degrees of freedom ue are confirmed, the displacements and rotations of every node along the centroid axis of the beam element are obtained by (4) and the deformed shape of the whole beam can be reconstructed by (1).


## 3. Minimum Number of the Section Strain Location Selection

Although the displacement shape function and the number of surface strain measurements in different loading cases have been discussed in [20, 21], n, the number of the sections where the section strains are distributed in, is not shown clearly. In (11), once the problem-dependent displacement boundary conditions are determined, the nonsingularity of the equation will depend on the selection of n. In our investigation, it is found that n can be determined by a specific loading case. More specifically, the orders of the section forces and the moments can be obtained by the loading case, which determines the order of the section strains. With the value of the section strain order, n is determined immediately.With the equilibrium equations (12), the section forces (N, Qy, and Qz) and moments (Mx, My, and Mz) can be obtained by the load concentrations qxx, qyx, and qzx along the x, y, and z directions (Figure 3).
(12)Nx+qxx=0,Qy,x+qyx=0,Qz,x+qzx=0,Mx,x=0,My,x−Qz=0,Mz,x−Qy=0.Figure 3
Beam section forces and moments.The relationship between section forces and moments and the section strainseii=1,2,…,6 can be interpreted as the following constitutive equations:
(13)N=Axe1,Mx=Jxe6,Qy=Gye5,My=Dye2,Qz=Gze4,Mz=Dze3.where Ax≡EA is the axial rigidity; Gy≡ky2GA and Gz≡kz2GA are the shear rigidities, with ky2 and kz2 denoting the shear correction factors; and Jx≡GIp and Dy≡EIy and Dz≡EIz are, respectively, the torsional rigidity and the bending rigidities. For the uniform section beam element, parameters mentioned above are constant. Then, the order of the section strains ei will be identical to that of the section forces and moments.Substituting (13) into (12) gives the relationship between the section strains and the loads as follows:
(14)e1=NAx=∫−qxxdxAx,e2=MyDy=∫QzdxDy=∫∫−qzxdxdxDy,e3=MzDz=∫QydxDz=∫∫−qyxdxdxDz,e4=QzGz=∫−qzxdxGz,e5=QyGy=∫−qyxdxGy,e6=MxJx=∫0dxJx=c1.In most cases, the load in thex direction is zero, that is, qxx=0; then,
(15)e1=∫−qxxdxAx=∫0dxAx=c2.The section strainse1 and e6 are constant in (14) and (15), and the order of the residual section strains, e2,e3,e4,ande5, will be discussed as the following two loading cases.Consider a beam element loaded by the end-node forces,qyx and qzx are zero along the x-axis. Then, e2,e3,e4,ande5 are determined by (14).
(16)e2=∫∫−qzxdxdxDy=∫∫0dxdxDy=a1x+b1,e3=∫∫−qyxdxdxDz=∫∫0dxdxDz=a2x+b2,e4=∫−qzxdxGz=∫0dxGz=a3,e5=∫−qyxdxGy=∫0dxGy=a4.where a1,a2,a3,a4,b1,andb2 are unknown constant parameters.As in (15), the highest order of the section strains eu is linear e2,e3, which means that the distribution of the bending moments is a skew line, that the corresponding errors in (7) can be obtained from the shaded area between two skew lines, and that each line can be confirmed by two different nodes, that is, two section strains in different sections (Figure 4(a)). Thus, n=2 for e2ande3 in (7). Meanwhile, the other section strains, e1,e4, e5,ande6, are constant, which means that their distributions are lines parallel to the x-axis, that the associative errors in (7) can be obtained from the shaded area between two parallel lines, and that every line can be confirmed by one node, that is, one section strain (Figure 4(b)). Thus, the number of the sections where section strains are distributed is n=1. Finally, the number of the sections, n, is 2 in (7) and (10) in the loading case of end-node force and moments.Figure 4
Distribution and error of section strains in different orders.eε is a section strain in eε computed from surface strain measurements by (8). e is a section strain in eu deduced by (3).
(a)
Constant section strain
(b)
Linear section strain
(c)
Quadratic section strainAnother case is the beam element loaded by uniformly distributed forces,qyx and qzx are constant along the x-axis, and the order of the section strains e2,e3,e4,ande5 can be deduced as follows.Assumption:qyx=d1, qzx=d2.Then,(17)e2=∫∫−qzxdxdxDy=∫∫−d2dxdxDy=−d22Dyx2+d3x+d4,e3=∫∫−qzxdxdxDz=∫∫−d1dxdxDz=−d12Dzx2+d5x+d6,e4=∫−qzxdxGz=∫−d2dxGz=−d2xGz+d7,e5=∫−qyxdxGy=∫−d1dxGy=−d1xGy+d8.Herein,dii=1,…,8 are unknown constant parameters. The highest order of the section strains eu is quadratic e2,e3, and then, the distribution of the bending moments is a parabola. The corresponding errors in (7) can be obtained from the shaded area between two parabolas which are confirmed by three different nodes (Figure 4(c)), that is, three section strains in three different sections. Thus, n=3 for e2ande3 in (7). Meanwhile, e4ande5 are linear, where the distribution is a skew line, and e1ande6 are constant. Similar to the loading case of end-node force, the corresponding numbers of sections are 1 and 2, respectively. Finally, the number of the sections, n, is 3 in (7) and (10) in the loading case of uniformly distributed forces.


## 4. Verification

A simple cantilevered solid beam and a three-dimensional frame structure were subjected to the end-node static loads to assess the iFEM potential for the flexible wing deformation estimation. The beam and frame structures were made of 6061-T6 aluminum alloy. The Young’s modulus isE = 73000 Mpa, the Poisson ratio is v=0.3, and its density is ρ=2712.63kg/m3. The frame was composed of three solid beams and middle plates; the span of each beam L = 640 mm and the radius R = 10 mm (see Figure 5). For the solid circular cross section, the shear correction factors are ky2=kz2=0.887 [21]. In our verification, every whole solid beam is regarded as one beam element. Accordingly, the principle axis of frame structure is divided into three elements.Figure 5
Loading case and the distributions of NDI identification points.
(a)
Clamped-free beam and loading
(b)
Identification points of NDI
(c)
Frame structure and loadingThe experimentally measured surface strains are obtained by fiber optic strain sensors. Displacement measurements are captured by 3D optical measurement instruction (NDI, Optotrak Certus), which determines the position of the identification point by using three CCD cameras to capture the infrared lights emitted by position sensors (Figure6). The instruction is also used to assess the iFEM-recovered deflections, where the accuracy of the NDI is 0.1 mm. A position sensor was placed as close as possible to the beam root to verify the effectiveness of the clamping arrangement. The force is achieved by placing several weights on a pothook.Figure 6
Static test setup. The single beam (A), frame structure (B), and NDI Optotrak Certus (C).For the frame structure, each beam is regarded as an element and the optical fiber sensor location scheme used for each element is identical. For the end-node static loads, the displacement field of the whole beam element is interpolated byC0 continuous shape function and the number of required strain sensors is 6 (see [21]). As the radius of the beam element is small (R = 10 mm) and the grid length of every fiber grating sensor is 10 mm, it is difficult to stick the six strain gauges on one section; six fiber optic strain gauges are placed on two different sections along the beam (Table 1).Table 1
Optical fiber sensor location in one element.
e
ε
1
x
,
θ
,
β
e
ε
2
x
,
θ
,
β
e
ε
3
x
,
θ
,
β
e
ε
4
x
,
θ
,
β
e
ε
5
x
,
θ
,
β
e
ε
6
x
,
θ
,
β
L
/
2
,
−
2
π
/
3
,
0
L
/
2
,
0
,
0
L
/
2
,
2
π
/
3
,
0
4
L
/
5
,
−
2
π
/
3
,
0
4
L
/
5
,
0
,
π
/
4
4
L
/
5
,
2
π
/
3
,
0The accuracy of the reconstitution is assessed by root-mean-square (RMS) and percentage difference (%Diff).(18)RMS=∑i=1mτiFEMxi−τNDIxi22,%Diffτ=100∗RMSmaxτNDIxi,where τ=v,w is the deformation displacement along the y- or z-axis and m is the number of deflection shape displacement measurement identification points of NDI (Figure 5(b)A) in the structure. For the beam, m=6 (Figure 5(a)), and for the frame, m=14 (Figure 5(c)). The superscript “iFEM” refers to the reconstitution by iFEM while “NDI” refers to the experimental measure from 3D optical measurement instruction; xi is the ith location along the axis where the displacement u is measured.At first, five different vertical force cases (Table2) were applied at the free end of the beam (Figure 5(a)). Figure 7 shows the comparison between iFEM reconstruction and NDI measurement along the beam axis in force direction, for the loading case 68 N. Symbol x/L indicates the location of identification point along the beam surface. Zero means identification point at the clamped end node, and 1 means identification point at the free end node. v(y) indicates the corresponding displacement in force direction.Table 2
Loading case for the clamped-free beam. Max deflections are captured from NDI.
Case 1
Case 2
Case 3
Case 4
Case 5
Loads
5 N
15 N
37 N
41 N
68 N
Max deflection inY
−0.79 mm
−2.5 mm
−6.11 mm
−6.82 mm
−10.96 mm
RMS inY
0.05 mm
0.16 mm
0.32 mm
0.33 mm
0.42 mm
%Diff(v) in Y
6%
6.4%
5.2%
4.8%
3.8%Figure 7
Comparison of iFEM reconstruction for a beam to NDI measurement inY.For the frame structure, five different vertical force cases (Table3) were applied at the free end of the frame (Figure 5(c)). The two deformation contour plots of the frame are plotted by high-fidelity direct FE analysis (ANSYS 12.0), which shows the comparison between iFEM reconstruction and NDI measurement along the beam axis in force direction, for the loading case 45.5 N (Figure 8). The deformation unit in the two figures is meter.Table 3
Loading case for the frame structure. Max deflections are captured from NDI.
Case 1
Case 2
Case 3
Case 4
Case 5
Loads
18 N
23 N
35 N
45.5 N
68 N
Max deflection inY
16.31 mm
20.5 mm
31.2 mm
41.3 mm
62.8 mm
RMS inY
0.22 mm
0.49 mm
2.0 mm
3.3 mm
6.8 mm
%Diff(v) in Y
1.3%
2.4%
6.4%
8.1%
10.9%
Max deflection inZ
2.53 mm
3.23 mm
5.1 mm
6.5 mm
10.2 mm
RMS inZ
0.32 mm
0.41 mm
0.73 mm
0.75 mm
1.1 mm
%Diff(w) in Z
12.6%
12.8%
14.6%
11.5%
10.8%Figure 8
Comparison of iFEM reconstruction for a frame to NDI measurement.
(a)
Contour plot of iFEM reconstruction
(b)
Contour plot of NDI measurementIt is seen from the preceding results that the iFEM methodology shows a good potential as a reliable estimation technique. The estimation results gained from the iFEM algorithm have a better approximation to the NDI measurement results. For the beam test, the accuracy of the reconstitution assessed by percentage difference does not exceed 6.4% and the error brings down with the loading increase. For the frame test in the loading case of 18 N, the percentage difference is 1.3%, but the difference grows up with the loading increase; especially from the loading case 23 N to 35 N, the difference increases 4%. The cause for this phenomenon is that the gripper of the frame structure is not very stable, which leads to the fact that the gap spacing between the frame and its carrier grows up with the loading increase that has a great influence on the boundary conditions.


## 5. Conclusion

This study investigates the application potential of the iFEM algorithm for the flexible wing and other frame structures. This method employs the in situ strain measurements and proper displacement shape function to estimate the deformation shape of the beam and frame structure, without the need to know the applied loads and material properties and to use modal shapes. When the boundary conditions of the reconstruction equation and the shape function are determined, the deformation shape will be accurately estimated in a nonsingular system. Although the boundary conditions are known in hypothesis, the number of location where the section strains are evaluated is not shown clearly; thereby, the singularity of the equation will be influenced.This paper discusses the minimum number of location where the section strains are evaluated with the loading case, using equilibrium equations in detail, and verifies the result that for the loading case of end-node force, the number is 2. The results of the static loads in the beam and the frame model show that the iFEM algorithm has a good potential for the flexible frame structure estimation; especially for a simply supported beam, the accuracy of the estimation by iFEM is around 6%, and for the frame structure in small load (18 N), the accuracy reaches up to 1.3%, while as the gripper of the frame is not firm enough, the estimation accuracy brings down with the loading increase. So, the next work is to explore the technology that reduces the influence of the gap spacing of the frame gripper on accuracy of the frame deformation estimation.



---

*Source: 1326309-2017-11-05.xml*