# Nonlinear Modal Analysis of a Full-Scale Aircraft

Gaëtan Kerschen, Maxime Peeters, Jean Claude Golinval, Cyrille Stéphan

► To cite this version:

Gaëtan Kerschen, Maxime Peeters, Jean Claude Golinval, Cyrille Stéphan. Nonlinear Modal Analysis of a Full-Scale Aircraft. Journal of Aircraft, 2013, 50 (5), pp.1409-1419. 10.2514/1.C031918. hal-01389708

HAL Id: hal-01389708

https://hal.science/hal-01389708

Submitted on 29 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 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

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 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

$\mathbf{F}$ = augmented two point boundary value problem, m and m/s

$\mathbf{f}_{nl}$ = restoring force vector, N

$f_{r,c}$ = restoring force at wing side, N

$\mathbf{g}$ = vector field, m/s and N/kg

$\mathbf{H}$ = shooting function, m and m/s

$h$ = phase condition

$I$ = identity matrix, 1/s

$K$ = stiffness matrix, N/m

$k$ = linear cubic spring, N/m

$k_{nl}^{-}$ = nonlinear cubic spring, N/m³

$\mathbf{M}$ = mass matrix, kg

$m_{c}$ = mass at wing side, kg

$n$ = number of degrees of freedom

$\mathbf{p}$ = tangent vector, m; m/s; and s

$\mathbf{T}$ = motion period, s

$\tilde{\mathbf{T}}$ = prediction of the motion period, s

$\mathbf{x}$ = displacement vector, m

$\tilde{\mathbf{x}}$ = velocity vector, m/s

$\tilde{\mathbf{x}}$ = acceleration vector, m/s²

$\tilde{x}_{c}$ = acceleration at wing side, m/s²

$x_{rel}$ = relative displacement across the nonlinear connection, m

$\mathbf{y}_{(o)}$ = linear mode shape of the full size finite element model

$\mathbf{y}_{(r)}$ = linear mode shape of the reduced finite element model

$\mathbf{z}$ = state space vector, m and m/s

$\tilde{\mathbf{z}}$ = prediction of the next state space vector, m and m/s

$\mathbf{z}_{0}$ = initial state space vector, m and m/s

small aircraft for which no finite element (FE) model is developed, GVT is thus the only means of providing a modal basis. This basis is used to form a mathematical model of the structural dynamic behavior that is used for flutter predictions. For larger aircraft, the main GVT objective is to update and validate the finite element model of the aircraft, which, in turn, will be used for making reliable flutter predictions.

Aircraft GVT has now become standard industrial practice [1,2]. However, nonlinearity is frequently encountered during these tests and represents a significant challenge for aerospace engineers [3,4]. For instance, besides the nonlinear fluid structure interaction, typical nonlinearities include backlash and friction in control surfaces and joints, hardening nonlinearities in engine to pylon connections, and saturation effects in hydraulic actuators. Although there has been some recent progress in nonlinearity detection and characterization in aerospace structures [5,6], further developments for accounting for nonlinear dynamic phenomena during aircraft certification are still needed [7]. They would clearly improve the confidence practicing engineers have in both aircraft testing and modeling.

In this context, the nonlinear normal mode (NNM) theory offers a solid theoretical tool for interpreting a wide class of nonlinear dynamical phenomena, yet NNMs have a clear conceptual relation to the linear normal modes (LNMIs) [8,9]. However, most engineers still view NNMs as a concept that is foreign to them. One reason suppor ting this statement is that most existing constructive techniques for computing NNMs are based on asymptotic approaches [10,11]. Because these techniques rely on fairly involved mathematical developments, they are usually limited to the analysis of systems with low dimensionality.

There have been few attempts to compute NNMs using numerical methods [12–15]. The objective of this paper is to demonstrate that 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 proposed in [16] is used in the present study to compute the NNMs of the airframe of the Morane Saulnier Paris aircraft, for which the ground vibration tests have exhibited some nonlinear structural behaviors. Very interesting nonlinear dynamics, including nonlinear modal interactions, will be revealed.

The paper is organized as follows. Section II briefly reviews the two main definitions of an NNM. Section III describes the algorithm used for NNM computation. Section IV introduces the mathematical modeling of the considered aircraft and discusses carefully the results of nonlinear modal analysis. The conclusions of the present study are summarized in Sec. V.

## I. Introduction

DURING 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

<footnote>*Professor, Department of Aerospace and Mechanical Engineering, 1, Chemin des Chevreuil, B 4000.</footnote>

<footnote>†Ph.D. Student, Department of Aerospace and Mechanical Engineering, 1, Chemin des Chevreuil, B 4000.</footnote>

<footnote>‡Research Scientist, 29 Avenue de la Division Leclerc, F 92320.</footnote>
## II. Brief Review of Nonlinear Normal Modes

A detailed description of NNMs and of their fundamental properties (e.g., frequency energy dependence, bifurcations, and stability) is given in [8,9]. For completeness, the two main definitions of an NNM are briefly reviewed in this section.

The free response of discrete conservative mechanical systems with n degrees of freedom (DOFs) is considered, assuming that continuous systems (e.g., beams, shells, or plates) have been spatially discretized using the finite element method. The governing equations of motion are

$$M\ddot{x}(t) + Kx(t) + f_{nl}\{x(t)\} = 0 \quad (1)$$

where M and K are the mass and stiffness matrices, respectively, x(t) and $\ddot{x}(t)$ are the displacement and acceleration vectors, respectively, and $f_{nl}$ 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]:

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 LNM 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

$$\ddot{x}_1 + (2x_1 - x_2) + 0.5x_1^3 = 0 \quad \ddot{x}_2 + (2x_2 - x_1) = 0 \quad (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

to the minimal period of the periodic motion and at an energy equal to the conserved total energy during the motion (i.e., the sum of kinetic and potential energies). A branch, represented by a solid line, is a family of NNM motions possessing the same qualitative features. The backbone of the plot is therefore formed by two branches, which represent in phase (S11+) and out of phase (S11−) NNMs. The FEP clearly shows that the nonlinear modal parameters have a strong dependence on the total energy in the system based on the following:

1) The frequency of both the in phase and out of phase NNMs increases with the energy level, which reveals the hardening characteristic of the system.

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 curve), whereas the out of phase NNM localizes to the first DOF (i.e., it resembles a horizontal curve).

<figure><img src="image_2.png" /><figcaption>Fig. 2 Algorithm for NNM computation.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3 Morane-Saulnier Paris aircraft.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 1 Frequency-energy plot of the two-DOF system.</figcaption></figure>
<table><caption>Table 1 Aircraft properties</caption><thead><tr><th>Length, m</th><th>Wingspan, m</th><th>Height, m</th><th>Wing area, m²</th><th>Weight, kg</th></tr></thead><tbody><tr><td>10.4</td><td>10.1</td><td>2.6</td><td>18</td><td>1945</td></tr></tbody></table>

Such a plot gives a very clear indication of the effects of nonlinearity on modal parameters.

### III. Numerical Computation of NNMs

The numerical method proposed here for the NNM computation relies on two main techniques, namely, shooting and pseudo arclength continuation. A detailed description of the algorithm is given in [16].

#### A. Shooting Method

The equations of motion of system (1) can be recast into state space form:

$$ \dot{\mathbf{z}} = \mathbf{g}(\mathbf{z}) \quad (3) $$

where $\mathbf{z} = [\mathbf{x}^*\dot{\mathbf{x}}^*]^*$ is the 2n dimensional state vector and star denotes the transpose operation, and

$$ \mathbf{g}(\mathbf{z}) = \begin{pmatrix} \dot{\mathbf{x}} \\ -M^{-1}[K\mathbf{x} + \mathbf{f}_{nl}(\mathbf{x}, \dot{\mathbf{x}})] \end{pmatrix} \quad (4) $$

is the vector field. The solution of this dynamical system for initial conditions $\mathbf{z}(0) = \mathbf{z}_0 = [\mathbf{x}_0^*\dot{\mathbf{x}}_0^*]^*$ is written as $\mathbf{z}(t) = \mathbf{z}(t, \mathbf{z}_0)$ in order to exhibit the dependence on the initial conditions, $\mathbf{z}(0, \mathbf{z}_0) = \mathbf{z}_0$. A solution $\mathbf{z}_p(t, \mathbf{z}_{p0})$ is a periodic solution of the autonomous system (3) if $\mathbf{z}_p(t, \mathbf{z}_{p0}) = \mathbf{z}_p(t + T, \mathbf{z}_{p0})$, where T is the minimal period.

The NNM computation is carried out by finding the periodic solutions of the governing nonlinear equations of motion (3). In this context, the shooting method is probably the most popular numerical technique. It solves numerically the two point boundary value problem defined by the periodicity condition

$$ \mathbf{H}(\mathbf{z}_{p0}, T) \equiv \mathbf{z}_p(T, \mathbf{z}_{p0}) - \mathbf{z}_{p0} = \mathbf{0} \quad (5) $$

$\mathbf{H}(\mathbf{z}_0, T) = \mathbf{z}(T, \mathbf{z}_0) - \mathbf{z}_0$ is called the shooting function and 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 initial conditions $\mathbf{z}_{p0}$ and the period T that realize a periodic motion.

Starting from some assumed initial conditions $\mathbf{z}_{p0}^{(0)}$, the motion $\dot{\mathbf{z}}_p^{(0)}(t, \mathbf{z}_{p0}^{(0)})$ at the assumed period $T^{(0)}$ can be obtained by numerical time integration methods (e.g., Runge Kutta or Newmark schemes). 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 $\mathbf{z}(t)$ is a solution of the autonomous system (3), then $\mathbf{z}(t + \Delta t)$ is geometrically the same solution in state space for any $\Delta 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 therefore computed by solving the augmented two point boundary value problem defined by

$$ \mathbf{F}(\mathbf{z}_{p0}, T) \equiv \begin{cases} \mathbf{H}(\mathbf{z}_{p0}, T) = 0 \\ h(\mathbf{z}_{p0}) = 0 \end{cases} \quad (6) $$

where $h(\mathbf{z}_{p0}) = 0$ is the phase condition.

#### B. Continuation of Periodic Solutions

Because of the frequency energy dependence, the modal parameters of an NNM vary with the total energy. An NNM family, governed by Eqs. (6), therefore traces a curve, termed an NNM branch, in the $(2n + 1)$ dimensional space of initial conditions and period $(\mathbf{z}_{p0}, T)$. Starting from the corresponding LNM at low energy, the computation is carried out by finding successive points $(\mathbf{z}_{p0}, T)$ of the NNM branch using methods for the numerical continuation of periodic motions [19,20]. The so called pseudo arclength continuation method is used herein.

Starting from a known solution $(\mathbf{z}_{p0}(j), T_{(j)})$, the next periodic solution $(\mathbf{z}_{p0(j+1)}, T_{(j+1)})$ on the branch is computed using a predictor step and a corrector step.

##### 1. Predictor Step

At step j, a prediction $(\tilde{\mathbf{z}}_{p0(j+1)}, \tilde{T}_{(j+1)})$ of the next solution $(\mathbf{z}_{p0(j+1)}, T_{(j+1)})$ is generated along the tangent vector to the branch at the current point $\mathbf{z}_{p0(j)}$:

$$ \begin{bmatrix} \tilde{\mathbf{z}}_{p0(j+1)} \\ \tilde{T}_{(j+1)} \end{bmatrix} = \begin{bmatrix} \mathbf{z}_{p0(j)} \\ T_{(j)} \end{bmatrix} + s_{(j)} \begin{bmatrix} \mathbf{p}_z(j) \\ p_{T(j)} \end{bmatrix} \quad (7) $$

where $s_{(j)}$ is the predictor stepsize. The tangent vector $\mathbf{p}_{(j)} = [\mathbf{p}_z(j) p_{T(j)}]^*$ to the branch defined by (6) is solution of the system

<figure><img src="image_5.png" /><figcaption>Fig. 4 a) Connection between external fuel tank and wing tip (top view). Close-up of b) front- and c) rear-bolted attachments.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5 Finite element model of the Morane–Saulnier Paris aircraft.</figcaption></figure>

$$ \left[ \frac{\partial \mathbf{H}}{\partial \mathbf{z}_{p0}} \begin{array}{c} (\mathbf{z}_{p0,(j)}, T_{(j)}) \\ \frac{\partial h}{\partial \mathbf{z}_{p0}} \end{array} \frac{\partial \mathbf{H}}{\partial T} \begin{array}{c} (\mathbf{z}_{p0,(j)}, T_{(j)}) \\ 0 \end{array} \right] \left[ \mathbf{p}_{z,(j)} \atop P_{T,(j)} \right] = \left[ \mathbf{0} \atop \mathbf{0} \right] \quad (8) $$

with the condition $||\mathbf{p}_{(j)}|| = 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 $\mathbf{z}_{p0,(j+1)}^{(k+1)} = \mathbf{z}_{p0,(j+1)}^{(k)} + \Delta \mathbf{z}_{p0,(j+1)}^{(k)}$ and $T_{(j+1)}^{(k+1)} = T_{(j+1)}^{(k)} + \Delta T_{(j+1)}^{(k)}$ are computed by solving the overdetermined linear system using the Moore–Penrose matrix inverse:

$$ \left[ \frac{\partial \mathbf{H}}{\partial \mathbf{z}_{p0}} \begin{array}{c} (\mathbf{z}_{p0,(j+1)}^{(k)}, T_{(j+1)}^{(k)}) \\ \frac{\partial h}{\partial \mathbf{z}_{p0}} \end{array} \frac{\partial \mathbf{H}}{\partial T} \begin{array}{c} (\mathbf{z}_{p0,(j+1)}^{(k)}, T_{(j+1)}^{(k)}) \\ 0 \end{array} \right] \left[ \mathbf{P}_{z,(j)}^{*} \atop P_{T,(j)} \right] = \left[ \begin{array}{c} -\mathbf{H}(\mathbf{z}_{p0,(j+1)}^{(k)}, T_{(j+1)}^{(k)}) \\ -h(\mathbf{z}_{p0,(j+1)}^{(k)}) \\ 0 \end{array} \right] $$

<table><caption>Table 2 Natural frequencies of the linear finite element model</caption><thead><tr><th>Mode number</th><th>Frequency, Hz</th><th>Modal shape</th></tr></thead><tbody><tr><td>1</td><td>0.0936</td><td>Rigid body mode</td></tr><tr><td>2</td><td>0.7260</td><td>Rigid body mode</td></tr><tr><td>3</td><td>0.9606</td><td>Rigid body mode</td></tr><tr><td>4</td><td>1.2118</td><td>Rigid body mode</td></tr><tr><td>5</td><td>1.2153</td><td>Rigid body mode</td></tr><tr><td>6</td><td>1.7951</td><td>Rigid body mode</td></tr><tr><td>7</td><td>2.1072</td><td>Rigid body mode</td></tr><tr><td>8</td><td>2.5157</td><td>Rigid body mode</td></tr><tr><td>9</td><td>3.5736</td><td>Rigid body mode</td></tr><tr><td>10</td><td>8.1913</td><td>Two node wing bending</td></tr><tr><td>11</td><td>9.8644</td><td>Fin bending</td></tr><tr><td>12</td><td>16.1790</td><td>Fin torsion</td></tr><tr><td>13</td><td>21.2193</td><td>First T tail symmetric bending</td></tr><tr><td>14</td><td>22.7619</td><td>Front fuselage torsion</td></tr><tr><td>15</td><td>23.6525</td><td>T tail torsion</td></tr><tr><td>16</td><td>25.8667</td><td>Antelevator torsion</td></tr><tr><td>17</td><td>28.2679</td><td>Two node vertical fuselage bending</td></tr><tr><td>18</td><td>29.3309</td><td>Three mode wing bending</td></tr><tr><td>19</td><td>31.0847</td><td>Symmetric wing torsion</td></tr><tr><td>20</td><td>34.9151</td><td>Antisymmetric wing torsion</td></tr><tr><td>21</td><td>39.5169</td><td>x (fore and aft) symmetric wing bending</td></tr><tr><td>22</td><td>40.8516</td><td>Antisymmetric wing torsion + front fuselage</td></tr><tr><td>23</td><td>47.3547</td><td>Three node vertical fuselage bending</td></tr><tr><td>24</td><td>52.1404</td><td>Three node T tail bending</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6 a) First wing bending mode, b) first wing torsional mode, and c) second wing torsional mode.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Instrumentation of the rear attachment of the right wing.</figcaption></figure>

where the prediction is used as initial guess. This iterative process is carried out until convergence is achieved. The convergence test is based on the relative error of the periodicity condition, i.e., $||\mathbf{H}(\mathbf{z}_{p0}, T)||/||\mathbf{z}_{p0}|| < \epsilon$ where $\epsilon$ is usually set to $10^{-6}$.

An important technical remark is that the method requires the evaluation of the $2n \times 2n$ Jacobian matrix:

$$ \frac{\partial \mathbf{H}}{\partial \mathbf{z}_0}(\mathbf{z}_0, T) = \left. \frac{\partial \mathbf{z}(t, \mathbf{z}_0)}{\partial \mathbf{z}_0} \right|_{t=T} - I \quad (10) $$

where $I$ is the $2n \times 2n$ identity matrix. The classical finite difference approach requires us to perturb successively each of the $2n$ initial
<figure><img src="image_9.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 FEP of an NNM involving T-tail symmetric bending.</figcaption></figure>
conditions and integrate the nonlinear governing equations of motion. This is extremely computational intensive for large scale finite element models such as the one considered in this study. Targeting a substantial reduction of the computational cost, sensitivity analysis is exploited for determining $\partial z(t, z_0)/\partial z_0$. It amounts to differentiating the equations of motion (3) with respect to the initial conditions $z_0$, which leads to

$$ \frac{d}{dt} \left[ \frac{\partial z(t, z_0)}{\partial z_0} \right] = \left. \frac{\partial g(z)}{\partial z} \right|_{z(t, z_0)} \left[ \frac{\partial z(t, z_0)}{\partial z_0} \right] \quad (11) $$

with

$$ \frac{\partial z(0, z_0)}{\partial z_0} = I \quad (12) $$

because $z(0, z_0) = z_0$. Hence, the matrix $\partial z(t, z_0)/\partial z_0$ at $t = T$ can be obtained by numerically integrating over $T$ the initial value problem defined by linear ordinary differential equations (11) with initial conditions (12).

The complete algorithm is shown in Fig. 2.

## IV. Nonlinear Modal Analysis of the Morane-Saulnier Paris Aircraft

The numerical computation of the NNMs of a complex, real world structure is addressed in this section. The structure is the airframe of the Morane Saulnier Paris aircraft, represented in Fig. 3. This French jet aircraft was built during the 1950s and was used as a trainer and

liaison aircraft. The structural configuration under consideration corresponds to the aircraft without its jet engines and standing on the ground through its three landing gears with deflated tires. General characteristics of the aircraft are listed in Table 1. A specimen of this airplane is present in ONERA's laboratory. It is used for training engineers and technicians for GVT and for research purposes. Ground vibration tests exhibited nonlinear behavior in the connection between the wings and external fuel tanks located at the wing tip. Figure 4 shows that this connection consists of bolted attachments.

### A. Aircraft Structural Model

#### 1. Finite Element Model of the Underlying Linear Structure

The FE model of the linear part of the full scale aircraft, illustrated in Fig. 5, was elaborated from drawings at ONERA. The wings, T tail, and fuselage are modeled by means of 2 D elements such as beams and shells. Three dimensional (3 D) spring elements, which take into account the structural flexibility of the tires and landing gears, are used as boundary conditions of the aircraft. At each wing tip, the front and rear connections between the wings and the fuel tanks are modeled using beam elements. The FE model, originally created in 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 [0 50 Hz] frequency range are given in Table 2. The first nine modes correspond to aircraft rigid body modes, i.e., six modes are landing gear suspension modes, whereas three modes involve rigid body motions of the control surfaces. The frequency range of the rigid body modes is comprised between 0.09 and 3.57 Hz, i.e., noticeably lower than the first elastic mode located at 8.19 Hz. The

<figure><img src="image_11.png" /><figcaption>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).</figcaption></figure>
modal shapes of different elastic normal modes of vibrations are depicted in Fig. 6. Figure 6a represents the first wing bending mode. The first and second wing torsional modes are depicted in Figs. 6b and 6c and correspond to symmetric and antisymmetric wing motions, respectively. These latter modes are of particular interest because there is a significant deformation of the nonlinear connections between the wings and fuel tanks. The other modes concern the aircraft tail and are consequently almost unaffected by the nonlinear connections.

## 2. Reduced Order Modeling

Because the nonlinearities are spatially localized, condensation of the linear components of the model can be achieved using the Craig Bampton reduction technique [21]. This will lead to a substantial decrease in the computational burden without degrading the computational accuracy, at least in the frequency range of interest. We stress that many practical applications possess spatially localized nonlinearities, and so the procedure developed herein is of wide applicability.

The Craig Bampton method expresses the complete set of initial DOFs in terms of retained DOFs and internal vibration modes of the primary structure clamped on the retained nodes. To introduce the nonlinear behavior of the connections between the wings and the tanks, the reduced order model of the aircraft is constructed by keeping one node on both sides of the attachments. For each wing, four nodes are retained, namely, two nodes for the front attachment and two nodes for the rear attachment. In total, eight nodes of the initial FE model possessing, each 6 DOFs and 500 internal modes of vibrations, are kept in the reduced model; the FE model is thus reduced to 548 DOFs.

To assess the accuracy of the reduced model, its LNM s are compared to those predicted by the initial FE model. The comparison

is performed in the space of the initial model after projecting the reduced modes back into the original space. The deviation between the mode shapes of the original $y_{(o)}$ and reduced $y_{(r)}$ models is determined using the modal assurance criterion (MAC):

$$MAC = \frac{|\mathbf{y}^*_{(o)}\mathbf{y}_{(r)}|^2}{|\mathbf{y}^*_{(o)}\mathbf{y}_{(o)}||\mathbf{y}^*_{(r)}\mathbf{y}_{(r)}|} \quad (13)$$

MAC values range from 0 in the absence of correlation to 1 for a complete correspondence. In the [0 100 Hz] range, MAC values between modes shapes are all greater than 0.999, and the maximum relative error on the natural frequencies is 0.2%. The accuracy of this linear reduced model is therefore excellent. We note that much less than 500 internal modes are sufficient to build a good reduced model in the [0 100 Hz] frequency range. However, a larger number of modes was deliberately chosen for two main reasons. On the one hand, it serves to illustrate that our NNM algorithm can deal with systems of relatively high dimensionality. On the other hand, due to the presence of harmonics, nonlinear modal interactions may occur between a mode in the frequency range of interest and another mode outside this range. More internal modes are therefore necessary to guarantee the accuracy of the reduced model in nonlinear regimes of motion.

## 3. Aircraft Nonlinearities

The existence of a softening nonlinear behavior was evidenced during ground vibration tests conducted at ONERA. In particular, frequency response function measurements revealed a decrease in resonant frequencies when the excitation level was increased. The connections between the wings and fuel tanks were assumed to cause this observed nonlinear effect. To confirm this hypothesis, the front

<figure><img src="image_12.png" /><figcaption>Fig. 11 FEP of an NNM involving symmetric wing torsion (linear frequency = 31.08 Hz). NNM s 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).</figcaption></figure>
and rear connections of each wing were instrumented, and measure
ments dedicated to nonlinearity characterization were carried out.
Specifically, accelerometers were positioned on both the wing and
tank sides of the connections, and two shakers were located at the
tanks. The instrumentation at one of the connection is shown in Fig. 7.

The dynamic behavior of these connections in the vertical direction
was investigated using the restoring force surface method [22].
Newton’s second law applied at the wing side of one connection
writes

$$
m_c \ddot{x}_c(t) + f_{r,c}\{x_c(t), \dot{x}_c(t)\} = 0 \quad (14)
$$

where $f_{r,c}$ is the restoring force at the considered DOF. The index $c$ is
related to the connection under consideration (i.e., either the rear or
front attachment of the left or right wing). From Eq. (14), the restoring
force is obtained by

$$
f_{r,c} = -m_c \ddot{x}_c \quad (15)
$$

Because the restoring force $f_{r,c}$ includes both the force related to the
connection of interest and the force generated by the elastic

deformation of the wing, a quantitative assessment of the nonlinear
effects cannot be achieved precisely. However, plotting the measured
acceleration signal $\ddot{x}_c(t)$ against the relative displacement and
velocity across the connection can provide meaningful qualitative
information about the type of nonlinearity present in the bolted
connections. To this end, the aircraft was excited close to the second
torsional mode, which is known to activate nonlinear behavior, using
a band limited swept sine excitation. Figure 8 presents the resulting
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
highlight elastic nonlinearities. A softening elastic nonlinearity with
a piecewise linear characteristic can clearly be observed. This result
seems to be compatible with what was previously reported for bolted
connections in the literature [23,24]. Similar nonlinear effects were
also observed for the front connections. They are not shown here,
because these connections participated in the aircraft response to a
lesser extent.

Because nonsmooth nonlinearities require dedicated and very
specific time integration methods, which are not yet available in our
NNM computation algorithm, the nonlinearities in the connections
were regularized using polynomial expansions. This is why the

<figure><img src="image_13.png" /><figcaption>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).</figcaption></figure>
stiffness curve was reconstructed by considering the mathematical model

$$f_{r,c} = kx_{rel} + k_{nl}x_{rel}^3 \quad (k_{nl} < 0) \quad (16)$$

A value of $-10^{13}$ N/m³ was found to be adequate for the cubic spring. Damping in the connections was not considered because the focus of this study is on the NNMs of the underlying Hamiltonian system.

The reduced order model, generated in the SAMCEF software, was exported to MATLAB® where nonlinearities were implemented. Four nonlinearities were added, each one being located between the two nodes defining a nonlinear connection.

## B. Nonlinear Normal Modes and Corresponding Oscillation Frequencies

The NNM computation was carried out in the MATLAB® environment using the nonlinear reduced order model built in the preceding sections. For conciseness, only the modes that are the most relevant ones for the purpose of this study are described in what follows.

The first mode examined herein is the nonlinear counterpart of the first T tail symmetric bending mode (i.e., mode 13 in Table 2). As in Sec. II, the computed backbone branch and related NNM motions are depicted in a FEP in Fig. 9. Because there is no visible wing deformation, the modal shape and the corresponding oscillation frequency remain practically unchanged when the total energy is increased in the system. Despite the fact that a nonlinear aircraft

model is considered, this mode is not affected by nonlinearity and can be considered as a purely linear mode.

The nonlinear extension of the two node wing bending mode (i.e., 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 compatible with physical constraints. The frequency of the NNM motions on the backbone slightly decreases with increasing energy levels, which results from the softening characteristic of the nonlinearity. A quantitatively similar decrease in frequency between low level and high level swept sine excitations was observed during the experimental tests. Even though the frequency decrease is small, we remark that an accuracy of 0.001 Hz is usually sought during 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 10c) is 0.99. Accordingly, the modal shapes do not change much over the considered energy range and resemble the corresponding LNM.

Figure 11 represents the FEP of the first symmetric wing torsional 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 of the connections. As a result, the oscillation frequency has a marked energy dependence along the backbone branch. This decay of the natural frequency is in agreement, at least qualitatively, with what observed experimentally and confirms the relevance of this nonlinear modal analysis. Conversely, the modal shapes on the backbone branch are only weakly altered by the nonlinearities; the MAC value between the NNMs at low and high energy levels is equal to 0.98.

In addition to the main backbone branch, three other NNM branches that are localized to a specific region of the FEP can be

<figure><img src="image_14.png" /><figcaption>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).</figcaption></figure>
noticed. These tongues bifurcate from the backbone branch of the
considered mode and bifurcate into the backbone branch of another
mode, thereby realizing an internal resonance between two modes.
Because of harmonics, modal interactions with distinct linear
frequencies can be generated in nonlinear regimes of motion.
Specifically, the 3:1, 5:1, and 9:1 internal resonances in Fig. 11
involve third, fifth, and ninth harmonics of the fundamental
frequency of the first symmetric wing torsional mode.

A close up of the 3:1 internal resonance is depicted in Fig. 12. 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 linear modal basis, but their mixing realized by the NNM in Fig. 12c is not. It is therefore a mode with no linear counterpart. These tongues of internal resonance represent an inherently nonlinear dynamical phenomenon. Focusing now on the frequency content along the branch, a third harmonic progressively appears for increasing energies. In the vicinity of the NNM in Fig. 12c, the aircraft undergoes a subharmonic motion with two dominant frequencies. The wing vibrates at a frequency around 31 Hz, whereas the tail vibrates at a frequency close to 93 Hz. Beyond this point, the relative importance of the third harmonic component grows, until a complete coalescence with the tail torsional mode at 93 Hz is realized.

A completely similar dynamical mechanism exists for the 5:1 and 9:1 internal resonances, but the wing torsional mode interacts with other modes of the linear basis at 153 and 275 Hz, respectively. These modes correspond to local modes of some aircraft parts and are not shown here. We note that such nonlinear modal interactions were also observed in [9] for the two DOF system (2). However, they appeared at very high energies, and this is why they were not present in Fig. 1.

Finally, the second (antisymmetric) wing torsional mode (i.e., mode 20 in Table 2) is plotted in Fig. 13. The oscillation frequency undergoes a decrease around 5% in the considered energy range. The MAC value between the modal shapes at low and high energy levels is 0.97. The modal shapes are therefore affected by nonlinearity but not to a large extent. As for the first wing torsion mode, one tongue of internal resonance is present in the FEP and indicates that nonlinear modal interactions are in fact generic in such a large scale structure.

The computation of the backbone branch up to the tongue in Fig. 13 requires 20 min of CPU time with a single 2.67 GHz processor. Because of the presence of bifurcations, the computation of the branch of internal resonance demands about 1 h. This shows that the computation of NNM s of real world structures is now within reach using a conventional laptop computer.

V. Conclusions

The objective of this paper was to compute the nonlinear normal modes (NNMs) of a full scale aircraft together with the corres ponding oscillation frequencies. NNM s were represented in a frequency energy plot, in which the impact that nonlinearities have on modal parameters can readily be seen.

Parts of aircraft, such as landing gears, connections between the
horizontal tail plane and vertical tail plane, and pylons, are known to
be nonlinear. When these components are excited during ground
vibration testing, complex nonlinear phenomena that cannot be
described using linear tools may be observed. Therefore, the first
practical implementation of this work is to provide a rigorous
framework for the understanding and interpretation of these
phenomena. For instance, the present study showed that only a few
aircraft modes were influenced by the nonlinear bolted connections.
Their frequencies were found to decrease by a few percent in the
considered energy range, whereas their mode shapes were practically
unaffected. All this information is relevant for the practicing
aerospace engineer. Because of mode bifurcations, supernumerary
nonlinear modes with no linear counterpart (i.e., that are not the direct
extension of linear normal modes) were also reported in this study. To
the authors' knowledge, this is the first time that such modes are
evidenced for a real life structure. The practical realization of these
nonlinear resonances is still an open question in view of the small

frequency region in which they live. This will be addressed in
subsequent research by applying, numerically and experimentally, a
slow sweep sine excitation and nonlinear force appropriation to this
aircraft.

Finally, this research paves the way for the constructive use of nonlinearity for design. Some nonlinear phenomena with no linear counterpart such as targeted energy transfer could be exploited for protecting more sensitive parts of aircraft by transferring or redistributing energy from one substructure to another substructure.

References

[1] Goge, D., Boswald, M., Fullekrug, U., and Lubrina, P., "Ground Vibration Testing of Large Aircraft State of the Art and Future Perspectives," *Proceedings of the 25th International Modal Analysis Conference*, Orlando, FL, Society of Experimental Mechanics, 2007.

[2] Peeters, B., Hendrix, W., Debille, J., and Climent, H., "Modern Solutions for Ground Vibration Testing of Large Aircraft," *Sound and Vibration*, Vol. 43, 2009, pp. 8-15.

[3] Rodriguez Ahlquist, J., Martinez Carreno, J., Climent, H., de Diego, R., and de Alba, J., "Assessment of Nonlinear Structural Response in A400M GVT," *Proceedings of the 28th International Modal Analysis Conference*, Jacksonville, FL, Society of Experimental Mechanics, 2010.

[4] Kehoe, M. W., and Freudinger, L. C., "Aircraft Ground Vibration Testing at the NASA Dryden Flight Research Facility," NASA TM 104275, 1993.

[5] Gloth, G., and Sinapius, M., "Influence and Characterisation of Weak Non Linearities in Swept Sine Modal Testing," *Aerospace Science and Technology*, Vol. 8, No. 2, 2004, pp. 111-120.

[6] Gloth, G., and Goge, D., "Handling of Non Linear Structural Characteristics in Ground Vibration Testing," *Proceedings of the International Conference on Noise and Vibration Engineering (ISMA)*, Leuven, Belgium, Katholieke Universiteit Leuven, 2004.

[7] Luber, W., "Application of Analytical Methods for Advanced Flight Flutter Test," *Proceedings of the 30th International Modal Analysis Conference*, Jacksonville, FL, Society of Experimental Mechanics, 2012.

[8] Vakakis, A. F., Manevitch, L. I., Mikhlin, Y. V., Pilipchuk, V. N., and Zevin, A. A., *Normal Modes and Localization in Nonlinear Systems*, Wiley, New York, 1996.

[9] Kerschen, G., Peeters, M., Golinval, J. C., and Vakakis, A. F., "Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist," *Mechanical Systems and Signal Processing*, Vol. 23, No. 1, 2009, pp. 170-194.
doi:10.1016/j.ymssp.2008.04.002

[10] King, M. E., and Vakakis, A. F., "An Energy Based Formulation for Computing Nonlinear Normal Modes in Undamped Continuous Systems," *Journal of Vibration and Acoustics*, Vol. 116, No. 3, 1994, pp. 332-340.
doi:10.1115/1.2930433

[11] Touzé, C., and Amabili, M., "Nonlinear Normal Modes for Damped Geometrically Nonlinear Systems: Application to Reduced Order Modelling of Harmonically Forced Structures," *Journal of Sound and Vibration*, Vol. 298, Nos. 4-5, 2006, pp. 958-981.
doi:10.1016/j.jsv.2006.06.032

[12] Slater, J. C., "A Numerical Method for Determining Nonlinear Normal Modes," *Nonlinear Dynamics*, Vol. 10, No. 1, 1996, pp. 19-30.
doi:10.1007/BF00114796

[13] Pesheck, E., Pierre, C., and Shaw, S. W., "A New Galerkin Based Approach for Accurate Non Linear Normal Modes Through Invariant Manifolds," *Journal of Sound and Vibration*, Vol. 249, No. 5, 2002, pp. 971-993.
doi:10.1006/jsvi.2001.3914

[14] Arquier, R., Bellizzi, S., Bouc, R., and Cochelin, B., "Two Methods for the Computation of Nonlinear Modes of Vibrating Systems at Large Amplitudes," *Computers and Structures*, Vol. 84, Nos. 24-25, 2006, pp. 1565-1576.
doi:10.1016/j.compstruc.2006.01.011

[15] Casini, P., and Vestroni, F., "Characterization of Bifurcating Non Linear Normal Modes in Piecewise Linear Mechanical Systems," *International Journal of Non Linear Mechanics*, Vol. 46, No. 1, 2011, pp. 142-150.
doi:10.1016/j.ijnonlinmec.2010.08.002

[16] Peeters, M., Viguié, R., Sérandour, G., Kerschen, G., and Golinval, J. C., "Nonlinear Normal Modes, Part II: Toward a Practical Computation Using Numerical Continuation Techniques," *Mechanical Systems and*
*Signal Processing*, Vol. 23, No. 1, 2009, pp. 195-216.
doi:10.1016/j.ymssp.2008.04.003

[17] Rosenberg, R. M., "Normal Modes of Nonlinear Dual Mode Systems," *Journal of Applied Mechanics*, Vol. 27, No. 2, 1960, pp. 263-268.
doi:10.1115/1.3643948

[18] Shaw, S. W., and Pierre, C., "Non Linear Normal Modes and Invariant Manifolds," *Journal of Sound and Vibration*, Vol. 150, No. 1, 1991, pp. 170-173.
doi:10.1016/0022-460X(91)90412-D

[19] Seydel, R., *Practical Bifurcation and Stability Analysis, from Equilibirum to Chaos*, 2nd ed., Springer Verlag, New York, 1994.

[20] Nayfeh, A. H., and Balachandran, B., *Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods*, Interscience, New York, 1995.

[21] Craig, R., and Bampton, M., "Coupling of Substructures for Dynamic Analysis," *AIAA Journal*, Vol. 6, No. 7, 1968, pp. 1313-1319.

[22] Masri, S. F., and Caughey, T. K., "A Nonparametric Identification Technique for Nonlinear Dynamic Problems," *Journal of Applied Mechanics*, Vol. 46, No. 2, 1979, pp. 433-447.
doi:10.1115/1.3424568

[23] Gaul, L., and Lenz, J., "Nonlinear Dynamics of Structures Assembled by Bolted Joints," *Acta Mechanica*, Vol. 125, Nos. 1-4, 1997, pp. 169-181.
doi:10.1007/BF01177306

[24] Hartwigsen, C. J., Song, Y., McFarland, D. M., Bergman, L. A., and Vakakis, A. F., "Experimental Study of Nonlinear Effects in a Typical Shear Lap Joint Configuration," *Journal of Sound and Vibration*, Vol. 277, Nos. 1-2, 2004, pp. 327-351.
doi:10.1016/j.jsv.2003.09.018