ebook img

NASA Technical Reports Server (NTRS) 20010066914: Time Integration Schemes for the Unsteady Navier-stokes Equations PDF

13 Pages·0.87 MB·English
Save to my drive
Quick download
Download
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview NASA Technical Reports Server (NTRS) 20010066914: Time Integration Schemes for the Unsteady Navier-stokes Equations

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

See more

The list of books you might like

Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.