![FIG. 9.2. The numerical solutions of ρ (“+”), m (“o”), and z (“*”) at t = 0.5 in x ∈ [−1, 1] for initial data (9.2) and (9.3) by NSP and NSPIF. ε = 10−8, ∆x = 0.01, and CFL = 0.5.](/figures/fig-9-2-the-numerical-solutions-of-r-m-o-and-z-at-t-0-5-in-x-8s0dp99y.webp)
![FIG. 9.5. The numerical solutions of ρ (“+”), m (“o”), and z (“*”) at t = 0.25 in x ∈ [0, 1] for initial data (9.4) and (9.5) by (from left to right, then top to bottom) SCNvL and SCNvLIF for CFL = 0.5 and SCNvL and SCNvLIF for CFL = 0.005. ε = 10−8, ∆x = 0.01.](/figures/fig-9-5-the-numerical-solutions-of-r-m-o-and-z-at-t-0-25-in-1y02pigw.webp)
![FIG. 9.4. The numerical solutions of ρ (“+”), m (“o”), and z (“*”) at t = 0.25 in x ∈ [0, 1] for initial data (9.4) and (9.5) by (from left to right, then top to bottom) SP1, SCN, SP1vL, SCNvL, NSP, NSPIF, and SCNvLIF. ε = 10−8, ∆x = 0.01. CFL = 1 for SP1 and SCN, CFL = 0.5 for others.](/figures/fig-9-4-the-numerical-solutions-of-r-m-o-and-z-at-t-0-25-in-23yw4f8f.webp)
![FIG. 9.3. The numerical solutions of ρ (“+”), m (“o”), and z (“*”) at t = 0.5 in x ∈ [−1, 1] for initial data (9.2) and (9.3) by (from left to right, then top to bottom) SP1, SCN, SP1vL, SCNvL, NSP, and NSPIF. ε = 0.02, ∆x = 0.02. CFL = 1 for SP1 and SCN, CFL = 0.5 for others.](/figures/fig-9-3-the-numerical-solutions-of-r-m-o-and-z-at-t-0-5-in-x-4iytzwsu.webp)
UNIFORMLY ACCURATE SCHEMES FOR HYPERBOLIC SYSTEMS
WITH RELAXATION
∗
RUSSEL E. CAFLISCH
†
, SHI JIN
‡
, AND GIOVANNI RUSSO
§
SIAM J. N
UMER.ANAL.
c
1997 Society for Industrial and Applied Mathematics
Vol. 34, No. 1, pp. 246–281, February 1997 011
Abstract. We develop high-resolution shock-capturing numerical schemes for hyperbolic sys-
tems with relaxation. In such systems the relaxation time may vary from order-1 to much less than
unity. When the relaxation time is small, the relaxation term becomes very strong and highly stiff,
and underresolved numerical schemes may produce spurious results. Usually one cannot decouple
the problem into separate regimes and handle different regimes with different methods. Thus it is
important to have a scheme that works uniformly with respect to the relaxation time. Using the
Broadwell model of the nonlinear Boltzmann equation we develop a second-order scheme that works
effectively, with a fixed spatial and temporal discretization, for all ranges of the mean free path.
Formal uniform consistency proof for a first-order scheme and numerical convergence proof for the
second-order scheme are also presented. We also make numerical comparisons of the new scheme with
some other schemes. This study is motivated by the reentry problem in hypersonic computations.
Key words. hyperbolic systems with relaxation, Broadwell model, stiff source, high-resolution
shock-capturing methods
AMS subject classifications. 35L65, 49M25, 34A65, 82B40
PII. S0036142994268090
1. Introduction. Hyperbolic systems with relaxation are used to describe many
physical problems that involve both convection and nonlinear interaction. In the
Boltzmann equation from the kinetic theory of rarefied gas dynamics, the collision (re-
laxation) terms describe the interaction of particles. In viscoelasticity, memory effects
are modeled as relaxation. Relaxation occurs in water waves when the gravitational
force balances the frictional force of the riverbed. For gas in thermal nonequilibrium
the internal state variable satisfies a rate equation that measures a departure of the
system from the local equilibrium. Relaxations also occur in other problems ranging
from magnetohydrodynamics to traffic flow.
In such systems, when the nonlinear interactions are strong, the relaxation rate
is large. In kinetic theory, for example, this occurs when the mean free path between
collisions is small (i.e., the Knudsen number is small). Within this regime, which
is referred to as the fluid dynamic limit, the gas flow is well described by the Euler
or Navier–Stokes equations of fluid mechanics, except in shock layers and boundary
layers. The characteristic length scale of the kinetic description of the gas is the micro-
scopic, collision distance; in the fluid dynamic limit it is replaced by the macroscopic
length scale of fluid dynamics. By analogy with the kinetic theory, we shall refer to
the limit of large relaxation rate (or small relaxation time) for a general hyperbolic
system with relaxation as the fluid dynamic limit.
The fluid dynamic limit is challenging for numerical methods, because in this
regime the relaxation terms become stiff. In particular, a standard numerical scheme
∗
Received by the editors May 18, 1994; accepted for publication (in revised form) April 5, 1995.
http://www.siam.org/journals/sinum/34-1/26809.html
†
Department of Mathematics, University of California, Los Angeles, CA 90024. The research of
this author was supported in part by ARO grant DAAL03-91-G-0162.
‡
School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332 (jin@math.
gatech.edu). The research of this author was supported in part by NSF grant DMS-9404157 and by
AFOSR grant F49620-92-J0098 as a visiting member at the Courant Institute.
§
Dipartimento di Matematica, Universit`a dell’Aquila, 67010 L’Aquila, Italy.
246
HYPERBOLIC SYSTEMS WITH RELAXATION 247
might fail to give physically correct solutions once the (microscopic) relaxation dis-
tance is much smaller than the spatial discretization. Although a full simulation of
the relaxation process would require a very fine (and expensive) discretization, it may
be possible to accurately compute the solution on a coarser fluid dynamic length scale.
The goal of this paper is to present a class of numerical methods that work with uni-
form accuracy from the rarefied regime to the fluid dynamic limit for the Broadwell
model of kinetic theory.
Numerical methods for hyperbolic systems with relaxation terms have attracted
a lot of attention in recent years [7], [23], [24], [16], [17]. Studying the numerical
behavior for these problems is important not only for the physical applications but
also for the development of new numerical methods for conservation laws, such as
kinetic schemes [15], [8], [25] and relaxation schemes [18]. Most kinetic or relaxation
schemes can be described as fractional step methods, in which the collision step is
just a projection of the system into a sort of discrete “local Maxwellian” or local
equilibrium. Although the goal of a kinetic scheme or relaxation scheme is different
from ours, nevertheless we use them as guidelines for the study of the properties of a
numerical scheme near the fluid dynamic regime.
In earlier works on the system of hyperbolic equations with relaxation, the goal
was to develop robust numerical schemes that handle the stiffness of the problem
effectively and avoid spurious numerical solutions when the grid spacing underresolves
the small relaxation time. In regions where the relaxation time is no longer small and
the problem becomes nonstiff, however, these schemes usually may not have high-order
accuracy uniformly with respect to the wide range of the relaxation time.
Our motivation differs from these earlier approaches in that we seek to develop
robust numerical schemes that work uniformly for a wide range of relaxation rates.
We consider a simpler model of the Boltzmann equation, and we derive a numerical
scheme which is of second order uniformly in the mean free path. This is motivated
by hypersonic computations of reentry problems. The new scheme we introduce here,
called NSPIF, is able to handle all different regimes from the rarefied gas to the fluid
limit (the stiff regime) with fixed spatial and temporal grids that are independent of
the mean free path. Although we develop our method based on the Broadwell model,
this scheme can be applied to a much wider class of hyperbolic systems with relaxation
terms and to other discrete velocity kinetic models. In particular, it applies to a class
of hyperbolic systems with relaxation characterized by Liu [22] and Chen, Levermore,
and Liu [5].
Probably the paper that is closest in spirit to our work is the one by Coron and
Perthame [7]. In that paper the authors derive a numerical scheme for solving the
BGK model of the Boltzmann equation under a wide range of mean free path. They
discretize velocity space and use a splitting scheme. The collision step is treated
by a semi-implicit method that guarantees positivity and entropy condition for the
time-discrete model. The scheme is first-order accurate in space and time.
Development of numerical methods for the problems considered here is consid-
erably aided by knowledge of the equations for the fluid dynamic limit. In other
stiff source problems, such as those arising in reacting flow computations, the corre-
sponding limit may be less well understood, and extra efforts are necessary to develop
underresolved numerical methods [6], [21], [9], [23].
The outline of the paper is as follows. In the next section we present the Broad-
well model and its fluid dynamic limit. Section 3 deals with the treatment of the
convective step. Upwind methods and flux limiters are briefly recalled. Sections 4–
7 are devoted to the development and analysis of a second-order scheme (in space
and time) called NSP. In section 4, the scheme is derived making use of truncation
248 R. E. CAFLISCH, S. JIN, AND G. RUSSO
analysis, and in the following section the linear stability of this scheme is studied.
In section 6 the Richardson extrapolation is used for the first time step in order to
guarantee second-order accuracy even in presence of an initial layer. The NSP com-
bined with this initial step is called the NSPIF. In section 7 the fluid dynamic limit
of the NSPIF is derived. In section 8 we study a first-order splitting scheme (called
the SP1). Important properties, such as positivity and entropy inequality, are proved
for the scheme. The main result in this section is a uniformly (in the relaxation time)
first-order consistency error estimate. Our analysis is based on the assumption of
smoothness of the solution. To the best of our knowledge this is the first estimate
on numerical methods for hyperbolic systems with stiff relaxation terms, even for
smooth solutions. Also, the uniformly first-order accuracy sharpens earlier work in
this direction [14]. Estimates on nonsmooth solutions or on higher order methods for
such inhomogeneous hyperbolic systems would require more advanced analytic tech-
niques and are beyond the scope of this paper. In the last section we present some
numerical results which show that NSPIF overperforms several other approaches and
always gives second-order physical solutions over a wide range of the mean free path.
2. The Broadwell model. A simple discrete velocity kinetic model for a gas
was introduced by Broadwell [2]. It describes a two-dimensional (2-D) (3-D) gas as
composed of particles of only four (six) velocities with a binary collision law and spatial
variation in only one direction. When looking for one-dimensional (1-D) solutions of
the 2-D gas, the evolution equations of the model are given by
∂
t
f + ∂
x
f =
1
ε
(h
2
− fg),
∂
t
h = −
1
ε
(h
2
−fg),(2.1)
∂
t
g − ∂
x
g =
1
ε
(h
2
− fg),
where ε is the mean free path, f, h,andgdenote the mass densities of gas particles
with speed 1, 0, and −1, respectively, in space x and time t. Similar equations are
obtained for 1-D solutions of a 3-D gas [3].
The fluid dynamic moment variables are density ρ, momentum m, and velocity u
defined by
ρ = f +2h+g, m = f − g, u =
m
ρ
.(2.2)
In addition, define
z = f + g.(2.3)
Then the Broadwell equations can be rewritten as
∂
t
ρ + ∂
x
m =0,(2.4)
∂
t
m + ∂
x
z =0,(2.5)
∂
t
z + ∂
x
m =
1
2ε
(ρ
2
+ m
2
− 2ρz).(2.6)
Note that if the fluid variables ρ, m,andzare known then f,g,andhcan be recovered
from (2.2) and (2.3) as
f =
1
2
(z + m),g=
1
2
(z−m),h=
1
2
(ρ−z).
HYPERBOLIC SYSTEMS WITH RELAXATION 249
A local Maxwellian is a density function that satisfies
Q(f, h, g)=h
2
−fg = ρ
2
+m
2
−2ρz =0,(2.7)
i.e.,
z = z
E
(ρ, m) ≡
1
2ρ
(ρ
2
+ m
2
)=
1
2
(ρ+ρu
2
).(2.8)
As ε → 0 equation (2.1) or (2.6) gives the local Maxwellian distribution (2.8). Apply-
ing (2.8) in (2.5), one gets the fluid dynamic limit described by the following model
Euler equations:
∂
t
ρ + ∂
x
(ρu)=0,
(2.9)
∂
t
(ρu)+∂
x
1
2
(ρ+ρu
2
)
=0.
To the next order, a model Navier–Stokes equation can be derived via the Chapman–
Enskog expansion [4]. For a description of the Broadwell model and its fluid dynamic
limit see, for example, [3] and [13].
Previously the numerical solution of the Broadwell equations has been considered
by several authors [10, 12, 1]. In these earlier works the attention was focused on
developing methods that work in the rarefied regime = O(1). Many of these methods
will have problems when → 0. Among the problems that arose in the fluid regimes
are numerical instability, poor shock, and rarefaction resolutions and even spurious
numerical solutions. See the numerical examples in section 9.
The Broadwell equations are a prototypical example for more general hyperbolic
systems with relaxations in the sense of Whitham [29] and Liu [22]. These problems
can be described mathematically by the system of evolutional equations
∂
t
U + ∂
x
F (U )=−
1
ε
R(U),U∈R
N
.(2.10)
We will call this system the relaxation system. Here we use the term relaxation in the
sense of Whitham [29] and Liu [22]. The relaxation term is endowed with an n × N
constant matrix Q with rank n<N such that
QR(U) = 0 for all U.(2.11)
This yields n independent conserved quantities v = QU . In addition we assume
that each such v uniquely determines a local equilibrium value U = E(v) satisfying
R(E(v)) = 0 and such that
QE(v)=v for all v.(2.12)
The image of E then constitutes the manifold of local equilibria of R.
Associated with Q are n local conservation laws satisfied by every solution of
(2.10) that take the form
∂
t
(QU)+∂
x
(QF(U))=0.(2.13)
These can be closed as a reduced system for v = QU if we take the local relaxation
approximation
U = E(v),(2.14)
250 R. E. CAFLISCH, S. JIN, AND G. RUSSO
∂
t
v + ∂
x
e(v)=0,(2.15)
where the reduced flux e is defined by
e(v) ≡QF(E(v)).(2.16)
A system of conservation laws with relaxation is stiff when ε is small compared with
the time scale determined by the characteristic speeds of the system and some ap-
propriate length scales. While we mainly concentrate on the Broadwell equation, the
analysis as well as the numerical schemes can certainly be applied to this class of
relaxation problems. In fact from time to time we will use the general equation (2.10)
to simplify the notation.
3. The discretizations of the convection terms. We introduce the spatial
grid points x
j+
1
2
,j=···,−1,0,1,··· with uniform mesh spacing ∆x = x
j+
1
2
− x
j−
1
2
for all j. The time level t
0
,t
1
,··· are also spaced uniformly with space step ∆t =
t
n+1
− t
n
for n =0,1,2,···. Here the assumption of a uniform grid is only for
simplicity. We use U
n
j
to denote the cell average of U in the cell [x
j−
1
2
,x
j+
1
2
]attime
t
n
,
U
n
j
=
1
∆x
Z
x
j+
1
2
x
j−
1
2
U(t
n
,x)dx.(3.1)
We use the method of lines, in which the time discretization and spatial discretiza-
tion are taken separately, for the Broadwell equations. In this section we shall discuss
the spatial discretization, which concerns the linear convection terms. Note that the
linear convection in the Broadwell equation is of hyperbolic type. Thus it is natural
to use upwind schemes.
A conservative spatial discretization to the Broadwell equations (2.4)–(2.6) is
∂
t
ρ
j
+
m
j+
1
2
− m
j−
1
2
∆x
=0,(3.2)
∂
t
m
j
+
z
j+
1
2
− z
j−
1
2
∆x
=0,(3.3)
∂
t
z
j
+
m
j+
1
2
− m
j−
1
2
∆x
=
1
2ε
(ρ
2
j
+ m
2
j
− 2ρ
j
z
j
).(3.4)
Equations (3.2)–(3.4) can be diagonalized into
∂
t
f
j
+
f
j+
1
2
− f
j−
1
2
∆x
=
1
ε
(h
2
j
− f
j
g
j
),
∂
t
h
j
= −
1
ε
(h
2
j
− f
j
g
j
),(3.5)
∂
t
g
j
−
g
j+
1
2
− g
j−
1
2
∆x
=
1
ε
(h
2
j
− f
j
g
j
),
where f, h,andgare exactly the original density functions for the Broadwell equa-
tion. The connection between (3.4) and (3.5) is established through the definition of
the fluid moment variables (2.2), (2.3). Since f and g travel along the constant char-
acteristics with speeds 1 and −1, respectively, upwind schemes can be easily applied
to them.