Numerical Simulations for Radiation Hydrodynamics. I. Diffusion Limit

JOURNAL OF COMPUTATIONAL PHYSICS ARTICLE NO.

142, 182–207 (1998)

CP985940

Numerical Simulations for Radiation Hydrodynamics. I. Diffusion Limit Wenlong Dai and Paul R. Woodward School of Physics and Astronomy/Laboratory for Computational Science and Engineering, University of Minnesota, 116 Church Street S. E., Minneapolis, Minnesota 55455 Received December 1, 1997

A second order accurate finite difference scheme is proposed for multidimensional radiation hydrodynamical equations in a diffusion limit. The radiation hydrodynamical equations in the limit constitute a hyperbolic system of conservation laws plus radiative heat transfer. A Godunov scheme including linear and nonlinear Riemann solvers is proposed for the set of conservation laws. The scheme with the linear Riemann solver works well for relatively strong shocks. The nonlinear Riemann solver is specially designed for flows involving strong shocks. The radiative heat conduction is treated implicitly. The treatment possesses a number of advantages over typical implicit methods. The most notable are the second order accuracy in both space and time, quick damping of numerical errors when the size of time steps is large, iterative solver and the fast convergence, the accurate treatment for the nonlinearity, and the energy conservation. Numerical examples are given to show the features of the schemes. °c 1998 Academic Press

1. INTRODUCTION

Radiation hydrodynamical equations play an important role in astrophysics. In the diffusion limit, the radiation hydrodynamical equations may be written as a hyperbolic system of conservation laws plus (nonlinear) radiative heat transfer. One of the major difficulties of standard numerical methods for the set of conservation laws is to resolve and keep track of shocks involved. A typical approach is to passively treat radiation pressure and energy, i.e., to solve the Euler equations first and then to update momentum and energy due to radiation. When the Euler equations are solved, the dynamical role of the radiation pressure and energy is temporarily ignored. For radiative heat transfer, both explicit and implicit methods are currently used. An explicit scheme, for example, the forward Euler scheme, is simple. But, the size of a time step is limited by a stability condition, and an explicit scheme is only first order accurate. At the other hand, implicit schemes normally involve solving a large set of nonlinear algebraic 182 0021-9991/98 $25.00 c 1998 by Academic Press Copyright ° All rights of reproduction in any form reserved.

RADIATION HYDRODYNAMICS

183

equations at each time step. Two typical implicit schemes are the backward Euler scheme and the Crank–Nicolson scheme [1]. The outstanding feature of the backward Euler scheme is that numerical errors undergo quick damping for a large time step. But the backward Euler scheme is only first order accurate. The Crank–Nicolson scheme is second order accurate, but numerical errors do not damp for a large time step. When a time step is large, the solution obtained through the Crank–Nicolson scheme is far away from an exact solution. The nonlinearity in the heat transfer is another headache in an implicit scheme. Typically, either linearization or the Newton iteration is used for the nonlinearity. The linearization makes an implicit scheme only first order accurate. The Newton iteration is very expensive in CPU in an implicit scheme because there is a large set of unknown variables. In this paper, we will develop a numerical method for radiation hydrodynamical equations in the diffusion limit. A second order accurate Godunov scheme will be proposed for the set of conservation laws, which includes linear and nonlinear Riemann solvers. An implicit scheme is proposed for the radiative heat transfer, which possesses a number of advantages over typical implicit methods: second order accurate in both space and time, quick damping in numerical errors for a large size of time steps, iterative and fast in convergence, and accurate in the treatment for nonlinearity. The plan of this paper is as follows. In the second section radiation hydrodynamical equations are given. A Godunov scheme for the set of conservation laws is described in the third section. In the fourth section we present our treatment for radiative heat transfer. Numerical examples are shown in the fifth section to demonstrate the features of the numerical schemes. The final section is the conclusion of this paper. 2. BASIC EQUATIONS

In radiating fluid, radiation often contains a fraction of momentum and energy. To describe the behavior of such flows we need conservation laws that account for both gas material and radiative contributions to the flow dynamics. If the material is extremely radiatively opaque, so that the mean-free-path of photons is much smaller than the typical length of flows. The radiation hydrodynamical equations [18] may be written as ∂ρ + ∇ · (ρu) = 0, ∂t µ ¶ 1 ∂ (ρu) + ∇ · (ρuu) + ∇ p + a R T 4 = 0, ∂t 3 µ ¶ ·µ ¶ ¸ 1 2 1 2 1 ∂ 4 4 4 ρe + ρu + a R T + ∇ · ρe + ρu + a R T + p + a R T u ∂t 2 2 3 = ∇ · [κ(T )∇T ].

(1) (2)

(3)

Here, ρ, p, u, e, and T are the mass density, thermal pressure, flow velocity, specific internal energy, and temperature, a R is a radiation constant, and κ(T ) is the nonlinear heat diffusivity. The radiation heat diffusivity 4a R cλ R T 3 /3 is included in κ(T ). Here c and λ R are the light speed and the Rosseland mean free path for photons. In Eqs. (1)–(3) we have assumed the radiative temperature equal to the fluid temperature, i.e., we are dealing with equilibrium diffusion. Like a R , λ R is assumed being a constant here. In Eqs. (2), (3), a R T 4 and a R T 4 /3 are radiative energy density and radiative pressure, respectively. The set of Eqs. (1)–(3) is

184

DAI AND WOODWARD

complete if an equation of state p = p(ρ, T ) is given. For the purpose of our test problems, we assume the γ -law for the equation of state, p = (γ − 1)ρe, and γ is the ratio of specific heats. Two kinds of physical processes are described in Eqs. (1)–(3), interaction between nonlinear waves, and heat transfer due to temperature gradient. We write the set of Eqs. (1)– (3) in two sets of equations: ∂U ∂F ∂G ∂H + + + = 0, ∂t ∂x ∂y ∂z

(4)

∂ (ρ E) = ∇ · [κ(T )∇T ]. ∂t

(5)

and

Here  ρ  ρu x     U≡  ρu y  ,  ρu z  ρE 



 ρu x  ρu x u x + p + 1 a R T 4    3   ρu x u y F(U) ≡  ,   ρu x u z   ¡ ¢ 1 4 ρu x E + p + 3 a R T u x

and G and H have similar expressions as F. In these definitions, E is the specific total energy, E ≡ e + ρu2 /2 + a R T 4 . Numerically, we may solve Eqs. (1)–(3) through consecutively solving Eq. (4) and Eq. (5). The mass density, momentum, and total energy are first updated through Eq. (4). The total energy obtained through solving Eq. (4) is only an intermediate result, which should be updated through Eq. (5). In this paper, Eq. (4) is solved through a second order Godunov scheme, and Eq. (5) is solved through an iterative implicit scheme. 3. SCHEMES FOR THE HYPERBOLIC SYSTEM OF CONSERVATION LAWS

In order to solve Eq. (4), it is natural to extend established methods for the Euler equations. During the last two decades, numerical methods for the Euler equations have been well developed (for example, see [2, 11, 13–15, 17, 19, 20, 22, 25–28], among which Godunov schemes are particularly efficient for shock problems. Godunov [2] supposed that initial data could be replaced by a piecewise constant data with discontinuities and used exact solutions of Riemann problems to advance piecewise constant data. A major extension of Godunov’s scheme was made by Van Leer in his MUSCL scheme [11] which used a Riemann solver to advance piecewise linear data. Roe developed an approximate Riemann solver suitable for the use in Godunov schemes. A nonlinear Riemann solver and a contact steepener were developed in the piecewise parabolic method (PPM) [20, 22]. Key points in Godunov schemes are the use of characteristic formulations, reconstruction of initial data, and an approximate Riemann solver which is suitable for computing a set of time-averaged fluxes at interfaces between grid cells. In this section, we will extend PPM to the system described by Eq. (4). Since we will employ a dimensionally split technique to treat the multi-dimensionality of Eq. (4), we may focus on the projection of Eq. (4) in the x-dimension. We write the

RADIATION HYDRODYNAMICS

185

one-dimensional version of Eq. (4) in a Lagrangian mass coordinate: ∂F ∂U + = 0. ∂t ∂m

(6)

Here  V  ux     U≡  uy  ,  uz  

 −u x  P     F(U) ≡   0 ,  0 

E



ux P

V is the specific volume, V ≡ 1/ρ, P is the total pressure, P ≡ p + 13 a R T 4 , and m is a mass coordinate, which is defined through dm ≡ ρ d x. Consider a numerical grid {xi } in the Lagrangian coordinate, integrate Eq. (6) in a grid cell xi ≤ x ≤ xi+1 and over a time step 0 ≤ t ≤ 1t, and we have Ui (1t) = Ui (0) +

1t ¯ (Fi − F¯ i+1 ). 1m i

(7)

Here 1m i is the mass contained in the cell, 1t is the time step, Ui (t) is a cell-averaged value of U at time t, F¯ i is a time-averaged flux F over the time step at the cell interface x = xi , and they are defined as Z xi+1 1 U(t, x)ρ(t, x) d x, (8) Ui (t) ≡ 1m i xi Z 1t ¯Fi ≡ 1 F[U(t, xi )] dt. (9) 1t 0 We should mention that Eq. (7) is exact. Therefore, the key point of the scheme is to approximately calculate the time-averaged fluxes at grid interfaces. A Riemann problem is an initial value problem, Eq. (6), subject to a specific initial condition: ½ L U if x < 0 U0 (x) = if x > 0. UR Here U L and U R are a pair of left and right states. For the purpose of the difference Eq. (7), we have to calculate only the time-averaged values for u x and P. We will use u ∗x and P ∗ to denote the time-averaged values. Riemann problems encountered in practical calculations may involve more than one strong shock, for example, for astrophysical problems. We propose two Riemann solvers here, a linear solver and a nonlinear solver, although other Riemann solvers, such as the Roe’s solver [13], may be applied here. Both Riemann solvers we will propose are based on characteristic formulations of Eq. (6); the linear Riemann solver is very simple, and the nonlinear Riemann solver assumes two shocks with arbitrary strengths present in a Riemann problem. The characteristic formulations may be obtained through the general outlines described in [6, 29]. In the mass coordinate, three characteristic speeds are found to be 0, Cs , and -Cs .

186

DAI AND WOODWARD

Here Cs is defined as · Cs ≡

P + (de/dT + 4a R V T 3 )(∂ T /∂ V ) P + a R T 4 (de/dT + 4a R V T 3 )(∂ T /∂ P)V

¸ 12

.

(10)

Here µ µ

∂T ∂V ∂T ∂P

¶ ¶

≡

(γ − 1)ρ 2 e , (γ − 1)ρ(de/dT ) + (4/3)a R T 3

≡

1 . (γ − 1)ρ(de/dT ) + (4/3)a R T 3

P

V

Along three characteristic curves defined by d x/dt = 0, d x/dt = cs , and d x/dt = −cs in the Lagrangian coordinate, the following differentials (of Riemann invariants) are vanishing, respectively, d R0 = 0,

du y = 0,

du z = 0 along d x/dt = 0,

d R+ = 0

along d x/dt = cs ,

(12)

d R− = 0

along d x/dt = −cs .

(13)

(11)

Here, cs ≡ Cs /ρ, and d R0 and d R± are defined as d R0 ≡ d P + Cs2 d V,

(14)

d R± ≡ d P ± Cs du x .

(15)

Our linear Riemann solver is based on the characteristic formulations Eq. (15). For a given pair of left and right states, the time-averaged values u ∗x and P ∗ are approximately obtained from the set of equations ¡ ¢ (P ∗ − P L ) + CsL u ∗x − u xL = 0, ¡ ¢ (P ∗ − P R ) − CsR u ∗x − u xR = 0. Here the subscript L (or R) denotes the evaluations at the left (or right) state. Like other linear Riemann solvers, this simple Riemann solver assumes the flow being smooth. But this Riemann solver works well even for relatively strong shocks if it is put into the conservation law Eq. (7). It is the conservation law which plays the crucial role in the numerical simulations for strong shocks. We will show the usefulness and limitation of the linear Riemann solver in numerical examples. For strong shocks involved in numerical simulations, a nonlinear Riemann solver is necessary. Our nonlinear Riemann solver is based on characteristic formulations as well as conservation laws across any discontinuous wave, i.e., jump conditions. From conservation laws, we may obtain a set of jump conditions across any discontinuity: W [V ] = −[u x ],

(16)

W [u x ] = [P],

(17)

W [E] = [u x P].

(18)

RADIATION HYDRODYNAMICS

187

FIG. 1. The relation between the pre- and post-shock states when α R 6= 0.

Here, W is the shock speed in the mass coordinate, and [U ] denotes the difference between the post-shock state U1 and the pre-shock state U0 . A typical relation between the pre-shock and post-shock states is shown in Fig. 1. If a Riemann problem involves only shocks and a contact discontinuity, then u ∗x and P ∗ are the values of u x and P at the intersection A, u x A , and PA , of two Rankine–Hugoniot curves L and R, as shown in Fig. 2, where the curve L is the one passing through the left

FIG. 2. An illustration for the solution of the time-averaged fluxes needed in a Godunov scheme.

188

DAI AND WOODWARD

FIG. 3. A contact discontinuity and two wave fronts divide the whole (x −t)-space into four possibly different regions.

state S L , and R is the one passing through the right state S R . If a Riemann problem involves rarefaction waves, then the u x A and PA in Fig. 2 will not be the solutions u ∗x and P ∗ . But, if rarefaction waves involved are weak, the point SL (or S R ) in Fig. 2 will be very close to the point A. The relation between S L (or S R ) and S A may be represented by a tangential line SL A1 (or S R A1 ), which is close to the Rankine–Hugoniot curve if the rarefaction wave is weak. Therefore, if rarefaction waves involved are weak, u x A and PA in Fig. 2 are the approximate solutions for u ∗x and P ∗ . Since rarefaction waves are smooth, the rarefaction waves involved in a Riemann problem of two adjacent cells is always weak except at the initial time. Therefore, we believe that in numerical simulations we may approximately use u x A and PA , the intersection between two Rankine–Hugoniot curves, as the time-averaged values u ∗x and P ∗ . In principle, u x A and PA in Fig. 2 may be iteratively obtained following a similar procedure described in [20, 22]. But each step in the iterative procedure needs to find an appropriate root of a polynomial, which makes the iterative procedure impractical. Therefore, we propose another iterative approach to find the intersection (u x A and PA ) in Fig. 2. The wave fronts of three possible waves involved in a Riemann problem divide the whole (x, t)-space into four possibly different regions, as shown in Fig. 3, where W− (or W+ ) is the wave speed propagating to the left (or right), and R2 and R3 are connected by a contact discontinuity. The flow velocity u x and total pressure P should be the same at regions R2 and R3. From the set of jump conditions, we may obtain the jump [u x ] from a jump [T ] across any discontinuity: [u x ]2 =

´ p 1 ³ B − A − C + (B − A − C)2 + 4β AB . 2−β

Here β, A, B, and C are defined as β≡

P0 − (1/3)a R T14 , P0 + a R T14

A ≡ (β + γ − 1)[e] +

µ

¶ 1 + β a R V0 [T 4 ], 3

B ≡ (2 − β)([e] + a R V0 [T 4 ]), ¢ ¡ C ≡ 2 P0 + aT T14 V0 .

(19)

RADIATION HYDRODYNAMICS

189

For a given pair of left and right states, our iterative procedure to find u x and P at the regions R2 and R3 is as follows. (1) Guess the temperatures T2 and T3 at the regions R2 and R3. (2) Calculate the states in R2 and R3 through Eq. (19) and jump conditions Eqs. (16)–(18), since [T ] is known. (3) Use the condition for a contact discontinuity to improve the initial guess on T2 and T3 : du x2 du x3 δT2 = u x3 (T3 ) + δT3 , (20) u x2 (T2 ) + dT2 dT3 P2 (T2 ) +

d P2 d P3 δT2 = P3 (T3 ) + δT3 . dT2 dT3

(21)

(4) Use the improved guess on the two temperatures, T2 = T2 + δT2 and T3 = T3 + δT3 , and go back to the second step above. In Eqs. (20), (21), we have used the Newton iteration. We find that a few iterations of the Riemann solver will give sufficiently accurate time-averaged fluxes needed in the scheme Eq. (7). Actually, we use only two iterations even for strong shocks, which will be shown in numerical examples. Our one-dimensional function code starts from a set of cell-averages of conserved quantities. Cubic polynomials are used to interpolate each independent variable a to find its values at grid interfaces. After we obtain the values at grid interfaces, the monotonicity constraint originally suggested by Van Leer is applied to these values at interfaces. As we know, interpolated structures are not always monotone increasing (decreasing) even though they have been constructed from monotone data. Over- and undershoots in interpolated internal zone structures eventually give rise to over- and undershoots in updated cellaveraged data. We assume that a(x) is continuous inside each zone but may have a big jump across a grid interface. After the monotonicity constraint is applied, we have three values of a for each grid cell: a cell-average ai and the values at two interfaces of the cell, ai L and ai R . For the calculation of the left and right states in the Riemann problem arising at the grid interface xi , a parabola defined by these three values is used to interpolate the structure of a(x) inside the cell. The effective left state is the domain-average of a over the domain xi − cs L 1t < x < xi , and the effective right state is the domain-average over the domain xi < x < xi + cs R 1t. The time-averaged fluxes at the interface are obtained through either the linear or the nonlinear Riemann solver described above for this set of left and right states. After obtaining the time-averaged fluxes, we update the conserved quantities in the Lagrangian coordinate according to the conservation law Eq. (7), by adding the flux advected into and subtracting the flux advected out from a cell during a time step. The cell-averaged values on a fixed Eulerian grid are obtained through a mapping procedure. The mapping procedure we used keeps the conservation laws exactly, and is second order accurate for each independent variable. Detail formulations of the interpolations and mapping used here may be found in [29]. The multi-dimensional Eq. (4) is solved through a dimensionally split technique [7]. Each time step of a multi-dimensional problem may be broken down into one-dimensional passes in which derivatives in other dimensions are set to zero. The dimensionally split technique may be written in the form y

y

x z z x L 1t L 1t L 1t L 1t L 1t Ui, j,k (0). Ui, j,k (21t) = L 1t

(22)

190

DAI AND WOODWARD

x Here L 1t is an one-dimensional operator in the x-direction with the time step 1t, which is described before. Since the one-dimensional operator is second order accurate, then the scheme constructed as Eq. (22) is second order accurate at every other time step.

4. NUMERICAL SCHEMES FOR HEAT TRANSFER

An extensive amount of literature exists on numerical methods for the solution of heat transfer (see, e.g., [12, 24] and references therein). As mentioned before, explicit schemes for heat transfer are only first order accurate, and the size of time steps is limited by numerical stability condition. On the other hand, in typical second order implicit schemes, numerical errors do not damp when the size of time steps is large. In this section, we will propose an iterative implicit scheme for nonlinear Eq. (5) which is second order accurate, and which gives the exact solution when the size of time steps is very large. Since we are using an operator split technique, the mass density, momentum, and energy are all updated through Eq. (4) solved, while through Eq. (5) we update only temperature and therefore the energy. Therefore, we may write Eq. (5) in the form ∂ (E T ) = ∇ · [κ(T )∇T ]. ∂t

(23)

Here E T is the part of total energy density depending on temperature, E T ≡ ρe + a R T 4 . For the simplicity in formulations, we consider only the two-dimensional situation here, and the extension to the three-dimensional situation is straightforward. Consider a grid cell xi < x < xi+1 and y j < y < y j+1 . Integrating Eq. (23) in the grid cell and over a time step (0, 1t) yields E TN = E T 0 +

1t 1t (q¯ x W − q¯ x E ) + (q¯ y S − q¯ x N ). 1x 1y

(24)

Here, E T 0 and E TN are cell-averaged values of E T at t = 0 and t = 1t, respectively, q ≡ −κ(T )∇T, q¯ x W , and q¯ x E (or q¯ y S and q¯ y N ) are the time-averaged values of heat fluxes at grid interfaces to the west and east (or to the south and north), respectively, and they are defined as Z y j+1 Z xi+1 1 E T (1t, x, y) d x d y, E TN ≡ 1x1y y j xi Z 1t Z y j+1 1 qx (t, xi , y) dy dt, q¯ x W ≡ 1t1y 0 yj Z 1t Z y j+1 1 q¯ x E ≡ qx (t, xi+1 , y) dy dt, (25) 1t1y 0 yj Z 1t Z xi+1 1 q y (t, x, y j ) d x dt, q¯ y S ≡ 1t1x 0 xi Z 1t Z xi+1 1 q¯ y N ≡ q y (t, x, y j+1 ) d x dt. (26) 1t1x 0 xi We should mention that the difference Eq. (24) is exact because no approximations have been involved yet.

191

RADIATION HYDRODYNAMICS

If the time-averaged fluxes in Eqs. (24) are replaced by a value at t = 0 (or at t = 1t), the result is the forward (or backward) Euler scheme which is first order accurate in time. If the time-averaged fluxes are replaced by their averaged values at t = 0 and t = 1t, the result is the Crank–Nicolson scheme, which is second order accurate in time. As stated before, the Crank–Nicolson scheme will give completely wrong solutions if the size of time steps is very large, because numerical errors in the scheme do not undergo damping for large time steps. When the size of time steps is very large, the solution should be independent of the initial condition, and the solution is determined only by the boundary condition. In order to introduce quick damping for numerical errors within the framework of the second order accuracy, we introduce an additional time level t = 1t/2. We approximately evaluate the time-averaged fluxes in Eq. (24) as their values at t = 1t/2, and Eq. (24) becomes E TN = E T 0 +

¢ 1t ¡ H ¢ 1t ¡ H qx W − qxHE + q y S − qxHN . 1x 1y

(27)

Here the superscript H denotes the evaluation at t = 1t/2, qxHW

1 ≡ 1y

q yHS ≡

1 1x

Z Z

y j+1

qx (1t/2, xi , y) dy,

qxHE

yj xi+1

q y (1t/2, x, y j ) d x,

1 ≡ 1y

q yHN ≡

xi

1 1x

Z

y j+1

qx (1t/2, xi+1 , y) dy,

(28)

yj

Z

xi+1

q y (1t/2, x, y j+1 ) d x. (29)

xi

For the first half time step, through the similar procedure for Eq. (27), we have E TH = E T 0 +

¢ ¢ 1t ¡ H 1t ¡ H q¯ − q¯ xHE + q¯ − q¯ xHN . 21x x W 21y y S

(30)

Here q¯ xHW and q¯ xHE (or q¯ yHS and q¯ yHN ) are the time-averaged values of heat flux over the first half time step, and they are defined as q¯ xHW ≡

2 1t1y

q¯ xHE

2 ≡ 1t1y

q¯ yHS

2 ≡ 1t1x

q¯ yHN

2 ≡ 1t1x

Z

1t/2

0

Z

1t/2

0

Z

1t/2

0

Z

0

1t/2

Z Z Z Z

y j+1

qx (t, xi , y) dy dt,

yj y j+1

qx (t, xi+1 , y) dy dt,

(31)

yj xi+1

q y (t, x, y j ) d x dt,

xi xi+1

q y (t, x, y j+1 ) d x dt.

(32)

xi

The time-averaged fluxes involved in Eq. (30) may be approximately calculated through an interpolation in time. As stated before, any approximate calculation for the time-averaged fluxes must not, even partially, depend on initial information when a time step is very large. Our interpolation for the time-averaged fluxes in Eq. (30) is linearly determined by their values at t = 1t and t = 1t/2. Therefore, the time-averaged fluxes in Eq. (30) are

192

DAI AND WOODWARD

approximately obtained 3 H 1 qx W − qxNW , 2 2 3 H 1 N ≈ qy S − qy S , 2 2

q¯ xHW ≈ q¯ yHS

3 H 1 qx E − qxNE , 2 2 3 H 1 N ≈ qy N − qy N . 2 2

q¯ xHE ≈ q¯ yHN

(33) (34)

Here, qxNW , qxNE , q yNS , and q yNN have the similar definitions as Eqs. (28), (29), except for that they are evaluated at t = 1t instead of 1t/2. In order to give specific forms of the fluxes at grid interfaces, we write q as q = −∇ K (T ). We use a center difference to approximately calculate the fluxes at four interfaces, for example, qxHW ≈

¤ 1 £ ¡ H¢ K TL − K (T H ) . 1x

Here, the subscript L refers to the left cell to the current cell; TLH is the cell-averaged value in the left cell at t = 1t/2. Writing K (T ) in a polynomial of T , we may have qxHW = ¢ ¡ Here, κ˜ TLH , T H is defined as κ˜

¡

¢ 1 ¡ H H ¢¡ H TL − T H . κ˜ TL , T 1x

TLH , T H

¢

(35)

¡ ¢ K TLH − K (T H ) . ≡ TLH − T H

˜ L , T ) = κ0 + κ3 (TL2 + T 2 )(TL − T )/4. For example, when κ(T ) = κ0 T + κ3 T 3 , then κ(T We should mention that we have not introduced any approximation in Eq. (35) for the nonlinearity. Similarly, we may write the fluxes at other three grid interfaces: ¢ 1 ¡ H H ¢¡ H TB − T H , κ˜ TB , T 1y ¢ 1 ¡ H H ¢¡ H =− κ˜ TR , T TR − T H , 1x ¢ 1 ¡ H H ¢¡ H =− κ˜ TT , T TT − T H . 1x

q yHS = qxHE q yHN

(36)

Here TRH , TBH , and TTH are the temperatures of the cell to the right, bottom and top, respectively, at t = 1t/2. We may similarly evaluate the fluxes qxNW , qxNE , qxNS , and qxNN . Applying these fluxes in Eqs. (27), (30), we obtain a set of nonlinear difference equations E TN = E T 0 + β H T H + Q H ,

(37)

3 1 E TH = E T 0 + [β H T H + Q H ] − [β N T N + Q N ]. 4 4

(38)

Here, β H and Q H are defined as ¤ ¤ 1t £ H 1t £ H κ˜ L + κ˜ RH + κ˜ B + κ˜ TH , 2 2 (1x) (1y) ¤ ¤ 1t £ H H 1t £ H H ≡ κ˜ T + κ˜ BH TBH + κ˜ T + κ˜ TH TTH , (1x)2 L L (1y)2 B B

βH ≡ QH

RADIATION HYDRODYNAMICS

193

and β N and Q N are exactly the same as the equations above except for the superscript H which should be replaced by N for β N and Q N . We should mention that the extra half time step and the interpolation based on the values at 1t/2 and 1t were first introduced in [23] for gas dynamics and are iteratively implemented in [30]. Compared to two typical implicit schemes, the backward Euler and the Crank–Nicolson schemes, the number of unknown variables has been doubled in Eqs. (37), (38). But, Eqs. (37), (38) have advantages we need, which is second order accurate and give the exact solution when the size of time steps is very large. Actually, three (or more) time-level implicit schemes for heat transfer have been developed in order to reach a high order accuracy. For example, the three time-level implicit Dupont scheme [9, 16] may produce excellent results for nonlinear problems. The threelevel scheme of Lees [5] doesn’t require iteration over a time step to handle nonlinearities. Since the flux calculation in the previous three (or more) time-level implicit schemes depends on initial data, numerical errors do not undergo damping when the size of time steps is large. The set of nonlinear algebraic Eqs. (37), (38) is what we want to solve. The linearized version of Eqs. (37), (38) (or the version for linear heat transfer) is a 2N × 2N matrix with ten non-zero bands. Here N is the number of the grid points in each direction. The nonlinearity in the set of Eqs. (37), (38) is another headache. One of the typical approach for nonlinearities is to use the Newton iteration. But the calculation for Jacobi coefficients is extremely time consuming if the number of unknown variables is large. We want to iteratively solve the nonlinear Eqs. (37), (38). A straightforward procedure is to evaluate the right hand sides (RHSs) of Eqs. (37), (38) using an initial guess for T H and T N . Then, from Eqs. (37), (38), we obtain improved E T at t = 1t/2 and t = 1t. From the improved E TH and E TN , we may obtain improved T H and T N . Unfortunately, this iterative procedure 2 is larger than unity, because through each iteration does not converge when κ1t/(1x) ˜ N numerical errors in E T are increased by a factor larger than unity. Our approach is as follows. We write Eqs. (37), (38) in the form δ E TN + β˜ H δ E TH = Q H − β˜ H E T 0 , ¶ µ 3 1 1 1 − β˜ N δ E TN + 1 + β˜ H δ E H = (3Q H − Q N ) + (β˜ N − 3β˜ H )E T 0 . 4 4 4 4

(39) (40)

Here δ E T ≡ E T − E T 0 , β˜ ≡ β/σ, and σ ≡ E/T . Equations (39), (40) may be further written in the form ¡ ¢ ¡ ¢ (1 + (3/4)β˜ H ) Q H + E˜ 0N − β˜ H (3/4)Q H − (1/4)Q N + E˜ 0H N (41) ET = ET 0 + 1 + (3/4)β˜ H + (1/4)β˜ H β˜ N ¡ ¢ (3/4)Q H − (1/4)Q N + E˜ 0H + (1/4)β˜ N Q H + E˜ 0N . (42) E TH = E T 0 + 1 + (3/4)β˜ H + (1/4)β˜ H β˜ N Here E˜ 0N ≡ −β˜ H E T 0 and E˜ 0H ≡ (β˜ N − 3β˜ H )E T 0 /4. Initially we guess the cell-averaged values of temperatures T N and T H , and evaluate the RHSs of Eqs. (41), (42). The improved T N and T H are obtained through E N and E H which are obtained through Eqs. (41), (42). If the improved solution doesn’t satisfied the accuracy requirement, we may consider the improved temperatures as the initial guess to continue the iteration. This primitive iterative approach is also called the Gauss–Seidel approach. Numerical experiments show that this iterative procedure converges even for nonlinear problems.

194

DAI AND WOODWARD

In the Gauss–Seidel approach, information is carried over only one grid cell through each iteration, and therefore the convergence is very slow for a large time step. In order to speed up the convergence, we divide all grid cells into two sets which are staggered with each other and called red and black sets. If we implement Eqs. (41), (42) for the red set first, then for the black set, in each iteration, the number of iterations may be reduced to half of that needed in the Gauss–Seidel approach for a given required accuracy, because information is carried over two cells through each iteration. From numerical heat transfer with a constant thermal diffusivity, it is known that numerical errors with high frequencies are efficiently killed in the first few iterations, and errors with low frequencies remain even after many iterations. These phenomena indicate that we may use a coarse grid to kill low frequency errors, and use a fine grid to kill high frequency errors. A typical algorithm of this kind of multigrid method starts from a certain level of grid. After a few iterations, if the error is tolerated, then the grid points in each direction are doubled, and another set of a few iterations are implemented in the finer grid with the initial guess obtained from the solution on the coarse grid. If the error is not tolerated, the grid points in each direction are reduced by a factor of two, and another set of a few iterations are implemented to kill the lower frequency errors. After the lower frequency errors are killed, then we go back to the original grid to kill the newly introduced high frequency errors. The multigrid method has been developed for many years. Fedorenko [3] and Backvalov [4] formulated multigrid algorithms for the standard five point finite difference discretization for the Poisson equation and the general linear elliptic partial differential equations. The paper by Brandt [8] is one of the earliest in which practical results were reported. At first there was much debate and skepticism about the true merits of the multigrid method. This led researchers to the development of more transparent convergence proofs (see, e.g., [10, 21] for a survey of theoretical development). Although the rate of convergence proofs of the multigrid method is complicated, their structure has now become more or less standard and transparent. Since the coarse grid is more efficient for a smoother solution, we may work on residue of the solution instead of the solution itself on the coarse grid. Suppose that C N (or C H ) and R N (or R H ) are the error and residue on a fine grid, i.e., C N ≡ T N − TeN ,

C H ≡ T H − TeH ,

R N ≡ E TN + β H T H − Q H − E T 0 ,

(43)

1 3 3 1 R H ≡ − β N T N + E TH + β H T H − Q H − Q N − E T 0 . 4 4 4 4

(44)

Here TeN and TeH denote the exact solutions. On the coarse grid, we iteratively solve the following two equations for C H and C N : σcN C N + β H C H = Q cH + R N , µ µ ¶ ¶ 3 3 H 1 N 1 Qc − Qc + RH . − β N C N + σcH + β H C H = 4 4 4 4

(45) (46)

If the thermal diffusivity κ(T ) is a constant, then σcN , σcH , βcH , βcN , Q cH , and Q cN in Eqs. (45), (46) are exactly the same as those in Eqs. (39), (40) except for the multipliers TLH , TRH , TBH , and TTH in Eqs. (39), (40), which are replaced by C LH , C RH , C BH , and C TH in Eqs. (45), (46). A few iterations of Eqs. (45), (46) on the coarse grid will give very

RADIATION HYDRODYNAMICS

195

FIG. 4. The temperature at t = 0 (the top) and at t = 1t (the bottom).

accurate solutions for C N and C H since C N and C H are smooth and the vanishing C N (or C H ) is a reasonable initial guess. After a few iterations for C N and C H on the coarse grid, we may correct the solution on the fine grid: T N = T N − CN,

(47)

T H = T H − CH.

(48)

196

DAI AND WOODWARD

This procedure is often called a “coarse grid correction.” We should mention that the initial guess obtained from the coarse grid correction often contains large high frequency errors in the fine grid, for which a few iterations on the fine grid are needed. Our problem is nonlinear, σcN , σcH , βcH , and βcN all depend on T N and TH , and Q N and Q H nonlinearly depend on T N and TH . Therefore, strictly speaking, Eqs. (45), (46) are no longer true for nonlinear problems. But, we should realize that a coarse grid is used only for killing low frequency errors, but not for the solution itself. After the coarse grid correction, we still have to go back to the fine grid for the solution. Therefore, equations we use in the coarse grid for the coarse grid correction do not necessarily have to be exact, as long as the coarse grid correction may speed up the convergence. Therefore, even for the nonlinear problem, we still use Eqs. (45), (46) for the coarse grid correction, in which σcN , σcH , βcH , and βcN are all evaluated at the solution obtained on the fine grid. Numerical experiments show that the coarse grid correction based on Eqs. (45), (46) works for nonlinear problems. We would like to demonstrate the convergence of the three kinds of iterative approach, Gauss–Seidel, red-black, and multigrid, in a calculation. The calculation is performed in a (0.7 − 0) × (0.7 − 0) domain with 256 × 256 grid cells. κ(T ) = 10−3 (1 + 10T 3 ). Initially, there is a circular band in which (ρ, T ) = (1, 1), and (ρ, T ) in the background is (0.25, 0.001), 1t = 0.2. The maximum value of the parameter 2κ1t/1x 2 , which should be less than unity in explicit schemes, is about 1.3 × 104 . Figure 4 shows the temperature at the initial time and after one time step. We run the problem using the three kinds of iterative approach mentioned above. The maximum value of the residues R H and R N versus either the number of iterations or the CPU time used is shown in Fig. 5, where broken lines,

FIG. 5. The convergence obtained from three kinds of iterative approach: Gauss–Seidel (broken lines), redblack (dashed lines), and multigrid (solid lines). The maximum value of the parameter 2κ1t/1x 2 is about 1.3 × 104 .

RADIATION HYDRODYNAMICS

197

dashed lines, and solid lines are respectively obtained from the Gauss–Seidel, red-black, and multigrid methods. Now we would like to discuss the conservative feature of the scheme for the heat transfer. After finding a sufficiently accurate solution for T H through the iterations described above, we substitute the solution into Eqs. (35), (36) to find the flux across each interface, qxHW and q yHS . Finally we follow the conservation law, Eq. (27), to update the energy density E T . 5. NUMERICAL EXAMPLES

The numerical schemes described in this paper have been tested for some one- and twodimensional problems for its correctness and robustness, a few of which are presented here to illustrate the features of the schemes. Uniform grids and 5/3 for γ are used in all the examples. Courant numbers which are defined as cs 1t/1x in the Godunov scheme are around 0.6. The linear Riemann solver is used in all these examples unless otherwise stated. In the first one-dimensional example, a wave is initially set up through nonvanishing d R+ , which is shown through the dashed lines in Fig. 6, a R = 1. A grid with 200 cells is used in the simulation domain 0 < x < 1. The solid lines in Fig. 6 show the simulation

FIG. 6. The propagation of a nonlinear wave. The dashed lines are the initial condition, and the solid lines are the profiles at t = 0.5, 1.0, 1.5, 2.0, and 2.5.

198

DAI AND WOODWARD

FIG. 7. A shock-tube problem. The left and right states for (ρ, T, u x ) are (1, 0.5, 50) and (2, 1, −40), which are shown by the dashed lines. The solid lines are the solution at t = 0.04, which is obtained using the linear Riemann solver. Two shocks with Mach numbers about 82 and 39 and a contact discontinuity are generated from the initial discontinuity.

results at t = 0.5, 1.0, 1.5, 2.0, and 2.5. Pr and cs in Fig. 6 denote the radiative pressure a R T 4 /3 and the wave speed. The smooth wave turns into a shock when propagating. The second one-dimensional problem is a shock-tube problem, in which the left and right states for (ρ, T, u x ) are (1, 0.5, 50) and (2, 1, −40), respectively, a R = 1. A grid with 200 cells is used in the domain 0 < x < 1. The dashed lines in Fig. 7 are the initial condition and the solid lines in Fig. 7 are the solution at t = 0.04. Two shocks with Mach numbers about 82 and 39 and a contact discontinuity are generated from the initial discontinuity. The exact solutions for (ρ, T, u x ) in the two middle regions R2 and R3 may be found through the nonlinear Riemann solver described in this paper, and they are (6.95456, 9.90081, −2.73959) and (13.8008, 9.86595, −2.73959), respectively. The numerical solution shown in Fig. 7 is in a good agreement with the exact solution. The scheme with the linear Riemann solver does not necessarily work for very strong shocks. For example, if we give the left and right states (1, 0.5, 150) and (2, 1, −100), respectively, for (ρ, T, u x ), then the scheme with the linear Riemann solver failed. But, the scheme with our nonlinear Riemann solver works well. The solid lines in Fig. 8 show the

RADIATION HYDRODYNAMICS

199

FIG. 8. A shock-tube problem. The left and right states for (ρ, T, u x ) are (1, 0.5, 150) and (2, 1, −100), which are shown by the dashed lines. The solid lines are the solution at one instant, which is obtained using the nonlinear Riemann solver. Two shocks with their Mach numbers about 227.4 and 107.7, and a contact discontinuity are involved in the problem.

solution obtained from our scheme with two iterations in the nonlinear Riemann solver. The exact solutions for (ρ, T, u x ) in the two middle regions R2 and R3 are (6.99036, 16.541, 3.54391) and (13.959, 16.5283, 3.54391). The Riemann problem involves two shocks with their Mach numbers about 227.4 and 107.7, and a contact discontinuity. Another one-dimensional example is a shock-tube problem involving rarefaction waves; a R = 1. The initial condition for (ρ, T, u x ) is (1, 1, −1) and (1, 1, 1) in the left and right states, respectively, as shown by the dashed lines in Fig. 9. The solid lines in Fig. 9 are our numerical solution at t = 0.2. We should point out that there is a “starting error” in this shock-tube problem, and the starting error results from the purely discontinuous initial condition. In order to show the correctness of our scheme with the nonlinear Riemann solver, which is based on shock jump conditions, we also simulated the problem using the scheme with the nonlinear Riemann solver. The result is also plotted in Fig. 9, which coincides with the solution obtained with the linear Riemann solver. The first two-dimensional example is the propagation of a one-dimensional wave in a two-dimensional domain. The initial wave is purely one-dimensional in a ξ -direction, which

200

DAI AND WOODWARD

FIG. 9. A shock-tube problem involving rarefaction waves. The initial condition for (ρ, T, u x ) is (1, 1, −1) and (1, 1, 1) in the left and right states, respectively, as shown by the dashed lines. The solutions at t = 0.2 obtained through nonlinear and linear Riemann solvers are plotted (solid lines) against each other. No difference may be observed between the two solutions.

makes an angle α (α = 30◦ ) with the x-axis. The simulation is performed in the domain (L x − 0) × (L y − 0), where L x = 1/ cos α and L y = 1/ sin α, so that periodic boundary conditions may be used. A grid with 128 × 128 cells is used in the simulation domain. The dashed lines in Fig. 10 are initial profiles along the line y = L y /2, and the solid lines are the profiles at t = 0.5, 1.0, 1.5, 2.0, and 2.5. The next two two-dimensional examples are shock-tube problems in a two-dimensional domain. The initial flows are purely one dimensional in the ξ -direction, which makes an angle α = 30◦ with the x-axis. The simulation domain is (L x −0)×(L y −0) with 128×128 grid cells. Here L x = L y = 1, and a R = 10−3 . Figure 11 shows the propagation of a Mach 61 shock. The dashed lines in Fig. 11 are the initial condition and the solid lines are the profiles at t = 0.1 along the line y = L y /2. Figure 12 shows a shock-tube problem involving two shocks. The dashed lines in Fig. 12 are the initial condition and the solid lines are the profiles at t = 0.1 along the line y = L y /2. Like other Godunov schemes for the Euler equations, there are starting errors in these two problems.

RADIATION HYDRODYNAMICS

201

FIG. 10. The propagation of a nonlinear wave at the ξ -direction in a two-dimensional domain (L x − 0) × (L y − 0). The ξ -direction is at an angle 30◦ with the x-axis. The dashed lines are the initial profiles along the line y = L y /2, and the solid lines are the results at t = 0.5, 1.0, 1.5, 2.0, and 2.5.

The next two-dimensional problem is the interaction between a shock and a denser cloud. Initially, a Mach 256.7 shock propagates in the x-direction. (ρ, T, u x , u y ) in the pre- and post-shock states are (1, 0.01, −22.9472, 0) and (6.57615, 20, 0, 0). In front of the shock, there is a circular cloud which is 100 times denser than the pre-shock state; a R = 0.01. The radius of cloud R = 0.15, and the distance between the shock front and the center of the cloud is 0.18. The simulation domain is 2 × 1 with 512 × 256 grid cells. Figure 13 shows the mass density at four instants. We have also simulated the problem with different grid resolutions. Figure 14 shows the mass densities at the same instant obtained from three simulations with different resolutions. In all the examples above, we have set the heat diffusivity κ(T ) equal to zero. Our final example is the interaction between a wind and denser cloud with a heat diffusivity κ(T ) = 10−3 (1 + 10T 3 ); a R = 1.0. The simulation is performed in a 2 × 1 domain with 512 × 256 grid cell. Initially, there is a circular cloud of radius R = 0.15 whose center is located at (0.3, 0.5). The cloud is 25 times denser than the ambient gas. The state for (ρ, T, u x , u y ) of the ambient gas is (1, 0.09, 0, 0). The temperature of the cloud is such

202

DAI AND WOODWARD

FIG. 11. The propagation of a Mach 61 shock propagating at the ξ -direction in a two-dimensional domain 1 × 1. The ξ -direction is at an angle 30◦ with the x-axis. The dashed lines are the profiles along the line y = 0.5, and the solid lines are the results at t = 0.1.

FIG. 12. A shock-tube problem along the ξ -direction in a two-dimensional domain 1×1. The initial condition is purely one-dimensional. The initial condition along the y = 0.5 is shown by the dashed lines and the solid lines are the profiles at t = 0.1.

RADIATION HYDRODYNAMICS

203

FIG. 13. The mass density during the interaction between a Mach 256 shock and a circular denser cloud.

204

DAI AND WOODWARD

FIG. 14. The solutions obtained through three different grid resolutions. From the top to bottom, 128 × 64, 256 × 128, and 512 × 256 grid cells are used.

RADIATION HYDRODYNAMICS

205

FIG. 15. The temperature at t = 1.2 obtained from three different grid resolutions: (a) 128×64, (b) 256×128, and (c) 512 × 256.

206

DAI AND WOODWARD

FIG. 16. The temperature at t = 1.2 for the interaction between a wind and a denser cloud when the heat diffusivity is set to zero.

that the cloud and ambient gas are mechanically in a equilibrium state. The wind is introduced through the left boundary, at which we always assign the state (1, 0.09, u 0 (t), 0) for (ρ, T, u x , u y ). Here u 0 (t) = 6(1 − e−10t ). Figure 15 shows the temperature at t = 1.2 obtained through three different grid resolutions: 128 × 64, 256 × 128, and 512 × 256. In order to compare the results between with and without the thermal diffusivity κ(T ), we have also simulated the problem without κ(T ) (i.e., κ(T ) = 0). Figure 16 shows the mass densities at the same instant as that in Fig. 15. 6. CONCLUSIONS

In this paper, we have proposed a numerical method for multidimensional radiation hydrodynamical equations in the diffusion limit, which include a second order accurate Godunov scheme and an implicit scheme for radiative heat transfer. The time-averaged fluxes in the Godunov scheme are calculated from approximate Riemann solvers. Both linear and nonlinear Riemann solvers have been proposed for the hyperbolic system of conservation laws. An iterative implicit scheme has been developed for the nonlinear heat transfer. The set of nonlinear algebraic equations arising from the implicit scheme is iteratively solved through a multigrid method. From numerical examples, it is shown that the proposed method keeps the principle advantages of Godunov schemes and the multigrid method. ACKNOWLEDGMENTS The work presented here has been supported by the Department of Energy through Grants DE-FG02-87ER25035 and DE-FG02-94ER25207, by the National Science Foundation through Grant ASC-9309829, by NASA through Grant USRA/5555-23/NASA, and by the University of Minnesota through its Minnesota Supercomputer Institute.

REFERENCES 1. J. Crank and P. Nicolson, A practical method for numerical evaluation of solution of partial differential equations of the heat conduction type, Proc. Camb. Philos. Soc. 43, 50 (1947). 2. S. K. Godunov, Difference methods for the numerical calculation of the equations of fluid dynamics, Math. Sb. 47, 271 (1959).

RADIATION HYDRODYNAMICS

207

3. R. P. Fedorenko, The speed of convergence of one iterative process, USSR Comput. Math. and Math. Phys. 4, 227 (1964). 4. N. S. Bachvalov, On the convergence of a relaxation method with natural constraints on the elliptic operator, USSR Comput. Math. and Math. Phys. 6, 101 (1966). 5. M. Lees, A linear three-level difference scheme for quasi-linear parabolic equations, Math. Comput. 20, 516 (1966). 6. R. Courant and K. O. Friedrichs, Supersonic Flow and Shock Wave (Interscience, New York, 1967), 5th ed. 7. G. Strang, On construction and comparison of differential schemes, SIAM J. Numer. Anal. 5, 506 (1968). 8. A. Brandt, Multi-level adaptive solution technique (MLAT) for fast numerical solution to boundary value problems, in Lecture Notes in Phys. (Springer-Verlag, Berlin, 1973), Vol. 18, p. 82. 9. T. Dupont, G. Fairweather, and J. P. Johnson, Three-level Galerkin methods for parabolic equations, SIAM J. Numer. Anal. 11, 392 (1974). 10. R. A. Nicolaides, On multiple grid and related techniques for solving discrete elliptic systems, J. Comput. Phys. 19, 418 (1975). 11. B. Van Leer, Toward the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, J. Comput. Phys. 32, 101 (1979). 12. S. V. Patankar, Numerical Heat Transfer and Fluid Flow (Hemisphere, New York, 1980). 13. P. L. Roe, Approximate Riemann solvers, parameters vectors, and difference schemes, J. Comput. Phys. 43, 357 (1981). 14. J. Glimm, E. Isaacson, D. Marchesin, and O. McBryan, Front tracking for hyperbolic system, Adv. in Appl. Math. 2, 91 (1981). 15. P. R. Woodward and P. Colella, High resolution difference schemes for compressible gas dynamics, in Lecture Notes in Phys. (Springer-Verlag, New York/Berlin, 1981), Vol. 141, p. 434. 16. M. A. Hogge, A comparison of two- and three-level integration schemes for non-linear heat conduction, in Numerical Methods in Heat Transfer, edited by R. W. Lewis, K. Morgan, and O. C. Zienkiewicz (Wiley, New York, 1981), pp. 75–90. 17. A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49, 357 (1983). 18. D. Mihalas and B. W. Mihalas, Foundations of Radiation Hydrodynamics (Oxford Univ. Press, New York, 1984). 19. P. R. Woodward and P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54, 115 (1984). 20. P. Colella and P. R. Woodward, The piecewise parabolic method (PPM) for gas-dynamical simulations, J. Comput. Phys. 54, 174 (1984). 21. W. Hackbusch, Mult-grid Methods and Applications (Springer-Verlag, Berlin, 1985). 22. P. R. Woodward, Piecewise-parabolic methods for astrophysical fluid dynamics, in Astrophysical Radiation Hydrodynamics, edited by K.-H. A. Winkler and M. L. Norman (1986), pp. 245–326. 23. B. A. Fryxell, P. R. Woodward, P. Colella, and K.-H. Winkler, An implicit-explicit hybrid method for Lagrangian hydrodynamics, J. Comput. Phys. 63, 283 (1986). 24. Y. Jaluria and K. E. Torrance, Computational Heat Transfer (Spring-Verlag, Berlin, 1986). 25. A. Harten and S. Osher, Uniformly high order accurate non-oscillatory schemes, I, SIAM J. Numer. Anal. 24, 279 (1987). 26. A. Harten, B. Enquist, S. Osher, and S. R. Chakravarthy, Uniformly high order accurate essentially nonoscillatory schemes, III, J. Comput. Phys. 71, 231 (1987). 27. J. B. Bell, P. Colella, and J. A. Trangenstein, Higher order Godunov methods for general systems of hyperbolic conservation laws, J. Comput. Phys. 82, 362 (1989). 28. R. J. LeVeque, Numerical Methods for Conservation Laws (Birkh¨auser, Basel, 1992). 29. W. Dai and P. R. Woodward, A simple Riemann solver and high-order Godunov schemes for hyperbolic systems of conservation laws, J. Comput. Phys. 121, 51 (1995). 30. W. Dai and P. R. Woodward, Iterative implementation of an implicit-explicit hybrid scheme for hydrodynamics, J. Comput. Phys. 124, 217 (1996).