The Origins and Further Development of theJamesonSchmidtTurkel (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 JamesonSchmidtTurkel scheme. After a description of the historical background and initial developmentof the scheme, the paper discusses three main developments: ﬁrst, 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 ﬂows.
I. Introduction
This paper is in response to an invitation to give a presentation on the srcins and subsequent development of theJamesonSchmidtTurkel (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 811259. 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 ﬁrst step in as yet incomplete research to ﬁnd a deﬁnitive 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 conﬁgurations 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 ﬂux 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 nonoscillatory shock capturingschemes. Thesepapersusheredintheeraofwhattheauthorlikestothinkofasthe“ﬂuxwars”(Godunov,StegerWarming,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 ﬂuxformulations 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 nonnegative, 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 ﬂux 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 ﬁrst introduced by Boris andBook in their ﬂux 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 ﬂux.While the JST scheme was srcinallydesignedfor steadystate solutions, it works quitewell for unsteadyﬂows 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 ﬂows, modeled by the Reynolds averaged NavierStokes (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 developcomputationalﬂuid dynamics(CFD) forthe purposeof airplanedesign. Starting with two dimensionaltransonicﬂow calculationsfor airfoils, the underlying strategy of this program was to make incremental advances in both the geometric complexity and the mathematicalmodels of the ﬂow, with the ultimate aim of calculating transonic and supersonic ﬂows over complete aircraft, using atleast the Euler equations of inviscid compressible ﬂow 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 sufﬁcient to calculate transonic potential ﬂow solutions for airfoils, using, for example,
FLO 6
10
, but in the author’s opinion not sufﬁcient to attempt Euler solutions for any geometry of signiﬁcant complexity.By 1975, it was feasible to calculate transonic potential ﬂow 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 transformation of the equations to a curvilinear coordinate system to represent the geometry. This approach could hardly beextended to any conﬁguration more complex than an isolated wing. We overcame this difﬁculty in
FLO 27
by developinga discretization scheme for arbitrary body ﬁtted hexahedral meshes
12
. Although we called this a ﬁnite volume scheme, itactually introduced isoparametric trilinear ﬁnite elements, stabilized by the addition of artiﬁcial diffusive terms to give anupwind bias. This method could have been used to calculate the ﬂow 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 ﬂow solution for a Falcon jet in 1982 usinga ﬁnite element method on a tetrahedral mesh.In the meanwhile, the present author had begun thinking seriously about ways to solve the Euler equations for steadyﬂows, and ﬁrst wrote an Euler solver,
EUL 1
‡
in 1976, which integrated the time dependent equations to reach a steadystate. This code used the Z scheme, deﬁned 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 semiimplicit 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: Zschemethe 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 ﬂow 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 Zscheme 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 Zscheme 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 ﬂow 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 ﬂows. The starting point wasan existing inhouse code to solve the two dimensional Euler equations which had been written by Rizzi and Schmidt.
¶
This used a ﬁnite volume implementation of the MacCormack scheme
16
. Because the MacCormack scheme producesoscillations at shock waves, it was modiﬁed by the addition of artiﬁcial diffusive terms with a coefﬁcient 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 ﬁrst step in developing the new code was to replace the MacCormack scheme by a semidiscrete 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 ﬁrst 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 beneﬁt 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 signiﬁcant change we implemented was in the construction of the artiﬁcial diffusive terms. Because theﬁnite volume equations are usually expressed as
V dwdt
+
R
(
w
) = 0
,
where
V
is the cell area or volume, and
R
(
w
)
is the ﬂux 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 ﬂux for the density
ρ
hadthe form
d
i
+
12
,j
=
ǫ
((
ρV
)
i
+1
,j
−
(
ρV
)
i,j
)
.
In the far ﬁeld where
ρ
should eventually be constant, these diffusive terms reduce to differences of the cell areas orvolumes, with the consequencethat uniformﬂow 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 deﬁne
qs
=
q q q
·
S S S, cs
=
c

S

.
Now in the case of a two dimensional ﬂow
∆
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 ﬂux 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 sufﬁcient 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 ﬁxed CFLnumber, such that waves should propagate from the body to the far ﬁeld in a ﬁxed number of time steps proportional tothe number of mesh intervals between the body and the far ﬁeld. 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,whileabandoninganyaspirationtotimeaccuracyforanunsteadyﬂowcalculation.A property of steady state Euler solutions is that the total enthalpy should be constant throughout the ﬂow ﬁeld. 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