Second-order-accurate discrete ordinates solutions of transient radiative transfer in a scattering slab with variable refractive index

International Communications in Heat and Mass Transfer 38 (2011) 1213–1218

Contents lists available at ScienceDirect

International Communications in Heat and Mass Transfer j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / i c h m t

Second-order-accurate discrete ordinates solutions of transient radiative transfer in a scattering slab with variable refractive index☆ Jia-Ming Wang, Chih-Yang Wu ⁎ Department of Mechanical Engineering, National Cheng Kung University, Tainan, 701, Taiwan, ROC

a r t i c l e

i n f o

Available online 16 July 2011 Keywords: Transient radiative transfer Anisotropical scattering Refractive media

a b s t r a c t The discrete ordinates method (DOM) with a second-order upwind interpolation scheme is applied to solve transient radiative transfer in a graded index slab suddenly exposed to a diffuse strong irradiation at one of its boundaries. The planar medium is absorbing and anisotropically scattering. From the comparison of the results obtained by the first-order DOM, the second-order DOM, the modified DOM and the Monte Carlo method, it can be seen that the numerical diffusion in the transient solutions obtained by the second-order DOM is less than that in the solutions obtained by the first-order DOM, but the numerical diffusion is still noticeable, especially for optically thin and moderate cases. By contrast, for optically thick cases the numerical diffusion due to the finite difference of the advection term of the transient radiative transfer equation is minor. In general, it is still necessary to adopt a DOM with a higher order scheme to capture the wave front of transient radiative transfer accurately. Besides, the influence of numerical diffusion is a little less noticeable for the case with a larger gradient of refractive index, and the distribution of direction-integrated intensity around the irradiation boundary decreases and that around the other boundary increases with the increase of the anisotropically scattering coefficient. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Over the past few decades a considerable number of studies have been made on transient radiative transfer in scattering media [1,2]. Recent studies on transient radiative transfer were reviewed by Kim and co-workers [2]. Only a few of them considered the radiative transfer problems incorporating the spatially continuous variation of refractive index in the medium. In these situations, because of the continuous variation of the refractive index, the speed of light is a function of the position and the radiation streams in curved paths rather than straight lines. The effects of the spatially continuous variation of refractive index on radiative transfer in a scattering medium are important in numerous applications, such as the propagation of laser pulse through a cloud of aerosols [3] and the thermal barrier coatings [4]. Several variations of the transient radiative transfer equation for a scattering medium with a spatially varying refractive index have been proposed by many authors [5–7]. Among those authors, Premaratne and co-workers [7] have shown that their result reduces to the standard radiative transfer equation when the refractive index is constant and to the geometrical optics equation ☆ Communicated by H. Yoshida and S. Maruyama. ⁎ Corresponding author. E-mail address: [email protected] (C.-Y. Wu). 0735-1933/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2011.06.015

when the medium is lossless and non-scattering. Wu [8] and Liu and Hsu [9] used the discrete ordinates method (DOM) and the discontinuous finite element method to solve transient radiative transfer in gradient index media, respectively. Various differential approximation methods for transient radiative transfer in refractive planar media were examined [10]. Recently, a Monte Carlo method (MCM) based on the path length and time of flight in closed form was presented [11] and a modified discrete ordinates method (MDOM) was developed [12] for transient radiative transfer in a refractive slab suddenly exposed to a diffuse irradiation. In this work, we adapt a second-order upwind scheme developed for solving Eulerian advection problems [13] to reduce numerical diffusion in the discrete ordinates solutions of transient radiative transfer in a refractive slab. The planar medium considered is anisotropically scattering and is suddenly exposed to a diffuse irradiation at one of its boundaries. Because the exact analytical path length in a slab with a linear refractive index is available [14], accurate solutions of transient radiative transfer in such a medium can be obtained by the MCM. Thus, for comparison purpose, we consider the case with a linear refractive index in this work. By comparing the present solutions with the results obtained by the MCM, the first-order DOM and the MDOM, we examine the performance of the second-orderaccurate DOM. Because of the steep features of the transient

1214

J.-M. Wang, C.-Y. Wu / International Communications in Heat and Mass Transfer 38 (2011) 1213–1218

0.6

Nomenclature

N S t tk t* wm z

anisotropically scattering coefficient coefficient in Eq. (11) coefficient in Eq. (11) speed of light in vacuum coefficient in Eq. (11) source term in Eq. (11) direction-integrated intensity Heaviside step function radiative intensity radiative intensity at the center of a directional grid intensity of incident diffuse radiation thickness of the slab total number of discrete ordinates refractive index refractive index of the medium at z = 0 refractive index of the medium at z = L gradient of the refractive index linear refractive index, defined as n(z) =n0 + (nL − n0) z/L total number of the slices source function time kth time level dimensionless time quadrature weight associated to μm geometrical coordinate

Greek symbols α coefficient; see Eq. (8) β extinction coefficient γ coefficient of the angular redistribution term; see Eq. (2) Δt size of a time level Δz L/N θ angle between the z-axis and the direction of radiation intensity σj coefficient defined in Eq. (13) μ μ = cos θ μm directional cosine of the discrete directions τL optical thickness Φ dimensionless direction-integrated intensity ω scattering albedo

N = 25 N = 50 N = 100

0.4

N = 200

Φ(z/L, t*)

a1 A B c0 C D G H I Im I0 L M n n0 nL n' n0 → nL

0.2

t* = 1.0

t* = 0.0563

0

0

0.2

t* = 0.125

0.4

0.6

0.8

1

z/L Fig. 1. Comparison of the transient solutions obtained by the second-order DOM using various total numbers of the slices for the case with τL = 0.1, n = 1 → 1.5, ω = 1 and a1 = 0.

radiative transfer, we compare the simulation results in the early stage. 2. Physical model and the second-order-accurate DOM Consider transient radiative transfer in an absorbing and anisotropically scattering slab with a continuously varying refractive index n(z). It is assumed that the slab is gray and the geometrical thickness L of the scattering slab is much greater than one wavelength. The transient radiative transfer equation obtained by simplifying the equation presented in [7] is expressed as h
 i 2 ∂ 1−μ I z; μ; t nðzÞ ∂I ðz; μ; t Þ ∂Iðz; μ; t Þ +μ + γ ð zÞ c0 ∂t ∂z ∂μ + βIðz; μ; t Þ = Sðz; μ; t Þ;

ð

Þ

ð1Þ

0.8 N = 25 N = 50

0.6

N = 100 N = 200

Φ(z/L, t*)

Subscripts 0 boundary at z = 0 e index of discrete ordinates method j jth spatial grid L location z = L m mth discrete ordinate t index of discrete ordinates method b index of discrete ordinates method w index of discrete ordinates method

0.4

0.2 t* = 10.0 t*=1.25

t* = 0.563

0 response to the step function irradiation, the present problem with a step function irradiation provides a stringent test of the schemes employed. Since the numerical diffusion of finitedifference solutions is noticeable in the early stage of transient

0

0.2

0.4

0.6

0.8

1

z/L Fig. 2. Comparison of the transient solutions obtained by the second-order DOM using various total numbers of the slices for the case with τL = 1, n = 1 → 1.5, ω = 1 and a1 = 0.

J.-M. Wang, C.-Y. Wu / International Communications in Heat and Mass Transfer 38 (2011) 1213–1218

1

with α1/2 = 0 and αM + 1/2 = 0, and the angular integration of the source function is replaced by a weighted sum. The finite-difference form of the spatial derivative in Eq. (7) is obtained by integrating Eq. (7) over a slice of width Δz = L/N, where N is the total number of the slices. The slice centers at zj and the edges of the slice are at zj ± 1/2. Moreover, applying the Euler implicit time differencing to the time derivative in Eq. (7), we obtain

N = 25 N = 50

0.8

Φ(z/L, t*)

N = 100 N = 200

0.6

k

t* = 100.0

0

0.2

k k + βIm; j = Sm; j ;

1≤ m ≤ M; 1≤ j ≤ N; k = 1; 2;:::;

0.4

0.6

ð9Þ

0.8

1

where the superscript k denotes the beginning of the time level (t k, t k + 1), cj =c0/n(zj) and

z/L Fig. 3. Comparison of the transient solutions obtained by the second-order DOM using various total numbers of the slices for the case with τL = 10, n = 1 → 1.5, ω = 1 and a1 = 0.

where c0 is the speed of light in vacuum, I the radiative intensity, μ = cos θ with θ denoting the angle between the positive z-axis and the direction of radiation intensity, t the time, γ defined as γ ð zÞ =

k

t* = 12.5

t* = 5.63

0

k

k−1

Im; j + 1=2 −Im; j−1=2 Im; j −Im; j + μm cj Δt Δz

 2 2 k k nj + 1=2 −nj−1=2 αm + 1=2 Im + 1=2; j −αm−1=2 Im−1=2; j + 2n2j Δzwm

0.4

0.2

1215

1 dnðzÞ ; nðzÞ dz

k

Sm; j =


 βω M k ∑ I 1+a1 μm μm′ wm′ 2 m′ = 1 m′ ; j

with wm denoting the quadrature weight associated to μm. Eq. (9) can be rewritten as k−1 wm Δz k ðI −Im; j Þ cj Δt m; j

k

k

k

k

k

k

+ Am ðIm;e −Im;w Þ + Bt; j It; j −Bb; j Ib; j + Cm Im; j = Dm; j ; 1≤ m ≤ M; 1≤ j ≤ N; k = 1; 2;:::; ð11Þ

ð2Þ

β the extinction coefficient and S(z, μ, t) the source function for linearanisotropic scattering expressed as

ð10Þ

where Am = |μm|wm, Bt, j = αm + 1/2|nj2+ 1/2 −nj2− 1/2|/2nj2, Bb, j = αm − 1/2 k k |nj2+1/2 − nj2− 1/2|/2nj2, Cm = βwmΔz, Dm, j =wmΔzSm, j, and the subscripts w, e, t and b denote w= j− sign(μm)(1/2), e= j+ sign(μm)(1/2), t= m+ sign(n′j)(1/2) and b= m− sign(n′j)(1/2), respectively, with n′j denoting

Sðz; μ; t Þ =

1

βω ∫ 2 −1

 0  1 + a1 μμ I z; μ ′ ; t dμ ′

ð3Þ

with ω denoting the scattering albedo and a1 denoting the anisotropically scattering coefficient. Here, the angular redistribution term with the angular derivative of the radiative intensity describes the curvature of the optical trajectory. Consider the surroundings and the medium at the boundaries have the same refractive index. A step function is adopted to describe the sudden irradiation. The boundary and initial conditions can be expressed as Ið0; μ; t Þ = Hðt ÞI0 ; μ > 0;

ð4Þ

IðL; μ; t Þ = 0; μ < 0;

ð5Þ

Iðz; μ; 0Þ = 0; −1 < μ < 1; 0≤ z ≤ L;

ð6Þ

where H denotes the Heaviside step function and I0 the intensity of diffuse irradiation. The discrete-ordinate representation of Eq. (1) can be written as αm + 1=2 Im + 1=2 −αm−1=2 Im−1=2 nðzÞ ∂Im ðz; t Þ ∂I ðz; t Þ + μm m + γ ð zÞ c0 wm ∂t ∂z + βIm ðz; t Þ = Sm ðz; t Þ; 1≤ m ≤ M;

ð7Þ

k

k

It;j = Im;j :

ð8Þ

ð12Þ

Because of the steep feature in the spatial distribution of radiation, we adopt the second-order-accurate scheme due to van Leer [13] to complete the spatial differencing procedure, k Im;e

!
 dI k n2e k Δt 2 m;j ; = 2 Im;j + signðμm Þne 1−σj with σj = cj Δz 2 nj

ð13Þ

where k k

2

ðiÞ dIm;j =

Im;j + 1 n2j + 1

Ik

− nm;j2

k

 k

Im;j n2j

j

Im;j

+ 1

n2j +

1

−

− 

k Im;j−1

k Im;j−1 n2j−1



; if

k Im;j

+ 1

n2j +

1

−

k Im;j

n2j

!

k Im;j

n2j

−

k Im;j−1

n2j−1

!

> 0; ð14Þ

n2j−1

k

ð15Þ

ðiiÞ dIm; j = 0; otherwise: When Eq. (14) holds, Eq. (13) can be expressed as k Im;e =

where Im is the radiative intensity evaluated at each of the discrete directions μm, the coefficients αm ± 1/2 are given by the recursion formula αm + 1=2 = αm−1=2 −2wm μm

the derivative of n(z) at zj. To complete the angular differencing procedure, we take the firstorder upwind method

n2e k 2 Im; j + signðμm Þne ð1−σj Þ n2j ( k ) 2 k 2 k 2 k 2 ðIm; j + 1 = nj + 1 −Im; j = nj ÞðIm; j = nj −Im; j−1 = nj−1 Þ k ðIm; j+1

k = n2j + 1 −Im; j−1

1≤ m ≤ M; 1≤ j ≤ N; k = 1; 2;:::;

= n2j−1 Þ

;

ð16Þ

1216

J.-M. Wang, C.-Y. Wu / International Communications in Heat and Mass Transfer 38 (2011) 1213–1218

0.6

k k When Eq. (14) does not hold, Im,e = (ne2/nj2)Im, j. Substituting k Eq. (12) and the expression of Im, e into Eq. (11), we can find

Monte Carlo method k k k−1 Am Im;w + Bb;j Ib;j + Dkm;j + wcmΔtΔz Im;j k j ; Im;j = wm Δz n2e c Δt + Am n2 + Bt;j + Cm

first-order DOM second-order DOM

0.4

j

Φ(z/L, t*)

MDOM

j

1≤ m ≤ M; 1≤ j ≤ N; k=1; 2; :::;

ð18Þ

The boundary and initial conditions can be expressed as k

k

Im;1 = 2 = I0 ; M = 2 + 1≤ m ≤ M; t > 0;

0.2

k

k

Im;N + 1 = 2 = 0; 1≤ m ≤ M = 2; t > 0; t* = 0.0563

0

ð19Þ

t* = 1.0

0

0.2

ð20Þ

t* = 0.125

0.4

0.6

0.8

0

Im; j = 0; 1≤ m ≤ M; 1 ≤ j ≤ N:

1

ð21Þ

z/L Fig. 4. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 0.1, n = 1 → 1.5, ω = 1 and a1 = 0.

Substituting Eqs. (12) and (16) into Eq. (11), we obtain (

)

Am signðμm Þn2e ð1−σj Þ

k

2

ðIm; j Þ k 2 k 2 n4j ðIm; j + 1 = nj + 1 −Im; j−1 = nj−1 Þ ( " 2 signðμm Þð1−σj Þ wm Δz n + Am e2 1 + k − k cj Δt nj Im;j + 1 = n2j + 1 −Im;j−1 = n2j−1 +

k Im;j

k k Am signðμm Þn2e ð1−σj Þ Im; j + 1 Im; j−1 2 2 2 k 2 + 1 = nj + 1 −Im;j−1 = nj−1 nj + 1 nj−1

k

Im; j + 1 n2j + 1

k

+

k

k

Im; j−1

!#

n2j−1

) k

+ Bt; j + Cm Im; j

k

+ Am Im;w + Bb; j Ib; j + Dm; j +

wm Δz k−1 = 0; I cj Δt m; j

1≤ m ≤ M; 1≤ j ≤ N; k = 1; 2;:::;

ð17Þ

k k In each time step, the intensities, Im, w and Ib, j, are known from a k previous spatial step or from boundary conditions. Since the Dm,j includes k k the unknown Im, ,j the formal solution of Im, obtained from Eq. (17) or (18) j is a system of equations for the intensities at time tk. The system of equations shall be solved by an iteration scheme similar to that presented in Ref. [12]. After solving the intensity field, one can obtain the directionintegrated intensity and the radiative flux by integrating the intensity over solid angle. The results in terms of dimensionless variables seem to be more general. For example, the dimensionless direction-integrated intensity is defined as



Φðz = L; t Þ =

Gðz = L; t  Þ = n2 ðz = LÞ ; 4πI0 = n20

ð22Þ

where t* = c0βt/n0 is the dimensionless time and the directionintegrated intensity G(z/L, t*) can be replaced by a weighted sum as M

k Gkj = 2π ∑ Im; j wm . m=1

k k If the Dm, j is assumed known, Eq. (17) is a quadratic equation of Im, j. It k k k is readily to find the formal solution of Im, in terms of I , I j m, j − 1 m, j + 1, k k k k−1 Im, w, Ib, j, Dm, j and Im, j from Eq. (17).

0.8 Monte Carlo method first-order DOM second-order DOM

0.6

Φ(z/L, t*)

MDOM

0.4

0.2 t* = 10.0 t*=1.25

t* = 0.563

0

0

0.2

0.4

0.6

0.8

1

z/L Fig. 5. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 1, n = 1 → 1.5, ω = 1 and a1 = 0.

3. Results and discussion First, a slice width (Δz =L/N) sensitivity study is performed to determine the width of slices used. The dimensionless directionintegrated intensity Φ(z/L, t*) in a slab with n=1→1.5 obtained by the DOM with second-order upwind scheme using Δt*=0.0001 for the case with τL =0.1, Δt*=0.001 for the case with τL =1 and Δt*=0.01 for the case with τL =10 are shown in Figs. 1–3, respectively. The notation n=n0 →nL represents a linear refractive index with n=n0 at z=0 and n=nL at z=L and the optical thickness is defined as τL =βL. The transient radiation field in an optically thin medium varies quickly. Thus, we use a smaller time step in the case with a smaller optical thickness. Besides, we have used a constant weight quadrature with 180 discrete values of μ similar to that used by Lemonnier and Le Dez [15]. From Figs. 1–3, we may find the following trends. (i) The influence of the slice width on the results around the incident plane at z=0 seems to be independent of the optical thickness. (ii) The influence of the slice width on the results around the wave front decreases with the increase of the optical thickness. (iii) The results obtained by using N=100 and N=200 are almost overlapping for all the cases considered. Thus, we adopt N=100 to obtain the discrete ordinates solutions shown in Figs. 4–10. The results obtained by using the DOM with the first-order upwind scheme, the DOM with second-order upwind scheme and the MDOM are shown in Figs. 4–10 for various combinations of parameters. Besides, to generate accurate results for comparison, we adapt the codes of the MCM developed by Wu [11] to solve the present problems and use large numbers of bundles, Nb =2×108 for the cases with τL =0.1 and Nb =109 for the cases with τL =1 and 10. Figs. 4–6 show the influence of the optical

J.-M. Wang, C.-Y. Wu / International Communications in Heat and Mass Transfer 38 (2011) 1213–1218

1

0.6 Monte Carlo method first-order DOM second-order DOM MDOM

0.8

Monte Carlo method first-order DOM second-order DOM

0.4

0.6

MDOM

Φ(z/L, t*)

Φ(z/L, t*)

1217

0.4

0.2 0.2 t* = 3.0

0

0

t* = 300.0

t* = 30.0

0.2

0.4

t* = 10.0 t* = 0.75

0.6

0.8

0

1

0

0.2

t* = 2.0

0.4

z/L

0.6

0.8

1

z/L

Fig. 6. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 10, n = 1 → 1.5, ω = 1 and a1 = 0.

Fig. 8. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 1, n = 1 → 3, ω = 1 and a1 = 0.

thickness, while Figs. 5, 7 and 8 show the influence of the gradient of refractive index. The curves of Φ(z/L, t*) for the cases with anisotropically scattering are shown in Figs. 9–10. Here, we lay some emphasis on radiative transfer in the early stage of a transient process. For the present cases the irradiation with μ = 1 on the interface z = 0 z takes t = ∫z = 0n(z′)dz′/c0 to arrive a location at z, which is the location of wave front in the early stage before radiation arriving the other boundary. The location of wave front can be expressed as z/L = {− n0 + [n02 + 2(nL − n0)t*/τL] 1/2}/(nL − n0). Thus, the location of wave front depends on the distribution of refractive index and the optical thickness. We can find the location of wave front z/L = 0.5 at t* = 0.0563 for the case with τL = 0.1 and n = 1 → 1.5, as shown in Fig. 4, z/L = 0.5 at t* = 0.563 for the cases with τL = 1 and n = 1 → 1.5, as shown in Figs. 5, 9 and 10, and z/L = 0.28 at t* = 3 for the case with τL = 10 and n = 1 → 1.5, as shown in Fig. 6. The wave front arrives z/L = 0.5 at t* = 0.625 for the case with τL = 1 and n = 1 → 2, while it arrives z/L = 0.5 at t* = 0.75 for the cases with τL = 1 and n = 1 → 3, as shown in Figs. 7 and 8. Figs. 4–10 show that the results by the MDOM and MCM are in excellent agreement for all the cases considered and the MDOM is capable of

capturing the wave front. However, the Φ(z/L, t*) curves of optically thin and moderate cases in the beginning of a transient process show that the first-order DOM may induce numerical diffusion, resulting in early transmitted radiation. The Φ(z/L, t*) curves obtained by the secondorder DOM show similar trend for optically thin and moderate cases, but the effect of numerical diffusion is less. By contrast, for the cases with an optically thick medium, the multiple scattering of radiation has more influence than the advection of radiation and so the numerical diffusion due to the finite difference of the advection term of the radiative transfer equation is minor, as shown in Fig. 6. Similarly, the influence of the multiple scattering also increases with the increase of the time. Thus, the Φ(z/L, t*) curves obtained by the four different methods are in excellent agreement after the beginning stage of a transient process. The above trends also hold for various gradients of refractive index and anisotropically scattering coefficients (a1), as shown in Figs. 7–10. From the comparison of Figs. 5, 7 and 8, we can see that the influence of numerical diffusion on the transient solutions is a little less noticeable for the case with a larger gradient of refractive index. This is because the larger gradient of refractive index causes more internal reflection and

0.8

0.8 Monte Carlo method

Monte Carlo method

first-order DOM

first-order DOM

second-order DOM

0.6

second-order DOM

0.6

MDOM

Φ(z/L, t*)

Φ(z/L, t*)

MDOM

0.4

0.2

0.4

0.2 t* = 10.0 t* = 0.625

0

0

0.2

t* = 10.0

t* = 1.5

0.4

0.6

t* = 1.25

t* = 0.563

0.8

1

z/L Fig. 7. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 1, n = 1 → 2, ω = 1 and a1 = 0.

0

0

0.2

0.4

0.6

0.8

1

z/L Fig. 9. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 1, n = 1 → 1.5, ω = 1 and a1 = 1.

1218

J.-M. Wang, C.-Y. Wu / International Communications in Heat and Mass Transfer 38 (2011) 1213–1218

0.8

present results shows that the influence of numerical diffusion is a little less noticeable for the case with a larger gradient of refractive index. Besides, the distribution of Φ(z/L, t*) around z = 0 decreases and that around z = L increases with the increase of anisotropically scattering coefficient a1.

Monte Carlo method first-order DOM second-order DOM

0.6

Φ(z/L, t*)

MDOM

Acknowledgments

0.4 This work was supported by the National Science Council of the Republic of China on Taiwan through Grant NSC 97-2221-E-006-225.

0.2 References

t* = 10.0 t* = 1.25

t* = 0.563

0

0

0.2

0.4

0.6

0.8

1

z/L Fig. 10. Comparison of the transient distributions of direction-integrated intensity obtained by the second-order DOM, the first-order DOM, the MDOM and the MCM for the case with τL = 1, n = 1 → 1.5, ω = 1 and a1 = − 1.

scattering, and so reduces the influence of numerical diffusion. Besides, comparison of Figs. 5, 9 and 10 shows the influence of the anisotropically scattering coefficient on the distribution of Φ(z/L, t*). The Φ(z/L, t*) around the irradiation boundary (z = 0) decreases and that around the other boundary (z = L) increases with the increase of a1. 4. Conclusion In this work, we applied the DOM with a second-order upwind scheme to solve transient radiative transfer in a graded index slab suddenly exposed to a diffuse strong irradiation at z = 0. The planar medium is absorbing and anisotropically scattering. From the comparison of the results obtained by the first-order DOM, the second-order DOM, the MDOM and the MCM, it may be concluded that the numerical diffusion in the transient solutions obtained by the second-order DOM is less than that appears in the solutions obtained by the first-order DOM, but the numerical diffusion is still noticeable and results in early transmitted radiation, especially for optically thin and moderate cases. By contrast, for optically thick cases, the multiple scattering of radiation has more influence than the advection of radiation and so the numerical diffusion due to the finite difference of the advection term of the radiative transfer equation is minor. The MDOM is an effective method to capture the wave front of radiation propagation, while it is still necessary to adopt a DOM with a higher order scheme to capture the wave front accurately. Comparison of the

[1] Y. Yamada, Light-tissue interaction and optical imaging in biomedicine, in: C.L. Tien (Ed.), Annual Review of Heat Transfer, vol. 6, Begell House, New York, 1995, pp. 1–59. [2] M.Y. Kim, S. Menon, S.W. Baek, On the transient radiative transfer in a onedimensional planar medium subjected to radiative equilibrium, International Journal of Heat and Mass Transfer 53 (2010) 5682–5691. [3] S.M. Chitanvis, A. Zardecki, Effect of thermal blooming on pulse propagation through vaporizing aerosols, Applied Optics 27 (1988) 2495–2501. [4] R. Siegel, C.M. Spückler, Variable refractive index effects on radiation in semitransparent scattering multilayered regions, Journal of Thermophysics and Heat Transfer 7 (1993) 624–630. [5] H.A. Ferwerda, The radiative transfer equation for scattering media with a spatially varying refractive index, Journal of Optics A: Pure and Applied Optics 1 (1999) L1–L2. [6] J.M. Tualle, E. Tinet, Derivation of the radiative transfer equation for scattering media with a spatially varying refractive index, Optical Communications 228 (2003) 33–38. [7] M. Premaratne, E. Premaratne, A.J. Lowery, The photon transport equation for turbid biological media with spatially varying isotropic refractive index, Optical Express 2 (2005) 389–399. [8] C.-Y. Wu, Discrete ordinates solution of transient radiative transfer in refractive planar media with pulse irradiation, in: G. de Vahl Davis (Ed.), The Annals of the Assembly for International Heat Transfer, Vol. 13, Begell House Inc., 2006, doi:10.1615/IHTC13.p4.120, Published online. [9] L.H. Liu, P.-F. Hsu, Analysis of transient radiative transfer in semitransparent graded index medium, Journal of Quantitative Spectroscopy and Radiative Transfer 105 (2007) 357–376. [10] M. Bhuvaneswari, C.-Y. Wu, Differential approximations for transient radiative transfer in refractive planar media with pulse irradiation, Journal of Quantitative Spectroscopy and Radiative Transfer 110 (2009) 389–401. [11] C.-Y. Wu, Monte Carlo simulation of transient radiative transfer in a medium with a variable refractive index, International Journal of Heat and Mass Transfer 52 (2009) 4151–4159. [12] J.-M. Wang, C.-Y. Wu, Transient radiative transfer in a scattering slab with variable refractive index and diffuse substrate, International Journal of Heat and Mass Transfer 53 (2010) 3709–3806. [13] B. van Leer, Towards the ultimate conservation difference scheme. IV. A new approach to numerical convection, Journal of Computational Physics 23 (1977) 276–299. [14] P. Ben Abdallah, V. Le Dez, Temperature field inside an absorbing–emitting semitransparent slab at radiative equilibrium with variable spatial refractive index, Journal of Quantitative Spectroscopy and Radiative Transfer 65 (2000) 595–608. [15] D. Lemonnier, V. Le Dez, Discrete ordinates solution of radiative transfer across a slab with variable refractive index, Journal of Quantitative Spectroscopy and Radiative Transfer 73 (2002) 195–204.