Continuation Methods for Nonlinear Flutter

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.


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

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: 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 : 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: and setting: where: is the Jacobian matrix of derivatives of f with respect to x.The tangent vector (t) is characterized by: where τ is the arc length along the curve.
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: known as a QR factorization.Substituting Equation (7) into Equation (3) yields the minimum-norm solution [7] (p.300): Corrector iterations continue using Equations ( 3)-( 8) until suitable convergence criteria are met, for example for some absolute tolerance a and relative tolerance r : 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, An arbitrary vector w is projected onto the null space with: where I 2 is the (n − m, n − m) identity and t = t/ t 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).

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

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: 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].
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.

Equations
Various formulations of the flutter equations have been proposed [9][10][11][12][13]; any of them can be used with the continuation method presented above.A general form of the frequency domain flutter characteristic equations is: 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: is a nonlinear eigenvalue problem with n s solution eigenpair (s i , y i ), i = 1, . . .n s , normalized with the 2-norm of ŷ 1.0 and the k-th component real.Contrast this with the generalized-coordinate vector q, which is a factor, η times ŷ: This distinction is important when the dynamic matrix is a function of generalized-coordinate amplitudes, particularly when η is small or zero.

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: so that: and the 2-norms of the eigenvectors and generalized-coordinate vectors are unchanged: where ( ) * is the conjugate transpose.

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.
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 | qj | is below δ and k 1 when it is above, is: 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].

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: If in addition η is to be held to a constant value η 0 , another equation is added: Therefore, to trace curves 2n s + 2, independent variables must be used if η is allowed to vary, or 2n s + 3 if it is held constant.

Continuation Formulations
The independent variables (x) comprise the 2n 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.

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.
The eigenvector is constrained to unit 2-norm with the addition of an equation like ( 21): y T y − 1 = 0, so there are m = 2n s + 2 equations, and three more parameters are necessary for n = m + 1 = 2n s + 3 independent variables.
For the demonstration parameters, the equations are: where I 2k: is row 2k of the order 2n s identity matrix.

Variable η
If η is allowed to vary, the additional normalization Equation ( 21) is not used, so there are m = 2n s + 1 equations: Equation ( 17) plus the phase normalization (Equation ( 20)), and two of the three parameters (V, σ, ω) are needed.Two combinations of the three parameters will be of interest, σ-ω-η at constant velocity (used to search for σ = 0) with: and V-ω-η at constant σ (used to trace LCO boundaries): If the first step in either of these processes starts from a known solution at η = 0 with: a modification in the tangent calculation is necessary, because the Jacobian and its null space vectors are then: The third column of Q 2 gives the desired greatest increase in η, so the first predicted solution for the σ-ω-η process is: In addition, an equation should be added to prevent converging to the trivial solution (q = 0), where η 0 is the 2-norm of the predicted q.For the first step, η 0 = α.

Demonstration Problem
A small, well-known problem illustrates the use of these continuation processes.The model is Example 9-1 in [14], implemented as Example HA145B in finite-element program MSC-NASTRAN [13], with a generalized-coordinate basis of n s = 10 free-vibration modes and doublet-lattice aerodynamics.Bilinear nonlinearities are applied to the first four diagonals of the (diagonal) stiffness matrix using the describing function of Equation ( 19) with δ = 0.05, r = 2 (hardening) for diagonals 1 and 2, r = 0.5 (softening) for diagonal 3 and r = 0.3 (softening) for diagonal 4. Figure 1 shows the model wing planform and unsteady aerodynamic grid.

Searching for Limit Cycle Oscillations
Limit cycle oscillation is a condition where the structure is oscillating with an amplitude that is neither growing nor decaying; it is in dynamic equilibrium and σ = 0. On the other hand, q = 0 is a static equilibrium with stability determined by the sign of σ.The goal then is to find all points within the region of interest where σ = 0.
The strategy is to find several points where σ = 0 using the formulations of Sections 3.6.1 and 3.6.2,then use these points to start V-ω-η processes to trace the variation of LCO amplitude with velocity.
An example of the results of an LCO analysis is shown in Figure 2 where η is plotted against velocity normalized with respect to the maximum velocity of interest, v = V/V max .Ordinarily, η would not be the plot variable of choice; instead, one or more of the nonlinear generalized-coordinate amplitudes would be plotted against velocity, but this choice is problem dependent, so η is used for generality.Of the two curves in Figure 2, only the bottom one extends to η = 0, so it is the only one that would be detected in a linear flutter analysis.Detecting the upper curve is why the processes of Sections 3.6.1 and 3.6.2are necessary.A systematic method for this search is to cover the range with a grid of V-σ-ω and σ-ω-η processes, any of which might encounter σ = 0 points.Figure 3 shows a minimal grid of six lines where horizontal lines are V-σ-ω processes (numbered 1, 3 and 4), and vertical lines are σ-ω-η processes (numbered 2, 5 and 6).Arrows indicate the direction of the continuation processes.Each of these lines represents one or more modes chosen to cover the frequency range of interest.This set of processes is started from a free-vibration solution, at the origin of Figure 3 (Point 0), with v = η = 0. Two processes are started from this free-vibration solution: Process 1, a V-σ-ω process along the v axis, and Process 2, a σ-ω-η process along the η axis.The remaining vertical and horizontal grid lines are started from these two by interpolating to the desired points.

Point 0: Free-Vibration
At V = 0, the aerodynamic term drops out, and with η = 0, q is zero for the purposes of evaluating the dynamic matrix D. The problem to be solved then is in general the complex nonlinear eigenvalue problem, Equation (14).Under certain conditions, Equation ( 14) reduces to a polynomial eigenvalue problem [15], or to a generalized eigenvalue problem [16].For example, in the absence of controls equations, viscous and structural damping and gyroscopics, the eigenvalue problem reduces to the real, generalized eigenvalue problem: Otherwise, it is usually necessary to resort to nonlinear eigenvalue solution techniques [15].Using the results of this eigensolution, representing the origin in Figure 2, the LCO search is begun with two continuation processes corresponding to the two axes of Figure 2: a V-σ-ω analysis at η = 0 (the horizontal or V axis) and a σ-ω-η analysis at V = 0 (the vertical or η axis).

Process 1:
V-σ-ω with η 0 = 0 A continuation process along the horizontal axis of Figure 2 is a V-σ-ω process at η 0 = 0, solving Equation ( 22), starting with select eigenvalues and eigenvectors from the free-vibration solution.This is equivalent to a linear flutter solution.Figure 4 shows results for the first four lowest frequency modes with green indicating a stable equilibrium and red unstable.
From these curves σ-ω-η Processes 5 and 6 in Figure 3 are started at the desired velocities by interpolating among solution points.Additionally, a V-ω-η processes could be started from the point where a curve crosses the σ = 0 axis, tracing the bottom curve in Figure 2.

Process 2: σ-ω-η at V = 0
In a σ-ω-η process η is allowed to vary; velocity is held constant; and the Equations are ( 23) and (23a).With V = 0, this process corresponds to the vertical axis of Figure 2. Start points for this process are again select eigenvalues and eigenvectors from the free-vibration solution.At zero velocity, the aerodynamic term drops out, and without viscous or structural damping, σ = 0 for each vibration mode.Only the frequency changes with increasing η, as shown in Figure 5.

Processes 3 and 4: V-σ-ω
Horizontal lines 3 and 4 in Figure 3 are V-σ-ω processes starting from Process 2 interpolated to the desired values of η. Figure 6 shows the result of these processes; as before, green indicates stability (negative σ), and red indicates instability.Three levels of η are highlighted for the only mode that crosses σ = 0.These crossings will be used in the V-ω-η processes of Section 3.6.2 to trace LCO boundaries.

Stability of Equilibria
At an LCO, if a small perturbation in the generalized coordinates decreases η and η continues to decrease with time, or if the perturbation increases η and η continues to increase with time, the equilibrium is unstable.On the other hand, if after the perturbation, the generalized coordinates return to the LCO amplitudes, the equilibrium is said to be stable.Mathematically, this means that stability is determined by the sign of the derivative of σ with respect to η: positive indicates instability, and negative indicates stability.This derivative can be computed from the components of the tangent vector corresponding to σ and q, recalling Equation ( 6), where τ again is the arc length along the curve.These tangent components are available at every step of a σ-ω-η process, but this information is only relevant at σ = 0.More important is the stability of each point of LCO boundaries, created with V-ω-η processes and Equation (23b).That Jacobian then must be converted to the Jacobian of Equation (23a) by replacing the first column with ∂ f /∂σ, where: It is not necessary to actually make this replacement and re-factorize the Jacobian; only the vectors u and v are of interest because with them, the QR factorization of the transposed Jacobian can be updated [7] (p.334): Updating the factorization requires a small fraction of the work required to factorize the Jacobian.From the updated Jacobian, the tangent used in Equation ( 29) is computed by projecting u onto the updated null space with Equation (11): As η increases or decreases from an unstable equilibrium, the next equilibrium may be stable or unstable.Often, an unstable equilibrium will be followed by a stable one; Figure 7 shows why this is true when the equilibriums are associated with the same aeroelastic mode.That this is not always true will be shown in Section 5.8, where an LCO is discovered associated with a different aeroelastic mode.5.7.Limit Cycles: V-ω-η Analysis at σ = 0 Among the six sets of continuation curves in Figures 6 and 7, there are five points where σ = 0, at the intersection of the green and red line segments in Figure 8.These points can be used to start V-ω-η processes, the goal in this exercise, shown as dotted curves in Figure 8.The upper dotted curve could be started from Point 1, 2 or 3; either start point will trace the same curve.Likewise, the lower dotted curve could be started from Point 4 or 5.

Optimal Path: V-σ-ω-η Analysis
Two of the grid lines in Figure 3 do not cross the σ = 0 axis (η = 0.5 in Figure 6 and v = 0.5 in Figure 7), but come close.If all four parameters were allowed to vary instead of holding one fixed and varying the other three, those two continuation processes might have also found limit cycle points.A continuation process where the number of unknowns is greater than the number of equations plus one requires only a minor modification of the continuation process and is in a sense an optimization procedure.
The idea is to allow any number of independent variables greater than the number of equations (m) and to trace a special curve on the surface (or volume, or higher-dimension manifold): the shortest path toward some goal, for example maximizing or minimizing some parameter β.The key to tracing this curve is the tangent to the curve, computed by projecting the vector: onto the null space of the Jacobian (Equation ( 11)).Aside from this modification, the continuation process is unchanged.
Here, the goal is σ = 0, a limit-cycle oscillation, and the equations are: Assuming the start point is at a negative σ, the tangent vector is computed by projecting: onto the null space of the Jacobian (Equation ( 11)).Figures 9 and 10 show the result of such a process, starting at a mode in Figure 4 that peaks around v = 0.9 and σ = −0.7 and continuing until σ = 0 is encountered.
Then, starting from the point where this curve meets σ = 0, a V-ω-η process traces a limit-cycle curve shown in Figure 11.The complete set of LCO curves, with stable and unstable LCO again shown in green and red, is shown in Figure 12.
Note that at v = 0.9 in Figure 12, as η increases from the unstable (static) equilibrium at zero, the next equilibrium encountered is also unstable, followed by two stable equilibriums and finally an unstable LCO.A σ-ω-η process would show continuity in two different modes, a three-Hertz mode and a 10-Hertz mode, and each of these would have the usual pattern of stable followed by unstable and vice versa.

Decreasing LCO Amplitudes
It is often desirable to reduce the amplitude of limit cycles, whether to avoid structural damage or to improve ride quality.Sometimes, this can be done with structural changes, such as reducing free play or increasing damping; in other situations, such changes might be impractical, and it may be necessary to resort to control systems.In either case, the technique of Section 5.8 can be used to decrease the amplitude of a generalized coordinate undergoing LCO.
For example, including a simple control system (matrix T in Equation ( 13)) can be used to reduce the amplitude of generalized coordinate 1 in the unstable LCO at normalized velocity 0.8 in Figure 12.Details of the control system are unimportant, as it is just for illustration; it simply feeds back displacement in generalized coordinates 1 and 2 to add stiffness to these coordinates with two complex gains: g 1 e iφ 1 and g 2 e iφ 2 .The equation to be solved is as before, Equation (23), with m = 2n s + 2, but now with three more unknowns than equations: A projection vector pointing towards decreasing | q1 | is: Beginning at v = 0.8 on the upper curve in Figure 12, where η = 1.18 and | q1 | = 1.16, the projection vector guides the solution toward | q1 | = 0, as shown in Figure 13.The stability of the LCO at each point on the curves is computed as in Section 5.6.

Bifurcation
When tracing a curve with a continuation method, it is possible to encounter a bifurcation, where two curves cross [1,17].Other characteristics of a bifurcation, useful in detecting their presence, are a rank-deficient Jacobian matrix and a change in the sign of the determinant of the Jacobian matrix augmented with the tangent vector.It can be shown [1] that a change in the sign of: indicates that a bifurcation has been passed.A method for computing µ is shown in Section 8 that is economical enough that it could (and should) be done at every step.One of the curves in the V-σ-ω process shown in Figure 6 encounters a bifurcation, although it is not evident at that scale.When plotted σ versus ω with an expanded scale, the bifurcation becomes more evident, as shown in Figure 14.Arrows indicate the direction of the positive determinant.Had the bifurcation not been detected, the process would have tracked only the green curve.Detecting the bifurcation and computing the bifurcation tangent allows the blue curve to be traced.
A brief summary of the steps necessary to the bifurcation branch is given here; for more detail, see [1,17].Branching from the bifurcation requires the tangent to the branch curve.For a simple bifurcation at x b , detected by a change in the sign of µ(x), the Jacobian will have one rank deficiency, and the Q 2 matrix has two columns, some combination of which gives the branch tangent; and another combination continues the original curve.
A singular value decomposition (SVD) [7] (p.76) of the Jacobian at x b , reveals the rank and is necessary for branching from the bifurcation.The last column of U (u m ) and the last two columns of V T (v m and v m+1 ) are left and right null vectors of the Jacobian, and the two tangents emanating from the bifurcation are combinations of these: where: One of these tangents is aligned with the original curve, the other with the bifurcation curve.Comparing with the previous tangent (t j ) determines which is the primary tangent (t c ) and the bifurcation tangent (t b ):

Programming Considerations
Much of the computational effort in the continuation method is in the factorization of the Jacobian matrix in Equation (7), where QR factorizations are used to solve the underdetermined system.Some continuation methods [1,18] add an equation in the corrector phase to constrain corrections, for example to hold one of the variables constant, creating a square set of linear equations that can be solved using an LU factorization instead of QR.Theoretically, LU is faster than QR; but QR is more stable, and modern computer architectures eliminate the speed advantage [7] (p.299).QR is also necessary for the optimal-path technique of Section 5.8.
The linear algebra library LAPACK [19] provides several variants of QR factorizations; the preferred routine for efficiency and stability is SGEQP3 [16] (Section 2.4.2.3).None of the QR routines in LAPACK return the Q and R matrices directly; instead, all of the information to create them is returned in the input matrix and a separate array named TAU.The R matrix is contained in the upper triangle of the input matrix, and subroutine STRSM solves a set of equations with it.Subroutine SORMQR is provided to multiply Q or Q T by a vector.
Projecting a vector onto the null space of the Jacobian can be done in LAPACK without actually forming the Q 2 matrix with two calls to subroutine SORMQR: that is, call SORMQR to multiply the vector by Q T , zero the first m elements of the result, then call it again to multiply by Q.
Updating the QR factorization to determine stability with Equation (32) can be done efficiently with the qrupdate package available at [20].
Evaluating the determinant in Equation ( 39) is simplified by recognizing that the tangent came from the QR factorization, that the determinant of an orthogonal matrix is ±1 and that the determinant of a triangular matrix is the product of the diagonals.Furthermore, the tangent and the QR factorization of the Jacobian it was computed from are available from the last corrector iteration, so the desired determinant is: where k is the number of non-zero elements in the TAU array [21] for LAPACK QR routines, and l is +1 if Q T t is positive, −1 otherwise.Derivatives are an essential part of continuation methods.Accurate derivatives of parameters and matrices with respect to the independent variables are necessary for computing the Jacobian matrix that is used both to iterate to solutions and predict the next solution.Deriving formulas for complicated expressions like Equation ( 19) can be tedious and error-prone.A technique that eliminates the need for closed expressions for the derivatives and can reduce coding errors is automatic differentiation [22].This programming technique gives exact derivatives at a cost on the order of evaluating the function itself, in contrast with difference approximations, which are not exact and much more expensive to compute.The autodiff.org[23] website is devoted to automatic differentiation resources.

Conclusions
A continuation method capable of treating frequency domain flutter equations with multiple nonlinearities approximated with describing functions has been developed.With it, a set of continuation processes can find with reasonable certainty all limit-cycle oscillations within a desired flight envelope, tracing curves of limit-cycle amplitude versus velocity and indicating stability with colors.
With a small modification to the continuation method, curves can be traced taking an optimal path with respect to any number of parameters toward some goal.This technique was used to seek out a limit cycle that was missed and to reduce limit-cycle amplitudes by adjusting control-system parameters.
Tracing curves with continuation methods, it is always possible to encounter bifurcation.Detecting bifurcation is inexpensive enough that it can be done at every step when tracing a curve; missing a bifurcation risks missing an important curve.When a bifurcation is detected, a straightforward calculation gives the two tangents, which are then used to start continuation processes to trace each curve from the bifurcation.

M
V/a, Mach number n s the number of generalized coordinates p s/V, complex reduced frequency x 2 2-norm of vector x = √ x T x q ρV 2 /2, dynamic pressure q, q complex and equivalent real generalized coordinates (g.c.) (Equations ( 12) and ( 16)) Q, R QR factors of the transposed Jacobian (Equation ( 7)) R n , R m×n the real n-vectors and m by n matrices (), () real and imaginary parts of a complex variable r, γ, δ bilinear stiffness describing-function variables (Equation ( 19)) s Laplace variable σ + iω t tangent vector V velocity (true airspeed) V max maximum velocity of interest (true airspeed) v normalized velocity V/V max w projection vector (Equation ( 11)) x independent variables x i j ith iteration at the jth continuation step (Equation ( 2)) x i ith element of vector x ŷ, y complex and equivalent real eigenvector z(t) motion of points of the structure (Equation ( 12))

Figure 3 .
Figure 3. Grid (blue) for finding limit cycles.Green means stable, red means unstable.
Figure 7  shows σ versus η at three velocities; for clarity, only one of the four modes traced in Process 1 is shown, the only one that encounters LCO.

Figure 8 .
Figure 8. Five start points for LCO curves (Green means stable, red means unstable).

Figure 11 .
Figure 11.LCO from the optimal path (black curve with arrow indicating tracking direction) (Green means stable, red means unstable).

1 |Figure 13 .
Figure 13.Reduction of an LCO amplitude with control system parameters (Green means stable, red means unstable).

Figure 14 .
Figure 14.Bifurcation in a V-σ-ω process.Arrows indicate the direction of positive µ, the blue curve arises from the bifurcation (Green means stable, red means unstable).

η 2 -
norm of the generalized coordinates (Equation (15)) ω oscillation frequency (imaginary part of the Laplace variable) Φ transformation from generalized to physical coordinates ρ fluid density σ growth factor (real part of the Laplace variable) τ arc length along a curve () T transpose