aerospace
Article
Continuation Methods for Nonlinear Flutter
Edward E. Meyer
Boreal Racing Shells, Seattle Rowing Center, 1116 West Ewing Street, Seattle, WA 98119, USA;
edward.e.meyer@borealracingshells.com; Tel.: +1-206-453-8778
Academic Editor: Konstantinos Kontis
Received: 27 October 2016; Accepted: 5 December 2016; Published: 9 December 2016
Abstract:
Continuation methods are presented that are capable of treating frequency domain flutter
equations, including multiple nonlinearities represented by describing functions. A small problem
demonstrates how a series of continuation processes can find all limit-cycle oscillations within
a specified region with a reasonable degree of confidence. Curves of the limit-cycle amplitude variation
with velocity, indicating regions of stability and instability with colors, give a compact view of the
nonlinear behavior throughout the flight regime. A continuation technique for reducing limit-cycle
amplitudes by adjusting various system parameters is presented. These processes are economical
enough to be a routine part of aircraft design and certification.
Keywords: aeroelasticity; multiple nonlinearity flutter; continuation methods; describing functions;
bifurcation; continuation optimization; controlling LCO amplitudes
1. Introduction
In the commercial airplane industry certification that an aircraft is free from flutter instabilities
is a major aspect of design and modification, requiring numerous analyses to cover the entire range
of flight conditions and design parameters. Because it is impossible to analyze all flight conditions,
such as altitude, speed, fuel loading and payload, limited sets of conditions are analyzed, often relying
on engineering judgment to pick the important conditions. To do the many analyses required, techniques
are typically limited to linear, frequency domain methods, using finite-element models reduced with
low-frequency free-vibration modes and linear unsteady aerodynamics from the doublet-lattice method.
Linear flutter equations assume infinitesimally small displacements; nonlinear flutter equations
are far more realistic, allowing for limit-cycle oscillation (LCO), self-excited, constant amplitude
oscillations that may be harmful or merely ride-quality problems depending on the amplitude.
The additional information provided by nonlinear analyses is important to the design of commercial
aircraft and should become a routine part of industrial flutter analyses.
Continuation methods, a class of methods for solving parameterized nonlinear equations [
1
],
when applied to linear flutter equations have proven advantageous for covering the range of flight
conditions by allowing for continuous variations in parameters, thereby reducing the number of
discrete conditions necessary [
2
]. For example, interpolating mass matrices at several fuel loadings
allows tracing the variation of flutter points with fuel loading in a continuous fashion. The success of
continuation methods in treating linear flutter equations is due to the fact that in spite of the term linear
flutter, the equations are nonlinear in several parameters, but linear in the generalized coordinates
(g.c.). Thus, it is a small step to extend continuation methods to treat nonlinear flutter equations,
identifying LCOs and tracing the variation of amplitude with parameters, such as velocity.
Single nonlinearities, for example free play or bilinear stiffness in one discrete degree-of-freedom,
have been studied for many years [
3
,
4
]. A single nonlinearity modeled with a describing function is
readily treated with linear flutter techniques, since the generalized coordinates can be normalized to
satisfy the amplitude of the single generalized coordinate. Multiple nonlinearities present a conceptually
Aerospace 2016, 3, 44; doi:10.3390/aerospace3040044 www.mdpi.com/journal/aerospace
Aerospace 2016, 3, 44 2 of 20
more difficult task because of the need to simultaneously adjust the generalized coordinates so that
they and the resulting nonlinearities satisfy the flutter equations. Various schemes have been proposed
to address this problem [
5
,
6
]. Here, it is shown how a continuation method designed for linear flutter
analysis can also be used with multiple nonlinearities with minor modifications. The result is a method
that is straightforward, efficient and familiar to flutter engineers.
Following a brief summary of the continuation method presented in [
2
] and its application to
nonlinear flutter equations, a small problem is used to illustrate the search for limit cycle oscillations.
With a series of continuation processes, the region of interest is covered with a grid of curves that find
almost all LCOs in the region. Because of the sparsity of the grid, an LCO is missed, but is discovered
using a continuation method that traces an optimal path toward a limit cycle. Starting from the limit
cycles encountered, curves of LCO amplitude versus velocity are traced with continuation. From one
of these curves at a specified velocity, the optimal-path continuation method is used to reduce the LCO
amplitude by adjusting gains and phases in a simple control system.
2. Continuation Method
Continuation methods solve systems of nonlinear equations that are functions of more variables
than equations, at discrete points over a range of the variables while maintaining continuity in the
nonlinear functions. Continuity in flutter equations is important because there are always multiple
solutions, known as aeroelastic modes, which occasionally become close and difficult to distinguish.
The problem to solve is:
f (x) = 0 ∈ R
m
, x ∈ R
n
, n > m (1)
starting from a known solution
x
0
and varying the independent variables through the desired range.
n −m
is the dimensionality of the solution: curves (1); surfaces (2); volumes (3), etc. The focus here is
on curves: n = m + 1.
The continuation method chosen to solve this is known as a pseudo-arc length method [
1
] with
a minimum-norm corrector. At step
j
, the next point is predicted using the tangent to the curve and
the known solution x
j
:
x
0
j+1
= x
j
+ αt
j
(2)
where the stepsize
α
is an approximate distance along the curve, hence the name pseudo-arc length.
Starting with the prediction
x
0
j+1
, the corrector iterates to find the solution
x
j+1
using a Newton-like
method, computing corrections h
i
satisfying:
J(x
i
)h
i
= −f (x
i
) (3)
and setting:
x
i+1
= x
i
+ h
i
, (4)
where:
J(x) = f
0
(x) =
"
∂ f
i
∂x
j
#
∈ R
m×n
(5)
is the Jacobian matrix of derivatives of f with respect to x. The tangent vector (t) is characterized by:
Jt = 0, ktk
2
=
p
t
T
t = 1, t =
∂x
∂τ
(6)
where τ is the arc length along the curve.
Aerospace 2016, 3, 44 3 of 20
Equation (3) is an underdetermined linear system that has an infinity of solutions; the smallest
such solution can be computed by factoring the transposed Jacobian into an
(n
,
n)
orthogonal matrix
Q times an (n, m) upper-triangular matrix R:
J
T
= QR =
[
Q
1
Q
2
]
"
R
1
0
#
, (7)
known as a QR factorization. Substituting Equation (7) into Equation (3) yields the minimum-norm
solution [7] (p. 300):
h
i
= −Q
1
R
−T
1
f (x
i
) (8)
Corrector iterations continue using Equations (3)–(8) until suitable convergence criteria are met,
for example for some absolute tolerance e
a
and relative tolerance e
r
:
kf (x
i+1
)k
2
< e
a
and kh
i
k
2
< e
a
+ e
r
kx
i+1
k
2
(9)
In preparation for the next predictor step, the tangent is computed from the (
n
,
n −m)
matrix
Q
2
in Equation (7), a by-product of the QR factorization.
Q
2
is an orthogonal matrix spanning the null
space of the Jacobian,
JQ
2
= 0 ∈ R
m×n−m
Q
T
2
Q
2
= I ∈ R
n−m×n−m
(10)
An arbitrary vector w is projected onto the null space with:
¯
t = Q
2
Q
T
2
w = Q
"
0 0
0 I
2
#
Q
T
w ∈ R
n
(11)
where
I
2
is the
(n −m
,
n −m)
identity and
t =
¯
t/k
¯
tk
2
to satisfy Equation (6). If
n = m +
1,
Q
2
consists
of one column, the projection only determines the sign of the tangent, and a logical choice for the
vector to project is the tangent at the previous continuation step to keep the curve tracking in the right
direction. In this case, it is not absolutely necessary to do the projection, but it avoids forming
Q
2
(see Section 8) and is necessary for the optimal-path technique (Sections 5.6, 5.8 and 6).
3. Flutter Equations
Treating the flutter equation in the frequency domain has computational advantages over the
time domain: the characteristic equation is a system of nonlinear algebraic equations instead of
a system of differential equations, which must be integrated at each flight condition. By contrast,
the algebraic equations can be solved over a continuous range of flight conditions. The disadvantage is
that nonlinearities must be approximated to conform to the harmonic motion assumption, Equation (12).
3.1. Assumed Motion
Linear, frequency domain flutter analyses are based on the assumption that the actual motion of
the structure (z) as a function of time is related to a set of complex generalized coordinates (ˆq) by:
z(t) = Φ<( ˆqe
st
) (12)
where
Φ
is typically a matrix of vibration modes,
t
is time,
s = σ + iω
is the Laplace variable,
ω
is the
frequency and oscillations are growing, neutral or decaying if the growth factor
σ
is positive, zero or
negative, respectively.
The actual motion in a nonlinear frequency domain flutter analysis does not in general
follow this assumption; retaining this assumption results in an a quasi-linear approximation [
8
].
Aerospace 2016, 3, 44 4 of 20
Describing functions and harmonic balance methods are examples of quasi-linear approximations.
Any approximation conforming to these assumptions should work with the methods presented here.
3.2. Equations
Various formulations of the flutter equations have been proposed [
9
–
13
]; any of them can be
used with the continuation method presented above. A general form of the frequency domain flutter
characteristic equations is:
h
s
2
M + sG + sV + (1 + id)K − qA(p, M) + T
i
ˆq =
ˆ
D ˆq = 0, (13)
where
M
,
K
,
G
,
V
,
A
and
T
are the (
n
s
,
n
s
) mass, stiffness, gyroscopic, viscous damping, unsteady
aerodynamic and control-system matrices, respectively,
d
is the structural damping coefficient,
q = ρV
2
/
2 is the dynamic pressure,
p = s/V
is the complex reduced frequency,
V
is the free-stream
velocity,
M = V/a
is the free-stream Mach number,
a
is the sonic velocity and
ˆ
D
is the complex
dynamic matrix.
The related equation:
ˆ
D(s) ˆy = 0, kˆyk
2
=
q
ˆy
∗
ˆy = 1, =(
ˆ
y
k
) = 0 (14)
is a nonlinear eigenvalue problem with
n
s
solution eigenpair
(s
i
,
y
i
)
,
i =
1,
. . . n
s
, normalized with the
2-norm of
ˆy
1.0 and the
k
-th component real. Contrast this with the generalized-coordinate vector
ˆq
,
which is a factor, η times ˆy:
ˆq = η ˆy, kˆqk
2
= η (15)
This distinction is important when the dynamic matrix is a function of generalized-coordinate
amplitudes, particularly when η is small or zero.
3.3. Conversion to Real
The continuation method presented is for real equations and variables; the flutter equations are
converted from a set of n
s
complex equations to 2n
s
real equations by setting:
D =
"
<(
ˆ
D
ij
) −=(
ˆ
D
ij
)
=(
ˆ
D
ij
) <(
ˆ
D
ij
)
#
y =
(
<(
ˆ
y
i
)
=(
ˆ
y
i
)
)
q = ηy, (16)
so that:
Dy = Dq = 0 ∈ R
2n
s
(17)
and the 2-norms of the eigenvectors and generalized-coordinate vectors are unchanged:
kˆyk
2
=
q
ˆy
∗
ˆy = kyk
2
=
q
y
T
y = 1
kˆqk
2
=
q
ˆq
∗
ˆq = kqk
2
=
q
q
T
q = η (18)
where
ˆ
( )
∗
is the conjugate transpose.
3.4. Matrix Parameterizations
The matrices in Equation (13) are, in general, not constant; rather, they are assumed to be
functions of problem-dependent parameters using various matrix parameterizations. For example,
aircraft mass matrices are often parameterized by payload and fuel loading by interpolating matrices
at several parameter values. More important for this study, the matrices might be parameterized by
generalized-coordinate amplitudes in various ways depending on the matrix and type of nonlinearity.
Aerospace 2016, 3, 44 5 of 20
The dynamic matrix is therefore a function of the generalized coordinates and a vector
p
of parameters,
such as
V
,
σ
, and
ω
:
D = D(p
,
q)
. Common examples of nonlinearities include structural free play,
nonlinear stiffnesses and control systems.
Nonlinear structural elements, such as free play and nonlinear stiffness, can be approximated
using describing functions, which conform to the harmonic motion assumption, Equation (12).
These stiffness elements are usually modeled with a single degree-of-freedom with no stiffness coupling
between it and other degrees-of-freedom [
4
]. A describing function for bilinear stiffness
k
jj
associated
with generalized coordinate
j
, with stiffness
k
0
when the displacement amplitude
|
ˆ
q
j
|
is below
δ
and
k
1
when it is above, is:
k
jj
(γ, r) = c(γ, r)k
0
c(γ, r) =
r +
2
π
(1 −r)
h
sin
−1
(γ) + γ
p
1 − γ
2
i
if γ ≤ 1
1 otherwise
r =
k
1
k
0
, γ =
δ
|
ˆ
q
j
|
(19)
Nonlinear controls equations can be treated by describing functions [
8
]; one of the first
multiple- nonlinearity flutter solutions was with structural nonlinearities in the control system [6].
3.5. Normalization
To make solutions unique, a phase normalization on the generalized coordinates is necessary,
for example by constraining the k-th component of ˆq to be real:
f
2n
s
+1
= =(
ˆ
q
k
) = q
2k
= 0 (20)
If in addition η is to be held to a constant value η
0
, another equation is added:
f
2n
s
+2
= q
T
q −η
2
0
= 0 (21)
Therefore, to trace curves 2
n
s
+
2, independent variables must be used if
η
is allowed to vary,
or 2n
s
+ 3 if it is held constant.
3.6. Continuation Formulations
The independent variables (
x
) comprise the 2
n
s
variables
y
or
q
, plus enough additional
parameters to give an underdetermined system. There are two cases to consider: constant
η
with
2n
s
+ 3 variables and variable η with 2n
s
+ 2 variables.
In what follows, various combinations of velocity, frequency,
σ
and
η
are used to demonstrate the
solution technique; other parameters are possible, but these three serve to illustrate typical analyses.
Each combination of parameters will be named according to which three of these parameters vary;
the fourth is implicitly constant. For example, a
V
-
σ
-
ω
analysis traces the variation in
σ
and
ω
with
velocity, holding
η
constant, as in a traditional linear flutter analysis, and a
V
-
ω
-
η
analysis traces LCO
boundaries, the goal of this study.
3.6.1. Constant η
Solving the flutter equations with the generalized-coordinate norm held to a constant value
η
0
presents a problem when
η
0
=
0 and possibly numerical difficulties when
η
0
is small. The vector
q = 0
is known as the trivial solution because even though it is technically a solution, it provides
no real information about the behavior of the structure. Yet solutions at
q = 0
are important for
determining the stability of the trivial solution, the basis of linear flutter solutions. For this reason,
when
η
is held constant, eigenvectors are used as variables instead of the generalized coordinates.