JOURNAL OF AIRCRAFT

Vol. 42, No. 3, May–June 2005

Direct Aeroelastic Bifurcation Analysis of a Symmetric

Wing Based on Euler Equations

K. J. Badcock,∗ M. A. Woodgate,† and B. E. Richards‡

University of Glasgow, Glasgow, Scotland G12 8QQ, United Kingdom

The application of a sparse matrix solver for the direct calculation of Hopf bifurcation points for the flexible
AGARD 445.6 wing in a transonic flow modeled using computational fluid dynamics is considered. The iteration
scheme for solving the Hopf equations is based on a modified Newton’s method. Direct solution of the linear system
for the updates has previously been restrictive for application of the method, and the sparse solver overcomes this
limitation. Previous work has demonstrated the scheme for aerofoil calculations. The current paper gives the first
three-dimensional results with the method, showing that an entire flutter boundary for the AGARD 445.6 wing can
be traced out in a time comparable to that required for a small number of time-marching calculations, yielding
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

two orders of magnitude improvement when compared to the time-marching approach.

Introduction expense of an outer iteration over the domains. This was tested on
the model problem in references in Beran and Carlson4 and Beran.5
T IME-DOMAIN analysis is the main solution method used in
computational aeroelasticity when the flow is modeled by the
Euler and Reynolds-averaged Navier–Stokes (RANS) equations.
The problems introduced by using a direct solver were resolved
in Badcock et al.,6 where a sparse matrix formulation was used
However, the need to execute searches over multiparameter space to make feasible the solution of the linear system for much larger
to identify stability behavior can lead to high computational cost grids. The Newton iteration was modified to enhance the efficiency
because of the need to do an unsteady coupled computational-fluid- of the scheme following work on approximate Jacobian matrices
dynamics (CFD) and computational-structural-dynamics (CSD) for CFD only problems reported in Badcock et al.7 The method was
calculation for each combination of parameters. This cost is not shown to be effective for tracing out flutter boundaries for symmetric
prohibitive when the intention is to examine behavior at previously aerofoils moving in pitch and plunge, with reductions of two orders
identified problem conditions, and there are several recent impres- of magnitude in the computational time required when compared
sive demonstrations of this kind for complete aircraft configurations with time marching.
(e.g., see Farhat et al.1 and Melville2 ). The current paper extends the method to calculate flutter bound-
A way of reducing the cost of parametric searches for stability aries for symmetric wings. The additional issues to be considered
behavior was proposed by Morton and Beran3 from the U.S. Air are the treatment of a moving grid around a deforming geometry (as
Force Laboratories. Their method uses dynamical systems theory opposed to rigid motions for the previous aerofoil cases), the use of a
to characterize the nature of the aeroelastic instability, with this addi- modal structural model (instead of the pitch-plunge equations), and
tional information concentrating the use of the CFD. In this way the the resulting requirement to pass information between nonmatch-
problem of locating a one-parameter Hopf bifurcation was reduced ing grids, and the larger problem size, and especially the impact of
from multiple time-marching calculations to a single steady-state this on the solution of the linear system. The formulation is con-
calculation of a modified system. This modified system calculates sidered in the following two sections and then results are presented
the value of the parameter for which an eigenvalue of the system for the AGARD 445.6 wing test case of Yates8 to demonstrate the
Jacobian matrix crosses the imaginary axis. A convection-diffusion feasibility of the method for three-dimensional problems.
problem was used to evaluate the approach in Beran and Carlson,4
and the method was then applied to an aeroelastic system consisting Aerodynamic and Structural Simulations
of an aerofoil moving in pitch and plunge in Morton and Beran.3 The Aerodynamics
linear system was solved using a direct method, and this motivated The three-dimensional Euler equations can be written in conser-
the use of an approximate Jacobian matrix to reduce the cost. Some vative form and Cartesian coordinates as
robustness problems were encountered when applying the method,
particularly at transonic Mach numbers. A complex variable for- ∂w f ∂Fi ∂Gi ∂Hi
mulation of the problem was introduced in Beran,5 which resolved + + + =0 (1)
∂t ∂x ∂y ∂z
some of these problems. An approach considered to reduce the dif-
ficulties of applying a direct solver to large linear systems was to
where w f = (ρ, ρu, ρv, ρw, ρ E)T denotes the vector of conserved
use domain decomposition to reduce the size of the system at the
variables. The flux vectors Fi , Gi and Hi are,
 
ρU ∗
Received 19 September 2003; revision received 20 February 2004; ac-
cepted for publication 27 February 2004. Copyright   ρuU ∗ + p 
c 2004 by University  
of Glasgow. Published by the American Institute of Aeronautics and Astro- F =
i
 ρvU ∗ 
 (2)
nautics, Inc., with permission. Copies of this paper may be made for personal  ρwU ∗ 
or internal use, on condition that the copier pay the $10.00 per-copy fee to
∗
the Copyright Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA U (ρ E + p) + ẋ
01923; include the code 0021-8669/05 $10.00 in correspondence with the
   
CCC.
∗ Reader, Computational Fluid Dynamics Laboratory, Department of
ρV ∗ ρW ∗
 ρuV ∗   ρuW ∗ 
Aerospace Engineering.    
G =  ρvV + p 
 ∗
H =  ρvW + p 
 ∗
† Research Assistant, Computational Fluid Dynamics Laboratory, Depart-
,
i i
 (3)
ment of Aerospace Engineering.
‡ Mechan Professor, Computational Fluid Dynamics Laboratory, Depart-  ρwV ∗   ρwW ∗ + p 
∗
ment of Aerospace Engineering. Associate Fellow AIAA. V (ρ E + p) + ẏ W ∗ (ρ E + p) + ż
731
 732 BADCOCK, WOODGATE, AND RICHARDS

In the preceding equations ρ, u, v, w, p, and E denote the density, for fsn + 1 . At convergence the fluid and structural solvers are prop-
the three Cartesian components of the velocity, the pressure, and the erly sequenced, at very little extra computational cost beyond what
specific total energy, respectively, and U ∗ , V ∗ , W ∗ the three Carte- is required for the aerodynamic solution.
sian components of the velocity relative to the moving coordinate The aerodynamic forces are calculated at face centers on the aero-
system, which has local velocity components ẋ, ẏ, and ż, that is, dynamic surface grid. The problem of communicating these forces
to the structural grid is complicated in the common situation that
U ∗ = u − ẋ (4) these grids not only do not match, but are also not even defined
on the same surface. This problem, and the influence it can have
V ∗ = v − ẏ (5) on the aeroelastic response, was considered in Goura13 and Goura
et al.,14 where a method was developed, called the constant-volume-
W ∗ = w − ż (6) tetrahedron (CVT) transformation. This method uses a combination
of projection of fluid points onto the structural grid, transformation
The variables here have been nondimensionalized with respect to of the projected point, and recovery of the out-of-plane component
the wing root chord c for x, y, and z; the freestream velocity U∞ to obtain a cheap but effective relation between deformations on the
for u, v, and w; the freestream density ρ∞ for ρ, U∞ /c for t and structural grid and those on the fluid grid. Denoting the fluid grid
ρ∞ U ∞2
for p. locations and aerodynamic forces as xa and fa , then
The time-marching results in the current work are obtained using
the Glasgow University code PMB (which stands for parallel multi- δxa = S(xa , xs , δxs )
block). A summary of some applications examined using the code
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

can be found in Badcock et al.7 A fully implicit steady solution of where S denotes the relationship defined by CVT. In practice this
the Euler equations is obtained by advancing the solution forward equation is linearized to give
in time by solving the discrete nonlinear system of equations
δxa = S(xa , xs )δxs
wnf + 1 − wnf  n+1

= Rf wf (7) and then by the principle of virtual work, fs = S T fa .
t The grid speeds on the wing surface are also needed, and these
The term on the right-hand side, called the residual, the discretiza- are approximated directly from the linearized transformation as
tion of the convective terms, given here by the approximate Riemann δ x˙a = S(xa , xs )δ ẋs
solver of Osher and Chakravarthy,9 MUSCL interpolation of Van
Leer,10 and Van Albada’s limiter. The sign of the definition of the where the structural grid speeds are given by
residual is opposite to convention in CFD, but this is to provide a
system of ordinary differential equations, which follows the con- δ x˙s =  α̇i φi . (11)
vention of dynamical systems theory, as will be discussed in the
next section. Equation (7) is a nonlinear system of algebraic equa- The geometries of interest deform during the motion. This means
tions, which is solved by an implicit method described in Cantariti that, unlike the rigid aerofoil problem, the aerodynamic mesh
et al.,11 the main features of which are an approximate linearization must be deformed rather than rigidly translated and rotated. This
to reduce the size and condition number of the linear system and the is achieved using transfinite interpolation of displacements (TFI)
use of a preconditioned Krylov subspace method to calculate the within the blocks containing the wing. More elaborate treatments
updates. that move blocks to maintain grid orthogonality are possible (see
The steady-state solver is applied to unsteady problems within a Tsai et al.15 ) but are not necessary here because only small wing
pseudo-time-stepping iteration as formulated by Jameson.12 deflections are encountered and the blocks in the mesh can be ex-
tended well away from the wing. The wing surface deflections are
Structural Dynamics, Intergrid Transformation, and Mesh Movement interpolated to the volume grid points xi jk as
The wing deflections δxs are defined at a set of points xs by
δxi jk = ψ 0j δxa,ik (12)
δxs = αi φi (8)
where ψ 0j are values of a blending function (see Gordon and Hall16 ),
where φi are the mode shapes calculated from a full finite element which varies between one at the wing surface (here j = 1) and zero
model of the structure and αi are the generalized coordinates. By at the block face opposite. The surface deflections xa,ik are obtained
projecting the finite element equations onto the mode shapes and from the transformation of the deflections on the structural grid and
assuming that the mode shapes have been scaled to give dimensional so ultimately depend on the values of αi . The grid speeds can be ob-
generalized masses m i = 1, the scalar equations tained by differentiating Eq. (12) to obtain their explicit dependence
on the values of α̇i .
d2 αi dαi
+ Di + ωi2 αi = µφiT fs (9) Formulation of Hopf Analysis
dt 2 dt
The semidiscrete form of the coupled CFD–CSD system is
are obtained where fs is the vector of aerodynamic forces at the
structural grid points and Di is the coefficient of structural damping. dw
= R(w, µ) (13)
Here a nondimensionalization consistent with the flow solver has dt
been used and the factor µ = ρ∞ /ρw is a density ratio, where ρw is
the density of the wing structure. These equations are rewritten as a where
system in the form w = [w f , ws ]T (14)
dws
= Rs (10) is a vector containing the fluid unknowns w f and the structural
dt unknowns ws and
where ws = (. . . , αi , α̇i , . . .)T and Rs = (. . . , α̇i , µφiT fs − ωi2 αi − R = [R f , Rs ]T (15)
Di α̇i , . . .)T . This equation can be solved by a two-stage Runge–
Kutta method, which requires a knowledge of fsn and fsn + 1 . To avoid is a vector containing the fluid residual R f and the structural resid-
introducing sequencing errors by approximating the value of fsn + 1 , ual Rs . The residual also depends on a parameter µ, which is
the Runge–Kutta solution is iterated in pseudotime along with the independent of w. An equilibrium of this system w0 (µ) satisfies
CFD solver, with the latest pseudoiterate being used to give a value R(w0 , µ) = 0.
 BADCOCK, WOODGATE, AND RICHARDS 733

Dynamical systems theory gives criteria for an equilibrium to be The linear system at each Newton step is solved using the sparse
stable (e.g., see Seydel17 ). In particular, all eigenvalues of the Jaco- matrix package Aztec (see Tuminaro et al.18 ). Although not opti-
bian matrix of Eq. (13), given by A = ∂R/∂w, must have negative mized for the current problem, the generality of the package has
real part. A Hopf bifurcation with respect to the parameter µ occurs allowed various experiments to be carried out. This package has
in the stability of the equilibrium at values of µ such that A(w0 , µ) three main solvers available, namely, GMRES, CGS, and TFQMR,
has one eigenvalue iω that crosses the imaginary axis. Denoting the although the differences in performance for the current problem
corresponding eigenvector by P = P1 + iP2 , a critical value of µ is were found to be small. The key issue for iterative linear solvers is
one at which there is an eigenpair ω and P such that usually the preconditioner. The incomplete LU factorization family
as described in Axelsson19 can be very effective at approximating the
AP = iωP. (16) inverse of the coefficient matrix with a small number of terms. For
CFD calculations, block incomplete upper lower factorizations with
This equation can be written in terms of real and imaginary parts as no fill in have proved very successful, as illustrated in Badcock et al.7
AP1 + ωP2 = 0 and AP2 − ωP1 = 0. A unique eigenvector is chosen One simplification arises if we are dealing with a symmetric prob-
by scaling against a constant real vector q to produce a fixed complex lem, in which case the equilibrium solution w0 is independent of µ
value, taken to be 0 + 1i. This yields two additional scalar equations and hence can be calculated from Eq. (13) independently of the
qT P1 = 0 and qT P2 − 1 = 0. other Hopf conditions in Eq. (17). Then, a smaller system can be
A bifurcation point can be calculated directly by solving the sys- solved for the bifurcation parameter. This algorithm has been im-
tem of equations plemented in a code, which will be referred to as the BIFOR code.
R A (w A ) = 0 The formulation of the coupled CFD–CSD problem follows that of
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

(17)
the PMB code already given.
where
  Calculation of the Jacobian Matrix
R
The difficult terms to form in the Jacobian matrix of the aug-
 AP1 + ωP2 
  mented system are A and Aµ . The calculation of A is most conve-
RA = 
 AP2 − ωP1 
 (18) niently done by partitioning the matrix as
 qT P1   
∂R f ∂R f
qT P2 − 1  
 ∂w f ∂ws 
A=  = Aff Afs
(21)
and w A = [w, P1 , P2 , µ, ω]T . If there are n components in w, then  ∂Rs ∂Rs  Asf Ass
w A has 3n + 2 components, as does R A , and hence Eq. (17) is closed. ∂w f ∂ws
The catch is that this is a large sparse system of nonlinear equations.
Newton’s method can be used to solve this type of problem. A The block Aff describes the influence of the fluid unknowns on
sequence of approximations wnA to a solution is generated by solving the fluid residual and has by far the largest number of nonzeros
the linear system when a modal structural model is used. The treatment of this term is
crucial to the efficiency of the scheme and is discussed in Badcock
∂R A et al.6 To drive the Newton iteration to convergence, the analytical
w A = −RnA (19)
∂w A Jacobian corresponding to the first-order spatial scheme is used.
This approach has proved successful for CFD-alone calculations
where w A = wnA+ 1 −wnA . The Jacobian matrix on the left-hand side shown in Badcock et al.7
of Eq. (19) is given in expanded form as The only consideration when choosing the Newton iteration ma-
  trix is that the scheme converges to the correct answer, the latter
A 0 0 Rµ 0
being determined by the calculation of the residual on the right-
(AP1 )w A Iω (AP1 )µ P2 
∂R A   hand side of the Newton iteration. Hence, the products AP1 and
=
(AP2 )w −I ω A (AP2 )µ −P1 
 (20) AP2 must be computed exactly. This can be done using a matrix-
∂w A  0 qT 0 0 0  free formulation of the form
0 0 qT 0 0 R(w + hx) − R(w − hx)
Ax ≈ (22)
2h
There are three key issues for the application of Eq. (19). First, a
good initial guess is required, or the iterations will diverge. Secondly, where x denotes the real or imaginary part of the critical eigenvalue
the Jacobian matrix ∂R A /∂w A is required. Thirdly, the large sparse and h is the increment applied. Computing this expression is not
linear system given in Eq. (19) must be solved. These points were costly as it requires only two residual evaluations. This gives an
considered for the aerofoil problem in Badcock et al.6 For the three- accurate approximation to the required product without having to
dimensional problem with a modal structural model, the Jacobian evaluate and store the exact A. Using the matrix-free evaluation of
calculation is the aspect which is different from the aerofoil case. the augmented residual reduces the memory requirements for the
This is therefore considered in the next section. The details of the scheme overall and simplifies the code considerably.
first and third points are as described earlier for the aerofoil case in The dependence of the fluid residual on the structural unknowns
Badcock et al.6 and are only summarized here. αi and α̇i is partially hidden by the notation used. The fluid residual
First, for the application of the scheme it is assumed that a good depends not only on the fluid cell values but also on the location
estimate of the flutter point and frequency is available from some of the grid points themselves and the cell volumes. The fluid and
other source, for example, linear theory, at the first Mach number structural unknowns are independent variables, and hence to calcu-
of interest. The inverse power method is then applied, again using late the term Afs the fluid unknowns are kept fixed. The influence
the sparse matrix formulation, to calculate the eigenvector corre- of the structural unknowns is felt through the moving grid. Using
sponding to the critical eigenvalue when the fluid Jacobian from the modal structural model, the updated grid locations and speeds
the first-order spatial scheme is used. This information is then used are calculated by moving the structural grid according to the values
as the starting solution for the Hopf calculation at the first Mach of the generalized coordinates and velocities, transferring these to
number, which is initially converged using the first-order fluid Ja- the fluid surface grid using the transformation and then applying
cobian before switching to the Jacobian of the full second-order TFI to transfer these boundary values to the volume grid. By using
spatial scheme. For subsequent Mach numbers, the converged solu- second-order finite differences, the terms in Afs can be calculated in
tion from the previous one provides an adequate starting solution. In 2n s evaluations of the aerodynamic residual if there are n s structural
this way the flutter boundary is traced for a range of Mach numbers. unknowns.
 734 BADCOCK, WOODGATE, AND RICHARDS

The term Asf involves calculating the dependence of the general-

ized fluid forces on the fluid unknowns. The surface forces on the
aerodynamic grid are calculated and then transferred to the struc-
tural grid using the transformation. The inner product is then formed
using the forces on the structural grid and the modal coefficients.
The Jacobian matrix for the forces on the structural grid with respect
to the fluid unknowns can be calculated first analytically. Then the
required terms for Asf can be calculated through matrix-vector mul-
tiplication.
Finally, the exact Jacobian matrix for the dependence of the struc-
tural equations on the structural unknowns is easy to calculate ana-
lytically. However, it is important to remember that the generalized
force will change with the structural unknowns because the sur-
face normals to the wing will change as the wing moves. A finite
difference calculation is used to calculate this effect.
The bifurcation parameter (µ in this case) only appears in the
structural equations. Therefore, for this case,
  a)
0 0
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

Aµ =  ∂ 2 Rs  (23)
0
∂µ∂ws

Because of the simple algebraic expression for ∂Rs /∂ws , it is

straightforward to calculate the required term analytically.

Results for Symmetric Problem

Test Case
An important problem with the development of aeroelastic simu-
lation tools is the lack of experimental data available for assessment.
The experiments are intrinsically destructive and require careful
model construction to ensure proper scaling, and hence the expense
is much higher than rigid model tests. A complete set of measure-
ments is available for the AGARD 445.6 wing, and results have
been included in most papers on computational aeroelasticity, giv-
ing a wide range of data with which to evaluate the current method.
However, a disadvantage of this test case is that it does not feature b)
significant nonlinear effects because the wing is thin. Despite this, it Fig. 1 Only the inner blocks above the wing are shown on the symme-
is commonly the first test case used to test time-marching codes and try plane: a) grid topology and b) medium surface mesh.
is suitable for the current work because it is symmetric. Previous
time-marching results are reviewed by Goura et al.20
The AGARD 445.6 wing is made of mahogany and has a 45-deg the published results is that in the majority of cases no structural
quarter-chord sweep, a root chord of 22.96 in. (0.583 m), and a damping was used. In the description of the experiment, a value of
constant NACA64A004 symmetric profile (see Yates8 ). A series of 2% is suggested by Yates,8 although it is not clear how certain this
flutter tests, which were carried out at the NASA Langley Transonic value is.
Dynamics Tunnel to determine stability characteristics, was reported
in 1963. Various wing models were tested (and broken). The case for Time-Marching Solutions
which most published results have appeared is the weakened wing An attempt was made to do a detailed grid-convergence study
(wing 3) in air. This wing had holes drilled out, which were filled within the limits of the computers available. All calculations re-
with plastic to maintain the aerodynamic shape while being struc- ported in this section were done with the PMB code. To optimize
turally weaker. Published experimental data include the dynamic the grid used, two requirements were identified. First, because the
conditions at which the wing was viewed to be unstable for Mach calculations are inviscid, and hence no wake needs to be preserved,
numbers in the range 0.338 to 1.141. The structural characteristics of an O grid was used, which helps to maximize the number of grid
the wing were provided in the form of measured natural frequencies points on the surface of the wing. A genuine multiblock topology
and mode shapes derived from a finite element model. Full details was used to allow a good quality mesh to be preserved in the tip
of the structural model used are given in Goura.13 Four modes are leading- and trailing-edge regions as shown in Fig. 1. Second, the
retained with the first two bending modes having frequencies 9.7 important quantities that must be well predicted for the flutter cal-
and 50.3 and the two torsional modes at 36.9 and 90.0. culations are the generalized modal forces. The pressure difference
A problem with the published results for the AGARD 445.6 wing between the upper and lower surface therefore must be predicted
is that most are of a demonstration nature in the sense that verifica- accurately. The flutter response is dominated by the first bending
tion is hardly ever shown. There is a significant spread of the results mode, which features some twist near the tip but is essentially a
obtained when using solutions of the Euler equations. The results plunging motion near the root. The significant pressure difference,
that cluster around the measured data in the region of the flutter and following from this the main contribution to the generalized
dip tend to be on coarser grids, with finer CFD grids generally giv- force, comes from the region towards the wing tip. The grid points
ing lower flutter speeds. The only published attempt at a systematic were therefore concentrated in this region. Most structured grids
grid refinement study was shown by Melville and Gordnier.21 In this shown in the literature have been of the C-H topology and are rea-
commendable study the fine- and medium-grid results were further sonably uniform in the spanwise direction.
apart than those on the medium and coarse grids, and hence grid The finest grid in the current work has 1.43 million points and
independence was not achieved, casting doubt on other published 17,700 on the wing surface. Medium and coarse grids that have
results on coarser grids. The main obstacle to a rigorous study is of 190,000 and 27,000 points, respectively, with 4453 and 1131 points
course the cost of the calculations. A second question mark against on the wing surface were extracted from this. The number of volume
 BADCOCK, WOODGATE, AND RICHARDS 735

points used in the refinement study of Gordnier and Melville,21

where the fine, medium and coarse grids have 2.1 million, 901,000,
and 274,000 points, respectively, are comparable to the current grids,
but significantly the number of surface points are less (9231, 5343,
and 2366, respectively). It is stressed however that the topology in
the current study would not be ideal for RANS calculations, which
were the main focus of Gordnier and Melville.21 Comparison is also
made with the structured22 and unstructured23 studies by Batina and
co-workers. The structured grid has 517,000 points with 5289 on
the wing surface. The unstructured grid has 129,000 tetrahedra, and,
although no information is given about the number of points on the
wing surface, the pictures shown in the paper suggest that the grid
points are more strongly clustered in the wing region than is possible
for structured grids.
A number of tests using the medium grid were made on the tempo-
ral parameters (time step and convergence level) at a Mach number
of 0.96, a freestream velocity of 308 m/s, a density of 0.08 kg/m3 ,
and structural damping of 0.5%. First, the responses when using 10
or 20 pseudosteps per real time step were identical, indicating that
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

10 steps was more than adequate. Secondly, using a reduced time

step of 0.01 and 0.02 also gave an identical response, and hence the Fig. 2 Variation of flutter speed index with the structural damping
larger time step is adequate. applied at Mach 0.96: ——, measurement.
The influence of grid resolution is harder to test because of the
calculation cost on the density of grids that are required. The three
grid levels were used to locate the flutter point for a range of Mach
numbers. Two calculations were run at different values of freestream
density for each Mach number and the growth of the response cal-
culated using the approach of Gordnier and Melville,21 where the
ratio of consecutive peaks was taken to define an amplification fac-
tor. Linear interpolation of the amplification was then used to esti-
mate the value of density at which a neutrally stable response would
be obtained. The medium grid calculations took about 5–6 h on a
2.5-GHz PC to calculate five periods of the flutter response.
The comparison at Mach 0.96, which is close to the bottom of
the flutter dip, between the predicted flutter speeds on the three
grid levels with other predictions is shown in Table 1. The results
of Gordnier and Melville21 show a downward (and accelerating)
trend in the flutter speed, with bigger differences between the fine
and medium than between the medium and coarse. The value from
Rausch et al.23 is lower still. The current results suggest that the
medium grid provides good spatial resolution, with the three grids
showing behavior consistent with spatial convergence. The grid-
converged value using no structural damping is much lower than Fig. 3 Comparison with measurements of the stability boundaries cal-
experiment. Using a value of structural damping of 2% shifts the culated on the medium grid using time-marching and the bifurcation
flutter speed index above the experimental value. solver.
The trends from these results suggest that the grid-converged so- bifurcation code has been written, and the main consequence of this
lution without structural damping is significantly below the exper- is that the medium grid is the largest problem that can be tackled
imental result. Adding structural damping brings the flutter speed on the computers that were available. This case requires 1.5 Gb of
back up into the range of the measurements, as shown in Fig. 2. memory. The PMB and BIFOR codes are based on similar formu-
The solution obtained using 0.5% structural damping is in good lations for the steady-state solutions. However, the PMB code has
agreement with the experimental values. been under development for 10 years and so has reached a level of
optimization well in advance of the newer BIFOR code. To quantify
Bifurcation Results this, an equivalent implicit step with the BIFOR code requires 2.35
The bifurcation solver was applied on the coarse and medium times the CPU time when compared with the PMB code. This will
grids using the BIFOR code. Currently no parallel version of the be accounted for in the comparisons between the bifurcation and
time-marching calculation costs by expressing these as multiples of
Table 1 Grid-refinement influence on flutter speed index
the CPU time required for a steady-state calculation with the same
at Mach 0.96
code. The timings are likely to be conservative when assessing the
Grid Grid Flutter speed performance of the bifurcation method because additional gains are
Reference volume surface Damping index likely from writing a dedicated linear solver (i.e., one that is not
Current Coarse 1,131 0 0.227 configured to handle general-sized matrix blocks).
Current Medium 4,453 0 0.192 Guided by the time-marching results, a value of structural damp-
Current Fine 17,700 0 0.175 ing of 0.5% was used. The comparison on the medium grid between
Current Coarse 1,131 2% 0.401 the measured, time-marching, and bifurcation results is shown for
Current Medium 4,453 2% 0.381 the dip region in Fig. 3. The bifurcation boundary was computed
Current Fine 17,700 2% 0.375 first for eight Mach numbers between 1.07 and 0.67, with a Mach
21 Coarse 2,366 0 0.314 number interval of 0.05, and subsequently for 12 Mach numbers in
21 Medium 5,340 0 0.304
the dip region, with an interval of 0.01. The frequency from the time-
21 Fine 9,231 0 0.285
23 Unstructured N/A 0 0.230 marching calculations was used with the inverse power method at the
22 Structured 5,289 0 0.294 largest Mach number to initiate the calculations. Good agreement
between the predictions of the two codes is observed as required.
 736 BADCOCK, WOODGATE, AND RICHARDS

Table 2 Average calculation cost using the PMB code for the first numbers is less than half the cost of a single time-marching calcu-
two rows and the BIFOR code for the bottom two rows in the table lation. Considering the number of time-marching calculations re-
(The relative costs have been scaled by the time for a steady-state quired to trace out a flutter boundary, the calculation cost can be re-
calculation with the appropriate code.) duced by two orders of magnitude by using the bifurcation method.
Calculation type CPU, s CPU (relative) One concern with the previous work on aerofoil problems was
that the performance of the linear solver would deteriorate for larger
Steady state 787 1 problems. On average for the aerofoil 30 iterations were required
Unsteady solution 19,810 45 to achieve a reduction of two orders in the residual. The number of
Steady calculation 1,767 1
Bifurcation calculation 3,304 1.87
linear solver iterations at each bifurcation iteration is shown in the
scatter plot in Fig. 5 along with the average number of iterations
required for the previous aerofoil calculations. The fact that the
number of linear solver iterations is spread about the average two-
dimensional cost indicates that the performance of the Krylov solver
has been maintained for the larger three-dimensional problems.

Conclusions
The performance of a bifurcation method for calculating flutter
boundaries has been evaluated for the AGARD 445.6 wing test case.
This is believed to be the first time that such a method has been used
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

to calculate a three-dimensional aeroelastic instability problem. The

good performance of the method previously observed for aerofoil
problems has been preserved. In particular the cost of the iterative
linear solver in terms of the number of iterations required has not
increased as the size of the matrix has increased. It is estimated that
the cost of tracing out a flutter boundary over 10 Mach numbers
has been reduced by two orders of magnitude compared with the
time-marching method.
The method presented relied on the system being symmetric,
which means that the equilibrium solution is independent of the bi-
furcation parameter. In the asymmetric case the equilibrium solution
has to be recalculated as the bifurcation parameter is updated during
Fig. 4 Convergence of the bifurcation parameter with the bifurcation the Newton iterations. Two approaches are being contemplated to
solver iteration number at Mach 0.96. The flutter speed index is shown realize this. First, the problem could be solved in a completely cou-
against the right-hand axes and the augmented residual against the left- pled way, with the Jacobian matrix driving the Newton convergence
hand axes.
given by Eq. (20). This has the likely advantage of optimizing the
convergence rate at a cost of having to store and use a matrix that
is around 50% larger. The alternative is to decouple the equilibrium
and bifurcation calculations, which in effect ignores the terms in the
first row, fourth and first columns, and third and fourth rows of the
matrix in Eq. (20). This approach will probably degrade the conver-
gence but has the benefit of reducing the storage requirements as the
equilibrium and bifurcation calculations can be sequenced. Further
work is required to evaluate the merits of these approaches.
The next stage of this work is to implement a parallel version of
the code with the linear solver tailored for the particular block sizes
encountered for this problem. This will allow the goal of assessing
the method for full aircraft flutter calculations to be realized.

Acknowledgments
This work was supported by Engineering and Physical Sciences
Research Council, Ministry of Defence, and BAE Systems and
forms part of the work programme of the Partnership for Unsteady
Methods Defence and Aerospace Research Partnership (PUMA
DARP).
Fig. 5 Number of linear solver steps per bifurcation solver iteration:
——, average number of steps for an aerofoil calculation.6
References
An assessment of the relative cost of the time-marching and bi- 1 Farhat, C., Geuzaine, P., and Brown, G., “Application of a Three-Field
furcation calculations can be made from the times in Table 2. Here Nonlinear Fluid-Structure Formulation to the Prediction of the Aeroelastic
the average CPU time of an unsteady calculation and the average Parameters of an F-16 Fighter,” Computers and Fluids, Vol. 32, No. 3, 2002,
bifurcation cost for each Mach number have been expressed in mul- pp. 3–29.
2 Melville, R., “Nonlinear Simulation of F-16 Aeroelastic Instability,”
tiples of the cost of a steady-state calculation using the respective
codes. The steady-state calculation in each case has been converged AIAA Paper 2001-0570, Jan. 2001.
3 Morton, S. A., and Beran, P. S., “Hopf-Bifurcation Analysis of Airfoil
five orders of magnitude. The bifurcation solver has been converged
to at least three significant figures for the flutter speed, as indicated Flutter at Transonic Speeds,” Journal of Aircraft, Vol. 36, 1999, pp. 421–429.
4 Beran, P. S., and Carlson, C. D., “Domain-Decomposition Methods for
in Fig. 4. The stopping criterion is based on reducing the magnitude
Bifurcation Analysis,” AIAA Paper 97-0518, 1997.
of the augmented residual, defined by Eq. (18), by one order from 5 Beran, P. S., “A Domain-Decomposition Method for Airfoil Flutter Anal-
the starting value. The time for the unsteady calculations is based ysis,” AIAA Paper 98-0098, 1998.
on 750 time steps resolving five cycles. 6 Badcock, K. J., Woodgate, M. A., and Richards, B. E., “Hopf Bifurcation
Similar conclusions to the previous aerofoil study can be drawn. Calculations for a Symmetric Airfoil in Transonic Flow,” AIAA Journal,
The time required to trace out the flutter boundary for eight Mach Vol. 42, No. 5, 2004, pp. 883–892.
 BADCOCK, WOODGATE, AND RICHARDS 737

7 Badcock, K. J., Richards, B. E., and Woodgate, M. A., “Elements of 15 Tsai, H. M., Wong, A. S. F., Cai, J., Zhu, Y., and Liu, F., “Unsteady
Computational Fluid Dynamics on Block Structured Grids Using Implicit Flow Calculations with a Parallel Multiblock Moving Mesh Algorithm,”
Solvers,” Progress in Aerospace Sciences, Vol. 36, 2000, pp. 351–392. AIAA Journal, Vol. 39, No. 6, 2001, pp. 1021–1029.
8 Yates, E. C., “AGARD Standard Aeroelastic Configurations for Dynamic 16 Gordon, W. J., and Hall, C. A., “Construction of Curvilinear Coordinate
Response 1: Wing 445,6,” AGARD, Rept. 765, 1988. Systems and Applications to Mesh Generation,” International Journal for
9 Osher, S., and Chakravarthy, S. R., “Upwind Schemes and Boundary Numerical Methods in Engineering, Vol. 7, 1973, pp. 461–477.
Conditions with Applications to Euler Equations in General Coordinates,” 17 Seydel, R., Practical Bifurcation Analysis and Stability Analysis, 2nd
Journal of Computational Physics, Vol. 50, 1983, pp. 447–481. ed., Springer-Verlag, Berlin, 1994.
10 Van Leer, B., “Towards the Ultimate Conservative Conservative Dif- 18 Tuminaro, R. S., Heroux, M., Hutchinson, S. A., and Shahid, J. N.,
ference Scheme II: Monotonicity and Conservation Combined in a Sec- “Official Aztec User’s Guide,” Ver. 2.1, Sandia Lab., SAND99-8801J, Al-
ond Order Scheme,” Journal of Computational Physics, Vol. 14, 1974, buquerque, NM, 1999.
pp. 361–374. 19 Axelsson, O., Iterative Solution Methods, Cambridge Univ. Press, Cam-
11 Cantariti, F., Dubuc, L., Gribben, B., Woodgate, M., Badcock, K., and bridge, UK, 1994.
Richards, B., “Approximate Jacobians for the Solution of the Euler and 20 Goura, G. S. L., Badcock, K. J., Woodgate, M. A., and Richards, B. E.,
Navier–Stokes Equations,” Univ. of Glasgow, Aerospace Engineering Rept. “Implicit Method for the Time Marching Analysis of Flutter,” Aeronautical
9704, Scotland, U.K., 1997. Journal, Vol. 105, No. 1046, April 2001, pp. 192–214.
12 Jameson, A., “Time Dependent Calculations Using Multigrid, with Ap- 21 Gordnier, R. E., and Melville, R. B., “Transonic Flutter Simulations
plications to Unsteady Flows Past Airfoils and Wings,” AIAA Paper 91-1596, Using an Implicit Aeroelastic Solver,” Journal of Aircraft, Vol. 37, No. 5,
1991. 2000, pp. 872–879.
13 Goura, G. S. L., “Time Marching Analysis of Flutter Using Computa- 22 Lee-Rausch, E. M., and Batina, J. T., “Wing Flutter Computations Using
tional Fluid Dynamics,” Ph.D. Dissertation, Dept. of Aerospace Engineering, an Aerodynamic Model Based on the Navier–Stokes Equations,” Journal of
Downloaded by UNIVERSITY OF OKLAHOMA on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/1.5323

Univ. of Glasgow, Scotland, U.K., Nov. 2001. Aircraft, Vol. 33, No. 6, 1996, pp. 1139–1147.
14 Goura, G. S. L., Badcock, K. J., Woodgate, M. A., and Richards, B. E., 23 Rausch, R. D., Batina, J. T., and Yang, H. T. Y., “Three-Dimensional
“Extrapolation Effects on Coupled CFD-CSD Simulations,” AIAA Journal, Time Marching Aeroelastic Analyses Using an Unstructured-Grid Euler
Vol. 41, No. 2, 2003, pp. 312–314. Method,” AIAA Journal, Vol. 31, No. 9, 1993, pp. 1626–1633.