jst_2015_updated_07_03_2015 | Computational Fluid Dynamics

Please download to get full document.

View again

of 41
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Information Report
Category:

Documents

Published:

Views: 4 | Pages: 41

Extension: PDF | Download: 0

Share
Related documents
Description
jst scheme
Transcript
  The Origins and Further Development of theJameson-Schmidt-Turkel (JST) Scheme Antony Jameson ∗  Department of Aeronautics and Astronautics, Stanford University, Stanford, CA, 94305 This paper is in response to an invitation to give a presentation on the srcins and subsequent developmentof the Jameson-Schmidt-Turkel scheme. After a description of the historical background and initial developmentof the scheme, the paper discusses three main developments: first, the development of convergence accelerationmethods, including residual averaging and multigrid; second, the extension to unstructured grids; and third thereformulation of the JST scheme as a total variation diminishing (TVD) scheme, and its relationship to symmetricTVD schemes. The paper concludes with a brief review of applications to unsteady and viscous flows. I. Introduction This paper is in response to an invitation to give a presentation on the srcins and subsequent development of theJameson-Schmidt-Turkel (JST) scheme in a session on history at the Aviation 2015 meeting in Dallas, Texas. The JSTscheme was presented at the AIAA 14 th Fluid and Plasma Dynamics Conference in Palo Alto, CA, June 1981 as AIAAPaper 81-1259. This paper has been widely cited, but it was never submitted for publication as a journal article, primarilybecause the present author considered it as only the first step in as yet incomplete research to find a definitive method tosolve the Euler equations.Section II. of the paper gives a synopsis of the srcins and goals of the JST scheme, which immediately proved to besuccessful in a wide variety of CFD simulations. The next steps were the further improvement of the rate of convergenceto a steady state by residual averaging, and the introduction of a full approximation agglomeration multigrid scheme,described in section III., and the extension of the method to treat arbitrarily complex configurations on unstructuredtetrahedral meshes, outlined in section IV..Soon after the srcinal presentation of the JST scheme in 1981, there appeared the seminal papers of Roe 2 , in whichhe proposed a characteristics based flux derived from an elegant linearization, and Harten 1 , in which he introduced theconcept of total variation diminishing (TVD) schemes for scalar conservation laws. Also, by this time the series of papersby Van Leer 3–7 had brought the pioneering work of Godunov 8 to the attention of the CFD community, and introducedthe concept of reconstruction with limiters † for the formulation of higher order accurate non-oscillatory shock capturingschemes. Thesepapersusheredintheeraofwhattheauthorlikestothinkofasthe“fluxwars”(Godunov,Steger-Warming,Roe, Osher, AUSM, HLL, HLLC, CUSP) and the “reconstruction wars” (MUSCL, ENO, WENO, SLIP, minmod limiter,superbee limiter, Van Leer limiter, Van Albada limiter, Venkatakrishnan limiter). Ultimately, one can mix and match fluxformulations and reconstruction procedures. In order to ensure that the discrete scheme is TVD for a scalar conservationlaw, the updated state should be a convex combination of the states at neighboring cells at the previous time step. It turnsout that this concept extends to the gas dynamics equations. If the discrete solution is physically realizable in the sensethat the pressure and density in each cell are non-negative, and the updated state in each cell is a convex combination of the states in neighboring cells, the updated state will also be physically realizable. Moreover, this property can be provenfor some flux vector splitting methods, and also for methods stabilized by scalar diffusion, similar to that used in the JST ∗ Professor of Engineering, Department of Aeronautics and Astronautics, Stanford University, AIAA Fellow † The idea of limiters to prevent the formation of new extrema in discreet solutions of the gas dynamics equations was first introduced by Boris andBook in their flux corrected transport (FCT) scheme 9 .1 of  41American Institute of Aeronautics and Astronautics  scheme. In section V., it is shown that when the JST scheme is applied to a scalar conservation law with an appropriate formulation of a switching function, it is actually a symmetric TVD scheme, which may also be labeled as a symmetriclimited positive (SLIP) scheme. It is also shown that the SLIP and JST schemes can be reformulated as second orderaccurate reconstruction schemes, which can then be combined with any preferred choice of the numerical flux.While the JST scheme was srcinallydesignedfor steadystate solutions, it works quitewell for unsteadyflows whenitis combined with an appropriate time stepping scheme. This is discussed in section VI., while section VII. presents some examples of the application of the JST scheme to compressible viscous flows, modeled by the Reynolds averaged Navier-Stokes (RANS) equations. Finally section VIII. presents some conclusions, together with some future requirements forCFD. II. Origins and goals of the JST scheme Thepresentauthorhad actuallyembarkedona longterm programsince 1970to developcomputationalfluid dynamics(CFD) forthe purposeof airplanedesign. Starting with two dimensionaltransonicflow calculationsfor airfoils, the under-lying strategy of this program was to make incremental advances in both the geometric complexity and the mathematicalmodels of the flow, with the ultimate aim of calculating transonic and supersonic flows over complete aircraft, using atleast the Euler equations of inviscid compressible flow as the mathematical model. The level of complexity that could beattempted at any stage of this program was paced by the available computing power, both processing speed and memory.In the early seventies, this was sufficient to calculate transonic potential flow solutions for airfoils, using, for example, FLO 6  10 , but in the author’s opinion not sufficient to attempt Euler solutions for any geometry of significant complexity.By 1975, it was feasible to calculate transonic potential flow solutions for swept wings using  FLO 22 11 , but the ControlData 6600, which represented the state of art in computing machinery, had only 131,000 words of core memory, with theconsequence that the solution had to be stored on the disk, and computed plane by plane.  FLO 22  used an analytic trans-formation of the equations to a curvilinear coordinate system to represent the geometry. This approach could hardly beextended to any configuration more complex than an isolated wing. We overcame this difficulty in  FLO 27   by developinga discretization scheme for arbitrary body fitted hexahedral meshes 12 . Although we called this a finite volume scheme, itactually introduced isoparametric trilinear finite elements, stabilized by the addition of artificial diffusive terms to give anupwind bias. This method could have been used to calculate the flow over a complete aircraft if it was provided with asuitable mesh, but at that time the required mesh generation technology had yet to be developed. Ultimately, the Frenchteam from INRIA and Dassault 13 managed to calculate a transonic potential flow solution for a Falcon jet in 1982 usinga finite element method on a tetrahedral mesh.In the meanwhile, the present author had begun thinking seriously about ways to solve the Euler equations for steadyflows, and first wrote an Euler solver,  EUL 1 ‡ in 1976, which integrated the time dependent equations to reach a steadystate. This code used the Z scheme, defined as follows:To solve the linear advection equation ∂u∂t  + a∂u∂x  = 0  (1)with  a >  0  (a right running waver), the discrete scheme on a uniform mesh with time and space intervals  ∆ t  and  ∆ x  is v n +1 j  =  v nj  −  λ 2  v n +1 j  − v n +1 j − 1  −  λ 2  v nj +1 − v nj  , where  v nj  is the discrete solution at time level  n  and meshpoint  j , and  λ  is the CFL number λ  =  a ∆ t ∆ x. The scheme is semi-implicit with the illustrated stencil, similar to a single step version of the MacCormack scheme 14 . Itis second order accurate and stable for an arbitrarily large positive CFL number λ . For a nonlinear problem such as ∂u∂t  +  ∂ ∂xf  ( u ) = 0 , ‡ The author still has a working copy of this code.2 of  41American Institute of Aeronautics and Astronautics  Figure 1: Z-schemethe scheme takes the form v n +1 j  =  v nj  −  ∆ t 2∆ x  f  ( v n +1 j  ) − f  ( v n +1 j − 1  )  −  ∆ t 2∆ x  f  ( v nj +1 ) − f  ( v nj  )  , which requires iterative solution for  v n +1 j  . At locally supersonic points the difference scheme was switched to the secondorder accurate backward difference formulaat both time levels n  and n +1 . This switch led to cleanly captured numericalshocks in transonic flow calculations. The resulting scheme reliably converged to a steady state and gave quite goodresults for simple geometries such as an axisymmetric boat tail. These studies of the Z-scheme were never published,and the scheme wound up being replaced by the LU implicit scheme 15 which could be regarded as forward and reversesweeps of the Z-scheme to give unconditional stability for waves traveling in both directions.In 1979 there was held a GAMM workshop in Stockholm on  “Numerical methods for the Computation of Inviscid Transonic Flow with Shock Waves” , organized by Rizzi and Viviand. The results of a variety of methods to calculateboth transonic potential flow and Euler solutions were presented and compared. The Euler methods generally relied ontime marching, and it appeared that none of them could reach a steady state. This was the general state of affairs whenthe author was approached by Wolfgang Schmidt § , then the leader of the CFD group at Dornier, with a proposal to visitDornier for three weeks during the summer of 1980, and embark on a collaboration to produce an improved Euler solver.Our primary goal from the outset was the calculation of fully convergedsteady state solutions, althoughsince we wereusing a time marching method, we also envisaged dual use for both steady and unsteady flows. The starting point wasan existing in-house code to solve the two dimensional Euler equations which had been written by Rizzi and Schmidt. ¶ This used a finite volume implementation of the MacCormack scheme 16 . Because the MacCormack scheme producesoscillations at shock waves, it was modified by the addition of artificial diffusive terms with a coefficient proportional toa sensor for sharp pressure gradients. In a one dimensional case, the form of the sensor is ǫ  =   p j +1 − 2  p j  +  p j − 1  p j +1  + 2  p j  +  p j − 1  ≤ 1 , where it is assumed that the pressure is positive. The code produced apparently reasonable answers but it would not fullyconverge to steady state. Moreover, the MacCormack scheme, like the Lax Wendroff scheme, bundles the space and timediscretizations, with the consequence that a steady state solution, if it could be reached, would depend on the time step ∆ t . The first step in developing the new code was to replace the MacCormack scheme by a semi-discrete scheme whichseparated the space and time discretizations (also called the method of lines). At the same time, we concluded that wecould use a simple central difference scheme for the spatial discretization and rely on added diffusive terms to preventoscillations in the vicinity of shock waves. The time stepping scheme we first used was based on the idea of iterating theresiduals of an implicit scheme. In order to solve dwdt  + R ( w ) = 0 , where  R ( w )  is the space residual for the solution  w , a second order accurate implicit scheme has the form w n +1 =  w n −  ∆ t 2  R ( w n +1 ) + R ( w n )  . § Sadly, Wolfgang Schmidt died of a sudden massive heart attack in 2008. ¶ This code was extremely cleanly written and caused the author to completely revise his programming style to emphasize clarity.3 of  41American Institute of Aeronautics and Astronautics  At any time step, denoting w (0) as the starting value, a simple iterative procedure is to set w (1) =  w (0) − ∆ tR ( w (0) ) w (2) =  w (0) −  ∆ t 2  R ( w (1) ) + R ( w (0) )  w (3) =  w (0) −  ∆ t 2  R ( w (2) ) + R ( w (0) )  ...  = ...It turns out from a linear analysis of the advection equation, that this iteration verges for CFL numbers ≤ 2, but there isno benefit in using more than 3 iterations, so we used a 3 stage scheme where  w (0) =  w n ,  w n +1 =  w (3) .The second significant change we implemented was in the construction of the artificial diffusive terms. Because thefinite volume equations are usually expressed as V  dwdt  + R ( w ) = 0 , where  V   is the cell area or volume, and  R ( w )  is the flux balance for the state vector, the existing Dornier code addeddiffusive terms which had the form D i,j  =  d i + 12 ,j − d i − 12 ,j  + d i,j + 12 − d i,j − 12 , where the indices  i  +  12  and  j  +  12  denote the cell interfaces, while, for example, the diffusive flux for the density  ρ  hadthe form d i + 12 ,j  =  ǫ (( ρV  ) i +1 ,j  − ( ρV  ) i,j ) . In the far field where  ρ  should eventually be constant, these diffusive terms reduce to differences of the cell areas orvolumes, with the consequencethat uniformflow would no longer satisfy the differenceequations. In order to rectify this,it was essential to base the diffusive terms on differences of the state variables without including cell metrics. It remainednecessary to determine an appropriate scaling. Since we wanted a steady state that would not depend on the time step  ∆ t ,we concluded that the scaling should depend on  ∆ t/ ∆ t ∗ , where  ∆ t ∗ is the time step for a CFL number of unity, whichmay be estimated as follows. On a face with vector area  S S S   with magnitude | S  | , denote the local velocity and speed of sound by q q q   and  c . Then define qs  =  q q q  · S S S, cs  =  c | S  | . Now in the case of a two dimensional flow ∆ t ∗ =  | S  | ( qs + cs ) i  + ( qs + cs ) j , where  ( qs + cs ) i  and  ( qs + cs ) j  are average values of   qs + cs  in each coordinate direction. The quantities  ( qs + cs ) i  and ( qs  + cs ) j  represent the spectral radii of the Jacobian matrices of the flux vectors in the  i  and  j  directions, so this scalingleads to diffusive terms scaled to the sum of the spectral radii, whereas it is now well known that it is sufficient to adddiffusive terms scaled to the separate spectral radii in each coordinate direction. The use of the sum of the spectral radii,however, results in a scheme that is quite robust for grids containing high aspect ratio cells.In order to accelerate convergence to a steady state, we used a variable local time step corresponding to a fixed CFLnumber, such that waves should propagate from the body to the far field in a fixed number of time steps proportional tothe number of mesh intervals between the body and the far field. The code also provided for the inclusion of enthalpydamping to accelerate convergence to a steady state. An analysis of these ideas 15 suggests that it might be possible toobtainexponentialconvergenceintime,whileabandoninganyaspirationtotimeaccuracyforanunsteadyflowcalculation.A property of steady state Euler solutions is that the total enthalpy should be constant throughout the flow field. In orderto satisfy this property, the diffusive terms for the energy equation were calculated from differences of the total enthalpy 4 of  41American Institute of Aeronautics and Astronautics
Recommended
View more...
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks