Graziani F. (editor) Computational Methods in Transport

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

Transport Approximations in Partially Diﬀusive Media 397

the case of a non-scattering tube, the limiting equation for U

as the thickness

of the tube tends to 0 is thus given by the local generalized diﬀusion model:

−∇ · D

∇U

+ σ

= q, Ω\Σ

∂U

∂n

=0,∂Ω

(x) lim

η→0



θ ·∇

⊥

(x + ηθ)dµ(θ)=−∇

· D

(x)∇

,Σ.

(92)

It may be easier to understand the solution of the above equation by recasting

it in a variational form: Find the unique solution U ∈ H

(Ω) such that for

all V ∈ H

(Ω), we have:



Ω

∇U ·∇V + σ

UV)dx +



∂Ω

UVdσ(x)



∇

U∇

Vdl(x)=



Ω

Vgdx . (93)

When n =3andd =2,thenl(x) is the (one-dimensional) arclength along

the curve Σ and ∇

is nothing but

∂

∂l

. The Hilbert space H

(Ω) is deﬁned

as the completion of the pre-Hilbert space of functions of class C

on Ω such

that U, U 

< ∞, where the inner product ·, ·

is deﬁned as:

U, V 



Ω

(∇U ·∇V + UV)dx +



∇

U∇

Vdl(x) . (94)

We know that H

(Ω) is a Hilbert space [28] and that thanks to the uniform

positiveness of D

, σ

,andD

on their domains of deﬁnition, the variational

formulation (93) admits a unique solution by the Lax Milgram theory.

Note that functions in H

(Ω), deﬁned in (31) as the natural space for

the classical diﬀusion approximation, do not necessarily admit traces on one-

dimensional curves (though functions in H

1+δ

(Ω) do for all δ>0). This

renders the use of the above completion argument necessary to construct

(Ω). This issue does not arise for (co-dimension one) surfaces as functions

in H

(Ω) indeed admit traces on surfaces. In any event, from the numerical

viewpoint, we observe that the non-scattering ﬁlament is simply modeled by

one additional integration over Σ in the variational formulation (93). This

renders its numerical simulation rather straightforward and much less costly

computationally than the full transport equation or the non-local diﬀusion

model (87).

Acknowledgments

I would like to thank Frank Graziani for inviting me to the very nice confer-

ence on computational methods in transport organized at the Granlibakken

398 G. Bal

conference center in September of 2004. I enjoyed many stimulating discus-

sions during my stay there. This research was funded in part by NSF Grant

DMS-0239097 and an Alfred P. Sloan fellowship.

A Local Second-Order Equation and Linear Corrector

This appendix provides some useful calculations related to the variational

formulations introduced in the text and sketches the derivation of the ex-

pression for the corrector of order ε.Weﬁrstverifythatbychoosingatest

function with support in Ω and by integrations by parts that the solution u

of (15) also solves

(I − ω ·∇G

−1

)(ω ·∇u + G

u − εq)=0. (95)

Let us introduce the shift operators

= F

−1

F, (

ˆu)

=ˆu

n±1

. (96)

We then verify that

ω ·∇= ∂S

−

∂S

Thanks to (40) and the above equalities, we deduce that (95) is equivalent

to:

−



∂

∂ˆu

n−2

− k

n−1

∂

∂ˆu

− k

n+1

+ ∂

∂ˆu

− k

n−1

∂

∂ˆu

n+2

− k

n+1



− k

ˆu

=ˆq

− ε



∂

ˆq

n−1

− k

n−1

∂

ˆq

n+1

− k

n+1



. (97)

Let u be the transport solution. The equation for its orthogonal projection

over linear functions in the angular variable Π

u is thus

−



∂

∂ˆu

−1

− k

∂

∂ˆu

− k

+ ∂

∂ˆu

− k



− k

ˆu

=ˆq

− ε



∂

ˆq

− k

∂

ˆq

− k



Up to terms of smaller order that can be estimated, this is

−



∂

∂ˆu

−1

− k

+ ∂

∂ˆu

− k



− k

ˆu

= −∂

ˆq

εσ

We check that this is also the equation veriﬁed by

= −G

−1

(ω ·∇U

) ,

up to smaller-order terms, which may be written in the Fourier domain as

Transport Approximations in Partially Diﬀusive Media 399



−1

(ω ·∇U

)

= ε

∂

− k

where U

is the diﬀusion approximation. Therefore, the expression for the

corrector is:

(x, ω)=−εH

(ω) ·∇U

(x) . (98)

This is the usual expression for the ﬁrst-order corrector to the transport

solution in the absence of boundaries.

References

1. R. T. Ackroyd, Completely boundary free minimum and maximum principles

for neutron transport and their least-squares and galerkin equivalents, Ann.

Nucl. Energy, 9 (1981), p. 95.

2. S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, The ﬁnite

element model for the propagation of light in scattering media: A direct method

for domains with nonscattering regions, Med. Phys., 27(1) (2000), pp. 252–264.

3. G. Bal, Particle transport through scattering regions with clear layers and in-

clusions, J. Comp. Phys., 180(2) (2002), pp. 659–685.

, Transport through diﬀusive and non-diﬀusive regions, embedded objects,

and clear layers, SIAM J. Appl. Math., 62(5) (2002), pp. 1677–1697.

5. G. Bal and Y. Maday, Coupling of transport and diﬀusion models in linear

transport theory, M2AN Math. Model. Numer. Anal., 36(1) (2002), pp. 69–86.

6. G. Bal and K. Ren, Generalized diﬀusion model in optical tomography with

clear layers, J. Opt. Soc. Amer. A, 20(12) (2003), pp. 2355–2364.

7. L. Bourhrara, New variational formulations for the neutron transport equa-

tion, Transport Theory Statist. Phys., 33(2) (2004), pp. 93–124.

8. P. S. Brantley and E. W. Larsen, The simpliﬁed P

approximation,Nucl.

Sci. Eng, 134 (2000), pp. 1–21.

9. S. C. Brenner and L. R. Scott, The mathematical theory of ﬁnite element

methods, Springer verlag, New York, 2002.

10. S. Chandrasekhar, Radiative Transfer, Dover Publications, New York, 1960.

11. M. Clark and G. C. Pomraning, The Variational Method Applied to the

Monoenergetic Boltzmann Equation, Part I, Nucl. Sci. Eng., 16 (1963), pp. 147–

154.

12. R. Dautray and J.-L. Lions, Mathematical Analysis and Numerical Methods

for Science and Technology. Vol.6, Springer Verlag, Berlin, 1993.

13. H. Dehghani, D. T. Delpy, and S. R. Arridge, Photon migration in non-

scattering tissue and the eﬀects on image reconstruction, Phys. Med. Biol., 44

(1999), pp. 2897–2906.

14. J. J. Duderstadt and W. R. Martin, Transport Theory, Wiley-Interscience,

New York, 1979.

15. F. Golse, S. Jin, and C. D. Levermore, The convergence of numerical

transfer schemes in diﬀusive regimes. I. Discrete-ordinate method.,SIAMJ.

Numer. Anal., 36 (1999), pp. 1333–1369.

16. F. Golse, S. Jin, and C. D. Levermore, A domain decomposition analysis

for a two-scale linear transport problem, M2AN Math. Model. Numer. Anal.,

37(6) (2003).

400 G. Bal

17. A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, Comparison of

ﬁnite-diﬀerence transport and diﬀusion calculations for photon migration in

homogeneous and heterogeneous tissues, Phys. Med. Biol., 43 (1998), pp. 1285–

1302.

18. A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE Press,

New York, 1997.

19. E. W. Larsen and J. B. Keller, Asymptotic solution of neutron transport

problems for small mean free paths, J. Math. Phys., 15 (1974), pp. 75–81.

20. E. W. Larsen and J. E. Morel, Asymptotic Solutions of Numerical Trans-

port Problems in Optically Thick Diﬀusive Regimes II,J.Comp.Phys.,83

(1989), p. 212.

21. E. W. Larsen, J. E. Morel, and W. F. Miller Jr., Asymptotic solutions

of numerical transport problems in optically thick, diﬀusive regimes,J.Comp.

Phys., 69 (1987), pp. 283–324.

22. C. D. Levermore, W. J. Morokoff, and B. T. Nadiga, Moment realizabil-

ity and the validity of the Navier-Stokes equations for rareﬁed gas dynamics,

Phys. Fluids, 10(12) (1998), pp. 3214–3226.

23. E. E. Lewis and W. F. Miller Jr., Computational Methods of Neutron

Transport, J. Wiley and sons, New York, 1984.

24. T. A. Manteuffel and K. Ressel, Least-squares ﬁnite-element solution of

the neutron transport in diﬀusive regimes, SIAM J. Numer. Anal., 35(2) (1998),

pp. 806–853.

25. T. A. Manteuffel, K. Ressel, and G. Starke, A boundary functional

for the least-squares ﬁnite-element solution of the neutron transport equation,

SIAM J. Numer. Anal., 37(2) (2000), pp. 556–586.

26. J. E. Morel and J. M. McGhee, A self-adjoint angular ﬂux equation,Nucl.

Sci. Eng., 132 (1999), pp. 312–325.

27. J. Planchard, M´ethodes math´ematiques en neutronique (in French), Collec-

tion de la Direction des Etudes et Recherches d’EDF, Eyrolles, 1995.

28. M. Reed and B. Simon, Methods of modern mathematical physics. I. Func-

tional analysis, Academic Press, Inc., New York, second ed., 1980.

29. J. Ripoll, M. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani,

Boundary conditions for light propagation in diﬀuse media with non-scattering

regions, J. Opt. Soc. Amer. A, 17(9) (2000), pp. 1671–1681.

30. R. P. Rulko, D. Tomasevic, and E. W. Larsen, Variational P1 Approxi-

mations of General-Geometry Multigroup Transport Problems, Nucl. Sci. Eng.,

121 (1995), p. 393.

31. H. Sato and M. C. Fehler, Seismic wave propagation and scattering in the

heterogeneous earth, AIP series in modern acoustics and signal processing, AIP

Press, Springer, New York, 1998.

32. W. M. Stacey, Variational Methods in Nuclear Reactor Physics, Academic

Press, New York, 1974.

33. M. D. Tidriri, Asymptotic analysis of a coupled system of kinetic equations,

C.R. Acad. Sci. Paris, t.328, S´erie I Math., (1999), pp. 637–642.

34.

, Rigorous derivation and analysis of coupling of kinetic equations and

their hydrodynamic limits for a simpliﬁed boltzmann model, J. Statist. Phys.,

104 (2001).

High Order Finite Volume Nonlinear Schemes

for the Boltzmann Transport Equation

Barna L. Bihari

and Peter N. Brown

Center for Applied Scientiﬁc Computing, Lawrence Livermore National

Laboratory, Livermore, California 94550

bihari@llnl.gov

Center for Applied Scientiﬁc Computing, Lawrence Livermore National

Laboratory, Livermore, California 94550

pnbrown@llnl.gov

Abstract. We apply the nonlinear WENO (Weighted Essentially Nonoscillatory)

scheme to the spatial discretization of the Boltzmann Transport Equation model-

ing linear particle transport. The method is a ﬁnite volume scheme which ensures

not only conservation, but also provides for a more natural handling of bound-

ary conditions, material properties and source terms, as well as an easier parallel

implementation and post processing. It is nonlinear in the sense that the stencil

depends on the solution at each time step or iteration level. By biasing the gradient

calculation towards the stencil with smaller derivatives, the scheme eliminates the

Gibb’s phenomenon with oscillations of size O(1) and recudes them to O(h

), where

h is the mesh size and r is the order of accuracy. Our current implementation is

three- dimensional, generalized for unequally spaced meshes, fully parallelized, and

up to ﬁfth order accurate (“WENO5”) in space. For unsteady problems, the result-

ing nonlinear spatial discretization yields a set of ODE’s in time, which in turn is

solved via high order implicit time-stepping with error control. For the steady-state

case, we need to solve the non-linear system, typically by Newton-Krylov iterations.

There are several numerical examples presented to demonstrate the accuracy, non-

oscillatory nature and eﬃciency of these high order methods, in comparison with

other ﬁxed-stencil schemes.

1 Introduction

In the deterministic modeling of uncharged particle transport, physical

processes are described by the Boltzmann Transport Equation (BTE) which

is a linear integro-diﬀerential equation to be solved for the scalar unknown Ψ,

usually called the particle ﬂux. During the solution process, iterative meth-

ods are routinely used to solve the large systems of equations resulting from

various discretizations, especially when it comes to solving the steady-state

problem. As noted in [FaMa89], such approaches can be viewed as iterative

solutions of a matrix equation

Ax = b (1)

This work was performed under the auspices of the U.S. Department of Energy

by the Lawrence Livermore National Laboratory under contract W-7405-Eng-48.

402 B.L. Bihari and P.N. Brown

for the unknown vector x representing scalar ﬂuxes Ψ , energy densities, or

other integrated quantities of interest.

The spatial discretization has been traditionally done by some ﬁxed stencil

ﬁnite diﬀerence or ﬁnite element method typically solving for point values of

the solution, thereby resulting in a “node-centered” scheme.

Finite volume schemes, on the other hand, were invented for the numerical

solution of high speed ﬂuid dynamics where conservation is crucial. Since they

solve for the cell average of the unknown instead of its point values, they are

sometimes refered to as “cell-centered” schemes as well. Application of these

schemes to the BTE ensures not only conservation, but also provides for a

more natural handling of boundary conditions, material properties and source

terms, as well as an easier parallel implementation and post processing. The

ﬁnite volume scheme also lends itself to an eﬃcient implementation of high

order spatial discretizations.

Material interfaces and time-dependent large source terms can introduce

severe oscillations even with second order ﬁxed stencil schemes. Slope limit-

ing, or essentially nonoscillatory (ENO) spatial interpolations eliminate these

oscillations, and make higher-than-second-order spatial accuracies possible.

A newer variation of these nonlinear schemes is the Weighted ENO (WENO)

scheme that makes the stencil transition less abrupt and boosts the accuracy

in smooth regions. For unsteady problems, the resulting nonlinear spatial

discretization yields a set of ODE’s in time, which in turn is solved via high

order implicit time-stepping with error control. For the steady-state case, we

need to solve the non-linear system, typically by Newton-Krylov iterations.

Both of these approaches would ideally expect a preconditioner in order to

obtain a reasonable rate of convergence.

We will discuss the advantages of using an ENO/WENO method, as well

as the various issues introduced by such nonlinear methods originally designed

for computing shocked ﬂuid ﬂows. There will be several 1-D, 2-D, and 3-D

numerical examples presented to demonstrate the accuracy, non-oscillatory

nature and eﬃciency of these high order methods, in comparison with other

ﬁxed-stencil schemes.

The paper is organized as follows. In Sect. 2, we describe the continu-

ous problem from which the discrete linear system (1) is derived. In Sect. 3,

the exact form of (1) is established via the introduction of the discrete ordi-

nate angular discretization together with the nonlinear spatial discretization.

Although these discretizations are described elsewhere in the literature, we

present them in some detail for the sake of clarity and completeness. In

Sect. 4 we show several numerical results, and close with concluding remarks

in Sect. 5.

High Order Schemes for BTE 403

2 Background

In this section we ﬁrst introduce the linear time-dependent BTE in three di-

mensional box geometry with general scattering [Pom73]. The spatial domain

is the box D≡{r =(x, y, z)|a

≤ x ≤ b

≤ y ≤ b

, and a

≤ z ≤ b

the direction variable is Ω ∈S

, the unit sphere in R

, the energy variable is

E ∈ (0, ∞), time is t, and the equation in the ﬂux ψ = ψ(r, Ω, E, t)isgiven

v(E)

∂ψ

∂t

(r, Ω, E, t)+Ω ·∇ψ(r, Ω, E, t)+σ(r, E)ψ(r, Ω, E, t)



∞



ψ(r, Ω



,t)σ

(r, Ω · Ω



→ E)dΩ



+ q(r, Ω, E, t) ,

(2)

with initial condition

ψ(r, Ω, E, t

)=ψ

(r, Ω, E)(3)

where ∇ψ ≡ (∂ψ/∂x,∂ψ/∂y,∂ψ/∂z), v(E) is the particle speed, and ψ

the initial state at time t = t

The code employs a general algorithm that solves for multiple energy

groups via a semi-discretization using a ﬁnite number of energy “bins”. How-

ever, for simplicity of notation and clarity of exposition we shall assume a

single energy group E for the rest of this paper, and thereby remove the

energy-dependency from ψ, v, σ and q in (2).

The scattering integral on the right hand side of (2) is handled by ex-

panding the ﬂux ψ(r, Ω, t) in surface harmonics according to

ψ(r, Ω)=

∞



n=0



m=−n

(r, t)Y

(Ω) .

Here, Y

(Ω) is a surface harmonic deﬁned by

(Ω)=a

|m|

(ξ)τ

(ϕ) ,

where Ω =(µ, η, ξ) = (sin θ cos ϕ, sin θ sin ϕ, cos θ), P

|m|

is an associated

Legendre polynomial [Lib80], and

(ϕ)=



cos mϕ,ifm ≥ 0, and

sin |m|ϕ,ifm<0 .

The constants a

are deﬁned by



(2n +1)(n −|m|)!

2(1 + δ

)π(n + |m|)!



1/2

404 B.L. Bihari and P.N. Brown

where δ

n,n



is the Kronecker delta,and

(r) ≡



ψ(r, Ω)Y

(Ω)dΩ,

is the (n, m)

moment of ψ. Similarly, the source q is expanded as

q(r, Ω, t)=

∞



n=0



m=−n

(r, t)Y

(Ω) ,

where

(r, t) ≡



q(r, Ω, t)Y

(Ω)dΩ .

For ease of exposition in what follows, we have elected to use real-valued

surface harmonics, all scaled to have unit norm in L

Given ψ in the above form, one is able to rewrite the scattering integral

in the form



(r, Ω · Ω



)ψ(r, Ω



,t)dΩ



∞



n=0

s,n

(r)



m=−n

(r, t)Y

(Ω),

where the σ

s,n

are given by

s,n

(r) ≡ 2π



−1

(r, µ

(µ

)dµ

and where µ

is the cosine of the scattering angle. The total cross section σ

is given by

σ(r) ≡ σ

(r)+2π



−1

(r, µ

)dµ

= σ

(r)+σ

s,0

(r) ,

where σ

is the absorption cross section.

Substituting the appropriate terms into (2), the equation to be solved

now becomes:

∂ψ

∂t

(r, Ω, t)+Ω ·∇ψ(r, Ω, t)+σ(r)ψ(r, Ω, t)

∞



n=0

s,n

(r)



m=−n

(r, t)Y

(Ω)+q(r, Ω, t) ,

(4)

Boundary conditions must also be speciﬁed so as to make (4) well-posed.

Various options include a reﬂecting condition on a face, or a Dirichlet condi-

tion in which the incident ﬂux is speciﬁed on a face. For simplicity, we will

consider only the latter case. Namely, we will consider boundary conditions

of the form

ψ(r, Ω)=g(r, Ω) for all r ∈ ∂D and Ω ∈S

with n(r) · Ω<0 , (5)

where n(r) is the outward pointing unit normal at r ∈ ∂D.

High Order Schemes for BTE 405

3 Discretization of the 3-D Problem

We now turn to the angular, spatial and temporal discretizations of (4).

In previous work [AsBr90], we derived a matrix version of the well-known

diamond diﬀerence discretization scheme for the 1-D slab problem analogous

to (4)–(5). We extend that development here to 3-D problems. The ﬁrst

subsection deals with the quadrature rules for approximating integrals on

, the next describes the spatial discretization, and the ﬁnal subsection

considers the discrete-ordinates method, i.e., the combined discretization of

problem (4)–(5).

3.1 Quadrature Rules

The quadrature rules to approximate integrals on S

use the standard sym-

metry assumptions. Following Carlson and Lathrop [CaLa68], the quadrature

rules we consider are of the form



ψ(Ω)dΩ ≈



=1



ψ(Ω



) , (6)

where Ω



≡ (µ



,η



,ξ



), for all  =1,...,L, with L = ν(ν +2) and ν is the

number of direction cosines (ν =2, 4, 6,...). Since Ω



∈S

for all ,wehave



+ η



+ ξ



= 1 for all . (7)

With full symmetry, latitudinal point arrangement, and ν direction cosines,

(7) becomes

+ ξ

ν/2+2−i−j

for i =1, 2,...,ν/2andj =1, 2,...,ν/2 − i + 1. This last equation can be

solved to give

= ξ

+(i − 1)

2(1 − 3ξ

)

ν − 2

, (8)

for i =1,...,ν/2and0<ξ

≤ 1/3.

For the weights, ﬁrst note that a constant function on S

must be inte-

grated exactly for all ν, and so we must have

4π =



1 · dΩ =



=1



For ν = 2 or 4, by requiring that weights be invariant under 90

◦

rotations

of the Ω coordinate system, it is easily seen that each weight must be the

same. For ν>4thereareν/2 distinct weights, and for these sets one can

determine the weights by requiring that

406 B.L. Bihari and P.N. Brown



=1



4π

2n +1

, (9)

for n =0,...,ν/2 −1. These conditions guarantee that as many even powers

of ξ as possible are integrated exactly by the quadrature rule. Note that

due to the symmetrical placement of the ξ



along the ξ axis, all odd powers

of ξ are integrated exactly. For ν ≥ 22, demanding that (9) holds for all

n =0,...,ν/2 − 1 leads to negative weights. As an alternative, one could

simply demand that the weights are all equal, in which case (9) always holds

for n =0and1.

For either type of quadrature rule, we only require that (9) holds with

n =0and 1, and that all the weights are positive. Finally, it also follows

from the symmetrical placement of the direction cosines that the following

additional results hold:



=1



=0,



=1



=0,and



=1



=0. (10)

3.2 Finite Volume Spatial Discretization

In the formulation of the ﬁnite volume method, we ﬁrst discretize D into cells

(also called “zones”). Introduce the spatial grids

≡ x

< ···<x

i−

< ···<x

M+

≡ b

≡ y

< ···<y

j−

< ···<y

J+

≡ b

,and

≡ z

< ···<z

k−

< ···<z

K+

≡ b

and deﬁne r

ijk

=(x

). Next, deﬁne

∆x

= x

− x

i−

for i =1,...,M

∆y

= y

− y

j−

for j =1,...,J,and

∆z

= z

− z

k−

for k =1,...,K .

Also deﬁne ∆r

ijk

≡ ∆x

∆y

∆z

.The{r

ijk

} are referred to as nodes (or grid

points), and function values at these points are called nodal values. Assume

that σ and σ

s,n

have constant values on each cell

ijk

≡



r|x

i−

<x<x

j−

<y<y

k−

<z<z



denoted by σ

ijk

and σ

ijk

s,n

, respectively. Function values that are constant on

cells will be referred to as cell-centered values. We use ψ

ijk

to denote the

approximation to ψ(r

ijk

), the true solution at r

ijk

To obtain the ﬁnite volume method, average equation (4) over cell Z

ijk

We then have a semi-discrete form