Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

68 A. Mezzacappa et al.

[29] Rampp M, Janka HT. Astronomy and Astrophysics 396:361 (2002)

[30] Mezzacappa A, Messer O. Astrophysical Journal submitted (2005)

[31] Burrows A, Thompson T. in Stellar Collapse, ed. C Fryer. Dordrecht: Kluwer

Academic Publishers (2004)

[32] Fuller GM, Mayle R, Wilson JR, Schramm D. Astrophysical Journal 322:795

(1987)

[33] Fuller GM, Mayle R, Meyer BS, Wilson JR. Astrophysical Journal 389:517

(1992)

[34] Qian YZ, Fuller G. Physical Review D 52:656 (1995)

[35] Fuller G, Qian YZ. Physical Review D submitted (2005)

[36] Pantaleone J. Physics Letters B 287:128 (1992)

[37] Pantaleone J. Physics Letters B 342:250 (1995)

[38] Sigl G, Raﬀelt G. Nuclear Physics B 406:423 (1993)

[39] Mezzacappa A, Bruenn SW. in The Identiﬁcation of Dark Matter,ed.

N Spooner, V Kudryavtsev. Singapore: World Scientiﬁc (1999)

Discrete-Ordinates Methods

for Radiative Transfer

in the Non-Relativistic Stellar Regime

Jim E. Morel

Mail Stop D413, Los Alamos National Laboratory, Los Alamos, New Mexico 87545

jim@lanl.gov

1 Introduction

The purpose of this paper is to brieﬂy describe the application of discrete-

ordinates methods to radiative transfer in the non-relativistic stellar regime.

We consider the non-relativistic regime to be characterized by material ve-

locities less than or equal to one percent of the speed of light. In most appli-

cations, the radiative transfer equations are coupled to the hydrodynamics

equations. We ﬁrst describe a radiation-hydrodynamics model that consists

of the Euler equations and the radiation transport equation together with

approximations for material-motion eﬀects based upon the assumption of

non-relativistic material velocities. For simplicity, we next focus on a model

of radiative transfer in a static medium that is adequate for the purpose of

describing numerical discretization and solution techniques for the transfer

equation that can also be applied in the more general context of radiation-

hydrodynamics.

2 The Approximate Radiation-Hydrodynamics Model

The hydrodynamics equations and the radiation moment equations can be

expressed to O(v/c) as follows assuming a grey radiation transport approxi-

mation with the transfer equation cast in the Eulerian frame:

• Conservation of ﬂuid mass

∂

ρ + ∂

(ρv

)=0, (1)

• Conservation of ﬂuid momentum

∂

(ρv

)+∂

(ρv

)+∂

p =

0,i

−



− E



, (2)

70 J.E. Morel

• Conservation of total ﬂuid energy:

∂



ρv

+ ρe



+ ∂



ρv

+ ρe + p





−cσ



− E



0,i

, (3)

• Conservation of radiation momentum

∂

+ ∂

= −

0,i



− E



, (4)

• Conservation of radiation energy

∂

E + ∂

= cσ



− E



−

0,i

, (5)

where ρ is the ﬂuid density, v

is component i of the ﬂuid velocity, c is the

speed of light, p is the ﬂuid pressure, e is the speciﬁc internal energy of the

ﬂuid, E is the radiation energy density, F

is component i of the radiation ﬂux,

is element ij of the radiation pressure tensor, E

is the comoving-frame

radiation energy density:

= E −

− v

E − v

) , (6)

0,i

is component i of the comoving-frame radiation ﬂux:

0,i

= F

− v

E − v

, (7)

is the macroscopic total cross section, σ

is the macroscopic absorption

cross section, and a is the radiation constant. The source terms in (2) through

(5) are expressed in terms of comoving-frame radiation variables. These

comoving-frame variables are deﬁned in terms of Eulerian-frame radiation

variables by (6) and (7). When the comoving-frame variables are eliminated

from (2) through (5) using (6) and (7), certain terms of O(v

) will ap-

pear. Even though (2) through (5) are only correct to O(v/c), these terms are

not to be neglected because they eliminate spurious equilibrium solutions [1]

and maintain consistency between the Eulerian-frame and comoving-frame

momentum and energy sources. We consider the nonrelativistic limit to be

characterized by v/c ≤ 0.01.

We next describe a very simple approximate model for radiation-hydro-

dynamics with multifrequency radiative transfer in the nonrelativistic limit.

The goal of this approximation is to maintain accuracy in an integral sense

while avoiding the complexities of the equations that are correct to O(v/c).

Since v/c is very small in the nonrelativistic limit, one might be tempted

to neglect material motion corrections to the radiative transfer equation en-

tirely. However, over the course of a calculation, this can result in a signif-

icant loss of energy conservation because the kinetic energy change in the

Discrete-Ordinates Methods for Radiative Transfer 71

ﬂuid due radiation momentum deposition is not removed from the radia-

tion energy ﬁeld. Thus we use a simple approximation for nonrelativistic

radiation-hydrodynamics that has the following properties:

• Total energy and momentum are conserved.

• Equilibrium solutions are preserved to O(

• The equilibrium diﬀusion limit is preserved to O(

In the nonrelativistic case, there is little diﬀerence between E and E

, but the

relative diﬀerence between F

0,i

and F

can be arbitrarily large. This relative

diﬀerence is essentially inﬁnite in the equilibrium state because F

0,i

while F

E to O(v/c). Since the system is locally equilibrated in the

equilibrium diﬀusion limit, it follows that one can expect a signiﬁcant relative

diﬀerence between F

0,i

and F

in this limit. Hence it is generally considered

important to preserve the equilibrium diﬀusion limit in any approximate

treatment of nonrelativistic radiation-hydrodynamics. See [1] for a formal

asymptotic derivation of the non-relativistic equilibrium diﬀusion limit.

The ﬁrst step in the derivation of our approximate model is modify the

hydrodynamics equations and the grey radiation moment equations. Speciﬁ-

cally, we make the following substitutions:

0,i

−



− E



⇒

0,i

, (8a)

cσ



− E



−

0,i

⇒ cσ



− E



−

0,i

, (8b)

0,i

= F

− v

E − v

⇒ F

−

E, (8c)

The term neglected in (8a) is a purely relativistic term. Although the radia-

tive transfer equation is always relativistic, the hydrodynamics equations are

purely classical to O(v/c). Furthermore, this term is zero in the equilibrium

diﬀusion limit. Thus it is reasonable to eliminate this term. In (8b) we replace

with E in the radiation energy source term. It can be seen from (6) that

this should be a good approximation when v/c ≤ 0.01. Furthermore, E

= E

in the equilibrium diﬀusion limit. Finally, we substitute E/3forP

i,j

in (8c).

This approximation is made simply to avoid calculating the radiation pressure

tensor. This tensor plays a small role in a small term, and the substitution is

exact in the equilibrium diﬀusion limit. The next step in the derivation is to

replace the modiﬁed grey expressions with frequency-dependent expressions

that are equivalent under the grey approximation. In particular,

0,i

⇒



∞



4π

(ν)

I(Ω, ν)



Ω

−



dΩ dν , (9a)

cσ



− E



⇒



∞



4π

(ν)



B(ν) − I



−→

Ω

,ν



dΩ dν , (9b)

0,i

⇒



∞



4π



−→

Ω

,ν





Ω

−



dΩ dν , (9c)

72 J.E. Morel

where

−→

Ω

is the photon direction vector, ν is the photon frequency, I



−→

Ω

,ν



is the angular intensity, B =

2hν



exp



hν



− 1



−1

is the Planck function

and k is the Boltzmann constant. Expression (9b) represents an exact substi-

tution, and (9c) represents a substitution that is exact to O(v/c). However,

(9a) is simply incorrect unless σ

is independent of frequency. The transfor-

mations between the Eulerian and co-moving frames being used here are only

correct for angle-energy-integrated quantities [2]. The transformation relat-

ing the co-moving frame ﬂux and the Eulerian frame ﬂux is eﬀectively being

used for frequency-dependent ﬂuxes in (9a).

Accounting for all of these substitutions, our approximate hydrodynamic

equations and multifrequency radiation moment equations can be expressed

as follows:

∂

ρ + ∂

(ρv

)=0 , (10)

∂

(ρv

)+∂

(ρv

)+∂

p =



∞



4π

(ν)



−→

Ω

,ν





Ω

−



dΩ dν , (11)

∂



ρv

+ ρe



+ ∂



ρv

+ ρe + p





−



∞



4π

(ν)



B(ν) − I



−→

Ω

,ν



dΩ dν



∞



4π

(ν)



−→

Ω

,ν





Ω

−



dΩ dν , (12)

∂

+ ∂

−



∞



4π

(ν)



−→

Ω

,ν





Ω

−



dΩ dν , (13)

∂

E + ∂



∞



4π

(ν)



B(ν) − I



−→

Ω

,ν



dΩ dν

−



∞



4π

(ν)



−→

Ω

,ν





Ω

−



dΩ dν . (14)

The radiative transfer equation for static media can be made consistent with

these model equations by adding a P

-type correction term to the right side

of the equation:

Discrete-Ordinates Methods for Radiative Transfer 73

∂I

∂t

−→

Ω

−→

∇

I + σ

I =

4π

φ + σ

B(T )+

4π

f,i

Ω

, (15)

where

φ =



4π



−→

Ω



,ν



dΩ



, (16)

= −



4π

(ν)



−→

Ω



,ν



Ω

−



dΩ



, (17)

f,i



4π

(ν)I



−→

Ω



,ν



dΩ



. (18)

Note that if we multiply (15) by Ω

/c and integrate over all directions and

frequencies, we obtain (13), and if simply we integrate (15) over all direc-

tions and frequencies, we obtain (14). We stress that (15) is an approximate

Eulerian frame transfer equation that is consistent with the other approx-

imations that we have made in the hydrodynamic and radiation moment

equations. The Eulerian frame transfer equation that is correct to O(v/c)

is considerably more complicated than (15) [2]. Our approximate radiation-

hydrodynamics model is completely deﬁned by (10) through (18).

3 Discretization and Solution Techniques

In this section, we discuss discretization and solution techniques for the radia-

tive transfer equation. For simplicity, we consider only the radiative transfer

equation for a static medium and neglect the hydrodynamics equations. In

general, the techniques we discuss here can be applied within a radiation-

hydrodynamics calculation. A discussion of discretization and solution tech-

niques for the full radiation-hydrodynamics equations is beyond the scope of

this paper. The equations we consider consist of a radiative transfer equation:

∂I

∂t

−→

Ω

−→

∇

I + σ

I =

4π

φ + σ

B(T ) , (19)

and an equation for the material temperature T (r, t):

∂T

∂t



∞

[ φ − 4πB(T )]dν , (20)

where C

is the heat capacity and φ is deﬁned by (16).

Equations (19) and (20) are discretized in direction using the standard S

approximation [2,3]. This is basically collocation at a set of discrete directions

that correspond to quadrature points on the unit sphere. The quadrature set

74 J.E. Morel

associated with those discrete directions is also used to perform directional

integrations. For instance, let {

−→

Ω

}

n=1

denote an M-point quadrature

set, then (19) and (20) are directionally discretized as follows:

∂I

∂t

−→

Ω

−→

∇

+ σ

4π

φ + σ

B(T ) , (21)

∂T

∂t



∞

[ φ − 4πB(T )]dν , (22)

where

= I



−→

Ω



, (23)

and

φ =



n=1

. (24)

Equations (19) and (20) are discretized in frequency using the stan-

dard multigroup method [2, 3]. The ﬁrst step in this method is parti-

tion the frequency domain into G disjoint contiguous intervals or groups:



[ν

,ν

),...,[ν

G−

,ν

G+

)



where the minimum frequency is ν

and the

maximum frequency is ν

G+

.Theg’th group refers to the interval,

[ν

g−

,ν

). Next (19) and (20) are discretized as follows:

∂I

∂t

−→

Ω

−→

∇

+ σ

t,g

4π

s,g

+ σ

a,g

(T ) ,g =1,G, (25)

∂T

∂t



a,g



[ φ



− 4πB



(T )] , (26)

where



g−

I(ν



) dν



, (27)



g−

B(ν



) dν



, (28)

and σ

denotes any desired average of σ(ν)overgroupg:

= σ(ν)

ν∈[ν

g−

,ν

)

. (29)

The physical characteristics of radiative transfer problems place severe

demands upon spatial discretization schemes. In particular, such schemes

must:

Discrete-Ordinates Methods for Radiative Transfer 75

• be at least second-order accurate,

• be highly damped in strongly absorbing media,

• be accurate in optically-thin regions,

• possess the thick diﬀusion limit,

• behave well with unresolved spatial boundary layers.

In general, lumped discontinuous spatial discretization schemes are required

for radiative transfer. These schemes are usually either ﬁnite-element based

[4, 5] or ﬁnite-volume based [6]. While lumping of the removal and source

terms is well deﬁned, lumping of the gradient term is also required in multidi-

mensional calculations [7]. There is no general procedure for gradient lumping

and it must generally be performed on a case-by-case basis. Gradient lumping

remains a research topic. The need to possess the thick diﬀusion limit places

additional demands upon spatial discretization schemes. It is important to

explain exactly what it means for a spatial discretization scheme to possess

the thick diﬀusion-limit. In this limit, the transport solution is diﬀusive and

the spatial cells can be arbitrarily thick with respect to a mean-free-path.

The diﬀusion length is the spatial scalelength for diﬀusive solutions, and the

diﬀusion length can be arbitrarily large with respect to a mean-free-path.

Hence, a spatial mesh with cells that are thin with respect to a diﬀusion

length but arbitrarily thick with respect to a mean-free-path might be intu-

itively expected to yield an accurate solution because the spatial variation of

the exact solution is well resolved by the mesh. Indeed, a scheme that yields

an accurate solution for highly diﬀusive problems whenever the variation of

the exact diﬀusive solution is well resolved by the mesh is said to possess the

thick diﬀusion limit. However, it is important to recognize that convergent

schemes, i.e. schemes with truncation errors that go to zero as the mesh is

reﬁned, do not necessarily possess the thick diﬀusion limit. A truncation error

analysis yields no information on the behavior of a scheme in this limit. Such

analyses indicate that accurate solutions will be obtained if the mesh cells

are thin with respect to a mean-free-path. Spatial resolution on the order of

a mean-free-path is completely impracical in highly diﬀusive problems. For

instance, as a problem is made increasingly diﬀusive, the characteristic width

of the system approaches a ﬁxed value when measured in diﬀusion lengths

and becomes inﬁnite when measured in mean-free-paths. To determine the

behavior of a spatial discretization scheme in the thick diﬀusion limit, one

must perform a discrete asymptotic analysis [4, 7–10].

Boundary layers that are a few mean-free-paths thick can exist at the

transitions between non-diﬀusive and diﬀusive regions. It can be very diﬃcult

to resolve such layers. In general, an adaptive approach is required because the

locations of these boundary layers usually vary dynamically. Furthermore, in

a multifrequency calculation one may not be able to resolve a boundary layer

on a mean-free-path basis for all frequencies. Thus it is desirable to have

schemes that behave well if not highly accurately with unresolved spatial

76 J.E. Morel

boundary layers. A discrete asymptotic analysis will reveal this behavior in

addition to the behavior of the scheme on the interior of a diﬀusive region.

The diﬀusion limit for radiative transfer is characterized by:

I = B(T ) , (30)

+4aT

)

∂T

∂t

−

−→

∇



4acT

3σ



−→

∇

T =0 , (31)

where σ

is the Rosseland-averaged total cross section:





∞

∗

(ν)

∂B(ν)

∂T

dν

 



∞

∂B(ν)

∂T

dν

!

−1

. (32)

Equations (19) and (20) are usually solved via a modiﬁed Newton’s

method. To illustrate this technique, we ﬁrst implicitly diﬀerence these equa-

tions in time:

c∆t



− I

k−



−→

Ω

−→

∇

+ σ

4π

+ σ

, (33)

∆t



− T

k−





∞



− 4πB



dν . (34)

We let T

∗

denote the latest Newton iterate for the temperature. Then the

linearized equations for the next Newton iteration are obtained by evaluating

the material properties at T

∗

and linearly expanding the Planck function

temperature dependence about T

∗

n+1

= B

∗

∂B

∗

∂T



− T

∗



, (35)

where a superscript “*” denotes a quantity evaluated at T

∗

.Thisisamodi-

ﬁed Newton’s method because the contribution of the material properties to

the Jacobian are neglected. With the expansion given in (35), the material

temperature can be eliminated from the transport equation. Suppressing the

temporal superscript “k +

”, the linearized temporally-diﬀerenced transport

equation can be expressed as:

−→

Ω

−→

∇

I + σ

∗

I =

4π

∗

φ +

4π

νχ



∞

∗

(ν



)φ(ν



) dν



+ ξ, (36)

where:

Discrete-Ordinates Methods for Radiative Transfer 77

= σ

+ τ, (37a)

τ =

c∆t

, (37b)

ν =



∞

∗

(ν)4π

∂B

∗

(ν)

∂T

dν

∗

∆t



∞

∗

(ν)4π

∂B

∗

(ν)

∂T

dν

, (37c)

χ(ν)=

∗

(ν)

∂B

∗

(ν)

∂T



∞

∗

(ν



)

∂B

∗

(ν



)

∂T

dν



, (37d)

ξ = σ

∗

+ τψ

k−

−

4π

νχ





∞

∗

(ν



)4πB

∗

(ν



) dν



∗

∆t



k−

− T

∗





, (37e)

and the material temperature is given by

= T

∗



∞

∗

(ν)[φ(ν) − 4πB

∗

(ν)] dν +

∗

∆t



k−

− T

∗



∗

∆t



∞

∗

(ν)4π

∂B

∗

(ν)

∂T

dν

. (38)

Equation (38) is used to calculate the temperatures after the linearized trans-

port equation has been solved. We note that (36) has the form of the steady-

state neutron transport equation with σ

playing the role of the ﬁssion cross

section, ν playing the role of the number of neutrons per ﬁssion, χ playing

the role of the ﬁssion spectrum, and ξ playing the role of a inhomogeneous

source.

The basic technique for solving (36) is the source iteration technique.

This method is based upon the fact that all coupling between direction and

frequency occurs on the right side of the equation. The operator on the left

side of (36) can be discretized in such a way that a block lower-triangular

coeﬃcient matrix is obtained for each direction and frequency group where

the unknowns in each block are those associated with a single spatial cell.

Such matrices are easily inverted in a direct manner. The inversion process

for each discrete direction and frequency starts where radiation enters the

mesh, proceeds across the mesh, and ends where radiation leaves the mesh.

Hence the inversion process is referred to as a sweep. The iterative process

can be represented in a simpliﬁed form as follows:

−→

Ω

−→

∇

+1

+ σ

∗

+1

4π

∗



4π

νχ



∞

∗

(ν



)φ



(ν



) dν



+ ξ,(39)

where  is the iteration index. The actual iterative process is nested such that

the inner iterations converge the scattering sources for each frequency group

while the absoprtion rate is held ﬁxed, and the outer iterations converge

the absorption rate. The inner source iteration process for group g can be

expressed as