# Geometrically Nonlinear Aeroelastic Stability Analysis and Wind Tunnel Test Validation of a Very Flexible Wing

Authors: Changchuan Xie; Yi Liu; Chao Yang; J. E. Cooper
Journal: Shock and Vibration
 (2016)

---

## Abstract

VFAs (very flexible aircraft) have begun to attract significant attention because of their good flight performances and significant application potentials; however, they also bring some challenges to researchers due to their unusual lightweight designs and large elastic deformations. A framework for the geometrically nonlinear aeroelastic stability analysis of very flexible wings is constructed in this paper to illustrate the unique aeroelastic characteristics and convenient use of these designs in engineering analysis. The nonlinear aeroelastic analysis model includes the geometrically nonlinear structure finite elements and steady and unsteady nonplanar aerodynamic computations (i.e., the nonplanar vortex lattice method and nonplanar doublet-lattice method). Fully nonlinear methods are used to analyse static aeroelastic features, and linearized structural dynamic equations are established at the structural nonlinear equilibrium state to estimate the stability of the system through the quasimode of the stressed and deformed structure. The exact flutter boundary is searched via an iterative procedure. A wind tunnel test is conducted to validate this theoretical analysis framework, and reasonable agreement is obtained. Both the analysis and test results indicate that the geometric nonlinearity of very flexible wings presents significantly different aeroelastic characteristics under different load cases with large deformations.


---

## Body

## 1. Introduction

Large-aspect-ratio wings produce a high lift-drag ratio and good flight performance for high-altitude long-endurance unmanned aerial vehicles (HALE UAVs). Advanced composite materials make wing structures lightweight but also introduce large elastic deformations during aerodynamic loads in flight, which induces significant geometric nonlinearity in aeroelasticity and flight dynamics. In particular, after the NASA Helios mishap, the shortcomings of traditional linear analysis were deemed unsuitable for VFAs with large deflections, and the following recommendation was made by the investigators of the accident [1]:[There is a need to] develop more advanced, multidisciplinary (structures, aeroelastic, aerodynamics, atmospheric, materials, propulsion, controls, etc.) time domain analysis methods appropriate to highly flexible, morphing vehicles.Structural geometric nonlinearity, a key feature of HALE UAVs, yield a nonlinear relationship between displacement and strain due to large elastic deformations. Because the structural stiffness and equilibrium equations all depend on instantaneous structural deflections and load condition, the aerodynamic calculations and dynamic equations should be established for a large deformed state [2]. Thus, geometrically nonlinear aeroelasticity could be defined as a subdiscipline of aeroelasticity that considers nonlinear large structural deformations and the aerodynamics on curved aerosurfaces simultaneously.In traditional linear analysis, aircraft were not particularly flexible, and the geometric nonlinearity were not significant; thus, linear structural finite elements based on the small displacement assumption combined with the planar doublet-lattice method [3] were widely used in engineering analyses and even imbedded in commercial software [4]. As aircraft have become more flexible, researchers have found linear aeroelastic analysis methods to be inaccurate and many new analysis methods for VFAs have been brought up. Hodges developed a nonlinear geometric exact beam element [5], reducing the order of structural nonlinearity, to analyse the nonlinear aeroelastic response of VFAs. He reported that structural geometric nonlinearity had a significant effect on the structural dynamics and dynamic aeroelastic characteristics of a high aspect ratio wing. The geometrically exact calculation of the angle of attack and aerodynamically consistent application of the air loads was also important for accurate aeroelastic characterisation [6]. Combined with Peter’s 2D inflow theory [7], an aeroelastic analysis toolbox named NATASHA was developed to analyse the nonlinear behaviour of VFAs. The nonlinear beam model and reduced order model (ROM) combined with 2D or quasi-3D aerodynamic theory help researchers to predict VFA performance [8–10] but prevent wide application in industry because a real structure is overly complex (i.e., it cannot be simply represented with a beam model); thus, it is necessary to find a 3D aerodynamic code to manage VFA aerodynamic computations. To obtain accurate aerodynamic loads, some researchers used CFD/CSD [11] methods to analyse the geometric nonlinearity of VFAs. Smith et al. [12] loosely coupled the Euler solver with geometric exact beam structural analysis and investigated the effects of adding aerodynamic nonlinearity to the elastic behaviour of high aspect ratio wings. For VFA stability analysis, considering the high computing costs of time domain aerodynamic computations and the importance of highlighting basic aeroelastic principles for unconventional wings with high aspect ratios [13], dynamic flexible motion of the system can be assumed with small amplitudes around a nonlinear static equilibrium state [14]; thus the linearized method and frequency domain solution are still a valuable approach in preliminary designs and even in the detailed design stage. Panel aerodynamic methods have been well understood and widely used in engineering design; extending the panel aerodynamic code into 3-Dimentional application can make geometrically nonlinear aeroelastic analysis easy to accept in industrial applications.Based on the above discussion, in this paper, a practicable geometrically nonlinear analysis system is established and a wind tunnel test is conducted to validate the analysis results. Nonlinear finite element method (FEM) is utilized so that the method could be easily used in industry. Actually there are some refined modelling methods like intrinsic beam [5], strain-based beam and plate elements [15], nonlinear substructure, ROM, and so forth. These methods are well developed theoretically and can illustrate the flexible structural characteristics to different degrees. Considering the practical problems we have met, the structures are often so complex that these modelling methods are not convenient to apply. Additionally, these simplified methods cannot well reflect the original structure especially in detailed parts. Moreover, FEM is the most often used modelling method in practical analysis and it will be very efficient and convenient if the FEM model can be directly used to analyse the VFAs structural geometric nonlinearity. The purpose of this paper is to establish an analysis framework easy to implement and able to reveal some nonlinear aeroelasticity of flexible wings. Based on the concerns above, the nonlinear FEM is responsible for structural nonlinear analysis. As to the aerodynamic computation, the importance of the nonplanar aerosurface effect and exact aerodynamic modelling consistent with structural deflection comes from our experience and many reference papers. Structural nonlinearity and nonplanar aerodynamic compotation must be considered simultaneously. Considering the easy programming and good inheritance of conventionally used linear method, NVLM (nonplanar vortex lattice method) and NDLM (nonplanar doublet-lattice method) are adopted and they can well present the nonplanar effect of aerodynamic for flexible wings in stability analysis. All of these methods together can describe the geometric nonlinearity of both the aerodynamics and structure of VFAs and can be conveniently used in industrial design. Wind tunnel test under different load cases indicates that different deformations under varying flight states result in different flutter speeds that may alter the aircraft’s flight envelope. Reasonable agreement between the analysis results and test results has been obtained, and all the results demonstrate that the static aeroelastic response and flutter characteristics are quite different from the results obtained from linear analysis.


## 2. Theory

### 2.1. Structure Geometric Nonlinearity

Large structural excursions induced by very flexible wings when undergoing aerodynamic loads in the air prevent the use of linear methods based on the small displacement assumption and call for geometrically nonlinear structural analysis. Geometric nonlinearity are based on the kinematic description of the body, and the strain on the wing should be defined in terms of local displacement of the wing for dynamic motion. These result in the nonlinear geometric equations including the quadric term of the displacement differential and require the nonlinear force equilibrium equation established on the deformed state of the structure. Structural geometrically nonlinear problems are often solved by the maturely developed nonlinear incremental finite element method [16] and two formulae called the total Lagrange formulation (TLF) and the updated Lagrange formulation (ULF) [17], which are well known. The ULF is used in this study, and the primary equations are presented briefly below. The core method of structural analysis has already imbedded in commercial software, so it is convenient to be used in engineering analysis.The relationship between the nonlinear Lagrange/Green strain and displacement is(1)
ε
i
j
t
=
1
2
u
t
i
,
j
+
u
t
j
,
i
+
u
k
,
i
t
u
k
,
j
t
,where u
i
,
j
t is the partial derivative of displacement component u
i to the coordinate x
j at time t.Despite a large elastic deformation, the material is still within the elastic limitation for a small strain, so the conjugate Kirchhoff stress tensorS
j
i at timet satisfies(2)
S
j
i
t
n
j
t
d
s
t
=
x
i
,
j
t
d
T
j
t
,where n
j
t is the direction cosine of a small area element d
s at time t and d
T
j is the corresponding surface force in which the follower force effect is considered. The linear elastic constitutive relation is given as follows:(3)
S
i
j
t
=
D
i
j
k
l
t
ε
k
l
t
,where D
i
j
k
l is the elastic tensor, which has a different form for isotropic or anisotropic material.The finite element method (FEM) based on energy principles is an effective approach to solve structural problems. For geometric nonlinear problems, considering the follower force effect, the incremental FEM is used. The strainε
i
j can be decomposed into a linear part e
i
j and a nonlinear part η
i
j of the current displacement:(4)
ε
i
j
t
=
e
i
j
t
+
η
i
j
t
.The stress can also be decomposed in increments, where S
i
j
t represents the equilibrium stress at time t and τ
i
j
t represents the incremental stress to be calculated at each time step:(5)
S
i
j
t
+
Δ
t
=
τ
i
j
t
+
S
i
j
t
.The integral equation is established by linearization in each incremental step:(6)
∫
V
t
D
i
j
k
l
t
e
k
l
δ
e
i
j
t
t
d
V
t
+
∫
V
t
τ
i
j
t
δ
η
i
j
t
d
V
t
=
Q
t
+
Δ
t
-
∫
V
t
τ
i
j
δ
e
i
j
t
t
d
V
t
,where Q
t
+
Δ
t is the incremental outer force including the aerodynamic force, engine thrust, and gravity, at the new time step. Considering a number of shape functions, the relationship between strain and deformation is presented as(7)
e
t
=
B
L
t
u
,
η
t
=
B
N
L
t
u
.Substituting these shape functions into (6) leads to the element governing equation [18]: (8)
K
L
t
+
K
N
L
t
u
=
Q
-
t
+
Δ
t
-
F
A
t
,where Q
-
t
+
Δ
t is the incremental outer force including the aerodynamic force, engine thrust, and gravity at the new time step. The stiffness matrix in (8) could be decomposed into a linear part K
L
t and nonlinear part K
N
L
t. The linear part is only related to the structure itself, whereas the nonlinear part is related to the deflected configuration and strain quality, which should be updated in each computation step.The corresponding dynamic equation can be expressed as(9)
M
u
¨
t
+
Δ
t
+
K
L
t
+
K
N
L
t
u
=
Q
-
t
+
Δ
t
-
F
A
t
,where u
¨
t
+
Δ
t is the structural displacement acceleration vector at new time step. The assumption of a small amplitude vibration around the static equilibrium state is suitable for many dynamic problems, including the flexible wing structural dynamic stability:(10)
u
=
u
-
+
x
,where u
- is the large deflecting equilibrium deformation from (8) and x is a small vibration deformation. According to (9) and the static equilibrium condition, the linearized structural quasimode can be obtained by generalized diagonalization, and the vibration equation of the system under steady forces reduces to(11)
M
T
x
¨
+
K
T
x
=
0
,where M
T is the inertial matrix of the structure at the nonlinear static equilibrium configuration and K
T is the corresponding stiffness matrix. Both of these parameters are nonlinear functions of u
- and vary under different equilibrium states, which is a key feature of geometric nonlinear structures. The mode shapes and frequencies can be deduced from (11).Introducing the harmonic oscillating assumptionx
=
ϕ
e
i
ω
t, the vibration equation can be written as (12)
K
T
-
ω
2
M
T
ϕ
=
0
,where ω is the vibration circular frequency and ϕ is the vibration mode matrix. If (12) has all-nonzero solutions, that demands (13)
det
⁡
K
T
-
ω
2
M
T
=
0
.That is a generalized eigenvalue problem about K
T and M
T. Solving (13), the vibration circular frequency ω
i should be obtained. Substituting ω
i into (12), structural eigenvector ϕ
i (structural mode shape) can be obtained.


### 2.2. Nonplanar Vortex Lattice Method

The modelling of flexible aircraft with significant structural deformations requires the incorporation of structural dynamics and aerodynamics in a unified framework. The aerodynamic loads in (9) should be computed on a deformed configuration. This section summarizes the primary characteristics of the NVLM and its application to a curved aerosurface. The NVLM is derived from three-dimensional potential flow theory and is suitable for most normal situations that VAFs may encounter. Because the simple programming effort is required, NVLM is easy to combine with structural dynamic computations to obtain the response results for aeroelastic structures. Additionally, the exact boundary condition is satisfied on the real wing surface in the NVLM, which can have camber and various platform shapes. Thus, the NVLM is convenient for use with very flexible wings, whose aerodynamic surfaces are subjected to large structural deformations [19]. The NVLM is implemented using vortex ring quadrilateral elements to discretize the curved lifting surface along with the wing’s deformation, as shown in Figure 1.Figure 1
Nonplanar vortex lattice model for a thin lifting surface.A Cartesian coordinate system is shown in Figure1, whose x
y-plane represents the undeformed wing plane, and the x-axis points in the flow direction when the angle of attack is 0. Both the structural displacement and aerodynamic computation will be computed and described in this unified coordinate system. The leading segment of the vortex ring is placed on the panel’s quarter chord line; the aeroforce of the panels acts on the midpoint of the segment, which is represented by “○” in Figure 2, and the collocation point (represented by “×” in Figure 2) is at the centre of the three-quarter chord line, where the nonpenetration boundary condition will be implemented. Aerofoil warp can be represented by its middle camber surface, and some typical panel elements are shown in Figure 2.Figure 2
Arrangement of the vortex ring elements.The velocity induced at an arbitrary point by a typical vortex ring can be calculated by applying the Biot-Savart law [20] to the ring’s four segments, but attention should be paid to the rings located at the trailing edge where two semi-infinite trailing vortex lines that model the wake flow are shed along the free stream direction, as shown in Figure 3. Then, the induced velocities at all collocation points could be represented as(14)
V
x
=
W
C
x
Γ
V
y
=
W
C
y
Γ
V
z
=
W
C
z
Γ
,where V
x, V
y, and V
z are the induced velocity components along the x-axis, y-axis, and z-axis, respectively, and W
C
x, W
C
y, and W
C
z are their influence coefficient matrices related to the current large deformed configuration.Figure 3
Typical vortex elements.
(a)
Common vortex ring element
(b)
Vortex element at trailing edgeThe Neumann boundary condition is implemented to obtain the vortex distribution. The nonpenetration boundary condition is applied on the aerodynamic lattice of the current configuration, which implies that the method is also geometrically nonlinear and can account for large wing deformations. For the collocation point of theith panel element,(15)
V
∞
+
V
i
·
n
i
=
0
,where V
∞ is the velocity of the free stream; V
i is the induced velocity at the ith collocation point by all vortex elements, V
i
=
V
x
i
V
y
i
V
z
i
T; and n
i is the normal vector of the ith panel element reflecting the local spatial deformation.The circulation of the bound vortices is obtained by solving (15). Thus, the pressure distribution over the lifting surface can be obtained through the Kutta-Joukowski lift theorem [21]. The aerodynamic force that acts on the ith panel element is (16)
f
A
i
=
ρ
V
∞
×
Γ
F
i
,where Γ
F
i is the total vortex strength at the ith panel element’s quarter chord line, Γ
F
i
=
l
i
Γ
F
i
′ has different values depending on whether the panel is the leading edge panel or not, and l
i is the vector describing the magnitude and direction of the ith panel element’s quarter chord line. Consider Γ
F
i
′
=
Γ
i
′ when the panel is located on the leading edge of the wing and Γ
F
i
′
=
Γ
i
′
-
Γ
i
-
1
′ when it is not.The curved surface distributed aerodynamic force obtained by the NVLM combined with the nonlinear FEM with the follower force effect can describe the large deformed static equilibrium state. Some extreme conditions, such as stall and hover, require further consideration by some possible modification or more accurate calculation method, such as CFD. Here, only a normal situation that NVLM can manage is considered, which is sufficient for normal stability analysis for flexible wings.


### 2.3. Nonplanar Doublet-Lattice Method

The calculation of unsteady loads plays an important role across much of the design and development of an aircraft and has significant impacts on structural design, aerodynamic characteristics, weight, flight control system design, control surface design, and performance [22]. In stability analysis, spatial aerodynamic modelling was found to perform well in situations dominated by small amplitude dynamics around large quasistatic wing deflections [23]; thus, the concept of frequency and modal analysis could also be inherited from vibration theory. Then, the frequency domain aerodynamic methods could also be suitable for VFAs in stability analysis, in which the planar DLM is often used in traditional aeroelastic analysis. In this section, the DLM code is extended into nonplanar cases to account for the 3D unsteady loads of large-aspect-ratio wings with large deflections and can be successively applied in engineering practice.To meet the demand of nonplanar aerodynamic computations, mesh dividing should be determined on the deformed surface and updated along with the structure deflection, as shown in Figure4. In addition to the spatial lattices, local coordinates should be established to reflect the exact nonplanar configuration of the wing. The nonplanar effect not only is reflected geometrically but also should be contained in the kernel function K.Figure 4
Typical nonplanar lattice on a curved lifting surface.The relationship between the pressure and normal wash distribution in unsteady potential subsonic flow can be written as [24] (17)
w
x
,
y
,
z
U
=
1
8
π
×
∬
Δ
p
ξ
,
η
,
ζ
K
x
,
y
,
z
,
ξ
,
η
,
ζ
;
k
,
M
q
r
2
d
ξ
d
η
.Rodemich demonstrated that the kernel function can be written in the following form; a detailed derivation and explanation can be found in [25]:(18)
K
=
e
-
i
ω
x
0
/
U
K
1
T
1
+
K
2
T
2
r
1
2
,where T
1 and T
2 are geometric relations that represent the local deformations. K
1 and K
2 could be evaluated by integrals or via other methods [24].The critical problem of NDLM is the implementation of exact geometric boundary conditions. The local normal wash velocity, shown in Figure5, is computed via (17); then, the boundary condition is expressed as(19)
U
n
i
=
U
n
b
i
.Figure 5
Typical doublet-lattice line and its normal vector.Due to large deformations, the wing cannot be treated as oscillating about thex
y-plane; thus, the real boundary condition is related to local geometric nonlinearity. The more generalized boundary condition written based on modes is(20)
w
=
i
k
b
f
+
∂
f
∂
x
,where w
=
U
n
b
/
V
∞ is the nondimensional induced velocity in the normal direction by the spatial lattice, k is the reduced frequency, and f is the modal vibration shape function of the wing structure in a local normal direction and varies under different equilibrium states; f already includes the large static deformation, which makes it different from the traditional DLM despite its similar format. Because of the nonlinear geometric stiffness effects discussed before, modal shape and structure stiffness should be updated according to the different deformed configurations when computing unsteady aerodynamic forces.Substituting (17) and (20) into (19), the boundary condition can be written as an aerodynamic influence coefficient (AIC) in matrix form, where D is the instantaneously geometry-related AIC matrix:(21)
w
=
D
Δ
c
p
.Although the nonplanar DLM has been developed and used for many years and it could not consider an excessive number of unsteady effects, it is still capable of use in flexible wing stability analysis because the structure can be linearized, and the hypothesis of small amplitude dynamics around large quasistatic wing deflections can be made.


### 2.4. Surface Spline Interpolation

The integration of the flexible structure motion and the corresponding aerodynamic computation relies on the coupling between aerodynamics and structure. This is a critical problem in aeroelastic analysis, particularly for VFAs, whose large 3D elastic deformations make conventional 1D or 2D interpolations invalid. The generalized surface spline interpolation [26] is used in this study to couple the aerodynamics and structure. This method can clearly describe the interpolations in arbitrary space dimensions, which makes it applicable for VFA aerodynamic/structure coupling, whose structural configuration is typically considered to be embedded in a 3D space. A brief introduction is presented here.Consider a given vector setX
i
=
{
x
i
1
,
…
,
x
i
N
}
(
i
=
1,2
,
…
,
n
) in N-dimension space and the corresponding image vectors W
i
=
{
w
i
1
,
…
,
w
i
M
}
(
i
=
1,2
,
…
,
n
) inM-dimension space. For the kth component of W, it can be expressed as(22)
w
k
X
=
c
1
k
+
∑
p
=
1
N
c
p
+
1
k
x
p
+
∑
i
=
1
n
c
N
+
1
+
i
k
r
i
2
ln
⁡
r
i
2
+
ε
,where the symbols in (22) are the same as those in [27] and the coefficient ε could be given different value with k. This can be used as the interpolation between structural grids and aerodynamic grids. Considering N given structural grids with coordinates X
S and a corresponding deformation vector U
S, the relationship between the coordinates and the deformation of the grids could be written in matrix form:(23)
A
S
C
=
W
S
,and according to (22) we know that A
S and W
S are constant matrices related to the coordinates X
S and corresponding deformation U
S. C is the coefficient matrix of the surface spline fitting function. When A
S is nonsingular, C can be solved as(24)
C
=
A
S
-
1
W
S
.Then, the deformation vectorU
A of M aerodynamic grids with the coordinates of X
A can be interpolated as(25)
U
A
=
A
A
A
S
-
1
W
S
,where A
A is the constant matrix related to the given coordinates of aerodynamic grids X
A. Equation (25) can be transformed into(26)
U
A
=
G
U
S
,where G is the spline matrix for displacement interpolation between the aerodynamic grids and structural grids. According to this relationship, the aerodynamic grid and meshes can be updated automatically, ensuring that the real deformed aerosurface is considered in every aerodynamic computation.The other transformation also occurs here between the aerodynamic and structural force systems to obtain an equivalent force on the structure for the structural deflection analysis. When the aerodynamic forcesF
A and their equivalent structure forces F
S perform the same virtual work on their respective virtual deflections, the structural equivalence of the two force systems is guaranteed by(27)
δ
U
A
T
F
A
=
δ
U
S
T
F
S
,where δ
U
A and δ
U
S are the arbitrary virtual deflections satisfying (26). The nonplanar aerodynamic force computed is F
A
=
F
A
x
F
A
y
F
A
z. Therefore,(28)
F
S
=
G
T
F
A
.The bidirectional coupling between aerodynamics and the structure helps to complete the iterative process to determine the equilibrium states, which is fundamental in nonlinear analysis.


### 2.5. Flutter Analysis

Nonlinear structural stiffness [28] and nonplanar aerodynamics due to large deformations have dramatic effects on the flutter response, and, thus, the conventional linear approach of aeroelastic stability analysis is not applicable. To determine the stability properties of the structure, a modal analysis is conducted at the nonlinear equilibrium state, and the aeroelastic equation can be derived as(29)
M
q
¨
+
K
q
=
Q
.Equation (29) has the same form as that of the linear case, but each item is obtained from the nonlinear analysis introduced previously. Q
=
(
1
/
2
)
ρ
V
2
A
q describes the generalized unsteady aerodynamics vector, and A is the unsteady nonplanar aerodynamic influence coefficient complex matrix related to the reduced frequency k and the deformed configuration given by the NDLM. Using the p
-
k method, (29) can be rewritten as(30)
p
2
M
-
p
b
2
k
ρ
V
Q
I
+
K
-
1
2
ρ
V
2
Q
R
q
=
0
k
=
b
V
Im
⁡
p
,where Q
R and Q
I are the generalized aerodynamic stiffness matrix and aerodynamic damping matrix, respectively; p
=
ω
(
γ
±
i
)
=
2
π
f
(
γ
±
i
) is the complex eigenvalue; γ is the decaying rate; and g
=
2
γ is the structural damping ratio.In detailed computations with a certain Mach number, air densityρ, and flying velocity V, solving (30) yields a series of ω/f and g. Drawing the V
-
g and V
-
f curves of each mode, when g
=
0, the corresponding flutter velocity and frequency can be obtained. However, in the nonlinear analysis, the flutter velocity obtained here is only a predicted value that may change under different static loads and deformation cases. Thus, the nonlinear static aeroelastic analysis and stability analysis should be considered jointly, not separately, as in the linear analysis. The exact flutter velocity can be searched iteratively to use the geometrically nonlinear aeroelastic analysis flowchart below.The geometric nonlinear aeroelastic analysis flowchart for a very flexible wing, as shown in Figure6, can be decomposed into two parts: a nonlinear static aeroelastic iterative analysis and a larger iterative stability analysis that searches for an exact flutter speed. The static aeroelastic analysis is implemented using an automatic aerosurface update via surface spline interpolation to acquire the accurate static aeroelastic results. The geometrically nonlinear static structural analysis is executed by the FEM, which could reduce the limitations of the simple beam model and can be easily used in engineering practice. Typically, after four to five times iterative calculation, the nonlinear static aeroelastic analysis can get a convergent result. Flutter analysis is implemented by coupling the NDLM computation and the quasimodes calculation of the wing at its nonlinear equilibrium state obtained by the static aeroelastic analysis. Once those two iterations reach the same flight speed, the aeroelastic flutter analysis is accomplished. Otherwise, the analysis changes the flight status, such as the Mach number or the dynamic pressure, and repeats the calculation. Commonly, dichotomy is used to search for the exact flutter speed. If the predicted flutter speed is higher than the current computed flight speed, then the intermediate speed can be selected as the next computed flight speed. Most often, after four times dichotomy computation, the exact flutter speed can be found.Figure 6
Geometric nonlinear aeroelastic analysis flowchart.


## 2.1. Structure Geometric Nonlinearity

Large structural excursions induced by very flexible wings when undergoing aerodynamic loads in the air prevent the use of linear methods based on the small displacement assumption and call for geometrically nonlinear structural analysis. Geometric nonlinearity are based on the kinematic description of the body, and the strain on the wing should be defined in terms of local displacement of the wing for dynamic motion. These result in the nonlinear geometric equations including the quadric term of the displacement differential and require the nonlinear force equilibrium equation established on the deformed state of the structure. Structural geometrically nonlinear problems are often solved by the maturely developed nonlinear incremental finite element method [16] and two formulae called the total Lagrange formulation (TLF) and the updated Lagrange formulation (ULF) [17], which are well known. The ULF is used in this study, and the primary equations are presented briefly below. The core method of structural analysis has already imbedded in commercial software, so it is convenient to be used in engineering analysis.The relationship between the nonlinear Lagrange/Green strain and displacement is(1)
ε
i
j
t
=
1
2
u
t
i
,
j
+
u
t
j
,
i
+
u
k
,
i
t
u
k
,
j
t
,where u
i
,
j
t is the partial derivative of displacement component u
i to the coordinate x
j at time t.Despite a large elastic deformation, the material is still within the elastic limitation for a small strain, so the conjugate Kirchhoff stress tensorS
j
i at timet satisfies(2)
S
j
i
t
n
j
t
d
s
t
=
x
i
,
j
t
d
T
j
t
,where n
j
t is the direction cosine of a small area element d
s at time t and d
T
j is the corresponding surface force in which the follower force effect is considered. The linear elastic constitutive relation is given as follows:(3)
S
i
j
t
=
D
i
j
k
l
t
ε
k
l
t
,where D
i
j
k
l is the elastic tensor, which has a different form for isotropic or anisotropic material.The finite element method (FEM) based on energy principles is an effective approach to solve structural problems. For geometric nonlinear problems, considering the follower force effect, the incremental FEM is used. The strainε
i
j can be decomposed into a linear part e
i
j and a nonlinear part η
i
j of the current displacement:(4)
ε
i
j
t
=
e
i
j
t
+
η
i
j
t
.The stress can also be decomposed in increments, where S
i
j
t represents the equilibrium stress at time t and τ
i
j
t represents the incremental stress to be calculated at each time step:(5)
S
i
j
t
+
Δ
t
=
τ
i
j
t
+
S
i
j
t
.The integral equation is established by linearization in each incremental step:(6)
∫
V
t
D
i
j
k
l
t
e
k
l
δ
e
i
j
t
t
d
V
t
+
∫
V
t
τ
i
j
t
δ
η
i
j
t
d
V
t
=
Q
t
+
Δ
t
-
∫
V
t
τ
i
j
δ
e
i
j
t
t
d
V
t
,where Q
t
+
Δ
t is the incremental outer force including the aerodynamic force, engine thrust, and gravity, at the new time step. Considering a number of shape functions, the relationship between strain and deformation is presented as(7)
e
t
=
B
L
t
u
,
η
t
=
B
N
L
t
u
.Substituting these shape functions into (6) leads to the element governing equation [18]: (8)
K
L
t
+
K
N
L
t
u
=
Q
-
t
+
Δ
t
-
F
A
t
,where Q
-
t
+
Δ
t is the incremental outer force including the aerodynamic force, engine thrust, and gravity at the new time step. The stiffness matrix in (8) could be decomposed into a linear part K
L
t and nonlinear part K
N
L
t. The linear part is only related to the structure itself, whereas the nonlinear part is related to the deflected configuration and strain quality, which should be updated in each computation step.The corresponding dynamic equation can be expressed as(9)
M
u
¨
t
+
Δ
t
+
K
L
t
+
K
N
L
t
u
=
Q
-
t
+
Δ
t
-
F
A
t
,where u
¨
t
+
Δ
t is the structural displacement acceleration vector at new time step. The assumption of a small amplitude vibration around the static equilibrium state is suitable for many dynamic problems, including the flexible wing structural dynamic stability:(10)
u
=
u
-
+
x
,where u
- is the large deflecting equilibrium deformation from (8) and x is a small vibration deformation. According to (9) and the static equilibrium condition, the linearized structural quasimode can be obtained by generalized diagonalization, and the vibration equation of the system under steady forces reduces to(11)
M
T
x
¨
+
K
T
x
=
0
,where M
T is the inertial matrix of the structure at the nonlinear static equilibrium configuration and K
T is the corresponding stiffness matrix. Both of these parameters are nonlinear functions of u
- and vary under different equilibrium states, which is a key feature of geometric nonlinear structures. The mode shapes and frequencies can be deduced from (11).Introducing the harmonic oscillating assumptionx
=
ϕ
e
i
ω
t, the vibration equation can be written as (12)
K
T
-
ω
2
M
T
ϕ
=
0
,where ω is the vibration circular frequency and ϕ is the vibration mode matrix. If (12) has all-nonzero solutions, that demands (13)
det
⁡
K
T
-
ω
2
M
T
=
0
.That is a generalized eigenvalue problem about K
T and M
T. Solving (13), the vibration circular frequency ω
i should be obtained. Substituting ω
i into (12), structural eigenvector ϕ
i (structural mode shape) can be obtained.


## 2.2. Nonplanar Vortex Lattice Method

The modelling of flexible aircraft with significant structural deformations requires the incorporation of structural dynamics and aerodynamics in a unified framework. The aerodynamic loads in (9) should be computed on a deformed configuration. This section summarizes the primary characteristics of the NVLM and its application to a curved aerosurface. The NVLM is derived from three-dimensional potential flow theory and is suitable for most normal situations that VAFs may encounter. Because the simple programming effort is required, NVLM is easy to combine with structural dynamic computations to obtain the response results for aeroelastic structures. Additionally, the exact boundary condition is satisfied on the real wing surface in the NVLM, which can have camber and various platform shapes. Thus, the NVLM is convenient for use with very flexible wings, whose aerodynamic surfaces are subjected to large structural deformations [19]. The NVLM is implemented using vortex ring quadrilateral elements to discretize the curved lifting surface along with the wing’s deformation, as shown in Figure 1.Figure 1
Nonplanar vortex lattice model for a thin lifting surface.A Cartesian coordinate system is shown in Figure1, whose x
y-plane represents the undeformed wing plane, and the x-axis points in the flow direction when the angle of attack is 0. Both the structural displacement and aerodynamic computation will be computed and described in this unified coordinate system. The leading segment of the vortex ring is placed on the panel’s quarter chord line; the aeroforce of the panels acts on the midpoint of the segment, which is represented by “○” in Figure 2, and the collocation point (represented by “×” in Figure 2) is at the centre of the three-quarter chord line, where the nonpenetration boundary condition will be implemented. Aerofoil warp can be represented by its middle camber surface, and some typical panel elements are shown in Figure 2.Figure 2
Arrangement of the vortex ring elements.The velocity induced at an arbitrary point by a typical vortex ring can be calculated by applying the Biot-Savart law [20] to the ring’s four segments, but attention should be paid to the rings located at the trailing edge where two semi-infinite trailing vortex lines that model the wake flow are shed along the free stream direction, as shown in Figure 3. Then, the induced velocities at all collocation points could be represented as(14)
V
x
=
W
C
x
Γ
V
y
=
W
C
y
Γ
V
z
=
W
C
z
Γ
,where V
x, V
y, and V
z are the induced velocity components along the x-axis, y-axis, and z-axis, respectively, and W
C
x, W
C
y, and W
C
z are their influence coefficient matrices related to the current large deformed configuration.Figure 3
Typical vortex elements.
(a)
Common vortex ring element
(b)
Vortex element at trailing edgeThe Neumann boundary condition is implemented to obtain the vortex distribution. The nonpenetration boundary condition is applied on the aerodynamic lattice of the current configuration, which implies that the method is also geometrically nonlinear and can account for large wing deformations. For the collocation point of theith panel element,(15)
V
∞
+
V
i
·
n
i
=
0
,where V
∞ is the velocity of the free stream; V
i is the induced velocity at the ith collocation point by all vortex elements, V
i
=
V
x
i
V
y
i
V
z
i
T; and n
i is the normal vector of the ith panel element reflecting the local spatial deformation.The circulation of the bound vortices is obtained by solving (15). Thus, the pressure distribution over the lifting surface can be obtained through the Kutta-Joukowski lift theorem [21]. The aerodynamic force that acts on the ith panel element is (16)
f
A
i
=
ρ
V
∞
×
Γ
F
i
,where Γ
F
i is the total vortex strength at the ith panel element’s quarter chord line, Γ
F
i
=
l
i
Γ
F
i
′ has different values depending on whether the panel is the leading edge panel or not, and l
i is the vector describing the magnitude and direction of the ith panel element’s quarter chord line. Consider Γ
F
i
′
=
Γ
i
′ when the panel is located on the leading edge of the wing and Γ
F
i
′
=
Γ
i
′
-
Γ
i
-
1
′ when it is not.The curved surface distributed aerodynamic force obtained by the NVLM combined with the nonlinear FEM with the follower force effect can describe the large deformed static equilibrium state. Some extreme conditions, such as stall and hover, require further consideration by some possible modification or more accurate calculation method, such as CFD. Here, only a normal situation that NVLM can manage is considered, which is sufficient for normal stability analysis for flexible wings.


## 2.3. Nonplanar Doublet-Lattice Method

The calculation of unsteady loads plays an important role across much of the design and development of an aircraft and has significant impacts on structural design, aerodynamic characteristics, weight, flight control system design, control surface design, and performance [22]. In stability analysis, spatial aerodynamic modelling was found to perform well in situations dominated by small amplitude dynamics around large quasistatic wing deflections [23]; thus, the concept of frequency and modal analysis could also be inherited from vibration theory. Then, the frequency domain aerodynamic methods could also be suitable for VFAs in stability analysis, in which the planar DLM is often used in traditional aeroelastic analysis. In this section, the DLM code is extended into nonplanar cases to account for the 3D unsteady loads of large-aspect-ratio wings with large deflections and can be successively applied in engineering practice.To meet the demand of nonplanar aerodynamic computations, mesh dividing should be determined on the deformed surface and updated along with the structure deflection, as shown in Figure4. In addition to the spatial lattices, local coordinates should be established to reflect the exact nonplanar configuration of the wing. The nonplanar effect not only is reflected geometrically but also should be contained in the kernel function K.Figure 4
Typical nonplanar lattice on a curved lifting surface.The relationship between the pressure and normal wash distribution in unsteady potential subsonic flow can be written as [24] (17)
w
x
,
y
,
z
U
=
1
8
π
×
∬
Δ
p
ξ
,
η
,
ζ
K
x
,
y
,
z
,
ξ
,
η
,
ζ
;
k
,
M
q
r
2
d
ξ
d
η
.Rodemich demonstrated that the kernel function can be written in the following form; a detailed derivation and explanation can be found in [25]:(18)
K
=
e
-
i
ω
x
0
/
U
K
1
T
1
+
K
2
T
2
r
1
2
,where T
1 and T
2 are geometric relations that represent the local deformations. K
1 and K
2 could be evaluated by integrals or via other methods [24].The critical problem of NDLM is the implementation of exact geometric boundary conditions. The local normal wash velocity, shown in Figure5, is computed via (17); then, the boundary condition is expressed as(19)
U
n
i
=
U
n
b
i
.Figure 5
Typical doublet-lattice line and its normal vector.Due to large deformations, the wing cannot be treated as oscillating about thex
y-plane; thus, the real boundary condition is related to local geometric nonlinearity. The more generalized boundary condition written based on modes is(20)
w
=
i
k
b
f
+
∂
f
∂
x
,where w
=
U
n
b
/
V
∞ is the nondimensional induced velocity in the normal direction by the spatial lattice, k is the reduced frequency, and f is the modal vibration shape function of the wing structure in a local normal direction and varies under different equilibrium states; f already includes the large static deformation, which makes it different from the traditional DLM despite its similar format. Because of the nonlinear geometric stiffness effects discussed before, modal shape and structure stiffness should be updated according to the different deformed configurations when computing unsteady aerodynamic forces.Substituting (17) and (20) into (19), the boundary condition can be written as an aerodynamic influence coefficient (AIC) in matrix form, where D is the instantaneously geometry-related AIC matrix:(21)
w
=
D
Δ
c
p
.Although the nonplanar DLM has been developed and used for many years and it could not consider an excessive number of unsteady effects, it is still capable of use in flexible wing stability analysis because the structure can be linearized, and the hypothesis of small amplitude dynamics around large quasistatic wing deflections can be made.


## 2.4. Surface Spline Interpolation

The integration of the flexible structure motion and the corresponding aerodynamic computation relies on the coupling between aerodynamics and structure. This is a critical problem in aeroelastic analysis, particularly for VFAs, whose large 3D elastic deformations make conventional 1D or 2D interpolations invalid. The generalized surface spline interpolation [26] is used in this study to couple the aerodynamics and structure. This method can clearly describe the interpolations in arbitrary space dimensions, which makes it applicable for VFA aerodynamic/structure coupling, whose structural configuration is typically considered to be embedded in a 3D space. A brief introduction is presented here.Consider a given vector setX
i
=
{
x
i
1
,
…
,
x
i
N
}
(
i
=
1,2
,
…
,
n
) in N-dimension space and the corresponding image vectors W
i
=
{
w
i
1
,
…
,
w
i
M
}
(
i
=
1,2
,
…
,
n
) inM-dimension space. For the kth component of W, it can be expressed as(22)
w
k
X
=
c
1
k
+
∑
p
=
1
N
c
p
+
1
k
x
p
+
∑
i
=
1
n
c
N
+
1
+
i
k
r
i
2
ln
⁡
r
i
2
+
ε
,where the symbols in (22) are the same as those in [27] and the coefficient ε could be given different value with k. This can be used as the interpolation between structural grids and aerodynamic grids. Considering N given structural grids with coordinates X
S and a corresponding deformation vector U
S, the relationship between the coordinates and the deformation of the grids could be written in matrix form:(23)
A
S
C
=
W
S
,and according to (22) we know that A
S and W
S are constant matrices related to the coordinates X
S and corresponding deformation U
S. C is the coefficient matrix of the surface spline fitting function. When A
S is nonsingular, C can be solved as(24)
C
=
A
S
-
1
W
S
.Then, the deformation vectorU
A of M aerodynamic grids with the coordinates of X
A can be interpolated as(25)
U
A
=
A
A
A
S
-
1
W
S
,where A
A is the constant matrix related to the given coordinates of aerodynamic grids X
A. Equation (25) can be transformed into(26)
U
A
=
G
U
S
,where G is the spline matrix for displacement interpolation between the aerodynamic grids and structural grids. According to this relationship, the aerodynamic grid and meshes can be updated automatically, ensuring that the real deformed aerosurface is considered in every aerodynamic computation.The other transformation also occurs here between the aerodynamic and structural force systems to obtain an equivalent force on the structure for the structural deflection analysis. When the aerodynamic forcesF
A and their equivalent structure forces F
S perform the same virtual work on their respective virtual deflections, the structural equivalence of the two force systems is guaranteed by(27)
δ
U
A
T
F
A
=
δ
U
S
T
F
S
,where δ
U
A and δ
U
S are the arbitrary virtual deflections satisfying (26). The nonplanar aerodynamic force computed is F
A
=
F
A
x
F
A
y
F
A
z. Therefore,(28)
F
S
=
G
T
F
A
.The bidirectional coupling between aerodynamics and the structure helps to complete the iterative process to determine the equilibrium states, which is fundamental in nonlinear analysis.


## 2.5. Flutter Analysis

Nonlinear structural stiffness [28] and nonplanar aerodynamics due to large deformations have dramatic effects on the flutter response, and, thus, the conventional linear approach of aeroelastic stability analysis is not applicable. To determine the stability properties of the structure, a modal analysis is conducted at the nonlinear equilibrium state, and the aeroelastic equation can be derived as(29)
M
q
¨
+
K
q
=
Q
.Equation (29) has the same form as that of the linear case, but each item is obtained from the nonlinear analysis introduced previously. Q
=
(
1
/
2
)
ρ
V
2
A
q describes the generalized unsteady aerodynamics vector, and A is the unsteady nonplanar aerodynamic influence coefficient complex matrix related to the reduced frequency k and the deformed configuration given by the NDLM. Using the p
-
k method, (29) can be rewritten as(30)
p
2
M
-
p
b
2
k
ρ
V
Q
I
+
K
-
1
2
ρ
V
2
Q
R
q
=
0
k
=
b
V
Im
⁡
p
,where Q
R and Q
I are the generalized aerodynamic stiffness matrix and aerodynamic damping matrix, respectively; p
=
ω
(
γ
±
i
)
=
2
π
f
(
γ
±
i
) is the complex eigenvalue; γ is the decaying rate; and g
=
2
γ is the structural damping ratio.In detailed computations with a certain Mach number, air densityρ, and flying velocity V, solving (30) yields a series of ω/f and g. Drawing the V
-
g and V
-
f curves of each mode, when g
=
0, the corresponding flutter velocity and frequency can be obtained. However, in the nonlinear analysis, the flutter velocity obtained here is only a predicted value that may change under different static loads and deformation cases. Thus, the nonlinear static aeroelastic analysis and stability analysis should be considered jointly, not separately, as in the linear analysis. The exact flutter velocity can be searched iteratively to use the geometrically nonlinear aeroelastic analysis flowchart below.The geometric nonlinear aeroelastic analysis flowchart for a very flexible wing, as shown in Figure6, can be decomposed into two parts: a nonlinear static aeroelastic iterative analysis and a larger iterative stability analysis that searches for an exact flutter speed. The static aeroelastic analysis is implemented using an automatic aerosurface update via surface spline interpolation to acquire the accurate static aeroelastic results. The geometrically nonlinear static structural analysis is executed by the FEM, which could reduce the limitations of the simple beam model and can be easily used in engineering practice. Typically, after four to five times iterative calculation, the nonlinear static aeroelastic analysis can get a convergent result. Flutter analysis is implemented by coupling the NDLM computation and the quasimodes calculation of the wing at its nonlinear equilibrium state obtained by the static aeroelastic analysis. Once those two iterations reach the same flight speed, the aeroelastic flutter analysis is accomplished. Otherwise, the analysis changes the flight status, such as the Mach number or the dynamic pressure, and repeats the calculation. Commonly, dichotomy is used to search for the exact flutter speed. If the predicted flutter speed is higher than the current computed flight speed, then the intermediate speed can be selected as the next computed flight speed. Most often, after four times dichotomy computation, the exact flutter speed can be found.Figure 6
Geometric nonlinear aeroelastic analysis flowchart.


## 3. Analysis Model

### 3.1. Wind Tunnel Test Model

The geometrically nonlinear characteristics of VFAs have attracted people’s attention for many years; however, there are still some basic principles that require a more in-depth description; this situation indicates the importance of wind tunnel tests. Thus, a wind tunnel test model was designed to validate the theoretical analysis established in this paper and to provide insight into VFA aeroelastic characteristics. Although the proposed theoretical analysis framework can be widely used for models of different complexities, a simple wing model here was used to perform the wind tunnel test and illustrate the nonlinearity of the aeroelasticity of very flexible wings. The tested semiwing, whose deformation is quite large during testing, is fixed at the root (right side of the64
×
8 lattices, as shown in Figure 8). The FEM model as shown in (Figure 7) has been modified by ground vibration tests, in which different weights were attached to test the static deformation and model frequency under horizontal suspended cases. The wing was also fixed upright to minimize the gravity effect and obtain its approximate natural frequency. The natural vibration information obtained by FEM after model calibration is shown in Table 1. Moreover, the GVT test results and nonlinear analysis results were also presented in the table. Although the wing had been fixed vertically, as shown in Figure 9, the nonlinear characteristics could still be observed. Since the gravity acted along spanwise, the stiffening effect can be noticed from the commonly higher frequency in the test results compared with linear natural frequency. The structural linearized frequency results shown in Table 1 were obtained via nonlinear static analysis and quasimode analysis under nonlinear equilibrium as introduced in Section 2.1. The acceptable small errors between test results and linearized analysis results indicate our structural modelling is credible and nonlinear analyses 64
×
8 were basically accurate and available to be used in the subsequent analysis.Table 1
Structure mode frequency.
Mode number
Mode name
Natural frequency (Hz)
Test frequency(Hz)
Linearized frequency (Hz)
The error between test and nonlinear analysis (%)
1
Vertical 1st bend
1.78
2.00
1.95
-
2.50
2
Horizontal 1st bend
10.60
10.56
10.63
-
0.66
3
1st twist
13.75
13.93
13.74
-
1.36
4
Vertical 2nd bend
17.36
18.10
17.60
-
2.76
5
Vertical 3rd bend
54.00
55.83
54.26
-
2.81Figure 7
Structure of the FEM model.Figure 8
Aerodynamic model.Figure 9
Ground vibration test (wing was fixed vertically).


## 3.1. Wind Tunnel Test Model

The geometrically nonlinear characteristics of VFAs have attracted people’s attention for many years; however, there are still some basic principles that require a more in-depth description; this situation indicates the importance of wind tunnel tests. Thus, a wind tunnel test model was designed to validate the theoretical analysis established in this paper and to provide insight into VFA aeroelastic characteristics. Although the proposed theoretical analysis framework can be widely used for models of different complexities, a simple wing model here was used to perform the wind tunnel test and illustrate the nonlinearity of the aeroelasticity of very flexible wings. The tested semiwing, whose deformation is quite large during testing, is fixed at the root (right side of the64
×
8 lattices, as shown in Figure 8). The FEM model as shown in (Figure 7) has been modified by ground vibration tests, in which different weights were attached to test the static deformation and model frequency under horizontal suspended cases. The wing was also fixed upright to minimize the gravity effect and obtain its approximate natural frequency. The natural vibration information obtained by FEM after model calibration is shown in Table 1. Moreover, the GVT test results and nonlinear analysis results were also presented in the table. Although the wing had been fixed vertically, as shown in Figure 9, the nonlinear characteristics could still be observed. Since the gravity acted along spanwise, the stiffening effect can be noticed from the commonly higher frequency in the test results compared with linear natural frequency. The structural linearized frequency results shown in Table 1 were obtained via nonlinear static analysis and quasimode analysis under nonlinear equilibrium as introduced in Section 2.1. The acceptable small errors between test results and linearized analysis results indicate our structural modelling is credible and nonlinear analyses 64
×
8 were basically accurate and available to be used in the subsequent analysis.Table 1
Structure mode frequency.
Mode number
Mode name
Natural frequency (Hz)
Test frequency(Hz)
Linearized frequency (Hz)
The error between test and nonlinear analysis (%)
1
Vertical 1st bend
1.78
2.00
1.95
-
2.50
2
Horizontal 1st bend
10.60
10.56
10.63
-
0.66
3
1st twist
13.75
13.93
13.74
-
1.36
4
Vertical 2nd bend
17.36
18.10
17.60
-
2.76
5
Vertical 3rd bend
54.00
55.83
54.26
-
2.81Figure 7
Structure of the FEM model.Figure 8
Aerodynamic model.Figure 9
Ground vibration test (wing was fixed vertically).


## 4. Aeroelastic Stability Analysis

### 4.1. Static Aeroelastic Analysis

For comparison purposes, the static aeroelasticity of the wing model is calculated by the traditional linear method (represented by the label “linear” in Figure10) and the nonlinear method developed in this paper (represented by the label “nonlinear” in Figure 10). In the linear method, the assumption of infinitesimal deformation is adopted, and the aerodynamic loads are computed by the steady planar vortex lattice method, which do not consider the wing’s deflection. Although the conventional planar vortex lattice method is not accurate and even erroneous in this case, it has been adopted for more than thirty years and still has been widely used in engineering analysis. However, that is our paper’s aim to show the inaccuracy of planar aerodynamic computation and establish an easily acceptable nonplanar aerodynamic computation to deal with the nonplanar lifting surfaces. Figure 10 shows the vertical displacements of the wing tip that are obtained by these two methods under a 0.5° angle of attack at the root; the results of the wind tunnel test are also presented.Figure 10
Vertical displacement at the wing tip versus airspeed.The test results in Figure10 indicate that vertical displacement of the wing tip increases with increasing wind speed, as does the curve slope of the test result until the velocity reaches approximately 33.5 m/s. When the velocity is further increased, the rate of wing tip displacement increase is slowed down and converges to a limit value due to geometric stiffness effect. The results obtained by the nonlinear method are consistent with the test results despite some small differences that may be caused by systematic measurement errors. However, the wing tip deflection obtained by the linear method increases with wind speed and reaches approximately 456.9 mm at 34.0 m/s, which is far beyond reality and is no longer significant.According to the test results, the vertical displacement at the wing tip is approximately 37.5% of the span under a 0.5° angle of attack when the wind speed is 34.0 m/s. The initial and final aerosurfaces in the nonlinear analysis are presented in Figure11. It is a typical nonlinear case of the flexible wing with a large deformation, and the nonplanar effect of the lifting surface is quite significant. Detailed displacements in three directions along the span in this case are presented in Figure 13. Figure 12 shows the linear planar case and nonlinear nonplanar case in static analysis. In the traditional linear analysis, the aerosurface remains planar, and the aerodynamic force F
A acts only in thez direction; considering the linear relation between displacement and force, deflections in thex-axisandy-axis are small, shown in Figure 13. However, in the nonlinear analysis, the aerosurface is automatically updated with structural deflections, and the follower force effect is included so that the displacements in thex-axis andy-axis can be considered. As shown in Figure 12, the aerodynamic force F
A acts vertically to the aerosurface. Thus, besides the lift, a large component side force is significant in nonlinear cases, which may induce large deflections in the y-axis (v in Figure 12) and reduce the effective lift. There is induced drag in x-axis in both linear case and nonlinear case, but it is only quite small part of force in x-axis. The deflection in x-axis under nonlinear case is mainly caused by the follower force effect other than the induced drag force. Once there is a little torsion, there is a component force of F
A in x-axis that may cause deflections in x-axis. In linear case, the aerosurface keeps planar and the aero force only acts in z
-axis, and the structural deflection in x-axis is rather small.Figure 11
Initial and final aerodynamic model of the flexible wing.Figure 12
Aerodynamic force in linear and nonlinear case.Figure 13
Displacement of the wing spar spanwise (0.5° AOA, 34 m/s).
(a)
T
x
(b)
T
y
(c)
T
zAdditionally, in the linear analysis, the effective aerodynamic load reduction due to the shortening of aerosurface and structural stiffening effect due to elastic deflections cannot be considered; as a result, the displacements along thez-axis (Figure 13) and the torsion angle (Figure 14) of the nonlinear analysis are smaller and more reliable than the linear results. Because the wings are fixed at their roots, the torsion angle increases from zero along the span. The tip torsion angle is still within the limitation of potential flow, so the NVLM still can be applied.Figure 14
Torsional angle of the wing spar spanwise (0.5° AOA, 34 m/s).The aerodynamic loads along spanwise, including the lift, drag, and side force, are shown in Figure15. These loads are the directional component forces of the total aerodynamic load obtained from the NVLM. Since the aerosurface keeps planar, the effective lifting surface area is much bigger than the nonlinear case. Additionally, the torsion angle is bigger than nonlinear case, so the lift in linear case is bigger than that in nonlinear case. Because of the follower force effect and the updated aerosurface, aerodynamic force may have significant component force in x-axis, so the drag in nonlinear case is bigger than that in linear case. Since the aerosurface keeps planar, there is no side force in linear case. In nonlinear case, the aerodynamic force acts vertically to the large deformed aerosurface so there is quite large component force acting along the span as the side force. The side force is comparable to the lift and contributes significantly to the root bend moment; thus, it cannot be ignored in the nonlinear analysis. The elastic torsional angle increases the effective angle of attack; thus, the maximum aerodynamic loads along the span move to the outer wing.Figure 15
Drag, side force, and lift along spanwise.
(a)
Lift
(b)
Drag
(c)
Side force


### 4.2. Flutter Analysis

The structural dynamics are different under various deflections and force cases due to geometrically nonlinear effects. Figure16 shows the structural dynamic frequency under different static load cases. The reduction of the horizontal bend frequency and the increase in the twist frequency are distinct. An iterative process is used to obtain the exact flutter boundary under certain flight conditions, and the effect of gravity is also considered. The results of different cases are compared with the experiment results to illustrate the nonlinear flutter characteristics and validate the theoretical analysis.Figure 16
Structural frequency under nonlinear static deflections.Figure17 shows the V
-
F and V
-
g curves of the traditional linear results. A classical bend/twist coupling flutter pattern results in a flutter speed at approximately 33.5 m/s, and the horizontal bend mode does not contribute to the flutter; its frequency remains unchanged. The 1st vertical mode frequency reduces to zero when the speed is over 34 m/s but the double calculation and the wind tunnel test indicate that there is no divergence occurring around 34 m/s.Figure 17
V
-
F and V
-
g curves of the linear flutter analysis.The nonlinear stability analysis results are presented below. AllV
-
g curves are found to be significantly different from the linear results with an unstable mode switching to a horizontal bend in the nonlinear analysis. This is caused by the geometrically nonlinear effect, which makes the horizontal bending mode contain large torsional components and becomes unstable when the structure experiences a large deformation, as shown in Figure 18. It is a remarkable characteristic of nonlinear stability analysis of VFAs that the horizontal bend mode plays an important role in flutter characteristics.Figure 18
Horizontal bending mode in a nonlinear equilibrium state.Figure19 shows the stability analysis results under the equilibrium states of 20 and 25 m/s. The corresponding predicted flutter speeds are approximately 31 and 33 m/s. Typically, the flutter speed decreases with increasing airspeed in the nonlinear analysis, but the test wing induces a negative deformation when the airspeed is below 28 m/s. Thus, the predicted flutter speed of 25 m/s is larger than that of 20 m/s. Because of the disagreement between the static equilibrium flight speed and predicted flutter speed, only the data near the static state marked red in Figure 19 are reliable. An iterative process is performed to narrow the disagreement and search for the exact flutter speed. Figure 20 shows the V
-
g and V
-
F curves from the dynamic aeroelastic analysis under the static equilibrium flight speed of 31.5 m/s. The predicted critical speed is 31.5 m/s, which is the theoretical exact nonlinear flutter speed of the analysis model under 0.5° angle of attack. Again, the data in Figure 20 is only accurate near 31.5 m/s, and if completely accurate V
-
g and V
-
F curves were desired, the divisional accurate data should be pieced together, as is in Figure 21, which requires more data to make the curve smoother and more reliable.Figure 19
Predicted flutter speed under two flight cases.Figure 20
V
-
g and V
-
F curves under 31.5 m/s.Figure 21
Accurate data pieced together roughly.The calculated cases above are under a 0.5° angle of attack. Due to geometric nonlinearity, cases under different angles of attack lead to different flutter speeds. The iterative processes searching for exact flutter boundary under different angles of attack have been implemented and the final results are shown in Figure22. The flutter results are summarized in Table 2. The aerodynamic loads and structural deflections increase as the angle of attack increases. Due to geometric nonlinearity, large deflections and load conditions easily cause the horizontal bending to become unstable. Thus, both the test and nonlinear analysis results indicate that the flutter speed decreases with an increasing angle of attack, whereas no changes are noted in the linear analysis under different cases. The first three cases present nearly identical deviations between the nonlinear analysis results and test results (in Figure 23). This effect may be caused by the conservative aerodynamic computation that can be validated in the deflection curve in Figure 10, where the analysed displacement is larger than that found in the test when the displacement is positive. Thus, it is understandable that the analysed flutter speed is smaller than the test results. The deviation narrows in the last case because the test result changes to a stall flutter depending on the flutter phenomena observed in the test; the stall model is not considered in theoretical aerodynamic computation, but the decline in flutter speed can still be reflected. In our test cases, the linear analysis results sometimes are less conservative than nonlinear results. This is not a referable conclusion. For most cases, the nonlinear analysis results are more accurate and reliable than linear results, especially when the nonlinear flutter speed is lower than the linear results.Table 2
Flutter speeds under different angles of attack.
Angle of attack
0°
0.5°
1.0°
1.5°
Linear
Computation
Experiment
Computation
Experiment
Computation
Experiment
Computation
Experiment
Computation
Flutter speed (m/s)
33.0
35.7
31.5
34.4
28.0
32.2
26.5
27.0
33.5
Flutter frequency (Hz)
7.60
6.94
7.82
6.88
7.34
7.00
7.53
13.38
7.29
Unstable mode
First horizontal bend in plane
First twistFigure 22
Nonlinear analysed flutter speed under different angles of attack.Figure 23
Flutter limit speed at different angles of attack.


## 4.1. Static Aeroelastic Analysis

For comparison purposes, the static aeroelasticity of the wing model is calculated by the traditional linear method (represented by the label “linear” in Figure10) and the nonlinear method developed in this paper (represented by the label “nonlinear” in Figure 10). In the linear method, the assumption of infinitesimal deformation is adopted, and the aerodynamic loads are computed by the steady planar vortex lattice method, which do not consider the wing’s deflection. Although the conventional planar vortex lattice method is not accurate and even erroneous in this case, it has been adopted for more than thirty years and still has been widely used in engineering analysis. However, that is our paper’s aim to show the inaccuracy of planar aerodynamic computation and establish an easily acceptable nonplanar aerodynamic computation to deal with the nonplanar lifting surfaces. Figure 10 shows the vertical displacements of the wing tip that are obtained by these two methods under a 0.5° angle of attack at the root; the results of the wind tunnel test are also presented.Figure 10
Vertical displacement at the wing tip versus airspeed.The test results in Figure10 indicate that vertical displacement of the wing tip increases with increasing wind speed, as does the curve slope of the test result until the velocity reaches approximately 33.5 m/s. When the velocity is further increased, the rate of wing tip displacement increase is slowed down and converges to a limit value due to geometric stiffness effect. The results obtained by the nonlinear method are consistent with the test results despite some small differences that may be caused by systematic measurement errors. However, the wing tip deflection obtained by the linear method increases with wind speed and reaches approximately 456.9 mm at 34.0 m/s, which is far beyond reality and is no longer significant.According to the test results, the vertical displacement at the wing tip is approximately 37.5% of the span under a 0.5° angle of attack when the wind speed is 34.0 m/s. The initial and final aerosurfaces in the nonlinear analysis are presented in Figure11. It is a typical nonlinear case of the flexible wing with a large deformation, and the nonplanar effect of the lifting surface is quite significant. Detailed displacements in three directions along the span in this case are presented in Figure 13. Figure 12 shows the linear planar case and nonlinear nonplanar case in static analysis. In the traditional linear analysis, the aerosurface remains planar, and the aerodynamic force F
A acts only in thez direction; considering the linear relation between displacement and force, deflections in thex-axisandy-axis are small, shown in Figure 13. However, in the nonlinear analysis, the aerosurface is automatically updated with structural deflections, and the follower force effect is included so that the displacements in thex-axis andy-axis can be considered. As shown in Figure 12, the aerodynamic force F
A acts vertically to the aerosurface. Thus, besides the lift, a large component side force is significant in nonlinear cases, which may induce large deflections in the y-axis (v in Figure 12) and reduce the effective lift. There is induced drag in x-axis in both linear case and nonlinear case, but it is only quite small part of force in x-axis. The deflection in x-axis under nonlinear case is mainly caused by the follower force effect other than the induced drag force. Once there is a little torsion, there is a component force of F
A in x-axis that may cause deflections in x-axis. In linear case, the aerosurface keeps planar and the aero force only acts in z
-axis, and the structural deflection in x-axis is rather small.Figure 11
Initial and final aerodynamic model of the flexible wing.Figure 12
Aerodynamic force in linear and nonlinear case.Figure 13
Displacement of the wing spar spanwise (0.5° AOA, 34 m/s).
(a)
T
x
(b)
T
y
(c)
T
zAdditionally, in the linear analysis, the effective aerodynamic load reduction due to the shortening of aerosurface and structural stiffening effect due to elastic deflections cannot be considered; as a result, the displacements along thez-axis (Figure 13) and the torsion angle (Figure 14) of the nonlinear analysis are smaller and more reliable than the linear results. Because the wings are fixed at their roots, the torsion angle increases from zero along the span. The tip torsion angle is still within the limitation of potential flow, so the NVLM still can be applied.Figure 14
Torsional angle of the wing spar spanwise (0.5° AOA, 34 m/s).The aerodynamic loads along spanwise, including the lift, drag, and side force, are shown in Figure15. These loads are the directional component forces of the total aerodynamic load obtained from the NVLM. Since the aerosurface keeps planar, the effective lifting surface area is much bigger than the nonlinear case. Additionally, the torsion angle is bigger than nonlinear case, so the lift in linear case is bigger than that in nonlinear case. Because of the follower force effect and the updated aerosurface, aerodynamic force may have significant component force in x-axis, so the drag in nonlinear case is bigger than that in linear case. Since the aerosurface keeps planar, there is no side force in linear case. In nonlinear case, the aerodynamic force acts vertically to the large deformed aerosurface so there is quite large component force acting along the span as the side force. The side force is comparable to the lift and contributes significantly to the root bend moment; thus, it cannot be ignored in the nonlinear analysis. The elastic torsional angle increases the effective angle of attack; thus, the maximum aerodynamic loads along the span move to the outer wing.Figure 15
Drag, side force, and lift along spanwise.
(a)
Lift
(b)
Drag
(c)
Side force


## 4.2. Flutter Analysis

The structural dynamics are different under various deflections and force cases due to geometrically nonlinear effects. Figure16 shows the structural dynamic frequency under different static load cases. The reduction of the horizontal bend frequency and the increase in the twist frequency are distinct. An iterative process is used to obtain the exact flutter boundary under certain flight conditions, and the effect of gravity is also considered. The results of different cases are compared with the experiment results to illustrate the nonlinear flutter characteristics and validate the theoretical analysis.Figure 16
Structural frequency under nonlinear static deflections.Figure17 shows the V
-
F and V
-
g curves of the traditional linear results. A classical bend/twist coupling flutter pattern results in a flutter speed at approximately 33.5 m/s, and the horizontal bend mode does not contribute to the flutter; its frequency remains unchanged. The 1st vertical mode frequency reduces to zero when the speed is over 34 m/s but the double calculation and the wind tunnel test indicate that there is no divergence occurring around 34 m/s.Figure 17
V
-
F and V
-
g curves of the linear flutter analysis.The nonlinear stability analysis results are presented below. AllV
-
g curves are found to be significantly different from the linear results with an unstable mode switching to a horizontal bend in the nonlinear analysis. This is caused by the geometrically nonlinear effect, which makes the horizontal bending mode contain large torsional components and becomes unstable when the structure experiences a large deformation, as shown in Figure 18. It is a remarkable characteristic of nonlinear stability analysis of VFAs that the horizontal bend mode plays an important role in flutter characteristics.Figure 18
Horizontal bending mode in a nonlinear equilibrium state.Figure19 shows the stability analysis results under the equilibrium states of 20 and 25 m/s. The corresponding predicted flutter speeds are approximately 31 and 33 m/s. Typically, the flutter speed decreases with increasing airspeed in the nonlinear analysis, but the test wing induces a negative deformation when the airspeed is below 28 m/s. Thus, the predicted flutter speed of 25 m/s is larger than that of 20 m/s. Because of the disagreement between the static equilibrium flight speed and predicted flutter speed, only the data near the static state marked red in Figure 19 are reliable. An iterative process is performed to narrow the disagreement and search for the exact flutter speed. Figure 20 shows the V
-
g and V
-
F curves from the dynamic aeroelastic analysis under the static equilibrium flight speed of 31.5 m/s. The predicted critical speed is 31.5 m/s, which is the theoretical exact nonlinear flutter speed of the analysis model under 0.5° angle of attack. Again, the data in Figure 20 is only accurate near 31.5 m/s, and if completely accurate V
-
g and V
-
F curves were desired, the divisional accurate data should be pieced together, as is in Figure 21, which requires more data to make the curve smoother and more reliable.Figure 19
Predicted flutter speed under two flight cases.Figure 20
V
-
g and V
-
F curves under 31.5 m/s.Figure 21
Accurate data pieced together roughly.The calculated cases above are under a 0.5° angle of attack. Due to geometric nonlinearity, cases under different angles of attack lead to different flutter speeds. The iterative processes searching for exact flutter boundary under different angles of attack have been implemented and the final results are shown in Figure22. The flutter results are summarized in Table 2. The aerodynamic loads and structural deflections increase as the angle of attack increases. Due to geometric nonlinearity, large deflections and load conditions easily cause the horizontal bending to become unstable. Thus, both the test and nonlinear analysis results indicate that the flutter speed decreases with an increasing angle of attack, whereas no changes are noted in the linear analysis under different cases. The first three cases present nearly identical deviations between the nonlinear analysis results and test results (in Figure 23). This effect may be caused by the conservative aerodynamic computation that can be validated in the deflection curve in Figure 10, where the analysed displacement is larger than that found in the test when the displacement is positive. Thus, it is understandable that the analysed flutter speed is smaller than the test results. The deviation narrows in the last case because the test result changes to a stall flutter depending on the flutter phenomena observed in the test; the stall model is not considered in theoretical aerodynamic computation, but the decline in flutter speed can still be reflected. In our test cases, the linear analysis results sometimes are less conservative than nonlinear results. This is not a referable conclusion. For most cases, the nonlinear analysis results are more accurate and reliable than linear results, especially when the nonlinear flutter speed is lower than the linear results.Table 2
Flutter speeds under different angles of attack.
Angle of attack
0°
0.5°
1.0°
1.5°
Linear
Computation
Experiment
Computation
Experiment
Computation
Experiment
Computation
Experiment
Computation
Flutter speed (m/s)
33.0
35.7
31.5
34.4
28.0
32.2
26.5
27.0
33.5
Flutter frequency (Hz)
7.60
6.94
7.82
6.88
7.34
7.00
7.53
13.38
7.29
Unstable mode
First horizontal bend in plane
First twistFigure 22
Nonlinear analysed flutter speed under different angles of attack.Figure 23
Flutter limit speed at different angles of attack.


## 5. Conclusion

A theoretical analysis framework has been introduced for very flexible wing structures that present notable geometric nonlinearity. The NVLM and NDLM are combined, and the proposed code has been developed to obtain the steady and unsteady aerodynamic loads for aeroelastic stability analysis. Wind tunnel testing is used to demonstrate the geometrically nonlinear characteristics of a large-aspect-ratio wing and validate the theoretical analysis results. Reasonable agreement between the computational results and test results has been obtained; all the results indicated that the static aeroelastic responses and flutter characteristics of very flexible wings are significantly different compared with the traditional linear analysis results. The geometrically nonlinear aeroelastic stability is related to certain flight conditions, and large elastic deformations make the horizontal bend mode unstable and may decrease the flutter speed, which have a significant influence on the flight envelope. The NVLM/NDLM and nonlinear FEM with the quasimodal dynamic approach make the nonlinear aeroelastic stability analysis feasible and more accurate. The theoretical analysis process established in this paper follows conventional linear analysis ways with some significant modifications considering VFAs’ geometrical nonlinearity both theoretically and practically and can thus be easily used to estimate the occurrence of critical aeroelastic behaviours of VFAs in engineering analysis. Efforts to apply the methodology to a flexible aircraft including control surfaces and coupled flight dynamics are ongoing; some aerodynamic derivatives may also be discussed in future work.



---

*Source: 5090719-2016-05-04.xml*