Application of central schemes for solving radiation hydrodynamical models

Computer Physics Communications 184 (2013) 1349–1363

Contents lists available at SciVerse ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

Application of central schemes for solving radiation hydrodynamical models Shamsul Qamar a,b,∗ , Waqas Ashraf a a

Department of Mathematics, COMSATS Institute of Information Technology, Islamabad, Pakistan

b

Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany

article

info

Article history: Received 29 July 2012 Received in revised form 21 November 2012 Accepted 19 December 2012 Available online 16 January 2013 Keywords: Radiation hydrodynamics Central schemes Nonlinear conservation laws Shock solutions Convection–diffusion problems

abstract This paper is concerned with the numerical investigation of radiation hydrodynamical models in one and two space dimensions. The flow equations are the set of nonlinear mixed type partial differential equations. The semi-discrete second order central upwind scheme is applied to solve the models. The proposed numerical scheme uses precise information of the local speeds of propagation and, thus, avoids the excessive numerical diffusion which is normally observed in the staggered central scheme. The second order accuracy of the scheme is achieved by using MUSCL-type reconstruction and Runge–Kutta time stepping method. Several case studies are carried out. For validation, the numerical results of the suggested scheme are qualitatively and quantitatively compared with the staggered central (NT) and kinetic flux-vector splitting (KFVS) schemes. The accuracy, efficiency and simplicity of the method demonstrate its potential to solve radiation hydrodynamical models. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Radiation hydrodynamical (RHD) models have wide range applications in different scientific and engineering disciplines [1–4]. The application areas include high-temperature hydrodynamics, modeling of gaseous stars in astrophysics, accretion disks, radiatively driven outflows, supernovas, laser fusion physics, combustion phenomena, reentry vehicles fusion physics, stellar convection and inertial confinement fusion [5,6]. The branch of fluid mechanics that deals with the study of moving fluid and changes in its state under diverse circumstances (e.g. internal and external forces) is called hydrodynamics. The absorption or emission of radiations through matter produces heating and cooling in a system, respectively. The existence of considerable radiation transport intimate the presence of temperature gradients and energy density, indicating the existence of pressure gradients as well. When sufficient time is available, pressure gradients generate considerable flow of fluid (radiation hydrodynamics) and changes in the density. In zero diffusion limit, the radiation hydrodynamical equations (RHEs) could form a system of hyperbolic conservation laws. However, they are different from the inviscid compressible Euler

∗ Corresponding author at: Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany. Tel.: +92 51 9235952; fax: +92 51 9247006. E-mail addresses: [email protected], [email protected] (S. Qamar), [email protected] (W. Ashraf). 0010-4655/$ – see front matter © 2013 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2012.12.021

equations of gas dynamics. One of the major difficulties of standard numerical methods for solving RHEs is to resolve and determine the exact path of strong shocks. Some pioneering work on the problems of radiation hydrodynamics can be found in [6–10]. Dai and Woodward [3] proposed the Godunov scheme including linear and nonlinear Riemann solvers for the solution of the RHEs. The numerical results showed that their methods have kept the key features of Godunov schemes. However, they were found to be computationally expensive. For astrophysical problems, the method proposed in [11] was found to be more successful one where operator splitting technique was combined with the Crank–Nicolson method. Further work on this algorithm and several advantages and disadvantages of Godunov type schemes for solving RHEs are exploited in [2,10,12–18]. Alhumaizi [19] compared different numerical schemes for solving RHEs and showed that flux-corrected transport (FCT), weighted essentially non-oscillatory (WENO) scheme and monotone upstream scheme for conservation laws (MUSCL) are accurate for various parameters of radiation hydrodynamics. Moreover, the work of Tang and Wu [4] on radiation hydrodynamics is a remarkable addition to the applications of kinetic flux-vector splitting (KFVS) methods. In this article, the same KFVS scheme is used for the validation of our numerical results. Central schemes are the simple ones in the sense that they avoid the eigenstructure of the problem. These schemes have been successfully applied in various areas including fluid mechanics, astrophysics, metrology, semiconductors, shallow flow and multicomponent flows [20,21]. The first-order Lax–Friedrichs

1350

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 1. Problem 1: The central upwind scheme results for different values of K on 50 grid points at t = 0.2.

scheme is the forerunner for such central schemes. The central Nessyahu–Tadmor (NT) scheme [22] offers higher resolution while retaining the simplicity of the Riemann-solver-free approach. The numerical viscosity present in the central scheme is of the order O(∆x2r /∆t ). In the convective regime where ∆t ∼ ∆x, the improved resolution of the NT scheme and its generalizations is achieved by lowering the amount of numerical viscosity with increasing r. At the same time, this family of central schemes suffers from excessive numerical viscosity when a sufficiently small time step is enforced, e.g., due to the presence of degenerate diffusion terms. Kurganov and Tadmor [23,24] improved these schemes by using the correct information of local propagation speeds and obtained the semi-discrete central upwind scheme. Analogously to the staggered non-oscillatory central (NT) scheme, it enjoys the benefits of high resolution, simplicity and robustness. However, the central upwind scheme reduces large amount of numerical dissipation present in the NT central scheme. The KFVS scheme is based on the splitting of macroscopic flux functions of the system of equations of the RHD model [4]. The upwinding bias in numerical flux function can be naturally obtained by considering fluid as a collection of particles. The movements of particles in the forward or backward directions automatically split the fluxes of mass, momentum and energy into forward and backward fluxes across the cell interface. In this scheme, we start with a cell averaged initial data of conservative variables and get back the cell averaged values of the conservative variables in the same cell at the next time step. In the two-dimensional case, the flux splitting is done in a usual dimensionally split manner, that is, formula for

the fluxes can be used along each coordinate direction. In order to get second-order accuracy, the MUSCL-type initial reconstruction and the Runge–Kutta time stepping method are employed. In this manuscript, the central upwind scheme is proposed for solving the system of RHEs. The NT central scheme is also applied for the first time to such models. Moreover, the results of the KFVS scheme are also presented for validation and comparison [4]. Several numerical case studies are considered in one and two space dimensions. It was found that central upwind scheme gives comparable solutions to the KFVS scheme. This article is organized as follows. In Section 2, the onedimensional RHD model is introduced. Afterwards, the onedimensional central upwind scheme is derived for the numerical approximation of these equations. We also present the onedimensional staggered central (NT) scheme very briefly. In Section 3, the RHD model and the corresponding numerical schemes are extended to the two-dimensional case. In Section 4, numerical test problems are presented. Finally, Section 5 gives conclusions and remarks. 2. One-dimensional RHD model In this section, the one-dimensional RHD model is presented. Afterwards, the model equations are expressed in analogous form of the compressible Euler equations. A radiatively opaque gas of same radiative and fluid temperatures is considered and the mean free path of photons is assumed narrower compared to the length of flow. The one-dimensional radiation hydrodynamical model is

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

1351

Fig. 2. Problem 1: L1 -error plots of density for different values of K .

given as [1–4]

∂ρ ∂(ρ u) + = 0, (1) ∂t ∂x  2  ∂(ρ u) ∂ ρ u + p + 13 ar T 4 + = 0, (2) ∂t ∂x    ∂(E + ar T 4 ) ∂ u E + p + 34 ar T 4 + = ∇ · [κ(T )∇ T ], (3) ∂t ∂x where ρ , p, uT , E, ar are the mass density, thermal pressure, flow velocity, temperature, total energy and a radiation constant, respectively. Moreover, p = (γ − 1)ρ e, E = ρ e + 21 ρ u2 , e is the specific internal energy and κ(T ) is the nonlinear heat diffusivity. Note that, κ(T ) includes two distinct physical effects. One is the ordinary (i.e. non-radiative) thermal conductivity and the other one is radiation diffusion which is defined as 4ar c λr T 3 /3. Here, c and λr denote the speed of light and the Rosseland mean free path for photons. The above equation can be re-written as

∂ρ ∂(ρ u) + = 0, (4) ∂t ∂x ∂(ρ u) ∂(ρ u2 + p∗ ) + = 0, (5) ∂t ∂x ∂ E∗ ∂[u(E ∗ + p∗ )] + = ∇ · [κ(T )∇ T ], (6) ∂t ∂x where p∗ = p + 31 ar T 4 and E ∗ = E + ar T 4 . Here, e = T being the specific internal energy. In the above model κ(T ) = 0

corresponds to a restrictive case where photon diffusion and energy exchange driven by differences in temperature between the gas and the radiation field are negligible in comparison to radiation work and advection of radiation. This is known as the dynamic diffusion limit [6]. Furthermore, the above system reduces to the inviscid Euler equations of gas dynamics if κ = 0 and ar = 0. The eigenvalues have a similar structure to the Euler equations, i.e.

λ1 = u − c ∗ ,

λ 2 = u,

λ3 = u + c ∗ ,

(7)

where

 ∗

c =

γ ∗ p∗ , ρ

(8)

is the speed of sound in the radiation hydrodynamics case. Let

µ=

1 ar T 4 3 p

and Γ = γ − 1, then

γ + [4Γ (2 − 3Γ ) + 12γ Γ ]µ + 16Γ µ2 (1 + µ)(1 + 12Γ µ)   γ + 20µ + 16µ2 γ −1  =  . 1 + 12µ (1 + µ) γ −1

γ∗ =

Note that, the sound speed can also be re-written as

(c ∗ )2 = c 2 + 4µ

p Γ (2 − 3Γ ) + 4µΓ

ρ

1 + 12µΓ

,

(9)

1352

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 3. Problem 1: L1 -error plots of density for different values of K .



where, c =

γp . ρ

If ar = 0 (i.e. if µ = 0) then c ∗ = c. To

express the temperature in terms of conserved variables, we can write E = ρ T + 12 ρ u2 + aR T 4 as T 4 + β T + η = 0, where β = then

 y=

q

ρ aR

− + 2

(10)

and η = − a1 (E − 12 ρ u2 ). Let p = −η and q = − R

β2 8

1 2

2

[0, xmax ]. For each i = 1, 2, . . . , N, ∆x is a constant width of each

mesh interval, xi denote the cell centers, and xi± 1 refer to the cell



q2 4

+

p3 27

 13

 −

q 2

 +

q2 4

+

p3 27

 13



  − 2y +

2β −2y + √

.

(11)

2y

(12)

After calculating T , the pressure can be obtained from p = (γ − 1)ρ T . 2.1. One-dimensional central upwind scheme In this section, the semi-discrete central-upwind scheme is derived [24]. The above radiation hydrodynamics model can be viewed as a system of convection–diffusion radiation system of the form Wt + F(W)x = Rx ,

x1/2 = 0,

xN +1/2 = xmax ,

(13)

xi+1/2 = i · ∆x,

for i = 1, 2, . . . , N .

(14)

Moreover, xi = (xi−1/2 + xi+1/2 )/2

 .

2

boundaries. We assign,

The quartic equation (10) has four roots from which the physically acceptable root is T =

where, W = (ρ, ρ u, E ∗ )T , F(W) = (ρ u, ρ u2 + p∗ , u(E ∗ + p∗ ))T and R = (0, 0, κ(T )Tx )T . The radiation effect is incorporated in convective term with eigenvalues given by Eq. (7). Before applying the scheme, it is required to discretize the computational domain. Let N be the number of discretization points and (xi− 1 )i∈{1,...,N +1} are partitions of the given interval

and xmax

(15)

∆x = xi+1/2 − xi−1/2 = . N +1   Let Ωi := xi−1/2 , xi+1/2 for i ≥ 1. In each Ωi , the averaged values of the conservative variable W(t ) are given as  1 W(x, t ) dx. (16) Wi := Wi (t ) = ∆x Ωi   By integrating Eq. (13) over the interval Ωi = xi−1/2 , xi+1/2 , we obtain dWi dt

=−

Si+ 1 (t ) − Si− 1 (t ) 2

2

∆x

+

Ri+ 1 (t ) − Ri− 1 (t ) 2

2

∆x

.

(17)

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

1353

Fig. 4. Problem 1: L1 -error plots of density for different values of K .

The numerical fluxes are defined as Si+ 1 =

F(W+ 1 ) i+ 2

+

F(W− 1 ) i+ 2

2

2

−

Here, MM denotes the min-mod nonlinear limiter

ai + 1  2

2

W+ 1 i+ 2

−

W− 1 i+ 2



,

(18)

and for R = (R1 , R2 , R3 )T , we have

(Rk )i± 1 = 0, 2

(R3 )i± 1 = 2

for k = 1, 2,

(19)

κ(Ti−± 1 ) + κ(Ti+± 1 ) Ti+± 1 − Ti−± 1 2 2 2 2 · . 2 ∆x

(20)

Here, W+ and W− are the point values of the piecewise linear ˜ = (ρ, reconstruction W ˜ ρ˜u, E˜∗ ) for W, namely: W+

i+ 12

= Wi+1 −

∆x 2

Wxi+1 ,

W−

i+ 21

= Wi +

∆x 2

Wxi .

(21)

The numerical derivatives Wxi are at least first-order approximations of Wx (xi , t ) and are computed using a nonlinear limiter that would ensure a non-oscillatory nature of the reconstruction (21). A possible computation of these slopes is given by the family of discrete derivatives parameterized with 1 ≤ θ ≤ 2, for example



Wxi = MM θ ∆Wi+ 1 , 2

where 1 ≤ differencing,

2

2

2

(23)

Further, the local one sided speed at xi+ 1 is given as: 2

     ∂F ∂F − + ai+ 1 (t ) = max ρ (W (t )) , ρ (W (t )) . (24) 2 ∂ W i+ 12 ∂ W i+ 12 To obtain the second order accuracy in time, we use a second order TVD Runge–Kutta scheme to solve Eq. (17). Denoting the righthand side of Eq. (17) as L(W), a second order TVD Runge–Kutta scheme update W through the following two stages W(1) = Wn + ∆tL(Wn ), Wn+1 =

(25a)

 Wn + W(1) + ∆t L(W(1) ) ,

 1 2

(25b)

n

n

n+1

where W is a solution at previous time step t and W is updated solution at next time step t n+1 . Moreover, ∆t represents the time step which is calculated under the following Courant–Friedrichs–Lewy (CFL) condition

(22)

2

θ ≤ 2 is a parameter and ∆ denotes central

∆Wi+ 1 = Wi+1 − Wi . 2

  θ ∆Wi+ 1 + ∆Wi− 1 , θ ∆Wi− 1 ,

 {xi } if xi > 0 ∀i,  min i MM {x1 , x2 , . . .} = max{xi } if xi < 0 ∀i,   i 0 otherwise.

∆t ≤ 0.5 min



∆x ∆ x2 , max(|u| + c ∗ ) max(2κ(T ))



,

(26)

where c ∗ is given by Eq. (8). The same CFL condition is also used for the other considered schemes.

1354

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 5. Problem 1: Results on 400 mesh cells at t = 0.2 for κ(T ) = 0.

2.2. One-dimensional central schemes

as [25,26]

Her, the high-resolution non-oscillatory central scheme of Nessyahu and Tadmor [22] is briefly presented. These predictor– corrector type methods are applied in two steps. In the predictor step, the midpoint values are predicted by using the nonoscillatory piecewise-linear reconstructions of the cell averages. In the second corrector step, staggered averaging, together with the predicted mid-values, are used to obtain the updated cell averaged solution. In summary, the scheme can be presented as Predictor:

n+ 21

Wi

= Wni −

2

Fx (Wni ),

1

(27) 1

(Wni + Wni+1 ) + (Wxi − Wxi+1 ) 8    n+ 12 n+ 12 n+ 21 − Fi + ξ Ri+1 − Ri ,

Corrector: Wn+11 = i+ 2

 1 n+ − ξ Fi+12

ξ

2

(28)

where, ξ = ∆t /∆x. Moreover, ∆1x Fx (Wi ) stands for an approximate numerical derivatives of the flux F(t , x = xi )

∂ F(w(t , x = xi )) + O(∆x). (29) ∆x ∂x The fluxes Fx (Wi ) are computed by the same manner as discussed 1

Fx (Wi ) =

for Wx in Eq. (22).

3. Two-dimensional RHD model In this section, we describe the RHEs in two space dimensions and the corresponding Euler-like form. The model can be written

∂ρ ∂(ρ u) ∂(ρv) + + = 0, (30) ∂t ∂x ∂y  2  ∂(ρ uv) ∂(ρ u) ∂ ρ u + p + 13 ar T 4 + + = 0, (31) ∂t ∂x ∂y   ∂(ρv) ∂(ρ uv) ∂ ρv 2 + p + 13 ar T 4 + + = 0, (32) ∂t ∂x ∂y       ∂ v E + p + 43 ar T 4 ∂(E + ar T 4 ) ∂ u E + p + 43 ar T 4 + + ∂t ∂x ∂y = ∇ · [κ(T )∇ T ], (33) where p = (γ − 1)ρ e, E = ρ e + 12 ρ(u2 + v 2 ) and u and v are fluid velocities in x and y-directions. The above system can be re-written as

∂ρ ∂(ρ u) ∂(ρv) + + = 0, ∂t ∂x ∂y

(34)

∂(ρ u) ∂(ρ u2 + p∗ ) ∂(ρ uv) + + = 0, ∂t ∂x ∂y

(35)

∂(ρv) ∂(ρ uv) ∂(ρv 2 + p∗ ) + + = 0, (36) ∂t ∂x ∂y ∂(E ∗ ) ∂[u(E ∗ + p∗ )] ∂[v(E ∗ + p∗ )] + + = ∇ · [κ(T )∇ T ], (37) ∂t ∂x ∂y where p∗ = p + 13 ar T 4 and E ∗ = E + ar T 4 .

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

1355

Fig. 6. Problem 2: Results on 400 mesh cells at t = 0.04.

3.1. Central upwind scheme

Now construct a piecewise linear interplant

Here, the semi-discrete central upwind scheme is presented in two space dimensions [24]. The above model can be expressed as: Wt + F(W)x + G(W)y = Rx + Ry ,

[x0 , xmax ] × [y0 , ymax ] which is covered by cells Cij ≡ xi− 1 , xi+ 1 2 2   × yj− 1 , yj+ 1 for 1 ≤ i ≤ Nx and 1 ≤ j ≤ Ny . The representative 2 2 coordinates in the cell Cij are denoted by (xi , yj ). Here

yj =

xi =

xi−1/2 + xi+1/2 , 2

yj−1/2 + yj+1/2

(39)

2

and

∆xi = xi+1/2 − xi−1/2 ,

∆yj = yj+1/2 − yj−1/2 .

(40)

The cell averaged values of Wi,j (t ) at any time t are given as Wi,j := Wi,j (t ) =

1

∆xi ∆yj

 Cij

W(x, y, t ) dydx.



Wi,j + (Wx )i,j (x − xi )

i,j

 + (Wy )i,j (y − yj ) χi,j ,

(38)

where, W = (ρ, ρ u, ρv, E ∗ )T , F(W) = (ρ u, ρ u2 + p∗ , ρ uv, u(E ∗ + p∗ ))T , G(W) = (ρ u, ρ uv, ρ u2 + p∗ , v(E ∗ + p∗ ))T , Rx = (0, 0, 0, κ(T )Tx )T and Ry = (0, 0, 0, κ(T )Ty )T . The eigenvalues in two space dimension is straightforward extension of Eq. (7). Let Nx and Ny be the large integers in x and y-directions, respectively. We assume a Cartesian grid with a rectangular  domain

(x1/2 , x1/2 ) = (0, 0),

W(x, y, t ) =

where χi,j is the characteristic function for the corresponding cell (xi− 1 , xi+ 1 ) × (yj− 1 , yj+ 1 ), (Wx )i,j and (Wy )i,j are the approxima2

2

2

2

tions of x and y-derivatives of W at the cell centers (xi , yj ). The generalized MM limiter is used for the computation of these partial derivatives to avoid oscillations

(Wx )ni,j   Wi+1,j − Wi,j Wi+1,j − Wi−1,j Wi,j − Wi−1,j , ,θ , = MM θ ∆x 2∆x ∆x (Wy )i,j   Wi,j+1 − Wi,j Wi,j+1 − Wi,j−1 Wi,j − Wi,j−1 = MM θ , ,θ , (43) ∆y 2∆y ∆y where 1 ≤ θ ≤ 2 and MM is given by Eq. (22). After integrating Eq. (38) over the control volume Cij , the two-dimensional extension of the scheme can be scripted as dWi,j dt

Sx

=−

i+ 12 ,j

+

− Sxi− 1 ,j 2

∆x Rx

(41)

(42)

i+ 21 ,j

− Rxi− 1 ,j 2

∆x

y

S

−

i,j+ 21

− Syi,j− 1 2

∆y y

R

+

i,j+ 21

− Ryi,j− 1 2

∆y

.

(44)

1356

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 7. Problem 2: L1 -error plots at κ(T ) = 0.

Here, ax

Here Sx 1 i+ 2 ,j S

F(W−

i+ 21 ,j

=

) + F(W+ ) i + 1 ,j 2

2 G(W−

y i,j+ 12

a

−

i+ 12 ,j

2



W+ 1 i+ 2 ,j

−

W− 1 i+ 2 ,j



,

) + G(W+ ) i,j+ 1

i,j+ 12

=

as

ax

−

2



i,j+ 21

W+

i,j+ 12

2

− W− i,j+ 1



.

(45)

2

y

Similarly, Rxk = 0 = Rk for k = 1, 2, 3 and

( )

Rx4 i± 1 ,j 2

(Ry4 )i,j± 1 2

∆x

= Wi,j +

W−

= Wi,j +

W+

= Wi,j+1 −

i,j+ 21

i,j+ 21

2 ∆y 2

(Wx )i,j ,

W+

i+ 12 ,j

= Wi+1,j −

(46)

∆x 2

(Wx )i+1,j

(Wy )i,j ,

∆y 2

(Wy )i,j+1 .

y i,j+ 21

y i,j+ 12

are the local speeds which can be calculated

 ∂F (W± , 1 ) ± ∂ W i+ 2 ,j   ∂G = max ρ (W± ) . ± ∂ W i,j+ 12 = max ρ



(48)

For complete derivation of the scheme the reader is referred to [24].

The two-dimensional central scheme was proposed by Jaing and Tadmor [27]. The scheme has again a two-step predictor–corrector form. Starting with the cell averages, Wni,j , we use the first-order predictor step for the evolution of the midpoint val-

The intermediate values are expressed as

i+ 12 ,j

i+ 12 ,j

and a

3.2. Two-dimensional central scheme

κ(Ti−± 1 ,j ) + κ(Ti+± 1 ,j ) Ti+± 1 ,j − Ti−± 1 ,j 2 2 2 2 = · , 2 ∆x κ(Ti−,j± 1 ) + κ(Ti+,j± 1 ) Ti+,j± 1 − Ti−,j± 1 2 2 2 2 = · . 2 ∆y

W−

ax a

2

y

i+ 21 ,j

(47)

n+ 1

ues, Wi,j 2 , followed by the second-order corrector step for com1 putation of the new cell averages Wni,+ j . Like the one-dimensional case, no exact (approximate) Riemann solvers are needed. The nonoscillatory behavior of the scheme is dependent on the reconstructed discrete slopes, Wx , Wy , Fx (W), and Gy (W). At each time step the grid is staggered to avoid the flux calculation at the cell interfaces. The scheme is summarized below. In the predictor step one has to calculate the midpoint values n+ 12

Wi,j

= Wni,j −

ξ 2

Fx (Wni,j ) −

η 2

Gy (Wni,j ),

(49)

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

1357

Fig. 8. Problem 3: Results on 400 mesh cells at t = 0.02.

where, ξ = ∆t /∆x and η = ∆t /∆y. This step is followed by a corrector step to get the updated values at the next time step Wn+11

i+ 2 ,j+ 21

=

1

(Wni,j + Wni+1,j + Wni,j+1 + Wni+1,j+1 )  1  ξ 1 n+ 2 n+ 21 x x Fi+1,j − Fi,j + (Wi,j − Wi+1,j ) − 16 2  1  1 ξ n+ n+ 1 + (Wxi,j+1 − Wxi+1,j+1 ) − Fi+12,j+1 − Fi,j+21 16 2   η 1 n+ 1 n+ 1 + (Wyi,j − Wyi,j+1 ) − Gi,j+21 − Gi,j 2 16 2   1 η n+ 21 n+ 12 y y + (Wi+1,j − Wi+1,j+1 ) − Gi+1,j+1 − Gi+1,j 16 2   ξ n+ 1 n+ 1 + (Rx )i+12,j − (Rx )i,j 2 4

2

+

η 2



n+ 1 Ry i,j+21

( )

n+ 1 Ry i,j 2

−( )



.

(50)

This completes the derivation of numerical schemes. 4. Numerical case studies In this section, six test problems are considered to validate the accuracy and performance of the proposed central schemes. For comparison, the KFVS scheme is also applied to this model [4].

Problem 1. This is a one-dimensional shock-tube problem involving two rarefaction waves moving in the opposite directions [25]. The diaphragm is placed at x = 0.5. The initial data are given as

 (1, 1, −1), x ≤ 0.5, (ρ, T , u) = (1, 1, 1) x ≥ 0.5.

(51)

The computational domain is [0, 1] which is subdivided into 50 cells and the final simulation time is t = 0.2. The nonlinear heat diffusivity is taken as κ(T ) = K (1 + 10T 3 ), where K is a scaling factor. To analyze the performance of numerical schemes for transport and diffusion dominated cases, different values of K are considered from the interval [0, 0.5]. Fig. 1 show the numerical results of central upwind scheme for K = 0, 10−3 , 10−2 , 10−1 , 0.5. One can observe that for K ≥ 10−2 the effect of diffusion term becomes visible. Tables 1–4 show the L1 -errors between the reference and numerical solutions of the central and KFVS schemes for different values of the scaling factor K . The reference solution was obtained at 2000 mesh points. The error plots of density, pressure and velocity for different values of the scaling factor are given in Figs. 2–4. Moreover, Fig. 5 shows the comparison of results for K = 0. It can be observed that the results of central-upwind scheme are in good agreement with KFVS and NT central schemes. All figures and tables show that the central upwind and the KFVS schemes give comparable results for the whole range of K while the staggered central (NT) scheme produces large errors in the solutions. In overall, the KFVS scheme has little edge over the central schemes. From these results it can be concluded that the proposed schemes have uniform behavior over the whole range of scaling factor K . Moreover,

1358

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 9. Problem 4: Results on 400 mesh cells at t = 0.18.

density, ρ

temperature, T 1

y–axis

y–axis

1

0.5

0

0.5

1 x–axis

1.5

0.5

2

0

0.5

Pressure, p

1.5

2

1.5

2

velocity, u 1

y–axis

1

y–axis

1 x–axis

0.5

0

0.5

1 x–axis

1.5

2

0.5

0

0.5

1 x–axis

Fig. 10. Problem 5: 2D results on 256 × 128 mesh cells at t = 0.6, κ(T ) = 0.

minor changes can be seen in the magnitudes of errors over the whole range of scaling factor K . Thus, the considered numerical

schemes behave uniformly in both transport and diffusion dominated limits.

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363 density, ρ

temperature, T 1

y–axis

y–axis

1

0.5

0

0.5

1 x–axis

1.5

0.5

0

2

0.5

1 x–axis

Pressure, p

1.5

2

1.5

2

velocity, u

1

1

y–axis

y–axis

1359

0.5

0

0.5

1 x–axis

1.5

0.5

0

2

0.5

1 x–axis

Fig. 11. Problem 5: 2D results on 256 × 128 mesh cells at t = 0.6, κ(T ) = 10−3 (1 + 10T 3 ). Table 1 Problem 1: Comparison of numerical errors at different grid points for K = 10−3 . Methods

Central upwind KFVS NT central

N = 80

N = 160

N = 320

ρ

T

p

ρ

T

p

ρ

T

p

0.013 0.012 0.020

0.007 0.007 0.013

0.014 0.013 0.027

0.005 0.006 0.014

0.003 0.003 0.008

0.006 0.006 0.015

0.002 0.002 0.009

0.001 0.001 0.005

0.002 0.002 0.009

Table 2 Problem 1: Comparison of numerical errors at different grid points for K = 10−2 . Methods

Central upwind KFVS NT central

N = 80

N = 160

N = 320

ρ

T

p

ρ

T

p

ρ

T

p

0.007 0.007 0.009

0.005 0.005 0.009

0.009 0.009 0.014

0.003 0.003 0.004

0.002 0.002 0.004

0.004 0.004 0.006

0.001 0.001 0.002

0.001 0.001 0.002

0.002 0.002 0.004

Table 3 Problem 1: Comparison of numerical errors at different grid points for K = 10−1 . Methods

Central upwind KFVS NT central

N = 80

N = 160

N = 320

ρ

T

p

ρ

T

p

ρ

T

p

0.015 0.014 0.030

0.003 0.003 0.005

0.011 0.010 0.019

0.007 0.006 0.016

0.001 0.001 0.002

0.005 0.005 0.010

0.003 0.003 0.012

0.001 0.001 0.001

0.002 0.002 0.007

Table 4 Problem 1: Comparison of numerical errors at different grid points for K = 0.5. Methods

Central upwind KFVS NT central

N = 80

N = 160

N = 320

ρ

T

p

ρ

T

p

ρ

T

p

0.014 0.016 0.050

0.002 0.002 0.003

0.012 0.011 0.030

0.008 0.008 0.036

0.001 0.001 0.002

0.005 0.005 0.021

0.003 0.003 0.021

0.000 0.000 0.001

0.002 0.002 0.012

Table 5 Problem 2: Comparison of numerical errors for K = 0 at different grid points and CPU times. Methods

Central upwind KFVS NT central

N = 200

N = 800

N = 1600

CPU (s)

ρ

T

p

ρ

T

p

ρ

T

p

N = 200

0.21 0.27 0.40

0.13 0.12 0.23

23.16 23.50 43.90

0.06 0.09 0.14

0.02 0.02 0.05

4.16 4.24 8.84

0.02 0.05 0.07

0.01 0.01 0.02

1.32 1.49 3.66

1.7 3.6 2.3

1360

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 12. Problem 5: 1D results along x = 0.25 at t = 0.6, κ(T ) = 0.

Problem 2. This one dimensional shock-tube problem was considered in [25]. The initial data are given as

 (1, 0.5, 50), x ≤ 0.6, (ρ, T , u) = (2, 1, −40) x ≥ 0.6.

(52)

The computational domain is [0, 1] which is partitioned into 400 grid cells. The simulation results at t = 0.04 for the density, velocity, pressure and temperature are shown in Fig. 6. Two shocks of Mach 82 and 39 and a contact discontinuity are produced from the initial discontinuities. The numerical results of central upwind and KFVS scheme are better than the staggered central (NT) scheme. The central upwind scheme seems to be superior in this case. To

Fig. 13. Problem 5: 1D results along y = 0.5 at t = 0.6, κ(T ) = 0.

quantitatively analyze the accuracy of proposed scheme, we have calculated the L1 -errors at different grid points given in Table 5. The reference solution was obtained at 2000 grid points. It is clear from the results that the central-upwind scheme produces less errors in the solution as compared to the KFVS and NT central schemes. Moreover, Fig. 7 shows the plots of L1 -errors which justify the results of Table 5. Problem 3. This problem is given by Jiang and Sun [25] that involves a strong shock in one-dimensional Riemann problem. Initially the values of (ρ, T , u) for x0 ≤ 0.5 are (1, 0.5, 150) and for

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363 density, ρ at t=0.25

temperature, T at t=0.25 1

y–axis

y–axis

1

0.5

0

0.5

1 x–axis

1.5

0.5

2

0

0.5

density, ρ at t=0.5

1.5

2

1.5

2

1.5

2

1

y–axis

y–axis

1 x–axis temperature, T at 0.5

1

0.5

0

0.5

1 x–axis

1.5

0

1 x–axis

1 x–axis

temperature, T at t=0.75

y–axis 0.5

0.5

1

0.5

0

0.5

2

density, ρ at t=0.75

1

y–axis

1361

1.5

2

0.5

0

0.5

1 x–axis

Fig. 14. Problem 5: Density and temperature contour at different time steps, κ(T ) = 0.

Fig. 15. Problem 5: 1D results along y = 0.55 at t = 0.5, κ(T ) = 0.

x0 > 0.5 are (2, 1, −100). The number of cells are 400 and the computational domain is [0, 1]. The simulation results at t = 0.018 for the density, velocity, pressure and temperature are shown in Fig. 8. The solution consists of two shocks with Mach numbers of about 227.4 and 107.7 and a constant discontinuity. Once again the results of KFVS and central upwind scheme are comparable and the NT scheme gives diffusive results. However, the central upwind scheme is superior in all three schemes.

(5, 1.5, −4). The number of cells are 400 and the domain is [0, 1]. The simulation results at t = 0.18 for the density, velocity, pressure and temperature are shown in Fig. 9. The solution consists of two shocks with Mach numbers of about 227.4 and 107.7 and a constant discontinuity. Once again the results of KFVS and central upwind schemes are comparable and NT scheme gives diffusive results. However, the central upwind scheme is superior in all three schemes.

Problem 4. We consider another test problem where the initial data (ρ, T , u) for x0 ≤ 0.5 are (5, 1.5, 4) and for x0 > 0.5 are

Problem 5. This is a two-dimensional problem describing the interaction of a wind and a denser circular cloud. In the rectangular

1362

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

Fig. 16. Problem 5: 1D results along x = 0.30 at t = 0.5, κ(T ) = 0.

Fig. 17. Problem 6: 2D results on 128 × 128 mesh cells at t = 0.5, κ(T ) = 0.

domain [0, 2] × [0, 1], there is a 25 times denser cylindrical bubble at (0.3, 0.5) whose radius is r = 0.15 and the number of grids are 256 × 128. The state of ambient gas for (ρ, T , u, v) is (1, 0.09, 0, 0) and the wind state is (1, 0.09, 6(1 − e−10t ), 0) which is introduced through the left boundary and ar = 1. The outflow boundary conditions are applied at the right, lower and upper boundaries of the domain. We computed the solution without and with diffusion limit, where, κ(T ) = 10−3 (1 + 10T 3 ) [3]. It is observed that the incoming shock produces shock in the bubble and a re-

fracted shock. The contours for ρ, T and u at t = 0.6 are shown in Figs. 10–14. The comparisons of results in Figs. 15 and 16 show good agreements between the proposed schemes and the results available in the literature. It can be observed that KFVS and central upwind scheme produce comparable solutions and the NT scheme is diffusive. However, The KFVS scheme produces a more resolved solution (see Figs. 15 and 16). Problem 6. Here, a square box of length 2.0 is considered. Further, a circular bubble of radius 0.15 is placed at the center of the

S. Qamar, W. Ashraf / Computer Physics Communications 184 (2013) 1349–1363

1363

central and KFVS schemes. Once again, a good performance of the central-upwind scheme can be seen. However, KFVS scheme shows better performance in all three schemes. 5. Conclusions We focused on the numerical solution of one- and twodimensional radiation hydrodynamical equations (RHEs). Two different types of central finite volume schemes were applied to solve these equations, the central-upwind and the staggered central schemes. The proposed numerical schemes preserves monotonicity due to using MUSCL-type reconstruction and also avoid the detailed knowledge of a complicated exact/approximate Riemann solver. For validation, the KFVS scheme was also applied to solve this model. A number of case studies were carried out and the accuracies of the schemes were analyzed quantitatively and qualitatively. It was found that central-upwind schemes have better resolved discontinuous profiles as compared to the NT schemes. Further, it was concluded that the proposed numerical scheme gives comparable results with the KFVS scheme at low CPU time. The suggested method gives comparable accuracy to the scheme that uses Riemann solvers, e.g. Sekora and Stone [2], and is computationally less expensive due to its Riemann-solver free algorithm. Acknowledgment A partial support by the Higher Education Commission (HEC) of Pakistan is gratefully acknowledged. References

Fig. 18. Problem 6: 1D results along y = 1.0 at t = 0.5, κ(T ) = 0.

box. The values for (ρ, T , u, v) are (1, 0.9, 0, 0) inside the circle and (25, 0.9, 0, 0) outside the circle. The simulation results are obtained at t = 0.5 on 128 × 128 grid points. Figs. 17 and 18 show the contour plots and comparison among the central-upwind, NT

[1] R.B. Lowrie, J.E. Morel, J. Quantitative Spectroscopy & Radiative Transfer 69 (2001) 475. [2] M.D. Sekora, J.M. Stone, J. Comput. Phys. 229 (2010) 6819. [3] W. Dai, P.R. Woodward, J. Comput. Phys. 142 (1998) 182. [4] H.Z. Tang, H.M. Wu, J. Comput. Fluids 29 (2000) 917. [5] G. Cox, Combustion Fundamentals of Fires, Academic Press, NewYork, 1995. [6] D. Mihalas, B. Weibel-Mihalas, Foundations of Radiation Hydrodynamics, University Press, Oxford, 1984. [7] J.I. Castor, Astrophys. J. 172 (1972) 779. [8] C.D. Levermore, G.C. Pomraning, Astrophs. J. 248 (1981) 321. [9] D. Mihalas, R. Klein, J. Comput. Phys. 46 (1982) 97. [10] G.C. Pomraning, The Equations of Radiation Hydrodynamics, Pergamon Press, Oxford, 1973, pp. 241–282. [11] J.M. Stone, D. Mihalas, M.L. Norman, Astrophys. J. Suppl. S. 3 (1996) 903. [12] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhauser Verlag, 1992. [13] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002. [14] S. Yip, Handbook of Materials Modeling, Vol. 1, Springer, 2005. [15] R.B. Lowrie, J.E. Morel, J. Quant. Spectrosc. Radiat. Transfer. 69 (2001) 475. [16] C. Buet, B. Despres, J. Comput. Phys. 215 (2006) 717. [17] S. Jin, C.D. Levermore, J. Comput. Phys. 126 (1996) 449. [18] G.A. Sod, J. Comput. Phys. 27 (1978) 1. [19] K. Alhumaizi, J. Comput. Chem. Eng. 28 (2004) 1759. [20] J. Balbas, S. Karni, Math. Model. Num 43 (2009) 333. [21] S. Qamar, G. Warnecke, Appl. Numer. Math. 50 (2004) 183. [22] H. Nessyahu, E. Tadmor, J. Comput. Phys. 87 (1990) 408. [23] A. Kurganov, C. Lin, J. Commun. Comput. Phys. 2 (2007) 141. [24] A. Kurganov, E. Tadmor, J. Comput. Phys. 160 (2000) 241. [25] S. Jiang, W. Sun, Int. J. Numer. Methods 53 (2007) 391. [26] W. Wenlong, R. Woodward, J.Comput. Phys. 157 (2000) 99. [27] G.-S. Jaing, E. Tadmor, SIAM J. Sci. Comput. 19 (1998) 1892.