International Forum on Aeroelasticity and Structural Dynamics

IFASD 2017
25-28 June 2017 Como, Italy

A FEM-BASED APPROACH FOR NONLINEAR AEROELASTIC

TRIM OF HIGHLY FLEXIBLE AIRCRAFT
Cristina Riso1 , Fausto G. Di Vincenzo2 , Markus Ritter3 , Carlos E. S. Cesnik4 , and
Franco Mastroddi1
1 Department of Mechanical and Aerospace Engineering
Sapienza University of Rome
via Eudossiana, 18 Roma 00184 Italy
[Link]@[Link]
2 [Link] Toulouse
Rue du Prof. Pierre Vellas, 4 Toulouse 31300 France
3 Institute
of Aeroelasticity
German Aerospace Center (DLR)
Bunsenstraße, 10 Göttingen 37073 Germany
4 Departmentof Aerospace Engineering
University of Michigan
Beal Avenue, 1320 Ann Arbor 48109 Michigan (USA)

Keywords: Highly Flexible Aircraft, Nonlinear Aeroelastic Trim, Inertia Relief, Finite Ele-
ment Method, Fluid-Structure Interaction.

Abstract: A novel approach is presented to compute nonlinear aeroelastic trim of highly flexi-
ble aircraft by coupling off-the-shelf solvers for structures and aerodynamics. The methodology
includes an inertia relief algorithm for large displacements to perform nonlinear static analysis
of unrestrained bodies, which is currently not possible in commercially available finite element
codes. This feature allows to simulate free-free boundary conditions at each step of the non-
linear aeroelastic trim iterative solution sequence. The proposed approach is implemented by
coupling the [Link] nonlinear structural solver SOL 400 with a three-dimensional vortex
lattice method code. Loads and displacements interpolations between aerodynamic and struc-
tural grids are performed using [Link] six degree-of-freedom splines. The computational
environment is tested on a beam-type model of the University of Michigan’s X-HALE, a highly
flexible experimental unmanned aerial vehicle for nonlinear aeroelastic flight tests. Numeri-
cal results are compared with the University of Michigan’s Nonlinear Aeroelastic Simulation
Toolbox and with an aeroelastic solver developed at the German Aerospace Center.

1 INTRODUCTION

Highly flexible aircraft exhibit low elastic frequencies and experience large trim deflections un-
der typical load conditions. Rigid-body and structural dynamics are thus inherently coupled,
and unique flight dynamic/aeroelastic stability, response, and control should be analyzed in
statically deformed configuration rather than for the baseline structure [1]. In the past decades,
unique formulations of flight dynamics and nonlinear aeroelasticity have been developed by

1
 IFASD-2017-170

coupling geometrically exact beam models [2, 3] with two-dimensional potential-flow aerody-
namics [4, 5]. As a result, low-order simulation tools are nowadays available to study very
flexible vehicles modeled as beam-type structures [6, 7]. However, high-fidelity modeling is
necessary to analyze complex configurations with sufficient detail since the early design phase.
In the limit of validity of linear aeroelastic theory, finite element method (FEM) solvers are
commercially available for trim, stability, and response analysis of flexible aircraft [8]. These
tools are standardly used in aeroelastic design and optimization of conventional, relatively stiff
configurations, for which it is adequate to assume linear aeroelastic behavior and neglect cou-
pling between rigid-body and elastic degrees-of-freedom. Conversely, high-fidelity simulation
of very flexible vehicles is still a challenge for computational aeroelasticity. Despite devel-
opment of state-of-the-art methodologies for structures, aerodynamics, and flight dynamics as
separate disciplines, there is still a lack of well-established multidisciplinary approaches and
computational environments for integrated analysis of complex highly flexible configurations.

This paper presents a novel approach to compute nonlinear aeroelastic trim of very flexible
aircraft by coupling off-the-shelf solvers for high-fidelity modeling of structures and aerody-
namics. The work is a first step to extend a FEM-based simulation environment for coupled
flight dynamics and aeroelasticity [9] in order to study very flexible vehicles, for which stability
and response should be addressed in statically deformed configuration. The main novelties in
the proposed FEM-based nonlinear aeroelastic trim methodology are 1) the use of six degree-of-
freedom splines for loads and displacements interpolations between aerodynamic and structural
grids; 2) a load-stepping relaxation procedure to implement a close fluid-structure coupling;
and 3) an inertia relief algorithm [10–12] for large displacements to simulate free-free boundary
conditions at each step of the iterative solution sequence. Indeed, clamped boundary conditions
are typically assumed in the loop when using commercial displacement-based FEM solvers to
compute structural deflection, since no inertia relief feature is currently available for nonlin-
ear analyses [13]. The proposed methodology is implemented for the [Link] nonlinear
structural solver SOL 400 [13] coupled with a three-dimensional vortex lattice method (VLM)
code [14] used in previous studies [15, 16]. [Link] six degree-of-freedom splines gen-
erated by means of the [Link] Hybrid Static Aeroelasticity (HSA) Toolkit [17, 18] are
used to transfer loads and interpolate displacements between structural and aerodynamic grids.
The developed computational environment is applied to a beam-type model of the University
of Michigan’s X-HALE [19], a highly flexible experimental unmanned aerial vehicle (UAV)
for nonlinear aeroelastic flight tests. Numerical results are compared with solutions from the
University of Michigan’s Nonlinear Aeroelastic Simulation Toolbox (UM/NAST) [7] and from
a FEM-based aeroelastic solver developed at the DLR [15, 16].

The paper is organized as follows. The proposed methodology for nonlinear aeroelastic trim
is presented in Sec. 2. Numerical results are discussed in Sec. 3. These include preliminary
test cases to assess the core [Link] SOL 400/VLM aeroelastostatic solver and a com-
plete nonlinear aeroelastic trim analysis of the X-HALE Reduced Risk Vehicle (RRV) [19, 20].
Concluding remarks and future work are outlined in Sec. 4.

2 NONLINEAR AEROELASTIC TRIM OF HIGHLY FLEXIBLE AIRCRAFT

Solving an aeroelastic trim problem consists of computing the combination of aerodynamic

angles, control inputs, and structural deflection for which a flexible vehicle is balanced in a
given equilibrium flight condition. Assuming linear elastic theory and linearized potential-
flow aerodynamics, aeroelastic trim analysis of complex configurations is standardly performed

2
 IFASD-2017-170

via commercial solvers that integrate displacement-based FEM structural modeling and steady
Doublet Lattice Method (DLM) aerodynamics [8]. Unfortunately, no integrated software is
commercially available to solve nonlinear aeroelastic trim problems. Hence, separate solvers
for structures and aerodynamics must be conveniently coupled for trim analysis of complex
very flexible configurations.

This section presents a novel methodology for nonlinear aeroelastic trim which is based on the
coupling of off-the-shelf solvers for structures and aerodynamics to allow high-fidelity model-
ing of both disciplines. The proposed methodology is compared to the algorithms implemented
in UM/NAST [7] and in a DLR FEM-based aeroelastic toolbox [15, 16].

2.1 Proposed algorithm

The proposed architecture for trim analysis of very flexible vehicles is shown in Fig. 1. The
solution sequence consists of the Newton-Raphson loop in Fig. 1(a) and the nested nonlinear
aeroelastostatic loop in Fig. 1(b). In the present work, the latter is implemented by coupling
the [Link] nonlinear structural solver SOL 400 with a three-dimensional VLM code [14]
developed for previous studies [15, 16]. However, the algorithm in Fig. 1(b) is applicable to
generic off-the-shelf structural and aerodynamic solvers. At each iteration of the outer Newton-
Raphson loop in Fig. 1(a) the nonlinear aeroelastostatic solver in Fig. 1(b) is called to obtain
the structural deflection and steady aerodynamic load for the current approximation to the trim
variables. These outputs are used to verify the vehicle equilibrium and, in the negative case, to
update the trim variables estimates by numerically computing the system Jacobian matrix. The
nonlinear aeroelastostatic solution algorithm in Fig. 1(b) is the core of the proposed methodol-
ogy and consists of the following steps:

1. The aerodynamic solver is called to compute the aerodynamic load distribution on the
current aerodynamic model configuration;
2. The aerodynamic load is transferred to the structural grids by means of six degree-of-
freedom splines;
3. The gravity field is added to compute the residual non-balanced load resultant;
4. The resulting inertia relief load is added to aerodynamics and gravity to obtain a self-
equilibrating load distribution;
5. The overall load is multiplied by a scaling factor for relaxation;
6. The structural solver is called to perform an incremental nonlinear static analysis that
starts from the deformed configuration achieved at the preceding iteration and takes into
account differential stiffness;
7. The incremental displacement field is interpolated by means of six degree-of-freedom
splines to update the aerodynamic model;
8. The preceding steps are repeated until aeroelastostatic deflection convergence is reached.

The main novelties of the proposed algorithm are 1) the use of six degree-of-freedom splines
to transfer loads and interpolate displacements between aerodynamic and structural grids; 2)
the load-stepping relaxation procedure that implements a close fluid-structure coupling; and
3) the large-amplitude inertia relief algorithm to simulate free-free boundary conditions in the
loop. These features are discussed in the following subsections. Note that the nonlinear aeroe-
lastostatic solution algorithm in Fig. 1(b) can be also used as a stand-alone solver. Clamped
structures can be simulated by skipping the inertia relief computations.

3
 IFASD-2017-170

(a) Nonlinear trim solver (b) Nonlinear aeroelastostatic solver

Figure 1: Proposed architecture for trim analysis of very flexible aircraft.

2.1.1 Fluid-structure coupling

Splines standardly used in [Link] aeroelastic solvers are thought for coupling of FEM
structural models with the linearized potential-flow aerodynamic models generated internally
within [Link], typically DLM models [8]. Consistently with linear aeroelastic theory,
the degrees-of-freedom utilized at the coupled grids include only normal displacements for
surface splines, and additionally torsional rotations for linear splines [8]. Hence, this approach
may be not suitable for fluid-structure coupling with external VLM or CFD models. Indeed,
the aerodynamic load defined on each body panel may be not limited to a lift force and pitching
moment. In addition, all six degrees-of-freedom should be taken into account in displacements
interpolation in presence of geometrically nonlinear behavior and/or to couple with viscous
aerodynamic models.

To improve coupling of FEM models with external aerodynamics, an enhanced spline tech-
nology has been developed in [Link] [18]. The new splines allows the aerodynamic
grids to have six degrees-of-freedom, so that all components of load and displacement fields
are taken into account in the fluid-structure interaction. Both surface and linear six degree-of-
freedom splines are available [17, 18], which can be used in place of the standard ones used for
[Link] FEM/DLM coupling [8]. Since beam-type structural models are considered in all
the present numerical studies, linear six degree-of-freedom splines are used in this work. These
enable to directly interface one-dimensional structures and aerodynamic surfaces with no need
to define additional leading-edge and trailing-edge grids [8]. The six degree-of-freedom splines
are generated by means of the [Link] Hybrid Static Aeroelasticity (HSA) Toolkit [18],
which enables to automatically subdivide structural and aerodynamic components in multiple
spline patches with suitable overlapping in order to improve the speed and smoothness of fluid-
structure coupling. The spline matrix may be updated in the loop (see Fig. 1(b)) to account for
the aerodynamic and structural models deflection. The deformation of the splines axes is taken
into account by defining a local coordinate system attached to a suitable structural grid for each
spline patch.

4
 IFASD-2017-170

2.1.2 Relaxation procedure

In case geometric nonlinearities are present in a static analysis, equilibrium of applied loads and
internal stresses must be satisfied in deformed configuration, which requires an iterative solution
approach. For aeroelastic analyses this implies repeated calls to the aerodynamic and structural
solvers as shown in Fig. 1(b). Note that the aerodynamic load considered in each nonlinear
structural analysis is evaluated on the deformed configuration achieved at the previous iteration.
Therefore, the load distribution may be not consistent with the updated model in presence of
large incremental deflections. In the proposed algorithm, this remark also applies to the inertia
relief loads that are included in the loop to simulate free-free boundary conditions. Indeed,
the inertia relief load distribution depends on the non-balanced force and moment resultant of
aerodynamics and gravity and on the model rigid-body mass matrix. All these quantities vary
in the loop due to structural deflection. Moreover, follower-force effects must be taken into
account in a nonlinear solution. In [Link] SOL 400, follower forces and moments can
be implemented by means of specific entries [13, 17]. However, these require to locally define
loads orientation in terms of suitably located grids [17], which may make implementation of
follower forces cumbersome for complex distributions.

In order to overcome these problems, a close fluid-structure coupling is implemented in the

nonlinear aeroelastostatic solver based on the load-stepping relaxation procedure shown in
Fig. 1(b). The structure is loaded incrementally in ns fluid-structure iterations, corresponding to
a load step size 1/ns . Each nonlinear structural analysis starts from the deformed configuration
achieved at the preceding step. The load field applied to the structure in the first ns iterations
is multiplied by an increasing scaling factor j/ns (j = 1, . . . , ns ). Hence, ns iterations of the
aeroelastic static loop are necessary to take the structure from an unloaded to the fully loaded
condition. The aerodynamic and inertia relief loads are recomputed at each iteration on the
updated model to replace the preceding distributions, and the scaling factor for the overall ap-
plied load increases until it becomes 1. At this point ns iterations of the aeroelastostatic loop
in Fig. 1(b) have been performed, so that the structure is in a fully loaded condition. However,
further iterations may be still necessary for convergence. Additional na fluid-structure iterations
are thus performed until the aeroelastostatic deflection converges using a load scaling factor al-
ways equal to 1. The total number of fluid-structure iterations in the loop is thus ni = ns + na ,
obtained by summing ns incremental load steps and na additional iterations performed after the
structure is fully loaded.

The proposed load-stepping relaxation procedure improves stability and accuracy of the aeroe-
lastostatic solution algorithm by implementing a close fluid-structure coupling, which may also
save computational time (see Subsecs. 3.1.2 and 3.2). Indeed, the step size can be tailored to
each test case to keep incremental deflections moderate between subsequent fluid-structure it-
erations, which improves the algorithm robustness and speed. This also allows to assume the
aerodynamic load as non-follower within each iteration, since incremental deflections are lim-
ited by the load stepping and follower-force effects can be thus accurately simulated by updating
the load at each call to the structural solver. Hence, standard non-follower force and moment
entries [17] can be used to input aerodynamic loads in the [Link] SOL 400.

2.1.3 Inertia relief algorithm

Conventional displacement-based FEM static analysis cannot be performed on an unrestrained

body, since the model stiffness matrix becomes singular in this condition. However, commercial

5
 IFASD-2017-170

displacement-based FEM solvers use the inertia relief concept [10–12] to simulate free-free
structures in linear structural and aeroelastic analyses [8, 12]. If a free structure is subjected to
a non-balanced load field, a rigid-body acceleration occurs. When using inertia relief, apparent
inertial loads are added in the analysis to obtain an overall self-balancing load distribution.
This allows to constrain the model at an arbitrarily selected grid for the purpose of stiffness
matrix inversion while still simulating free-free boundary conditions, as the constraint gives
zero reaction in presence of a self-balancing load field. In this case, elastic deflections obtained
with different choices of the constrained grid can be superimposed by means of a rigid-body
displacement, so that they represent the unique deformed state of the unrestrained structure
under the original non-balanced load [21].

Unfortunately, no inertia relief mechanism is currently implemented for nonlinear static and
dynamic analyses in commercially available FEM solvers [13]. As a consequence, clamped
boundary conditions must be assumed. In the case of nonlinear aeroelastic trim analyses, the
applied load distribution is by definition self-balancing when the trim condition is achieved,
so that the converged aeroelastostatic deflection is representative of an unrestrained structure
even if the model is clamped at some grid for stiffness matrix inversion. However, free-free
boundary conditions are not simulated in the nonlinear structural analyses performed in the trim
loop unless apparent inertial loads are taken into account. Indeed, aerodynamics and gravity do
not give a self-balancing load field until the trim condition is verified, which gives a non-zero
reaction experienced at the constrained grid if the inertia relief concept is not applied. This
may lead to larger deflections within the loop while converging to the trim solution, which may
negatively affect algorithm stability, accuracy, and computational time (see Subsec. 3.2).

In order to simulate free-free boundary conditions at each step of the nonlinear aeroelastic trim
solution sequence, an inertial relief algorithm for large displacements is implemented externally
to [Link] SOL 400 as shown in Fig. 1(b). The algorithm extends the inertial relief feature
implemented in [Link] linear aeroelastic solvers [8, 11, 12] by iteratively recomputing
and updating inertia relief loads in the loop. This is necessary for nonlinear problems since 1)
aerodynamic loads are updated in the loop; 2) the application points of aerodynamic and gravity
loads move; and 3) the model mass distribution varies with structural deflection.

The proposed inertia relief algorithm is implemented as follows. A fixity structural grid is
chosen that experiences zero displacement in the nonlinear aeroelastic trim solution, so that
the relative displacement field is computed. The fixity grid position in the FEM model in the
assumed coordinate system is given by x c and remains constant in the iterative solution process.
At the jth iteration of the nonlinear aeroelastostatic loop, the force and moment resultants of
aerodynamics and gravity with respect to the fixity grid are given by:
N
X N
X
F jtot = (f j
i,a +f j
i,g ) M jtot = [ (x ji − x c ) × ( f ji,a + f ji,g ) + m ji,a ] (1)
i=1 i=1

In the above equation, the quantities f ji,a , f ji,g , and m ji,a are respectively the pin aerodynamic
force, gravity force, and aerodynamic moment applied at the jth iteration on the ith grid of
current position x ji . The apparent rigid-body translational and rotational accelerations due to
the net force and moment resultants in Eq. (1) are recovered from:
( j ) ( )
j
a rb,t F tot
M jrb =− (2)
a jrb,r M jtot

6
 IFASD-2017-170

where M jrb is the 6 × 6 model rigid-body mass matrix with respect to the fixity grid evaluated
at the jth iteration. The acceleration loads due to the apparent accelerations in Eq. (2) give the
inertia relief load distribution that is added to aerodynamics and gravity as shown in Fig. 1(b).
With this approach, the overall load distribution applied to the structure is self-balancing at
each iteration even if the trim condition is not satisfied yet. Consequently, free-free boundary
conditions are simulated in the loop. Indeed, although the fixity grid is clamped for stiffness
matrix inversion, no constraint is actually experienced by the model since there is no reaction
in presence of a self-balancing load distribution. Besides the present use in the nonlinear aeroe-
lastic trim architecture shown in Fig. 1, note that the proposed inertia relief algorithm can be
also applied to simulate free-free boundary conditions in generic nonlinear static and dynamic
analyses.

2.2 Comparison with existing algorithms

In the present work, numerical results obtained with the proposed methodology are compared
with predictions from the nonlinear aeroelastic trim solvers implemented in UM/NAST [7] and
in an aeroelastic toolbox developed at the DLR [15, 16]. The main features of the examined
approaches are summarized below and in Tab. 1.

Nonlinear structures are modeled in UM/NAST using a strain-based geometrically exact beam
formulation [3]. Nonlinear aeroelastic trim is performed by either solving the force/moment
equilibrium or forcing zero linear and angular accelerations at the origin of the body reference
system. This is accomplished by means of an outer Newton-Raphson loop as the one shown
in Fig. 1(a) and an inner nonlinear aeroelastostatic solver that couples the strain-based beam
formulation with strip theory or, more recently, with a three-dimensional VLM code [14] along
with gravity loads and load factors. The aerodynamic model and loads are updated at each
aeroelastostatic iteration by recovering structural displacements from strains. Fluid-structure
coupling is performed by transferring each aerodynamic chordwise load distribution to the cor-
responding beam axis location and by recovering the rigid-body cross-sectional displacement
field from the beam axis deflection. Inertia relief is not necessary in UM/NAST since the model
stiffness matrix is always definite positive in the strain-based formulation.

The DLR aeroelastic toolbox is a FEM-based simulation environment that allows analysis
of complex configurations using linear, reduced-order, or fully nonlinear structural formula-
tions [15,16]. The nonlinear aeroelastic trim solver consists of a main Newton-Raphson loop as
the one in Fig. 1(a) and a nested nonlinear aeroelastostatic loop based on [Link] SOL 400
coupled with the VLM code [14] considered in Ref. [15, 16]. The same code is used to imple-
ment the proposed methodology and it has been also interfaced with UM/NAST to allow three-
dimensional aerodynamic modeling. Fluid-structure coupling is performed using Radial Basis
Functions (RBF) and, for beam-type FEMs, the coupling model approach [15, 16] to expand
one-dimensional structures. No load-stepping relaxation procedure is currently implemented
in the loop. Each [Link] SOL 400 nonlinear static analysis starts from the undeformed
configuration, but considers the aerodynamic load distribution computed on the updated VLM
model. A structural grid is clamped in the loop to invert the FEM stiffness matrix.

3 NUMERICAL RESULTS
Preliminary analyses to assess the core nonlinear aeroelastostatic solver are discussed in Sub-
sec. 3.1. The proposed nonlinear aeroelastic trim methodology is applied to the X-HALE
RRV [19, 20] in Subsec. 3.2.

7
 IFASD-2017-170

Table 1: Main features of the implemented nonlinear aeroelastic trim solver based on the solution algorithm pro-
posed in Fig. 1 and comparison with the UM/NAST and DLR solvers.

Feature Present solver UM/NAST solver DLR solver

Structural modeling Complex FEM Strain-based beam Complex FEM
Aerodynamic modeling VLM VLM VLM
Coupling approach 6 DOF splines Section-based RBF
Relaxation Yes No No
Inertia relief Yes Not applicable No

Figure 2: Highly flexible 16-meter wing aeroelastic model [15].

3.1 Assessment of the nonlinear aeroelastostatic solver

The nonlinear aeroelastostatic solution algorithm in Fig. 1(b) is implemented for [Link]
SOL 400 coupled with the VLM code used in Refs. [15, 16]. Fluid-structure coupling is
performed via [Link] six degree-of-freedom splines [17] generated by means of the
[Link] HSA Toolkit [18]. The developed nonlinear aeroelastostatic solver is assessed
by computing the response of clamped structures at prescribed free-stream velocity and angle
of attack. The inertia relief algorithm is switched off in the solution sequence in Fig. 1(b) to
implement clamped boundary conditions. Two test cases of increasing complexity are con-
sidered: 1) the highly flexible 16-meter wing studied in Ref. [15]; and 2) the X-HALE RRV
analyzed in Ref. [16]. Beam-type FEM models are considered for both applications to directly
compare the deflections with UM/NAST. The results are also compared with the [Link]
SOL 400/VLM nonlinear aeroelastostatic solver of the DLR toolbox.

3.1.1 Test case 1: 16-meter beam

The first test case is the highly flexible 16-meter wing studied in Ref. [15]. Nonlinear aeroelas-
tostatic analyses are carried out using 1) the present [Link] SOL 400/VLM solver; 2) the
UM/NAST solver based on strip-theory aerodynamics corrected with weighting factors [15]; 3)
the UM/NAST solver based on VLM aerodynamics; and 4) the [Link] SOL 400/VLM
solver of the DLR toolbox. The aeroelastic model used for the present computations is shown in

8
 IFASD-2017-170

Table 2: Nonlinear aeroelastostatic response of the 16-meter wing at sea level, U = 40 m/s, and α = 5◦ with no
gravity: tip vertical displacement over semispan, computational time, and number of iterations using the
proposed approach with different step size.

Step size (%) uztip (%) ∆uztip (%) Time (s) ns na ni

100 20.3053 – 125.52 1 4 5
50 20.2622 -0.2124 176.44 2 4 6
20 20.2147 -0.2343 206.60 5 4 9
10 20.1959 -0.0931 264.00 10 3 13
5 20.1862 -0.0481 512.84 20 3 23

Table 3: Nonlinear aeroelastostatic response of the 16-meter wing at sea level, U = 40 m/s, and α = 3◦ , 4◦ , 5◦ with
no gravity: tip vertical displacement over semispan (%) obtained using the proposed approach, UM/NAST
with strip theory [13], UM/NAST with VLM, and the DLR toolbox.

α = 3◦ α = 4◦ α = 5◦
Proposed approach 12.5996 16.5067 20.1959
UM/NAST with VLM 12.3633 16.1815 19.7798
UM/NAST with strip theory 12.5964 16.4977 20.1787
DLR toolbox 12.6233 16.5733 20.3294

Fig. 2 and consists of a single-member beam-type structure with quadratically varying stiffness
properties coupled with a unit-chord rectangular lifting surface [15].

The aeroelastostatic response of the 16-meter wing is studied at sea level for fixed free-stream
velocity U = 40 m/s and increasing angle of attack α = 3◦ , 4◦ , 5◦ . Gravity is not taken into
account. In order to choose the number of load steps for the relaxation procedure in Fig. 1(b), a
preliminary sensitivity analysis is carried out for α = 5◦ by performing multiple computations
with different step size. The variation of tip vertical displacement over semispan, computational
time, and number of iterations with the step size is summarized in Tab. 2. The tip vertical
displacement decreases by reducing the step size, but no practical change is found for step size
below 10%. Hence, this value is assumed for all the analyses. For this test case, the total
number of iterations ni increases by reducing the step size due to the consequent increasing

Figure 3: Nonlinear aeroelastostatic response of the 16-meter wing at sea level, U = 40 m/s and α = 3◦ , 4◦ , 5◦
with no gravity: deformed configuration in the vertical plane obtained using the proposed approach,
UM/NAST with strip theory [13], UM/NAST with VLM, and the DLR toolbox.

9
 IFASD-2017-170

Figure 4: X-HALE RRV aeroelastic model [16].

number of load steps. However, the na additional iterations for convergence after the structure
is fully loaded decrease from 4 to 3 for smaller step size. This behavior is expected since a
closer fluid-structure coupling is implemented by reducing the step size, which gives a smoother
approaching of the converged solution. This does not result in smaller computational time
for this simple test case since the analysis converges fast even without load-stepping, and the
number of iterations necessary after the structure is fully loaded shows just a slight dependency
on the step size (see Tab. 2). However, the load-stepping relaxation procedure may allow to
save computational time for more complex configurations as shown in Subsec. 3.1.2.

The aeroelastostatic deflection in the vertical plane obtained for α = 3◦ , 4◦ , 5◦ with different
approaches are compared in Fig. 3. The tip vertical deflections over semispan are also compared
in Tab. 3. The solutions predicted with the proposed approach match the UM/NAST results
based on VLM aerodynamics at all the angles of attack. The DLR toolbox solution practically
lies on top of the previous ones, just showing a 0.65% larger tip vertical displacement for α =
5◦ . The UM/NAST deflection based on strip-theory aerodynamic is always smaller than with
VLM. The relative difference in the tip vertical displacement increases from 1.85% to 2% with
the angle of attack. This is justified by the larger deflection experienced at higher angle of
attack, since the weighting factors used to tune the strip-theory aerodynamic load distribution
are evaluated for the undeformed wing [15].

3.1.2 Test case 2: X-HALE RRV

The second test case is the X-HALE RRV studied in Ref. [16] and also considered for the non-
linear aeroelastic trim analyses. This case study enables to assess the nonlinear aeroelastostatic
solver for a realistic configuration composed by multiple structural members and interacting
aerodynamic surfaces. Nonlinear aeroelastostatic analyses are carried out using 1) the present
[Link] SOL 400/VLM solver; 2) the UM/NAST solver based on VLM aerodynamics;
and 3) the [Link] SOL 400/VLM solver of the DLR toolbox. The UM/NAST results
based on strip-theory aerodynamics in Ref. [16] are not considered for comparison since strip
theory does not take into account mutual interactions between aerodynamic surfaces, which are
significant for the X-HALE RRV. The aeroelastic model used for the present computations is

10
 IFASD-2017-170

Table 4: Nonlinear aeroelastostatic response of the X-HALE RRV for ρ = 1.22161 kg/m3 , U = 16 m/s, and
α = 1◦ with gravity: right-half wing tip vertical displacement over semispan, computational time, and
number of iterations using the proposed approach with different step size.

Step size (%) uztip (%) ∆uztip (%) Time (s) ns na ni

100 17.4326 – 1093.96 1 28 29
50 16.7863 -0.6463 727.21 2 20 22
20 16.4141 -0.3722 808.01 5 19 24
10 16.2952 -0.1189 931.56 10 17 27
5 16.2368 -0.0584 1227.02 20 15 35

Table 5: Nonlinear aeroelastostatic response of the X-HALE RRV for ρ = 1.22161 kg/m3 , U = 16 m/s, and
α = 0◦ , 0.5◦ , 1◦ with gravity: right half-wing tip vertical displacement over semispan (%) obtained using
the proposed approach, UM/NAST with VLM, and the DLR toolbox.

α = 0◦ α = 0.5◦ α = 1◦
Proposed approach 11.7419 14.0734 16.2952
UM/NAST with VLM 11.8773 14.3971 16.8117
DLR toolbox 11.6878 14.0520 16.3138

shown in Fig. 4 and consists of a flexible beam-type member to model the wing elastic axis and
rigid bars to represent booms, tails, pods, and fins. The wing is divided in six unit segments,
with the outer ones presenting a 10◦ positive dihedral. The wing aerodynamic model features
a reverse-camber EMX07 airfoil and has an incidence of 5◦ with respect to the horizontal axis.
Flat-plate surfaces are used to model pods, fins, and tails. Further details on the X-HALE RRV
configuration and model are found in Refs. [16, 19, 20]

The aeroelastostatic response of the X-HALE RRV is studied for ρ = 1.22161 kg/m3 , fixed
free-stream velocity U = 16 m/s and increasing angle of attack α = 0◦ , 0.5◦ , 1◦ . Gravity is
included in the analysis (g = 9.80665 m/s2 ). Solution sensitivity to the step size is preliminarily
analyzed for α = 1◦ as summarized in Tab. 4. The tip vertical displacement decreases by
reducing the step size as found for the 16-meter wing, although a higher sensitivity to the load
path is observed for this more complex test case. In contrast with what observed for the 16-meter

Figure 5: Nonlinear aeroelastostatic response of the X-HALE RRV for ρ = 1.22161 kg/m3 , U = 16 m/s and
α = 0◦ , 0.5◦ , 1◦ with gravity: right half-wing deformed configuration in the vertical obtained using the
proposed approach, UM/NAST with VLM, and the DLR toolbox.

11
 IFASD-2017-170

wing, note that computational time does not monotonically increase by reducing the step size for
the X-HALE RRV. Indeed, it first decrease by reducing the step size from 100% (no relaxation)
to 50%, then it increases for further reductions. However, the time necessary to complete the
analysis with load steps of 10% is still smaller than with no relaxation. This occurs because
the number of load-step iterations ns increases by reducing the step size, but the additional
na = ni − ns iterations performed after the structure is fully loaded significantly decrease from
28 to 15 (see Tab. 4). Since incremental deflections are very small after taking the structure to
a fully loaded condition, the last na iterations are faster than the first ns ones associated with
load steps. As a consequence, the computational time necessary with relaxation is smaller than
with no relaxation for step size not below 10%, since in this case most of the total iterations
are performed after the structure is fully loaded (ns < na ). Conversely, computational time
increases with respect to the case of no relaxation for step size below 10%, since in this case
most of the total iterations are load steps (ns > na ). This behavior is not observed for the 16-
meter wing, since the nonlinear aeroelastic analyses are very fast due to the model simplicity
and the number of iterations performed after the structure is fully loaded shows just a slight
dependency on the step size (see Tab. 2). Taking into account the above remarks and that the
X-HALE aeroelastostatic response practically does not change for step size below 10% (see
Tab. 4), this value is assumed for all the analyses. For α = 1◦ , the analysis with step size equal
to 10% is 15% faster than with no relaxation (see Tab. 4).

The aeroelastostatic results obtained for α = 0◦ , 0.5◦ , 1◦ with different solvers are compared in
Fig. 5 and Tab. 5. The responses predicted with the proposed approach and the DLR toolbox
practically overlap for all the angles of attack, whereas UM/NAST predicts larger displace-
ments. The difference at tip increases from 1.1% to 3% with the angle of attack. This fact may
be justified by considering that a section-based fluid-structure coupling is used in UM/NAST
(see Subsec. 2.2), so that the load transferred to each structural grid is only influenced by the
aerodynamic grids in a close neighborhood. Conversely, the fluid-structure coupling approaches
used in the present solver and in the DLR toolbox involve regions of the aerodynamic and
structural models, so that mutual influence between cross-sections is taken into account in load
transferal.

Additional analyses performed on the X-HALE RRV for α = 1◦ and increasing free-stream
velocity have pointed out that for larger nonlinearities the proposed solution algorithm is not
stable in the absence of relaxation (step size equal to 100%), and that the DLR solver cannot
achieve a convergent solution as well. This issue may be critical in presence of very large
nonlinearities. Further investigation on the solution sensitivity to the load path in terms of
stability, accuracy, and computational time for highly geometrically nonlinear behavior will be
carried out in a future work.

3.2 Nonlinear aeroelastic trim of the X-HALE

A complete longitudinal nonlinear aeroelastic trim analysis of the X-HALE RRV is performed
for level flight at ρ = 1.22161 kg/m3 and U = 16 m/s, which is a typical operative condition
for the vehicle [19, 20]. The aeroelastic model considered for trim is the same as described in
Subsec. 3.1.2, with the horizontal tails modeled as control surfaces to satisfy pitch equilibrium.
The nonlinear aeroelastic trim analysis is carried out using 1) the present solver composed by
the outer Newton-Raphson loop in Fig. 1(a) and the nonlinear aeroelastostatic solver assessed in
Subsec. 3.1; 2) the UM/NAST nonlinear aeroelastic trim based on VLM aerodynamics; and 3)
the nonlinear aeroelastic trim solver of the DLR toolbox. The computation of inertia relief loads

12
 IFASD-2017-170

is activated in the present solver to simulate free-free boundary conditions at each iteration. The
fixity grid for the inertia relief algorithm in Subsec. 2.1.3 is chosen at wing beam axis centerline,
so that the aeroelastic trim deflection is evaluated relatively to this grid. With this choice the
result can be directly compared with the DLR toolbox solution, which is computed by clamping
the same grid in the loop in order to invert the model stiffness matrix. The structural deflection
given by UM/NAST is also referred to this grid, but this is computed with no need to either
impose a constraint or use inertia relief since the strain-based stiffness matrix is always definite
positive. Since the X-HALE RRV experiences tip deflections of about 5% semispan in the
examined flight condition [16, 19], a larger step size equal to 20% is assumed in the present
solver for the trim analysis in place of the smaller value 10% used in Subsec. 3.1.2. Indeed,
the latter does not give any change in the solution for deflections of this order of magnitude.
With this choice, the trim solution converges 10% faster than with no relaxation (step size
equal to 100%). Linear aeroelastic trim predictions are also obtained for comparison using
the [Link] SOL 144 solver [8], which performs a one-shot fully linear aeroelastic trim
analysis based on steady DLM aerodynamics. Since only flat-plate surfaces can be modeled
via DLM, a correction to account for the X-HALE RRV wing incidence and camber must be
included in the analysis [8]. Two approaches are considered for this purpose. The first one
is to input a downwash distribution to take into account rotation of the DLM panels normal
vectors due to incidence and camber in the aerodynamic boundary conditions [8, 17]. The
second approach is to include an external load field to take into account the aerodynamic loads
experienced at zero angle of attack, which are not captured by DLM. In this work, the supplied
external load is given by the VLM aerodynamic load distribution on the undeformed vehicle at
the examined trim velocity, zero angle of attack, and zero tails deflection.

The aeroelastic trim results obtained with different approaches are summarized in Tab. 6, while
the trim right half-wing deflections in the vertical plane are compared in Fig. 6. The true-scale
trim deformed configuration and the VLM force distribution are shown in Fig. 7. A better agree-
ment in terms of structural deflection and tails rotation is found between the results from the
present solver and the DLR toolbox, whereas the UM/NAST and DLR toolbox trim angles of
attack are slightly closer. However, differences between UM/NAST and the other approaches
could be expected from the ones discussed in Subsec. 3.1.2. Since the X-HALE RRV ex-
periences a moderate deflection in the examined flight condition [16, 19], a good agreement
between nonlinear and linear predictions is found. A better matching of the VLM-based results
is achieved by accounting for wing incidence and camber using the external aerodynamic load
input than the downwash input [8,17]. This is justified by considering that only the aerodynamic
boundary conditions are influenced by downwash input, which does not alter the DLM model
geometry and thus the computation of the Aerodynamic Influence Coefficient (AIC) matrix.
Conversely, the external load included in [Link] SOL 144 to account for aerodynamic
loads at zero angle of attack is based on the actual VLM model geometry, and thus provides a
more accurate modeling of incidence and camber effects.

In order to point out the advantage of including inertia relief in the trim loop when using a com-
mercial (displacement-based) FEM solver compute structural deflection, an additional analysis
is carried out with the present methodology by assuming the same parameters as before but
switching the inertia relief algorithm off (see Fig. 1). The result is compared with the one com-
puted using inertia relief to focus only on the effect of this feature, which cannot be isolated by
comparing the present solution with the DLR toolbox one. Free-free boundary conditions are
simulated at each iteration in the present solver when the inertia relief algorithm is activated,
whereas clamped boundary conditions are used in the loop when inertia relief is off. Note that

13
 IFASD-2017-170

Table 6: Aeroelastic trim of the X-HALE RRV for level flight at ρ = 1.22161 kg/m3 and U = 16 m/s: right half-
wing tip displacement over semispan, angle of attack, and tails deflection obtained using the proposed
approach (inertia relief in the loop), UM/NAST with VLM, the DLR toolbox, and [Link] SOL
144 [8].

uztip (%) α (deg) δe (deg)

Proposed approach 5.0696 0.7503 2.4162
DLR toolbox 4.9992 0.7853 2.4716
UM/NAST 4.7010 0.7922 2.6162
[Link] SOL 144 with external load input 4.8966 0.7330 2.5962
[Link] SOL 144 with downwash input 4.6096 0.7155 2.6099

Figure 6: Aeroelastic trim of the X-HALE RRV for level flight at ρ = 1.22161 kg/m3 and U = 16 m/s: right half-
wing deformed configuration in the vertical plane obtained using the proposed approach (inertia relief in
the loop), UM/NAST with VLM, the DLR toolbox, and [Link] SOL 144 [8].

Figure 7: Aeroelastic trim of the X-HALE RRV for level flight at ρ = 1.22161 kg/m3 and U = 16 m/s: vehi-
cle deformed configuration in the vertical plane in true scale and VLM aerodynamic force distribution
obtained with the proposed approach (inertia relief in the loop).

14
 IFASD-2017-170

(a) Inertia relief off (b) Inertia relief on

Figure 8: Aeroelastic trim of the X-HALE RRV for level flight at ρ = 1.22161 kg/m3 and U = 16 m/s: right half-
wing deformed configuration in the vertical plane at each iteration of the trim loop (blue) and converged
solution (red) obtained using the proposed approach with inertia relief off/on.

Table 7: Aeroelastic trim of the X-HALE RRV for level flight at ρ = 1.22161 kg/m3 and U = 16 m/s: right
half-wing tip vertical displacement over semispan, angle of attack, tails deflection, computational time,
and number of iterations using the proposed approach with inertia relief off/on.

uztip (%) α (deg) δe (deg) Time (s) ni

Inrel off 5.1898 0.7203 2.4414 2511.21 87
Inrel on 5.0696 0.7503 2.4162 1111.14 32
∆% -0.12 4.16 -1.03 -55.75 -63.22

the converged solution is in both cases representative of an unrestrained vehicle, since the aero-
dynamic and gravity load fields give a self-balancing distribution at trim and no reactive load
is thus experienced at the clamped grid. However, the path followed to achieve the converged
solution is different in the two cases, since aerodynamics and gravity are not self-balancing in
the loop.

The intermediate right half-wing deflections in the loop (blue) and the converged solution (red)
obtained with inertia relief off/on are shown in Fig. 8. The trim right half-wing tip vertical
displacement over semispan, angle of attack, tails deflection, computational time, and number of
fluid-structure iterations are compared in Tab. 7. Figure 8 shows that a more regular convergence
to the trim solution is achieved using inertia relief, which results in a 56% computational time
saving with respect to the case of inertia relief off. This behavior is justified by considering
that a physics-based load path is followed using inertia relief, as the structure is simulated as
unrestrained while it incrementally deflects toward the trim solution. Indeed, aerodynamics
and gravity do not give a self-balancing load field in the loop, since the trim condition is not
satisfied yet. Therefore, the model experiences a non-zero reaction at the clamped grid when the
inertia relief algorithm is switched off. Conversely, although a grid is still kept fixed when the
inertia relief algorithm is activated, in this case the model does not experience any constraint.
The apparent inertial loads are added to aerodynamics and gravity in each nonlinear structural
analysis, so that the structure is always subjected to a self-balancing load distribution and no
reactive load is consequently experienced at the fixed grid. This results in smaller deflections
in the loop as shown in Fig. 8, which saves fluid-structure iterations and thus computational
time as shown in Tab. 7. The trim solutions obtained with inertia relief off/on are moderately

15
 IFASD-2017-170

different, which is justified by considering that the stiffness matrix is updated in the loop and it
is thus sensitive to the load path. The main difference is observed in the trim angle of attack,
whereas the obtained aeroelastostatic deflections practically overlap. Note that including an
inertia relief algorithm in the loop to simulate free-free boundary conditions at any iteration is
of concern only when using a displacement-based FEM solver (like standard commercial FEM
solvers) to compute structural deflection. Indeed, the model stiffness matrix is always definite
positive when the independent variables are element strains rather than nodal displacements.
For this reason, simulating unrestrained bodies is not critical in the case of a strain-based finite
element analysis.

4 CONCLUDING REMARKS
A novel approach has been presented to compute nonlinear aeroelastic trim of highly flexible
aircraft by coupling off-the-shelf solvers for high-fidelity modeling of structures and aerody-
namics. The main features of the proposed trim methodology are 1) the use of six degree-of-
freedom splines for loads and displacements interpolation between aerodynamic and structural
grids; 2) a load-stepping relaxation procedure to implement a close fluid-structural coupling;
and 3) a large-amplitude inertia relief algorithm to simulate free-free boundary conditions in the
loop. The proposed approach has been implemented for [Link] SOL 400 coupled with
a three-dimensional VLM code using six degree-of-freedom splines. The developed computa-
tional environment has been assessed by performing both nonlinear aeroelastic static analyses
of a 16-meter wing and of the University of Michigan’s X-HALE RRV at different prescribed
flight conditions and a complete nonlinear aeroelastic trim analysis of the latter configuration.
Beam-type FEM models have been considered in all the numerical studies in order to directly
compare the results with UM/NAST. The results have been also compared with a [Link]
SOL 400/VLM nonlinear aeroelastic toolbox developed at the DLR and, in the case of trim,
with linear predictions obtained via [Link] SOL 144. For the 16-meter wing, match-
ing results in terms of nonlinear aeroelastic static response have been achieved with all the
methodologies. For the X-HALE, a close agreement has been found between the present re-
sults and the DLR toolbox ones, while slight differences have been observed with UM/NAST.
A good agreement between nonlinear and linear trim predictions has been also pointed out due
to the moderate deflections experienced by the X-HALE RRV in the examined flight condition.
Solution sensitivity to the step size used in the relaxation procedure and to the inertia relief
algorithm have been analyzed, showing that computational time can be saved by conveniently
setting the number of fluid-structure load steps and by simulating free-free boundary conditions
at each iteration of the nonlinear aeroelastic trim loop. Future work will address application
the proposed approach to complex FEM models and further investigation on the influence of
six degree-of-freedom coupling, load-stepping relaxation, and large-amplitude inertia relief in
nonlinear aeroelastic analyses. Comparison with experimental results from the X-HALE RRV
flight tests will be carried out to further assess the methodology. The developed nonlinear aeroe-
lastic trim solver will be integrated in an existing FEM-based computational environment for
coupled flight dynamics and aeroelastic simulation to extend the toolbox for analysis of very
flexible vehicles.

ACKNOWLEDGMENTS
The authors would like to thank the PhD students Patricia Teixeira and Ryan Kitson of the
University of Michigan’s Active Aeroelasticity and Structures Research Laboratory (A2 SRL)
for providing the UM/NAST results.

16
 IFASD-2017-170

5 REFERENCES

[1] Noll, T.E., Brown, J.M., Perez-Davis, M.E., Ishmael, S.D., Tiffany, G.C., and Gaier,
M., “Investigation of the Helios Prototype Aircraft Mishap,” Volume I, January 2004,
NASA Mishap Report.

[2] Hodges, D. H., “Geometrically Exact, Intrinsic Theory for Dynamics of Curved and
Twisted Anisotropic Beams,” AIAA Journal, Vol. 41, No. 6, 2003, pp. 1131–1137.

[3] Su, W., and Cesnik, C. E. S., “Strain-based Geometrically Nonlinear Beam Formulation
for Modeling Very Flexible Aircraft,” International Journal of Solids and Structures,
Vol. 48, Nos. 16–17, 2011, pp. 2349–2360.

[4] Peters, D. A., and Johnson, M. J., “Finite-State Airloads for Deformable Airfoils on
Fixed and Rotating Wings,” Symposium on Aeroelasticity and Fluid/Structure Inter-
action, ASME Winter Annual Meeting, edited by P. P. Friedmann and J. C. I. Chang,
American Society of Mechanical Engineers, New York, Nov. 1994, pp. 128.

[5] Peters, D. A., Karunamoorthy, S., and Cao, W. M., “Finite State Induced Flow Mod-
els Part I: Two-Dimensional Thin Airfoil,” Journal of Aircraft, Vol. 32, No. 2, 1995,
pp. 313–322.

[6] Patil, M.J., and Hodges, D.H., “Flight Dynamics of Highly Flexible Flying Wings,”
Journal of Aircraft, Vol. 43, No. 6, 2006, pp. 1790–1799.

[7] Su, W., and Cesnik, C.E.S., “Nonlinear Aeroelasticity of a Very Flexible Blended-Wing-
Body Aircraft,” Journal of Aircraft, Vol. 47, No. 5, 2010, pp. 1539–1553.

[8] Rodden, W. P., Johnson, E. H., [Link] Aeroelastic Analysis User’s Guide,
[Link] Corporation, Santa Ana, CA, USA, 2010.

[9] Saltari, F., Riso, C., Mastroddi, F., and De Matteis, G., “Finite-Element Based Modeling
for Flight Dynamics and Aeroelasticity of Flexible Aircraft,” Journal of Aircraft, 2017,
in publication.

[10] Cavin, R. K. III, and Dusto, A. R., “Hamilton’s Principle: Finite-Element Methods and
Flexible Body Dynamics,” AIAA Journal Vol. 15, No. 12, 1977, pp. 1684–1690.

[11] Rodden, W. P., and Love, R. J., “Equations of Motion of a Quasisteady Flight Vehicle
Utilizing Restrained Static Aeroelastic Characteristics,” Journal of Aircraft, Vol. 22, No.
9, 1985, pp. 802–809.

[12] Anon., [Link] Linear Statics User’s Guide, [Link] Corporation, Santa
Ana, CA, USA, 2012.

[13] Anon., [Link] Nonlinear Users Guide SOL 400, [Link] Corporation,
Santa Ana, CA, USA, 2009.

[14] Katz, J., and Plotkin, A., Low-Speed Aerodynamics, Cambridge Aerospace Series, Cam-
bridge University Press, 2001.

17
 IFASD-2017-170

[15] Ritter, M., Cesnik, C. E. S., and Kruger, W. R., “An Enhanced Modal Approach
for Large Deformation Modeling of Wing-Like Structures,” AIAA SciTech, 56th
AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference,
Kissimmee, Florida, 5–9 January 2015.

[16] Ritter, M., Jones, J. R., Cesnik, C. E. S., “Enhanced Modal Approach for Free-Flight
Nonlinear Aeroelastic Simulation of Very Flexible Aircraft,” AIAA SciTech, 15th Dy-
namics Specialists Conference, San Diego, California, 4–8 January 2016.

[17] Anon., [Link] Quick Reference Guide Guide, [Link] Corporation, 2012.

[18] Di Vincenzo, F. G., and Castrichini, A., Hybrid Static Aeroelasticity Toolkit (HSA)
User’s Guide, [Link] Toulouse, Toulouse, France, 2013.

[19] Cesnik, C. E. S., Senatore, P. J., Su, W., Atkins, E. M., and Shearer, C. M., “X-HALE:
A Very Flexible UAV for Nonlinear Aeroelastic Tests,” AIAA Journal, Vol. 50, No. 12,
2012, pp. 2820–2833.

[20] Jones, J. R., and Cesnik, C. E. S., “Preliminary Flight Test Correlations of the X-HALE
Aeroelastic Experiment,” Aeronautical Journal of The Royal Aeronautical Society, Feb.
2015.

[21] Neto, A. B. G., Silva, R. G. A., Paglione P., and Silvestre F. J., “Formulation of the
Flight Dynamics of Flexible Aircraft Using General Body Axes,” AIAA Journal, Vol. 54,
No. 11, 2016, pp. 3516–3534.

COPYRIGHT STATEMENT

The authors confirm that they, and/or their company or organization, hold copyright on all of
the original material included in this paper. The authors also confirm that they have obtained
permission, from the copyright holder of any third party material included in this paper, to
publish it as part of their paper. The authors confirm that they give permission, or have obtained
permission from the copyright holder of this paper, for the publication and distribution of this
paper as part of the IFASD-2017 proceedings or as individual off-prints from the proceedings.

18