Time Integration Schemes for the Unsteady Navier-Stokes Equations

13 Pages·2001·0.87 MB·English

by  Bijl, Hester

j
W •
AIAA-2001-2612
Time Integration Schemes for the Unsteady
Navier-Stokes Equations
HesterBijl
Delft Universityof Technology
Delft,The Netherlands
Mark H.Carpenter
NASA LangleyResearchCenter
Hampton,VA
Veer N. Vatsa
NASA LangleyResearchCenter
Hampton,VA
15th Computational Fluid Dynamics Conference
June 11-14, 2001 / Anaheim, CA
For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics
1801 Alexander Bell Drive, Suite 500, Reston, VA 20191-4344
AIAA 2001-2612
TIME INTEGRATION SCHEMES FOR THE UNSTEADY NAVIER-STOKES
EQUATIONS
Hester Bijl, *Mark H. Carpenter t and Veer N. Vatsa t
Abstract schemes ranging from first to fifth order. A similar
study was performed for only first and second order
The efficiency and accuracy of several time inte-
time integration methods by Marx u
gration schemes are investigated for the unsteady
Implicit in any comparison of efficiency is a precise
Navier-Stokes equations. This study focuses on the
error tolerance requirement. Unfortunately, seldom
efficiency of higher-order Runge-Kutta schemes in
does a single scheme prove to be optimal over a wide
comparison with the popular Backward Differenc-
range of solution error tolerances. It is well known
ing Formulations. For this comparison an unsteady
that high-order schemes (fourth-, fifth-, etc.) out-
two-dimensional laminar flow problem is chosen, i.e.
perform low order schemes for error tolerances that
flow around a circular cylinder at Re=1200. It is
are small (_ 10-7). For example, Kennedy s com-
concluded that for realistic error tolerances (smaller
pares explicit third-, fourth- and fifth-order Runge-
than 10-1) fourth- and fifth-order Runge Kutta
schemes are the most efficient. For reasons of robust- Kutta (RK) schemes on Direct Navier-Stokes simu-
lations (DNS), and determines the optimal order for
ness and computer storage, the fourth-order Runge-
a given temporal error tolerance. For DNS it was
Kutta method is recommended. The efficiency of the
found, the fourth-order methods are optimal over a
fourth-order Runge-Kutta scheme exceeds that of
surprisingly broad range of error tolerances, and are
second-order Backward Difference Formula (BDF2)
competitive at large error tolerances as well.
by a factor of 2.5 at engineering error tolerance lev-
The hallmark of large-scale aerodynamics calcula-
els (10-1-10-2). Efficiency gains are more dramatic
tions is that they seldom require small error toler-
at smaller tolerances.
ances. Calculations that are accurate to one or two
Introduction significant digits, which translates into an error tol-
erance of (10-L10-2), are frequently sufficient. For
Due to constraints of computing costs, the devel- these aerodynamics calculations, the second-order
opment of numerical techniques for fluid flow simula- accurate BDF2 scheme is currently the method of
tions in the past has focused mainly on steady state choice. There is little question that calculations re-
calculations. However, many physical phenomena quiring low error tolerances (10-4-10 -5) will be well
of interest are inherently unsteady; a %w examples suited for fourth-order RK formulations. The cen-
being separated flows, wake flows and buffet, fluid tral question of this study, however, is the feasibility
actuators and maneuvering. With the continuous re- of using fourth-order RK formulations for simula-
duction of computer costs recently more attention is tions requiring error tolerances of (10-L10-2).
devoted to the simulation of these unsteady flows 1,4. A production aerodynamics solver is needed for
However, the need for further reduction of computer this study. For this, the extensively tested and well
time for unsteady flow computations is still appar- documented solver of Vatsa (TLNS3D) is is chosen.
ent. Therefore, in this paper we investigate possible This multi-block structured grid solver is representa-
reductions in computer time due to the choice of tive of a broad class of commonly used solvers. The
an efficient time integration scheme from a series of TLNS3D (Thin-Layer Navier-Stokes 3-Dimensions)
code utilizes a special form of the unsteady thin-
*Pro%ssor, Department of Aerospace Engineering, Delft
layer Navier-Stokes equations. The spatial terms
University ofTechnology, Delft, The Netherlands
?Senior Research Scientist, Computational Methods and are discretized using a conventional cell-centered fi-
Simulation Branch, NASA Langley Research Center, Hamp- nite volume scheme with artificial dissipation added
ton, VA 23681-0001 for stability. Time is discretized in a fully implicit
Copyright @2001 bythe American Institute of Aeronau-
sense using both multistep BDF and multistage RK
tics and Astronautics, Inc. No copyright is asserted in the
United States under Title 17, U.S. Code. The U.S. Govern- schemes. The resultant nonlinear algebraic equa-
ment has aroyalty-free license to exercise all rights under the tions are solved iteratively in pseudo-time with a
copyright claimed herein for government purposes. Allother multi-grid acceleration used to speed up the conver-
rights are reserved by the copyright owner.
genceto(pseudo-timset)eady-state. and oscillations in the vicinity of shock waves and
Inthispaper,first,thegoverninfglowequations stagnation points 18. The viscous terms are cen-
andthespacediscretizatioanregiven.Thereafter, tral differenced with second-order formulas. The
the timediscretizationtechniqueesmployedi,.e.
turbulence models used are Baldwin-Lomax 2, and
severaml ultistepBackwardDifferencinFgormula-
Spalart-Allmaras 15
tionsandmultistageRunge-Kuttaschemeasr,eex-
tensiveldyiscussedT.hisisfollowedbyadescription Time Discretization
ofthesolutionalgorithmT. hatis,theimplicittime
Consider the integration of the system of ordi-
integrationoftheflowequationasndtheiterativeal-
nary differential equations (ODE's) represented by
gorithmforsolvingtheresultingnon-lineairmplicit
the equation
equationT.hroughapseudo-timsetabilityanalysis
theexactnatureofthepseudo-timseub-iterationiss dU
_ s(t, u(t)) .
analyzedA.ftervalidationofthespace-timmeethod dt
forunsteadylaminarflowaroundacircularcylin-
The vector S in our case, results from the semi-
der,accuracayndefficiencoyfthetimeintegration
schemeissdiscussed. discretization of the equations of fluid mechanics
plus a suitable turbulence model. Inclusion of the
turbulence model enables the simulation of high
Governing Equations
Reynolds number flows in excess of 107, but indi-
In the present work, a modified version of the thin- rectly introduces extremely fine length-scales in the
layer Navier-Stokes equations is used to model the near-wall regions. Near wall stiffness in the range of
flow. The equation set is obtained from the com- 103- 104 is not uncommon in practical engineering
plete Navier-Stokes equations by retaining only the problems, and increases with Reynolds number. It
viscous diffusion terms normal to the solid surfaces. is imperative that any efficient solver of these equa-
For a body-fitted coordinate system (_, _1,_) fixed in tions be able to maintain stability at arbitrarily large
time, these equations can be written in the conser- time steps, thus allowing the potential to "step over"
vative form as: unimportant boundary layer time scales.
The term "stiffness" is difficult to define. A prac-
O(U) c_(F - Fv)
+ + tical definition for the purpose of this work, com-
0t
pares an implicit timestep At I with a baseline ex-
0(G - Gv) c_(H - Hv) plicit timestep AtE. The implicit timestep is the
+ - 0, (1)
&J oc maximum allowed by accuracy considerations, and
the baseline explicit timestep is obtained from sta-
where U represents a combination of the transforma-
bility consideration. All other variables are identi-
tion Jacobian J and the conserved variable vector.
cal in the comparison. The stiffness is the ratio of
The vectors E, G, H, and Ev, Gv, Hv represent the the two timesteps. The implicit scheme is capable of
convective and diffusive fluxes in the three trans-
stepping over large negative eigenvalues that are not
formed coordinate directions, respectively. These important for solution accuracy, while the explicit
equations represent a generalized form of the classi- scheme must bound them in the finite stability en-
cal thin-Layer Navier-Stokes equations and include velop.
all normal components of the viscous stress terms. There are two mathematical properties that all
The TLNS3D computer code, is used in this study
candidate numerical integrators should possess. The
to solve equation (1). Many references exist detail-
first (and most important) is the "A-stability" prop-
ing the discretization and implementational issues of erty which guarantees that all eigenvalues lying in
TLNS3D. We include only a brief summary of the
the left half of the complex plane (LHP) will have
general features, and refer to the work of Vatsa 18
an amplification of no more than 1, independent of
for further details.
the chosen step size. The only restriction on the
time-step with an A-stable scheme is the consider-
Space Discretization
ation of solution accuracy. The second is the "L-
The spatial terms in Equation (1) are discretized stability" property which guarantees that eigenval-
using a standard cell-centered finite volume scheme. ues approaching -oc are damped in one time-step.
The convection terms are discretized with second- These spurious eigenvalues are generated by high fre-
order central differences with scalar/matrix artifi- quency information in the spatial discretization, and
cial dissipation (second- and fourth- difference dissi- by incomplete solution of the non-linear system at
pation) added to suppress the odd-even decoupling each time-step. (The nonlinear system is never con-
vergedtothelevelsofmachineprecision)I.f these
Method Steps Order L(c 0 flk
spuriousmodesarenotstronglydampedtheycan 1
BDF1 1 1 L(90 °)
buildupandcauseinstability.
I BDF2 I 2 I 2 I L(90°) I _- I
Mostnumericamlethodssuitableforthesecompu-
I BDF3 I 3 I 3 I L(86.03°) I _ I
tationscanbecategorizeindtotwobroadclasses1:)
multistepa,nd2)multistagee,achhavingitsadvan- I BDF4 I 4 I 4 I L(73.35o) I ,2
tagesanddisadvantageTsh.ecurrent"methodsof I BDF5 I 5 I 5 I L(51.84°) I _ I
choicei"nthecomputatioonflargescaleengineering I BDF6I 6 I 6 IL(lr.S4°)l_ I
flowsarethemultistepBDFformulasa,ndinpar-
ticulartheBDF2schemeT. heseschemeaschieve Table 1. Properties of the BDF methods.
greatefficiencybecausetheysolveonlyonenon- The BDF1 and BDF2 schemes are L-stable for eigen-
linearsetof equationspertime-step.Theysuf- values anywhere in the left half of the complex plane.
fer,howeverfr,omnotbeingself-startinga,redif- Beyond second order the BDF schemes are L(_)-
ficulttousewithvariabletimesteps,andarenot stable, with the _ being defined in Table (1). This
A-stablebeyondsecond-ordetermporalaccuracy. definition of stability implies that the schemes will
Multi-stageRunge-Kuttaschemersequiremultiple be stable for eigenvalues lying in the wedge bounded
nonlineasrolvespertime-stepb,ut areselfstart- above and below by the complex lines -t-_. For fur-
ing,areeasilyimplementeindvariabletime-stepping ther details see the work of Hairer 6
modea, ndcanbedesignewdithA-and The
L-stability
Runge-Kutta methods are multistage schemes and
properties for any temporal order.
are implemented as
Practical experience indicates that large scale en-
gineering computations are seldom stable if run with
k
BDF4 12. The BDF3 scheme is often stable, but di- U k = u° + (At)Za_js (uJ), k = 1,,
verges for certain problems and some spatial opera- j=,
tors. These failures result from portions of the Left
Half Plane that are not stable for arbitrarily large u_+_ = u_ + (At)_ bjs(uJ) (a)
timesteps. Thus, a reasonable practitioner uses the
BDF2 scheme exclusively for large scale computa-
tions. The essential question this paper seeks to 0_+1 = u_ + (At)_ bjs(uJ),
address is whether high-order RK schemes can be
j=l
designed which are more efficient than the BDF2
schemes, and if so, what is the optimal order of the where s is the number of stages and aii and b/ are
the Butcher coefficients of the scheme. The vectors
RK scheme.
The general formula for a k-step BDF scheme can U, and U are the pth-order and (p- 1)th;order so-
lutions at time level n + 1. The vector U is used
be written as
solely for estimating error and is obtained with little
additional computational overhead. In this work we
choose to focus on the ESDIRK class of RK schemes,
which refers to the Explicit first stage, Single diag-
k-1
onal coefficient, Diagonally Implicit Runge-Kutta.
u(_+_) = - Z _u(_+° + (At)V_s(_+_)(2)
The Butcher tableau for these schemes (here repre-
i=0
sented with s = 5) takes the form
0 0 0 0 0 0
C2 a21 akk 0 0 0
The BDF formulas involve at each time-step the
C3 a31 a32 akk 0 0
storage of k + 1levels of the solution vector U, and
C4 a41 a42 a43 akk 0
the implicit solution of one set of nonlinear equa- 1 bl b2 b3 b4 akk
tions. The BDF schemes are subject to the famous bl b2 b3 b4 akk
Dahlquist barrier, which proves that A-stability is
impossible for linear multistep schemes beyond sec-
ond order. Table (1) presents the stability prop- where ci denotes the point in the time interval
erties and diagonal contribution flk of several BDF t --+ t+At where the intermediate stage is evaluated,
schemes. and bj are the coefficients used for the embedded er-
ror predictor. The fully implicit RK schemes are not
pursuedowingtothecomplexityoftheirimplemen-
form. Thus, the pseudo-time term _ is added to
tationinourgeneraalerodynamiscoslveTr LNSaD.
each stage k, yielding the expression
TheESDIRKschemeasreusedratherthanthestan-
dardSDIRKscheme(sail= 0,seeHairer7)be- c_Ut: Ut: - U_ t:
causetheESDIRKschemecsanachievestageor- 0T + St +Zat:jsJ=0, k=l, (5)
der2(eachstageisatleastsecond-ordaecrcurate), j=l
andtheextracolumnofnonzeroButchecroefficients
In this form, equation (5) is amenable to all the non-
(a11)allowsmoreflexibilitywhendesigningnew
schemesF.orallschemetsh,e"stittty-accuratea"s- linear solving machinery available in TLNS3D. Each
nonlinear equation is marched in pseudo-time with
sumption(asj=
bj) is enforced, which automatically
multi-grid acceleration until a predetermined con-
extends A-stability into L-stability. In addition, it
vergence criterion is satisfied.
eliminates the potentially damaging explicit update
Equation (5) can be rewritten in the form
U n+l = U n +At _=1 bJSj from the algorithm, and
replaces it with the condition U"+1 = Us. OUt: Ut:
Numerous ESDIRK schemes were developed and o_ + X/+ at:t:st:+ E = 0, k = 1,, (6)
tested in this study. Three schemes, one each of
third-, fourth-, and fifth-order accuracy were chosen with
as being representative of each of these orders. No
k-1
claim is made as to their optimality. Henceforth, _U n
they are referred to as ESDIRK3, ESDIRK4, and E -- At q- E at:JSJ'
j=l
ESDIRK5, respectively. All three schemes are pre-
sented in Kennedy 9. The coefficients for the recom- where E includes all the iteration independent in-
mended ESDIRK4 scheme are included in the Ap- formation at stage k of the RK scheme. Discretiz-
pendix. ing equation (6) in the variable r using an Implicit-
Explicit Runge-Kutta (RK_) scheme, yields
Solution Algorithm
Discretizing equation (1) with an s-stage ESDIRK U kp _ U k° U kp (r)
+ %(A7- + at:t:St:('-_) + E) = 0
scheme represented in equation (4) yields for stage AT
k, the expression
where the p is the stage value of the /_K_ scheme.
Note that the contribution from the time term is
U k _ U n k
q-EakJSJ = 0, k = 1,s (4) treated implicitly, the inviscid and viscous flux terms
At
are treated explicitly, and that first-order temporal
j=l
accuracy is sufficient for this scheme. Adding, sub-
with tracting, rearranging terms and accounting for resid-
ual smoothing in equation (7) yields
sj _ c_(-F + Fv) + cq(-G q- Gv)
0_ 0rj c_pAT t:p c_pAT t:p-_
(1+ _7-)v = v t:°+ _7-v
c_(-H + Hv)
+
oC Ut:p-1
- apArLi-_ls[ A7 + at:t:St:'-_ + E], (8)
The vector Uk is the solution at stage k, U n is the
solution at the previous time level n, and akj are where Li_ls is the implicit residual smoothing matrix.
the Butcher parameters for the RK-method used. Note that the bracketed term in equation (8) is still
Again note that in our case the last stage gives the the residual of the physical time RK scheme.
solution at the new time level, that is U n+l = U s.
Pseudo-Time Stability Analysis
Advancing the solution vector Uk from time level
Each grid level of the multi-grid process yields an
n --+ n + 1, requires solving the sequential set of
equation similar to equation (8), with minor differ-
s nonlinear algebraic equations defined in equation
ences in the variable E accounting for grid trans-
(4).
fer contributions. The pseudo-time sub-iteration
Pseudo-Time Iterative Algorithm should converge rapidly, and any remaining residual
We follow here the work of Melson et al. 12 at the end of the iteration should be devoid of high-
which was originally developed for a BDF algorithm. frequency content. A Fourier stability analysis is
Equation (4) can be difficult to solve in its present used to analyze the exact nature of the pseudo-time
sub-iterationandtooptimizetheresiduaslmooth-
in high Reynolds number calculations.
ingparameters. An exhaustive study in the parameters 7, /_,
A necessarfyirst testfor thesub-iterational-
yields the following general conclusions. The resid-
gorithmis goodbehaviofrorthe1-Dwaveequa-
ual smoothing parameter 7 is most productive when
tion
Ut + aUx = 0. Convergence and stability are the effective stability limit of the sub-iteration is in-
strong functions of the discrete spatial operator Ux. creased by a factor 2 - 3. Parameter values in the
TLNS3D uses spatial discretizations that closely re- range 1.0 _< 7 _< 1.5 produce this increase. Figures
semble an approximation of the form (1)-(3) show a comparison of the damping charac-
teristics as a function of _, with the value 7 = 1.5.
OU Ui+l - Ui__
Shown are the cases /_ = 0, /_ = _, and /_ = 1,
Ox 2Ax
with 1 _< _ < 7 varying on each plot. Plots of
Ui+2 - 4Ui+_ + 6Ui - 4Ui__ + Ui-2
+ /_ < _ are virtually indistinguishable fl'om the case
32Ax
/_ = 0. The three plots reveal that the pseudo-time
The analysis begins by substituting this spatial op- sub-iteration is stable for _ _ 7 ,and has the best
erator into the 1-D wave equation on a periodic do- damping characteristics in the 7_r _< 0 _<rrrange for
main, discretizing time using the ESDIRK scheme, values of _ near the stability limit.
and assuming solutions of the form U(x,t) = The parametric study yields the following algo-
ii(t)ei 0o). Implementing the pseudo-time sub-
iteration algorithm described in equation (8) then
yields for the pth stage of the /_K_ operator, the
expression
(l+%R)g p = gO+%Rgp-1
p-1 + akk ¢gP- ], (9)
with
¢ = i sin(0) + _ sin(0/2) 4
R =0, CFL= 1.0
Li-_l_ = (1 + 47 sin(0/2)2) -1 , R =0, CFL= 2.0
R =0, CFL= 3.0
R =0, CFL= 4.0
and R =0, CFL= 5.0
R =0, CFL= 6.0
aAr R =0, CFL= 7.0
gP = i__;p R = AAtt. ' _ = akk--.;Ax i i i i II i i i i 2I i i i i 3I
Fourier Mode
The variable E is eliminated fl'om consideration in
the stability analysis because it is constant during Figure 1. Damping characteristics for the case R =
the pseudo-time sub-iteration and therefore has no a_ O, with "CFL" = _.
At --
influence on the stability. The diagonal coefficient
akk is the same on each stage of an ESDIRK, which rithm for the pseudo-time sub-iteration. The /_K_
yields an identical stability analysis for each stage. scheme is implemented in five stages with the co-
(A stability analysis for the BDF scheme given by efficients being defined by c; = [1/4,1/6,3/8,1/2,1].
equation (2) yields a similar result, with the variable Three evaluations of artificial dissipation terms
akk replaced by ilk.) (computed at the odd stages) are used to obtain a
Several observations about equation (9) are in- larger stability bound, which allows a higher CFL
structive. The independent parameters are 7: the number in the presence of physical diffusion terms.
level of implicit residual smoothing, R: the ratio of The stability limit of the numerical method is fur-
time-steps, and _: the CFL condition for the pseudo- ther increased with the use of the implicit residual
time-stepping scheme multiplied by the diagonal co- smoothing technique that employs grid aspect-ratio-
efficient from the ESDIRK scheme. In general the dependent coefficients is and local time-stepping is
parameter R satisfies the condition 0 < R < 1. The used in each cell. The efficiency of the solution pro-
limit R --+ 0 or At --+ cc reduces to the steady cess is significantly enhanced through the use of a
state formulation. The limit R --+ 1 corresponds multigrid acceleration technique. For a detailed de-
to a physical time-step for which an explicit method scription of the resulting iterative method see Vatsa
would be stable, a situation that is unlikely to occur et al. is The same concept for the computation of
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 .c£ 0.6
8
0.5 0.5
_= , = .
0.4 0.4
R =1/10, CFL =1.0
0.3
R =1/10, CFL =2.0
R =1/10, CFL =3.0 R =1,CFL= 3.0
0.2 R =1/10, CFL =4.0 0.2 _ R =1,CFL= 4.0
R =1/10, CFL =5.0 R =1,CFL= 5.0
0.1 R =1/10, CFL =6.0 0.1 + R =1,CFL= 6.0
R =1/10, CFL =7.0 0.3! + RR ==11,,CCFFLL== 72..00
% i i i i I i i i i I i i i i I 0 i i i i I i i i i I i i i i I
I 2 3 I 2 3
Fourier Mode Fourier Mode
Figure 2. Damping characteristics for the case R = Figure 3. Damping characteristics for the case R =
AAtr = --]10,with "CFL" = i. A_ _ l, with "CEL" = A.
At --
unsteady flows is used by other researchers, see for fifth-order accuracy. The ESDIRK3, ESDIRK4, and
example Martinelli 10 and Swanson and Turkel 16. ESDIRK5, used in this study (as well as the work of
Kennedy 9) are representative candidates from this
Diagonal Coetficients 1
study of ESDIRK schemes. The value akk =
The pseudo-time sub-iteration is strongly influ-
in the ESDIRK4 scheme is an example of five im-
enced by the diagonal coefficients akk and ilk. In
plicit stages producing a more efficient fourth-order
TLNS3D the pseudo-time sub-iteration is always ad- scheme than do four. This is consistent with the
vanced with the maximum allowable scaled pseudo-
findings of other investigators (see Hairer 7).
time-step t = akkt _ 7. The rate of relaxation in
non - scaled pseudo-time is therefore, inversely pro- Method Stgs Impl Stgs Order akk
portional to the diagonal coefficient akk: the smaller ESDIRK3 4 3 3 0.435
the value of akk the more rapidly the pseudo-time ESDIRK4 6 5 4 0.250
sub-iteration progresses. The values akk and flk vary ESDIRK5 7 6 5 0.184
significantly for different physical time-advancement
Table 2. Properties of the ESDIRI¢ methods.
schemes. Table (1) presents the diagonal coefficients
Validation of the Space-Time Method
flk for the BDF schemes, where the coefficients vary
by approximately a factor of 2. Table (2) presents The accuracy of the space-time integration meth-
the diagonal contribution of the ESDIRK schemes, ods is investigated for an unsteady laminar flow test
as well as the number of stages, the implicit stages case. The test problem is laminar flow around a two-
and the order. The akk coefficients vary by ap- dimensional circular cylinder at a Reynolds number
proximately a factor of two, and are generally much of 1200 and a Mach number of 0.3. The initial flow
smaller than the BDF flk values. is symmetric with zero lift. As the wake behind the
Unlike the BDF schemes, the ESDIRK schemes cylinder starts to grow, it becomes unstable and be-
can be optimized to improve their efficiency. The gins to shed from alternate sides of the cylinder. De-
important parameters are the diagonal coefficient tailed numerical and experimental investigations of
akk and the number of stages s. A fortuitous trend this flow have been performed by several authors
observed with the ESDIRK schemes is a general de- a, a, 14,16 The computational grid of 97 x 65 is
crease in the value of akk with increasing stage num- shown in Figure 4. The boundary is a distance of
ber, at a given accuracy and L-stability property. 20 times the diameter of the cylinder away from the
Increasing the number of implicit stages does not wall, while the distance between the wall and the
always decrease the efficiency of the scheme, and first grid point is 0.001 times the diameter of the
sometimes yields greater" efficiency. Fifteen ESDIRK cylinder. Grid points are clustered in the wake. A
schemes ranging from 3to 8 stages were investigated density contour plot is shown in Figure 5 as calcu-
to determine good choices of third-, fourth-, and lated on the 97 x 65 grid. Note that the near wall
vorticasltructureasppeatrobesufficientlryesolved,
butthatresolutioinslostasthegridexpandisnthe
far-field.Inthesepreliminarycalculationa, small
time-stepwaschosens,othatthedominanctompo-
nentoferroristhespatiacl ontribution.
_J
\\\\
Figure 5. The density contours as calculated on the
97 x 65 O-Grid.
Temporal Accuracy
A temporal refinement study is performed to as-
sess the accuracy of various ESDIRK and BDF
Figure 4. The 97 x 65 O-Grid used in the cireular schemes. The initial condition for the study was
cylinder study. obtained by simulating the limit cycle behavior of
the flow for approximately 20 shedding cycles, with
a relatively small time-step (At = ½). After 20
Figure 6shows a comparison of the lift history for cycles, the solution was stored in a restart file for
one shedding cycle as calculated on the 97 × 65 and use as the initial condition in the subsequent stud-
193 × 129 grids. The 193 × 129 data is shifted so ies. A classical temporal study was then performed
that the first zero in lift coincides on both curves. from this initial condition. A typical practitioner is
The Strouhal number in each case was 0.2489, and interested in lift, drag, pitching moment, skin fric-
0.2459 as calculated on the 97 × 65 and 193 × 129 tion, and the frequency spectrum over several cy-
grids, respectively. Small differences in the lift are cles. Thus, the time interval of the study includes
seen near the peaks of the cycle. These values are approximately 1_ shedding cycles. This interval is
larger than 0.21, reported from experiments 4,14 sufficient in length to allow accumulation of tempo-
As Mittal and Balachander report 13 this might be ral error during the shedding cycle. No exact solu-
caused by the onset of three- dimensional effects, tion is know for this problem, so a "numerical exact"
not captured in our two- dimensional computations. solution was obtained using a small time-step and
Two-dimensional computations performed by other a small iteration tolerance on the non-linear solves.
researchers 5,17resulted in larger Strouhal numbers The "exact" time-step used was At = 0.05, and a
too, in the range of 0.23-0.24. This study and a more stopping criterion for the iterative solve of the im-
exhaustive study of inviscid vortex propagation sup- plicit equations, was max(residual) < 10-6. The
port the conclusion that the space-time operator is "exact" solution is accurate to approximately 6 sig-
working properly, and that the spatial operator con- nificant digits in the lift.
verges at the design spatial accuracy. The coarse The lift on the body is used as the representative
grid (97 × 65) provides adequate spatial resolution measure of error in all calculations. Other integral
to capture the relevant large scales features in the measures including drag, skin friction, total drag,
shedding process. As such, the 97 × 65 grid is cho- pitching moment, as well as L2 and Loo norms over
sen as the basis of most of the temporal refinement the domain were studied. All integral measures yield
studies presented in this work. nearly the same quantitative conclusions, although
1.5 _ 193x129
97x65 -1.5
1 -2
-2.5
" -3
_'_-3.5
,0-I -4
-4.5
f
-5 _ Lift
DragT
Pitch
-5.5 _ _ DragV
45 50 55 -_6.751 I I I_0_.5 .... -0'25. .... 0_ .... 0.25' I
Time Log_0 (dt)
Figure 6. Comparison of lift over one shedding cycle Figure 7. Convergence behavior for lift, drag and
between coarse and fine grids pitching moment as calculated with a fourth-order
ES'DIRK scheme.
the qualitative errors are different for each case. A
nonlinear equations at each stage (step) are solved to
detailed investigation of the location of the maxi-
tolerances which are small compared with the abso-
mum temporal error revealed that the vortex gen-
lute error in the calculation. There is a dramatic in-
eration process in the near-wall regions is the most
temporal demanding portion of the cycle. It is not crease in accuracy in going from BDF1 to ESDIRK5.
For example, an error tolerance of 10-1 is achieve
surprising, therefore, that the lift integral is a good
measure of total error in the calculation. with a time-step of At = 1 for the ESDIRK4 and
Detailed results from the study are now presented. ESDIRK5 schemes, while the BDF1 and BDF2 re-
quire At = 10.2 and At = 10-1 , respectively.
Figure 7 shows a detailed refinement study with an
ESDIRK4 scheme. Shown are the solution errors in The BDF3 and the ESDIRK3 schemes have nearly
the lift, viscous drag, total drag, and pitching mo- the same absolute level of error, although the con-
ment, as a function of logarithm of time-step. In vergence behavior of the BDF3 scheme is sporadic.
No explanation of this behavior was identified, al-
all cases, the nonlinear system is solved to strict
tolerances to eliminate "iteration" error as a con- though a possible explanation is the schemes lack of
A-stability.
taminating variable in the study. At coarse time-
As mentioned previously in this work, the BDF2
steps the solution accuracy deteriorates away from
scheme is consistently used by practitioners because
design accuracy. For sufficiently small time-steps,
of its robustness, simplicity and efficiency. The re-
design accuracy is clearly demonstrated in all vari-
sults in Figure 8 clearly show that the ESDIRK4
ables. For example, a least-squares fit of the finest
four data points on each curve reveals that the con- scheme can be used at time-steps which are a fac-
tor of ten larger than those used in the BDF2, while
vergence rates for [lift, drag V, drag T, pitch] are
achieving similar accuracy. At fine tolerances, this
[3.9628, 4.0257, 4.0478, 4.0251]. The coarsest data
point corresponds to a time-step of At = 2, for which difference becomes even greater. Figure 8 still can
not be used to conclude that the ESDIRK4 scheme
twelve points resolve the shedding cycle. Note that
15 - 20 points are needed to resolve the cycle be- is more efficient than the BDF2 scheme. To do so
fore the fourth-order scheme converges with design requires a detailed accounting of the work involved
accuracy. This resolution is consistent with conven- in each algorithm, and is not a simple task.
tional estimates of "points-per-wavelength" needed
to resolve a periodic wave. Temporal efficiency
Figure 8 presents the error in lift versus the time For large computations, the work involved in ad-
step (log-log) for three BDF schemes and three ES- vancing the solution from t = To to t = Tf is
DIRK schemes. The ESDIRK schemes presented in proportional to the number of non-linear solves re-
this figure are summarized in table (2). Again, the quired over that interval, but also depends strongly
CIRCULAR CYLINDER CIRCULAR CYLINDER
-0.5
_1.5
£ _c .J
/
[_ -2 .c -2
=E_ /
,,_1
LU
#-3 _-2.5
o
,,_1 BDF1 ,o,_1
BDF2 -3 BDF1
BDF3 BDF2
ESDIRK3 BDF3
ESDIRK4 -3.5 ,_ ESDIRK3
ESDIRK5 ESDIRK4
ESDIRK5
i i i i iiii I I I I I IIII I I I I I I I IIIII I i i i iiiii
101 10U 104 10_
Timestep Work
Figure 8. Convergence as afunction of time-step for Figure 9. Convergence as a function of work for
BDF and ESDIRI< schemes. BDF and ESDIRK schemes.
on how quickly each non-linear solve converges• The order scheme at extremely coarse tolerances• For
number of non-hne•ar solves .is si.mply NL = Tj-TAt 0Is, example, the BDF2 is 2.5 times less efficient than
where Is is the number of implicit stages per time- the ESDIRK4 scheme at an error tolerance of 10%
step. (An explicit stage is virtually "free": hence the (one significant digit in the solution accuracy)• As
use of the ESDIRK instead of the SDIRK schemes.) the error tolerance becomes more strict, the high-
For the BDF schemes Is = 1. A much more dif- order schemes easily outperform the BDF1 and the
ficult aspect to predict, however, is how rapidly a BDF2 scheme. For a desired temporal error of 1% in
nonlinear solve will converge• Smaller physical time- lift the fourth-order integration method ESDIRK4,
steps provide a much better initial guess for the non- only requires 1.5% of the work required by BDF1,
linear iteration, which implies an advantage for the an efficiency increase with a factor 70. For a tem-
BDF schemes. The ESDIRK schemes, however, have poral error of 0.1% solution ESDIRK4 requires 10%
a much smaller diagonal coefficient akk which in- of the work of BDF2, an efficiency increase with a
creases the asymptotic convergence rate of the multi- factor 10. Note again that the BDF3 scheme shows
grid process• an irregular behavior• Perhaps this is due to the fact
Figure 9 presents the convergence of the six that the BDF3 solution oscillates around the exact
schemes as a function of the required work. Three solution for different time-steps, and sometimes by
accuracy levels are chosen, [10-1,10-2,10 .3] as coincidence yields unusually low levels of error•
representative of desired engineering accuracy lev- It is interesting to note that the fourth- and
els. The appropriate values of At needed for each fifth-order ESDIRK schemes have exactly the same
method are obtained from Figure 8. An iteration accuracy-work ratio, and that their convergence be-
tolerance a factor of 200 times smaller than the de- havior appears to be logarithmic in nature• Both
sired error level, is used for the nonlinear iteration, have different time-steps, stages, and diagonal coef-
with the error based on the Loo norm of the resid- ficients akk, and there is no reason why they should
ual. For example, if the desired error is 10.2 then at lie on the same line, or have logarithmic convergence
each stage (step) the nonlinear system is solved until behavior• Though both have the same efficiency, the
the maximum residual is 1°-_ The work from each ESDIRK5 scheme requires more storage and is less
W6_-"
method is measured as the total number of multi- robust than the ESDIRK4 scheme (internal stabil-
grid cycles used in the entire time interval• ity problems at huge time-steps). Therefore, for an
efficient and robust solution in time the ESDIRK4
An obvious conclusion from the study presented
scheme is recommended.
in Figure 9is that the BDF1 scheme (Euler Implicit)
will never compete with the higher-order schemes in A major contributor to the inefficiency of implicit
terms of efficiency. Similarly, the BDF2 scheme is methods is solving the nonlinear systems at each
only competitive with the third-, fourth- and fifth- stage (step) to inappropriate sub-iteration tolerance