JOURNAL OF AIRCRAFT

Vol. 50, No. 5, September–October 2013

Nonlinear Modal Analysis of a Full-Scale Aircraft

G. Kerschen,∗ M. Peeters,† and J. C. Golinval∗
University of Liège, 4000 Liège, Belgium
and
C. Stéphan‡
ONERA–The French Aerospace Lab, 92320 Châtillon, France
DOI: 10.2514/1.C031918
Nonlinear normal modes, which are defined as a nonlinear extension of the concept of linear normal modes, are a
rigorous tool for nonlinear modal analysis. The objective of this paper is to demonstrate that the computation of
nonlinear normal modes and of their oscillation frequencies can now be achieved for complex, real-world aerospace
structures. The application considered in this study is the airframe of the Morane–Saulnier Paris aircraft. Ground
vibration tests of this aircraft exhibited softening nonlinearities in the connection between the external fuel tanks and
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

the wing tips. The nonlinear normal modes of this aircraft are computed from a reduced-order nonlinear finite
element model using a numerical algorithm combining shooting and pseudo-arclength continuation. Several
nonlinear normal modes, involving, e.g., wing bending, wing torsion, and tail bending, are presented, which highlights
that the aircraft can exhibit very interesting nonlinear phenomena. Specifically, it is shown that modes with distinct
linear frequencies can interact and generate additional nonlinear modes with no linear counterpart.

Nomenclature small aircraft for which no finite element (FE) model is developed,
F = augmented two-point boundary-value problem, m and m∕s GVT is thus the only means of providing a modal basis. This basis is
fnl = restoring force vector, N used to form a mathematical model of the structural dynamic
fr;c = restoring force at wing side, N behavior that is used for flutter predictions. For larger aircraft, the
g = vector field, m∕s and N∕kg main GVTobjective is to update and validate the finite element model
H = shooting function, m and m∕s of the aircraft, which, in turn, will be used for making reliable flutter
h = phase condition predictions.
I = identity matrix, 1∕s Aircraft GVT has now become standard industrial practice [1,2].
K = stiffness matrix, N∕m However, nonlinearity is frequently encountered during these tests
k = linear cubic spring, N∕m and represents a significant challenge for aerospace engineers [3,4].
k−nl = nonlinear cubic spring, N∕m3 For instance, besides the nonlinear fluid–structure interaction, typical
M = mass matrix, kg nonlinearities include backlash and friction in control surfaces and
mc = mass at wing side, kg joints, hardening nonlinearities in engine-to-pylon connections, and
n = number of degrees of freedom saturation effects in hydraulic actuators. Although there has been
p = tangent vector, m; m∕s; and s some recent progress in nonlinearity detection and characterization in
T = motion period, s aerospace structures [5,6], further developments for accounting for
T~ = prediction of the motion period, s nonlinear dynamic phenomena during aircraft certification are still
x = displacement vector, m needed [7]. They would clearly improve the confidence practicing
x_ = velocity vector, m∕s engineers have in both aircraft testing and modeling.
x
= acceleration vector, m∕s2 In this context, the nonlinear normal mode (NNM) theory offers a
x
c = acceleration at wing side, m∕s2 solid theoretical tool for interpreting a wide class of nonlinear
xrel = relative displacement across the nonlinear connection, m dynamical phenomena, yet NNMs have a clear conceptual relation to
y
o = linear mode shape of the full-size finite element model the linear normal modes (LNMs) [8,9]. However, most engineers still
y
r = linear mode shape of the reduced finite element model view NNMs as a concept that is foreign to them. One reason suppor-
z = state space vector, m and m∕s ting this statement is that most existing constructive techniques for
z~ = prediction of the next state space vector, m and m∕s computing NNMs are based on asymptotic approaches [10,11].
z0 = initial state space vector, m and m∕s Because these techniques rely on fairly involved mathematical
developments, they are usually limited to the analysis of systems with
low dimensionality.
I. Introduction There have been few attempts to compute NNMs using numerical
methods [12–15]. The objective of this paper is to demonstrate that
D URING the development of a new aircraft, testing plays a key
role for flutter qualification. Before flight flutter tests, ground
vibration testing (GVT) is performed on the full-scale aircraft. For
the computation of NNMs of complex, real-world structures is now
within reach and that it can be of great help for understanding and
interpreting experimental observations. Specifically, the algorithm
Received 10 April 2012; revision received 19 March 2013; accepted for proposed in [16] is used in the present study to compute the NNMs of
publication 21 March 2013; published online 29 August 2013. Copyright © the airframe of the Morane–Saulnier Paris aircraft, for which the
2013 by the American Institute of Aeronautics and Astronautics, Inc. All ground vibration tests have exhibited some nonlinear structural
rights reserved. Copies of this paper may be made for personal or internal use, behaviors. Very interesting nonlinear dynamics, including nonlinear
on condition that the copier pay the $10.00 per-copy fee to the Copyright modal interactions, will be revealed.
Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA 01923; include The paper is organized as follows. Section II briefly reviews the
the code 1542-3868/13 and $10.00 in correspondence with the CCC.
two main definitions of an NNM. Section III describes the algorithm
*Professor, Department of Aerospace and Mechanical Engineering, 1,
Chemin des Chevreuils, B-4000. used for NNM computation. Section IV introduces the mathematical
†
Ph.D. Student, Department of Aerospace and Mechanical Engineering, 1, modeling of the considered aircraft and discusses carefully the results
Chemin des Chevreuils, B-4000. of nonlinear modal analysis. The conclusions of the present study are
‡
Research Scientist, 29 Avenue de la Division Leclerc, F-92320. summarized in Sec. V.
1409
 1410 KERSCHEN ET AL.

II. Brief Review of Nonlinear Normal Modes to the minimal period of the periodic motion and at an energy equal to
A detailed description of NNMs and of their fundamental the conserved total energy during the motion (i.e., the sum of kinetic
properties (e.g., frequency–energy dependence, bifurcations, and and potential energies). A branch, represented by a solid line, is a
stability) is given in [8,9]. For completeness, the two main definitions family of NNM motions possessing the same qualitative features.
of an NNM are briefly reviewed in this section. The backbone of the plot is therefore formed by two branches, which
The free response of discrete conservative mechanical systems represent in-phase (S11) and out-of-phase (S11−) NNMs. The FEP
with n degrees of freedom (DOFs) is considered, assuming that clearly shows that the nonlinear modal parameters have a strong
continuous systems (e.g., beams, shells, or plates) have been spatially dependence on the total energy in the system based on the following:
discretized using the finite element method. The governing equations 1) The frequency of both the in-phase and out-of-phase NNMs
of motion are increases with the energy level, which reveals the hardening
characteristic of the system.
Mx
t

 Kx
t  fnl fx
tg  0 (1) 2) The modal curves change for increasing energies. The in-phase
NNM tends to localize to the second DOF (i.e., it resembles a vertical
where M and K are the mass and stiffness matrices, respectively, x
t curve), whereas the out-of-phase NNM localizes to the first DOF
and x
t

are the displacement and acceleration vectors, respectively, (i.e., it resembles a horizontal curve).
and fnl is the nonlinear restoring force vector.
There exist two main definitions of an NNM in the literature due to
Rosenberg [17] and Shaw and Pierre [18]:
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

1) Targeting a straightforward nonlinear extension of the LNM

concept, Rosenberg defined an NNM motion as a vibration in unison
of the system (i.e., a synchronous periodic oscillation).
2) To provide an extension of the NNM concept to damped
systems, Shaw and Pierre defined an NNM as a two-dimensional
(2-D) invariant manifold in phase space. Such a manifold is invariant
under the flow (i.e., orbits that start out in the manifold remain in it for
all time), which generalizes the invariance property of LNMs to
nonlinear systems.
At first glance, Rosenberg’s definition may appear restrictive in
two cases. First, it cannot be easily extended to nonconservative
systems. However, the damped dynamics can often be interpreted
based on the topological structure of the NNMs of the underlying
conservative system [9]. Second, in the presence of internal
resonances, the NNM motion is no longer synchronous, but it is still
periodic.
In the present study, an NNM motion is therefore defined as a
(nonnecessarily synchronous) periodic motion of the conservative
mechanical system (1). As we will show, this extended definition is
particularly attractive when targeting a numerical computation of the
NNMs. It enables the nonlinear modes to be effectively computed
using algorithms for the continuation of periodic solutions.
For illustration, the NNMs of the two-DOF system

x
1 
2x1 − x2   0.5x31  0 x
2 
2x2 − x1   0 (2)

are depicted in Fig. 1. Because of the frequency–energy dependence

of nonlinear oscillations, an appropriate graphical depiction is to
represent NNMs in a frequency–energy plot (FEP). An NNM is
represented by a point, which is drawn at a frequency corresponding

Fig. 2 Algorithm for NNM computation.

Fig. 1 Frequency–energy plot of the two-DOF system. Fig. 3 Morane–Saulnier Paris aircraft.
 KERSCHEN ET AL. 1411

Table 1 Aircraft properties Starting from some assumed initial conditions z
0 p0 , the motion
z
0
0
p
t; zp0  at the assumed period T

0 can be obtained by numerical
Length, m Wingspan, m Height, m Wing area, m2 Weight, kg
time integration methods (e.g., Runge–Kutta or Newmark schemes).
10.4 10.1 2.6 18 1945
A Newton–Raphson iteration scheme is therefore used to correct the
initial guess and to converge to the actual solution.
The phase of the periodic solutions is not fixed. If z
t is a solution
Such a plot gives a very clear indication of the effects of of the autonomous system (3), then z
t  Δt is geometrically the
nonlinearity on modal parameters. same solution in state space for any Δt. Hence, an additional
condition, termed the phase condition, has to be specified in order to
remove the arbitrariness of the initial conditions. An isolated NNM is
III. Numerical Computation of NNMs therefore computed by solving the augmented two-point boundary-
The numerical method proposed here for the NNM computation value problem defined by
relies on two main techniques, namely, shooting and pseudo-
arclength continuation. A detailed description of the algorithm is 
given in [16]. H
zp0 ; T  0
F
zp0 ; T ≡ (6)
h
zp0   0
A. Shooting Method
The equations of motion of system (1) can be recast into state space
where h
zp0   0 is the phase condition.
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

form:

z_  g
z (3) B. Continuation of Periodic Solutions

Because of the frequency–energy dependence, the modal
where z  x x_   is the 2n-dimensional state vector and star parameters of an NNM vary with the total energy. An NNM family,
denotes the transpose operation, and governed by Eqs. (6), therefore traces a curve, termed an NNM

 branch, in the
2n  1-dimensional space of initial conditions and
x_ period
zp0 ; T. Starting from the corresponding LNM at low energy,
g
z  (4)
−M−1 Kx  fnl
x; x
_ the computation is carried out by finding successive points
zp0 ; T of
the NNM branch using methods for the numerical continuation of
is the vector field. The solution of this dynamical system for initial periodic motions [19,20]. The so-called pseudo-arclength continu-
conditions z
0  z0  x0 x_ 0  is written as z
t  z
t; z0  in order ation method is used herein.
to exhibit the dependence on the initial conditions, z
0; z0   z0 . A Starting from a known solution
zp0;
j ; T
j , the next periodic
solution zp
t; zp0  is a periodic solution of the autonomous system solution
zp0;
j1 ; T
j1  on the branch is computed using a
(3) if zp
t; zp0   zp
t  T; zp0 , where T is the minimal period. predictor step and a corrector step.
The NNM computation is carried out by finding the periodic
solutions of the governing nonlinear equations of motion (3). In this 1. Predictor Step
context, the shooting method is probably the most popular numerical At step j, a prediction
z~ p0;
j1 ; T~
j1  of the next solution
technique. It solves numerically the two-point boundary-value
zp0;
j1 ; T
j1  is generated along the tangent vector to the branch
problem defined by the periodicity condition at the current point zp0;
j :
H
zp0 ; T ≡ zp
T; zp0  − zp0  0 (5)      
z~ p0;
j1 zp0;
j pz;
j
  s
j (7)
H
z0 ; T  z
T; z0  − z0 is called the shooting function and T~
j1 T
j pT;
j
represents the difference between the initial conditions and the
system response at time T.
The shooting method consists of finding, in an iterative way, the where s
j is the predictor stepsize. The tangent vector p
j 
initial conditions zp0 and the period T that realize a periodic motion. pz;
j pT;
j  to the branch defined by (6) is solution of the system

Fig. 4 a) Connection between external fuel tank and wing tip (top view). Close-up of b) front- and c) rear-bolted attachments.
 1412 KERSCHEN ET AL.

Fig. 5 Finite element model of the Morane–Saulnier Paris aircraft.

2   3
∂H  ∂H    
5 pz;
j  0
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

∂zp0
z ∂T
zp0;
j ;T
j 

4  p0;
j ;T
j 
∂h   (8)
∂zp0 0 pT;
j 0

zp0;
j 

with the condition kp
j k  1.

2. Corrector Step
The prediction is corrected by a shooting procedure in order to
solve Eq. (6) in which the variations of the initial conditions and the
period are forced to be orthogonal to the predictor step. At iteration k,
the corrections z
k1
k
k
p0;
j1  zp0;
j1  Δzp0;
j1 and T
j1 

k1

k
k
T
j1  ΔT
j1 are computed by solving the overdetermined linear
system using the Moore–Penrose matrix inverse:
2   3
∂H  ∂H
∂zp0
z
k ;T
k
 ∂T
z
k
;T
k
 2 3
6 p0;
j1
j1 7
k
7 Δzp0;
j1
p0;
j1
j1
6 
6 ∂h   0 74 5
6 ∂zp0
z
k  7
4 p0;
j1 5 ΔT
k

j1
pz;
j pT;
j
2 3
−H
z
k
k
p0;
j1 ; T
j1 
6 7
6 −h
z
k 7
4 p0;
j1  5
0 Fig. 6 a) First wing bending mode, b) first wing torsional mode, and
c) second wing torsional mode.
Table 2 Natural frequencies of the linear finite element model
Mode Frequency, Modal shape
number Hz
1 0.0936 Rigid-body mode
2 0.7260 Rigid-body mode
3 0.9606 Rigid-body mode
Wing side

Tank side

4 1.2118 Rigid-body mode

5 1.2153 Rigid-body mode
6 1.7951 Rigid-body mode
7 2.1072 Rigid-body mode
8 2.5157 Rigid-body mode
9 3.5736 Rigid-body mode
10 8.1913 Two-node wing bending
11 9.8644 Fin bending Fig. 7 Instrumentation of the rear attachment of the right wing.
12 16.1790 Fin torsion
13 21.2193 First T-tail symmetric bending where the prediction is used as initial guess. This iterative process is
14 22.7619 Front fuselage torsion carried out until convergence is achieved. The convergence test is
15 23.6525 T-tail torsion based on the relative error of the periodicity condition, i.e.,
16 25.8667 Antielevator torsion
17 28.2679 Two-node vertical fuselage bending
kH
zp0 ; Tk∕kzp0 k < ϵ where ϵ is usually set to 10−6 .
18 29.3309 Three-mode wing bending An important technical remark is that the method requires the
19 31.0847 Symmetric wing torsion evaluation of the 2n × 2n Jacobian matrix:
20 34.9151 Antisymmetric wing torsion 
21 39.5169 x (fore and aft)-symmetric wing bending ∂H ∂z
t; z0 

z ; T  −I (10)
22 40.8516 Antisymmetric wing torsion  front ∂z0 0 ∂z0 tT
fuselage
23 47.3547 Three-node vertical fuselage bending
24 52.1404 Three-node T-tail bending where I is the 2n × 2n identity matrix. The classical finite-difference
approach requires us to perturb successively each of the 2n initial
 KERSCHEN ET AL. 1413
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

Fig. 8 Left: restoring force surface, right: stiffness curves (zero relative velocity). Top: rear tank of the left wing, bottom: rear tank of the right wing.

21.3

21.28

21.26

21.24
Frequency (Hz)

21.22
a b
21.2

21.18

21.16

21.14

21.12

21.1 −4 −3 −2 −1 0 1
10 10 10 10 10 10
Energy (log scale)
a)

b) c)
Fig. 9 FEP of an NNM involving T-tail symmetric bending.
 1414 KERSCHEN ET AL.

conditions and integrate the nonlinear governing equations of liaison aircraft. The structural configuration under consideration
motion. This is extremely computational intensive for large-scale corresponds to the aircraft without its jet engines and standing on the
finite element models such as the one considered in this study. ground through its three landing gears with deflated tires. General
Targeting a substantial reduction of the computational cost, characteristics of the aircraft are listed in Table 1. A specimen of this
sensitivity analysis is exploited for determining ∂z
t; z0 ∕∂z0 . It airplane is present in ONERA’s laboratory. It is used for training
amounts to differentiating the equations of motion (3) with respect to engineers and technicians for GVT and for research purposes.
the initial conditions z0 , which leads to Ground vibration tests exhibited nonlinear behavior in the connection
     between the wings and external fuel tanks located at the wing tip.
d ∂z
t; z0  ∂g
z ∂z
t; z0  Figure 4 shows that this connection consists of bolted attachments.
 (11)
dt ∂z0 ∂z z
t;z0  ∂z0
A. Aircraft Structural Model
with 1. Finite Element Model of the Underlying Linear Structure
The FE model of the linear part of the full-scale aircraft, illustrated
∂z
0; z0  in Fig. 5, was elaborated from drawings at ONERA. The wings, T tail,
I (12)
∂z0 and fuselage are modeled by means of 2-D elements such as beams
and shells. Three-dimensional (3-D) spring elements, which take into
because z
0; z0   z0 . Hence, the matrix ∂z
t; z0 ∕∂z0 at t  T can account the structural flexibility of the tires and landing gears, are
be obtained by numerically integrating over T the initial-value used as boundary conditions of the aircraft. At each wing tip, the front
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

problem defined by linear ordinary differential equations (11) with and rear connections between the wings and the fuel tanks are
initial conditions (12). modeled using beam elements. The FE model, originally created in
The complete algorithm is shown in Fig. 2. NASTRAN, was converted into the SAMCEF FE software. The
complete FE model has more than 80,000 DOFs.
The natural frequencies of the underlying linear system in the
IV. Nonlinear Modal Analysis of the Morane–Saulnier [0–50 Hz] frequency range are given in Table 2. The first nine modes
Paris Aircraft correspond to aircraft rigid-body modes, i.e., six modes are landing
The numerical computation of the NNMs of a complex, real-world gear suspension modes, whereas three modes involve rigid-body
structure is addressed in this section. The structure is the airframe of motions of the control surfaces. The frequency range of the
the Morane–Saulnier Paris aircraft, represented in Fig. 3. This French rigid-body modes is comprised between 0.09 and 3.57 Hz, i.e.,
jet aircraft was built during the 1950s and was used as a trainer and noticeably lower than the first elastic mode located at 8.19 Hz. The

8.3

8.25

8.2
Frequency (Hz)

a
8.15 b

8.1

8.05

8
−4 −3 −2 −1 0
10 10 10 10 10
Energy (log scale)
a)

b) c)
Fig. 10 FEP of an NNM involving wing bending (linear frequency
8.19 Hz). NNMs at energy levels marked in the FEP are inset; they are given in terms
of the initial displacements (m) that realize the periodic motion (with zero initial velocities assumed).
 KERSCHEN ET AL. 1415

modal shapes of different elastic normal modes of vibrations are is performed in the space of the initial model after projecting the
depicted in Fig. 6. Figure 6a represents the first wing bending mode. reduced modes back into the original space. The deviation between
The first and second wing torsional modes are depicted in Figs. 6b the mode shapes of the original y
o and reduced y
r models is
and 6c and correspond to symmetric and antisymmetric wing determined using the modal assurance criterion (MAC):
motions, respectively. These latter modes are of particular interest
because there is a significant deformation of the nonlinear connec- jy 
o y
r j2
tions between the wings and fuel tanks. The other modes concern the MAC  (13)
jy
o y
o jjy 
r y
r j

aircraft tail and are consequently almost unaffected by the nonlinear
connections.
MAC values range from 0 in the absence of correlation to 1 for a
2. Reduced-Order Modeling complete correspondence. In the [0–100 Hz] range, MAC values
Because the nonlinearities are spatially localized, condensation between modes shapes are all greater than 0.999, and the maximum
of the linear components of the model can be achieved using the relative error on the natural frequencies is 0.2%. The accuracy of this
Craig–Bampton reduction technique [21]. This will lead to a linear reduced model is therefore excellent. We note that much less
substantial decrease in the computational burden without degrading than 500 internal modes are sufficient to build a good reduced model
the computational accuracy, at least in the frequency range of interest. in the [0–100 Hz] frequency range. However, a larger number of
We stress that many practical applications possess spatially localized modes was deliberately chosen for two main reasons. On the one
nonlinearities, and so the procedure developed herein is of wide hand, it serves to illustrate that our NNM algorithm can deal with
systems of relatively high dimensionality. On the other hand, due to
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

applicability.
The Craig–Bampton method expresses the complete set of initial the presence of harmonics, nonlinear modal interactions may occur
DOFs in terms of retained DOFs and internal vibration modes of the between a mode in the frequency range of interest and another mode
primary structure clamped on the retained nodes. To introduce the outside this range. More internal modes are therefore necessary to
nonlinear behavior of the connections between the wings and the guarantee the accuracy of the reduced model in nonlinear regimes of
tanks, the reduced-order model of the aircraft is constructed by motion.
keeping one node on both sides of the attachments. For each wing,
four nodes are retained, namely, two nodes for the front attachment 3. Aircraft Nonlinearities
and two nodes for the rear attachment. In total, eight nodes of the The existence of a softening nonlinear behavior was evidenced
initial FE model possessing, each 6 DOFs and 500 internal modes of during ground vibration tests conducted at ONERA. In particular,
vibrations, are kept in the reduced model; the FE model is thus frequency response function measurements revealed a decrease in
reduced to 548 DOFs. resonant frequencies when the excitation level was increased. The
To assess the accuracy of the reduced model, its LNMs are connections between the wings and fuel tanks were assumed to cause
compared to those predicted by the initial FE model. The comparison this observed nonlinear effect. To confirm this hypothesis, the front

31.2 31.0360

31.1

31 a
tongue 3:1
Frequency (Hz)

30.9
31.0358
0
10
30.8

30.7 30.6

30.6
tongue 5:1
30.5
b 30.58
tongue 9:1
30.4
−4 −2 0 2 4 0
10 10 10 10 10 10
Energy (log scale)
a)

b) c)
Fig. 11 FEP of an NNM involving symmetric wing torsion (linear frequency
31.08 Hz). NNMs at energy levels marked in the FEP are inset; they are
given in terms of the initial displacements (m) that realize the periodic motion (with zero initial velocities assumed).
 1416 KERSCHEN ET AL.

and rear connections of each wing were instrumented, and measure- deformation of the wing, a quantitative assessment of the nonlinear
ments dedicated to nonlinearity characterization were carried out. effects cannot be achieved precisely. However, plotting the measured
Specifically, accelerometers were positioned on both the wing and acceleration signal x
c
t against the relative displacement and
tank sides of the connections, and two shakers were located at the velocity across the connection can provide meaningful qualitative
tanks. The instrumentation at one of the connection is shown in Fig. 7. information about the type of nonlinearity present in the bolted
The dynamic behavior of these connections in the vertical direction connections. To this end, the aircraft was excited close to the second
was investigated using the restoring force surface method [22]. torsional mode, which is known to activate nonlinear behavior, using
Newton’s second law applied at the wing side of one connection a band-limited swept sine excitation. Figure 8 presents the resulting
writes 3-D plots for the rear connections of the left and right wings. 2-D
sections corresponding to zero relative velocities are also shown to
mc x
c
t  fr;c fxc
t; x_c
tg  0 (14) highlight elastic nonlinearities. A softening elastic nonlinearity with
a piecewise linear characteristic can clearly be observed. This result
where fr;c is the restoring force at the considered DOF. The index c is seems to be compatible with what was previously reported for bolted
related to the connection under consideration (i.e., either the rear or connections in the literature [23,24]. Similar nonlinear effects were
front attachment of the left or right wing). From Eq. (14), the restoring also observed for the front connections. They are not shown here,
force is obtained by because these connections participated in the aircraft response to a
lesser extent.
fr;c  −mc x
c (15) Because nonsmooth nonlinearities require dedicated and very
specific time integration methods, which are not yet available in our
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

Because the restoring force fr;c includes both the force related to the NNM computation algorithm, the nonlinearities in the connections
connection of interest and the force generated by the elastic were regularized using polynomial expansions. This is why the

31.036
a
Frequency (Hz)

31.036
b c
31.0359

31.0358

31.0358
0 2
10 10
Energy
a)

−4 −4
x 10 x 10
4 5 0.01
Displacement (m)

Displacement (m)
Displacement (m)

2 0.005

0 0 0

−2 −0.005

−4 −5 −0.01
0 0.01 0.02 0.03 0 0.01 0.02 0.03 0 0.01 0.02 0.03
Time (s) Time (s) Time (s)
−4 −4 −3
x 10 x 10 x 10
4 6 10

4
2 5
Ai

Ai

Ai

2
0 0
0

−2 −2 −5
1 3 5 1 3 5 1 3 5
Harmonic i Harmonic i Harmonic i
b) c) d)
Fig. 12 3∶1 internal resonance between the first wing torsional mode and a higher tail mode. a) Close up in the FEP of the 3∶1 tongue of Fig. 11. Bottom
plots: NNM motions at b) beginning of the tongue (in the vicinity of the backbone of the first wing torsional mode), c) middle of the tongue, and d) extremity
of the tongue. From top to bottom: NNM shapes, time series of the vertical displacements at the rear tip of the right tank ( — — ) and at the right side of the
horizontal tail (–––), and Fourier coefficients of both displacements (in gray and black, respectively).
 KERSCHEN ET AL. 1417

stiffness curve was reconstructed by considering the mathematical model is considered, this mode is not affected by nonlinearity and can
model be considered as a purely linear mode.
The nonlinear extension of the two-node wing-bending mode (i.e.,
fr;c  kxrel  k−nl x3rel
k−nl < 0 (16) mode 10 in Table 2) is illustrated in Fig. 10. The FEP reveals that this
mode is weakly affected by the nonlinearities, at least until energies
A value of −1013 N∕m3 was found to be adequate for the cubic compatible with physical constraints. The frequency of the NNM
spring. Damping in the connections was not considered because the motions on the backbone slightly decreases with increasing energy
focus of this study is on the NNMs of the underlying Hamiltonian levels, which results from the softening characteristic of the
system. nonlinearity. A quantitatively similar decrease in frequency between
The reduced-order model, generated in the SAMCEF software, low-level and high-level swept sine excitations was observed during
was exported to MATLAB® where nonlinearities were imple- the experimental tests. Even though the frequency decrease is small,
mented. Four nonlinearities were added, each one being located we remark that an accuracy of 0.001 Hz is usually sought during
between the two nodes defining a nonlinear connection. GVTs; it is therefore important to account for it. The MAC value
between the NNMs at low and high energy levels (see Figs. 10b and
B. Nonlinear Normal Modes and Corresponding Oscillation 10c) is 0.99. Accordingly, the modal shapes do not change much over
Frequencies the considered energy range and resemble the corresponding LNM.
The NNM computation was carried out in the MATLAB® Figure 11 represents the FEP of the first symmetric wing torsional
environment using the nonlinear reduced-order model built in the mode (i.e., mode 19 in Table 2). For this mode, the relative motion of
the fuel tanks is more important, which enhances the nonlinear effect
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

preceding sections. For conciseness, only the modes that are the most
relevant ones for the purpose of this study are described in what of the connections. As a result, the oscillation frequency has a marked
follows. energy dependence along the backbone branch. This decay of the
The first mode examined herein is the nonlinear counterpart of the natural frequency is in agreement, at least qualitatively, with what
first T-tail symmetric bending mode (i.e., mode 13 in Table 2). As in observed experimentally and confirms the relevance of this nonlinear
Sec. II, the computed backbone branch and related NNM motions are modal analysis. Conversely, the modal shapes on the backbone
depicted in a FEP in Fig. 9. Because there is no visible wing branch are only weakly altered by the nonlinearities; the MAC value
deformation, the modal shape and the corresponding oscillation between the NNMs at low and high energy levels is equal to 0.98.
frequency remain practically unchanged when the total energy is In addition to the main backbone branch, three other NNM
increased in the system. Despite the fact that a nonlinear aircraft branches that are localized to a specific region of the FEP can be

35

33.5344
34.8 a
33.5343

33.5342
34.6
33.5341 tongue 9:1
33.534
Frequency (Hz)

34.4
33.5339
0 2
34.2 10 10

34

33.8
b

33.6

33.4
−4 −2 0 2 4
10 10 10 10 10
Energy (log scale)
a)

b) c)
Fig. 13 FEP of an NNM involving antisymmetric wing torsion. NNMs at energy levels marked in the FEP are inset; they are given in terms of the initial
displacements (m) that realize the periodic motion (with zero initial velocities assumed).
 1418 KERSCHEN ET AL.

noticed. These tongues bifurcate from the backbone branch of the frequency region in which they live. This will be addressed in
considered mode and bifurcate into the backbone branch of another subsequent research by applying, numerically and experimentally, a
mode, thereby realizing an internal resonance between two modes. slow sweep sine excitation and nonlinear force appropriation to this
Because of harmonics, modal interactions with distinct linear aircraft.
frequencies can be generated in nonlinear regimes of motion. Finally, this research paves the way for the constructive use of
Specifically, the 3∶1, 5∶1, and 9∶1 internal resonances in Fig. 11 nonlinearity for design. Some nonlinear phenomena with no linear
involve third, fifth, and ninth harmonics of the fundamental counterpart such as targeted energy transfer could be exploited for
frequency of the first symmetric wing torsional mode. protecting more sensitive parts of aircraft by transferring or
A close-up of the 3∶1 internal resonance is depicted in Fig. 12. redistributing energy from one substructure to another substructure.
Modal shapes at three different locations on the tongue are also
represented. When the energy gradually increases along the tongue, a
smooth transition from the first wing torsional mode to a higher-
frequency tail torsional mode occurs. Both modes are present in the References
linear modal basis, but their mixing realized by the NNM in Fig. 12c [1] Goge, D., Boswald, M., Fullekrug, U., and Lubrina, P., “Ground
is not. It is therefore a mode with no linear counterpart. These tongues Vibration Testing of Large Aircraft – State-of-the-Art and Future
Perspectives,” Proceedings of the 25th International Modal Analysis
of internal resonance represent an inherently nonlinear dynamical
Conference, Orlando, FL, Society of Experimental Mechanics, 2007.
phenomenon. Focusing now on the frequency content along the [2] Peeters, B., Hendricx, W., Debille, J., and Climent, H., “Modern
branch, a third harmonic progressively appears for increasing Solutions for Ground Vibration Testing of Large Aircraft,” Sound and
energies. In the vicinity of the NNM in Fig. 12c, the aircraft Vibration, Vol. 43, 2009, pp. 8–15.
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918

undergoes a subharmonic motion with two dominant frequencies. [3] Rodriguez Ahlquist, J., Martinez Carreno, J., Climent, H., de Diego, R.,
The wing vibrates at a frequency around 31 Hz, whereas the tail and de Alba, J., “Assessment of Nonlinear Structural Response in
vibrates at a frequency close to 93 Hz. Beyond this point, the relative A400M GVT,” Proceedings of the 28th International Modal Analysis
importance of the third harmonic component grows, until a complete Conference, Jacksonville, FL, Society of Experimental Mechanics,
coalescence with the tail torsional mode at 93 Hz is realized. 2010.
[4] Kehoe, M. W., and Freudinger, L. C., “Aircraft Ground Vibration
A completely similar dynamical mechanism exists for the 5∶1 and
Testing at the NASA Dryden Flight Research Facility,” NASA TM-
9∶1 internal resonances, but the wing torsional mode interacts with 104275, 1993.
other modes of the linear basis at 153 and 275 Hz, respectively. These [5] Gloth, G., and Sinapius, M., “Influence and Characterisation of Weak
modes correspond to local modes of some aircraft parts and are not Non-Linearities in Swept-Sine Modal Testing,” Aerospace Science and
shown here. We note that such nonlinear modal interactions were also Technology, Vol. 8, No. 2, 2004, pp. 111–120.
observed in [9] for the two-DOF system (2). However, they appeared [6] Gloth, G., and Goge, D., “Handling of Non-Linear Structural
at very high energies, and this is why they were not present in Fig. 1. Characteristics in Ground Vibration Testing,” Proceedings of the
Finally, the second (antisymmetric) wing torsional mode (i.e., International Conference on Noise and Vibration Engineering (ISMA),
mode 20 in Table 2) is plotted in Fig. 13. The oscillation frequency Leuven, Belgium, Katholieke Universiteit Leuven, 2004.
[7] Luber, W., “Application of Analytical Methods for Advanced Flight
undergoes a decrease around 5% in the considered energy range. The
Flutter Test,” Proceedings of the 30th International Modal Analysis
MAC value between the modal shapes at low and high energy levels Conference, Jacksonville, FL, Society of Experimental Mechanics,
is 0.97. The modal shapes are therefore affected by nonlinearity but 2012.
not to a large extent. As for the first wing torsion mode, one tongue of [8] Vakakis, A. F., Manevitch, L. I., Mikhlin, Y. V., Pilipchuk, V. N., and
internal resonance is present in the FEP and indicates that nonlinear Zevin, A. A., Normal Modes and Localization in Nonlinear Systems,
modal interactions are in fact generic in such a large-scale structure. Wiley, New York, 1996.
The computation of the backbone branch up to the tongue in [9] Kerschen, G., Peeters, M., Golinval, J. C., and Vakakis, A. F.,
Fig. 13 requires 20 min of CPU time with a single 2.67 GHz “Nonlinear Normal Modes, Part I: A Useful Framework for the
processor. Because of the presence of bifurcations, the computation Structural Dynamicist,” Mechanical Systems and Signal Processing,
Vol. 23, No. 1, 2009, pp. 170–194.
of the branch of internal resonance demands about 1 h. This shows
doi:10.1016/j.ymssp.2008.04.002
that the computation of NNMs of real-world structures is now within [10] King, M. E., and Vakakis, A. F., “An Energy-Based Formulation for
reach using a conventional laptop computer. Computing Nonlinear Normal Modes in Undamped Continuous
Systems,” Journal of Vibration and Acoustics, Vol. 116, No. 3, 1994,
pp. 332–340.
V. Conclusions doi:10.1115/1.2930433
[11] Touzé, C., and Amabili, M., “Nonlinear Normal Modes for Damped
The objective of this paper was to compute the nonlinear normal Geometrically Nonlinear Systems: Application to Reduced-Order
modes (NNMs) of a full-scale aircraft together with the corres- Modelling of Harmonically Forced Structures,” Journal of Sound and
ponding oscillation frequencies. NNMs were represented in a Vibration, Vol. 298, Nos. 4–5, 2006, pp. 958–981.
frequency–energy plot, in which the impact that nonlinearities have doi:10.1016/j.jsv.2006.06.032
on modal parameters can readily be seen. [12] Slater, J. C., “A Numerical Method for Determining Nonlinear Normal
Parts of aircraft, such as landing gears, connections between the Modes,” Nonlinear Dynamics, Vol. 10, No. 1, 1996, pp. 19–30.
horizontal tail plane and vertical tail plane, and pylons, are known to doi:10.1007/BF00114796
be nonlinear. When these components are excited during ground [13] Pesheck, E., Pierre, C., and Shaw, S. W., “A New Galerkin-Based
Approach for Accurate Non-Linear Normal Modes Through Invariant
vibration testing, complex nonlinear phenomena that cannot be
Manifolds,” Journal of Sound and Vibration, Vol. 249, No. 5, 2002,
described using linear tools may be observed. Therefore, the first pp. 971–993.
practical implementation of this work is to provide a rigorous doi:10.1006/jsvi.2001.3914
framework for the understanding and interpretation of these [14] Arquier, R., Bellizzi, S., Bouc, R., and Cochelin, B., “Two Methods for
phenomena. For instance, the present study showed that only a few the Computation of Nonlinear Modes of Vibrating Systems at Large
aircraft modes were influenced by the nonlinear bolted connections. Amplitudes,” Computers and Structures, Vol. 84, Nos. 24–25, 2006,
Their frequencies were found to decrease by a few percent in the pp. 1565–1576.
considered energy range, whereas their mode shapes were practically doi:10.1016/j.compstruc.2006.01.011
unaffected. All this information is relevant for the practicing [15] Casini, P., and Vestroni, F., “Characterization of Bifurcating Non-
aerospace engineer. Because of mode bifurcations, supernumerary Linear Normal Modes in Piecewise-Linear Mechanical Systems,”
International Journal of Non-Linear Mechanics, Vol. 46, No. 1, 2011,
nonlinear modes with no linear counterpart (i.e., that are not the direct pp. 142–150.
extension of linear normal modes) were also reported in this study. To doi:10.1016/j.ijnonlinmec.2010.08.002
the authors’ knowledge, this is the first time that such modes are [16] Peeters, M., Viguié, R., Sérandour, G., Kerschen, G., and Golinval, J. C.,
evidenced for a real-life structure. The practical realization of these “Nonlinear Normal Modes, Part II: Toward a Practical Computation
nonlinear resonances is still an open question in view of the small Using Numerical Continuation Techniques,” Mechanical Systems and
 KERSCHEN ET AL. 1419

Signal Processing, Vol. 23, No. 1, 2009, pp. 195–216. [21] Craig, R., and Bampton, M., “Coupling of Substructures for Dynamic
doi:10.1016/j.ymssp.2008.04.003 Analysis,” AIAA Journal, Vol. 6, No. 7, 1968, pp. 1313–1319.
[17] Rosenberg, R. M., “Normal Modes of Nonlinear Dual-Mode Systems,” [22] Masri, S. F., and Caughey, T. K., “A Nonparametric Identification
Journal of Applied Mechanics, Vol. 27, No. 2, 1960, pp. 263–268. Technique for Nonlinear Dynamic Problems,” Journal of Applied
doi:10.1115/1.3643948 Mechanics, Vol. 46, No. 2, 1979, pp. 433–447.
[18] Shaw, S. W., and Pierre, C., “Non-Linear Normal Modes and Invariant doi:10.1115/1.3424568
Manifolds,” Journal of Sound and Vibration, Vol. 150, No. 1, 1991, [23] Gaul, L., and Lenz, J., “Nonlinear Dynamics of Structures Assembled
pp. 170–173. by Bolted Joints,” Acta Mechanica, Vol. 125, Nos. 1–4, 1997, pp. 169–
doi:10.1016/0022-460X(91)90412-D 181.
[19] Seydel, R., Practical Bifurcation and Stability Analysis, from doi:10.1007/BF01177306
Equilibirum to Chaos, 2nd ed., Springer–Verlag, New York, [24] Hartwigsen, C. J., Song, Y., McFarland, D. M., Bergman, L. A., and
1994. Vakakis, A. F., “Experimental Study of Nonlinear Effects in a Typical
[20] Nayfeh, A. H., and Balachandran, B., Applied Nonlinear Dynamics: Shear Lap Joint Configuration,” Journal of Sound and Vibration,
Analytical, Computational, and Experimental Methods, Interscience, Vol. 277, Nos. 1–2, 2004, pp. 327–351.
New York, 1995. doi:10.1016/j.jsv.2003.09.018
Downloaded by MONASH UNIVERSITY on November 26, 2013 | http://arc.aiaa.org | DOI: 10.2514/1.C031918