Graziani F. (editor) Computational Methods in Transport

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

Regularized Scattering Operator 449

feature, but its accuracy is strictly limited to a small region around the pencil

beam axis. The error increases dramatically with increasing lateral distance

away from the axis and this deﬁciency in the solution can be traced to the

neglect of large angle scattering which is manifested through the neglect of

all momentum transfer moments higher than the ﬁrst. While this conclusion

suggests that greater accuracy may be achieved with the retention of addi-

tional terms in the expansion, more detailed theoretical analysis shows that in

fact all truncations beyond the ﬁrst term lead to unstable or unbounded solu-

tions, in then sense of transforms at least. A similar conclusion was reached

with the use of higher order Fokker-Planck expansions [LL01, PP01]. The

infeasibility of adopting higher order truncations eliminates the possibility

of systematically obtaining closed form generalizations of the Fermi solution

(however, see [PP98]). It is necessary to then consider approximations which

require numerical solution methods but prove more eﬃcient to solve than

the analog problem. The considerations of this section are still very relevant,

however, for the equivalence of the Boltzmann and diﬀerential forms of the

collision operator, namely (2) and (9), strongly suggests that approximate

models which correctly reproduce a certain number of the cross section an-

gular moments deﬁned in (11) are likely to be more accurate than those that

do not. This principle underlies the moment-preserving approach described

in the next section.

3 Regularized Collision Operator

Consider the following functional of the angular ﬂux

G[ψ]=A exp(βL) ψ − Aψ , (16)

where L is the 2D planar Laplacian deﬁned above and A and β are unspeciﬁed

constants. Expanding the exponential up to third order yields

G[ψ]=AβLψ +

Aβ

ψ + O(Aβ

) . (17)

Comparing (9) and (17), it follows that the J-andG-functionals can be made

formally equivalent through second order if A and β are chosen to satisfy

A =2

,β=

, (18)

so that we can write

J = G + O(σ

/32σ

+ σ

/288) . (19)

Thus, the two formulations of the scattering integral give identical ﬁrst and

second momentum transfer moments. The higher order terms will diﬀer, of

450 A.K. Prinja and B.C. Franke

course, but if scattering is suﬃciently forward peaked (σ

 σ

···)the

G-functional should represent an accurate approximation to the analog scat-

tering integral. In particular, under circumstances where the Fokker-Planck

approximation is accurate, the G-functional should also be accurate. More-

over, because the latter retains all higher order moments, although approxi-

mately, it has the potential of being more accurate than the Fokker-Planck

approximation. We now demonstrate the signiﬁcance of the G-functional for-

mulation. Consider the quantity

χ(s, η, ξ)=exp(s L) ψ(η, ξ) , (20)

where s is a continuous non-negative real variable. With this deﬁnition, (16)

can alternatively be written as

G[ψ]=Aχ(β,η,ξ) − Aψ . (21)

Now, diﬀerentiating (20) with respect to s it is not diﬃcult to show that

χ satisﬁes the following “time dependent, inﬁnite medium heat conduction”

equation

∂χ

∂s

∂

∂η

∂

∂ξ

,s≥ 0, −∞ <η,ξ<∞ , (22)

with auxiliary data

χ(0,η,ξ)=ψ(η, ξ), lim

|η|,|ξ|→∞

χ(s, η, ξ)=0. (23)

This initial-boundary value problem can be readily solved to obtain

χ(s, η, ξ)=

4πs



∞

−∞



∞

−∞

exp −



(η − η



)

+(ξ − ξ



)



ψ(η



,ξ



)dη



dξ



(24)

Finally, setting s = β in (24) and inserting the result into (21) yields the

following alternate representation of the G-functional

J[ψ] ≈ G[ψ]=

4πβ



∞

−∞



∞

−∞

exp −



(η − η



)

+(ξ − ξ



)



4β

ψ(η



,ξ



)dη



dξ



− Aψ(η,ξ) . (25)

For obvious reasons, we refer to the above representation of the scattering

integral as a moment-preserving approximation. It can form the basis for

generalizing Fermi-Eyges theory to account for the eﬀects of large angle scat-

tering on electron dose deposition, and we hope to report on this in a future

communication. Before we comment further on this result, however, we ﬁrst

note that (25) may be expressed in a form that is more practical for numer-

ical implementation by restoring the domain of the direction variable to the

unit sphere. This immediately follows upon using (3), getting

Regularized Scattering Operator 451

J[ψ]=

4πβ



4π

exp



−

(1 − µ

)

2β



ψ(Ω



)dΩ



− Aψ(Ω) . (26)

It should be noted that in so restricting the domain an error of O(exp(−1/β))

is introduced in the cross section moments. Although the coeﬃcients can be

rescaled to restore the exact moments, this is not necessary since β is typically

a small number. That is, the error is exponentially small and has no practical

consequence.

There are a number of attractive features about the approximate collision

integral obtained above and its generalization described below. First, we ob-

serve that the scattering operator has been regularized, i.e., the singularity at

forward directions has been eliminated, and there is a concomitant reduction

in the magnitude of the total scattering cross section. This smoothed approx-

imation of the transport process greatly improves prospects for application of

standard deterministic and Monte Carlo solution techniques. Also, the direct

use of existing neutral particle transport codes becomes a distinct possibility

since the algorithmic modiﬁcations necessary when diﬀerential Fokker-Planck

approximations are employed are not required with the present formulation.

Second, the approach presented here yields a univeral diﬀerential cross section

that depends only on cross section angular moments. That is, the ﬁnal result

is independent of the detailed form of the underlying analog cross section,

making it feasible to handle empirical or semi-empirical cross section repre-

senations with relative ease. Finally, we note that sampling deﬂection angles

from the exponential kernel in a Monte Carlo implementation is straight-

forward. This is in contrast to the inﬁnite spherical harmonics expansion

representation that results from the use of a generalized Fokker-Planck ex-

pansion [LL01].

Our approach can be readily extended to preserve higher moments of the

analog DCS by generalizing (16) as follows

G[ψ]=



n=1

exp(β

L) ψ −





n=1



ψ. (27)

The 2N free parameters {A

,β

,n=1, 2 ...N} can be selected by requiring

the consolidated generalized Fermi expansion of (27) to be identical to that of

the analog expansion up to the appropriate order, thereby ensuring that the

ﬁrst 2N momentum transfer moments of the DCS are preserved. This pro-

cedure will yield a nonlinear system of algebraic equations for the expansion

parameters which in general must be solved numerically. Retaining higher

order terms enables larger scattering angle eﬀects to be accommodated al-

though experience with other moment preserving strategies [FPKL1,FPKL2]

indicates that in practice two terms (four moments) is usually suﬃcient to

accurately capture angular distributions even for relatively thin targets. How-

ever, a more eﬀective way of incorporating large scattering angle eﬀects is to

ﬁrst extract from the analog diﬀerential cross section a component that is

452 A.K. Prinja and B.C. Franke

not forward peaked and is characterized by a long mean free path. By fur-

ther requiring the shape and amplitude of this smooth cross section to be

identical to or closely resemble the analog cross section, high order moments

can be accurately reproduced and large angle scattering eﬀectively captured.

The residual component will carry the singular part of the cross section and

is more suitable for the regularization procedure described above. We have

previously implemented this hybrid procedure with a discrete scattering an-

gle representation of the residual cross section and observed the approach to

be very accurate indeed [FPKL2].

4 Numerical Results

Monte Carlo simulations have been performed for 1 MeV electrons incident

on a gold target of thickness T =3.85 × 10

−3

cm [FPKL2], and we com-

pare transmitted and reﬂected angular distributions obtained using our new

regularized model, in pure and hybridized form, and our previous discrete

scattering angle representation [FPKL1] against corresponding benchmark

analog results. The screened Rutherford scattering cross section [FPKL1]

was used for sampling angular deﬂections in the analog model and for com-

puting the necessary momentum transfer moments. Energy dependence was

neglected in these simulations so that the electrons were either transmitted

or reﬂected from the inﬁnite slab.

In Fig. 1, the relative diﬀerences with analog results are shown for the 1-

and 2-discrete angle model and the exponential kernel model. The absence of

ray artifacts with our regularized model stands in stark contrast to the dis-

crete models. The exponential model is less accurate than the 2-discrete angle

model near grazing angles for reﬂection, but this can partly be attributed to

the discrete model preserving two additional moments. Because of the greater

expense of sampling from an exponential distribution, our new model is only

slightly faster than the 2-discrete angle method. The discrete model shows

ray eﬀects in the reﬂected angular distribution due to the presence of a dis-

crete scattering angle in the backward direction while the exponential model

yields a smooth distribution across all angles, as expected. Both methods

perform well in calculating the transmitted distribution, with the discrete

method performing better at large angles relative to the surface normal. The

2 discrete angle method accurately captures four angular scattering moments,

while the exponential model accurately captures only two moments. Because

of the relative size of the interaction cross section and the greater expense

of sampling from an exponential distribution, the exponential model is only

slightly faster than the 2 angle method.

In Fig. 2, the relative diﬀerences with analog results are shown for the

exponential model and the hybrid model. The hybrid model shows greatly

improved accuracy for both transmitted and reﬂected distributions, a con-

sequence of more accurately treating the higher order momentum transfer

Regularized Scattering Operator 453

Reflection Angle (Degrees)

(ψ

Analog

- ψ

Approx

)/ψ

Analog

0 30 60 90

-0.2

-0.1

0.1

1 Discrete Angle

2 Discrete Angles

Exponential

Transmission Angle (Degrees)

(ψ

Analog

- ψ

Approx

)/ψ

Analog

0 30 60 90

-0.5

-0.4

-0.3

-0.2

-0.1

0.1

0.2

1 Discrete Angle

2 Discrete Angles

Exponential

Fig. 1. Error in Angular Escape Distributions Relative to Outward Normal for

Electrons Incident on a Gold Slab for Discrete and Exponential Kernels

moments. However, this improved accuracy is achieved at a slightly increased

computational expense, with the runtime for the hybrid model approximately

two times greater than for the pure exponential model.

Conclusions

We have demonstrated that a closed form regularized scattering kernel for

electron transport, derived using an angular moment-preserving technique,

yields accurate escape distributions. Unlike our previous discrete scattering

454 A.K. Prinja and B.C. Franke

Transmission Angle (Degrees)

(ψ

Analog

- ψ

Approx

)/ψ

Analog

0 30 60 90

-0.2

-0.1

0.1

Exponential

Hybrid

Reflection Angle (Degrees)

(ψ

Analog

- ψ

Approx

)/ψ

Analog

0 30 60 90

-0.2

-0.1

0.1

Exponential

Hybrid

Fig. 2. Error in Angular Escape Distributions Relative to Outward Normal for

Electrons Incident on a Gold Slab for Exponential Kernels

angle representation, the present model is free of ray artifacts yet it is com-

putationally competitive. Moreover, the procedure described herein can be

readily and systematically generalized to incorporate more accurate physics,

withough increasing algorithmic complexity, by retaining successively higher

order angular moments of the analog diﬀerential cross section. Finally, nu-

merical testing shows that the best performance, from accuracy and compu-

tational eﬃciency considerations, is realized when the cross section is ﬁrst

decomposed into smooth, analog-like and singular components and the reg-

ularization procedure applied to the latter.

Regularized Scattering Operator 455

Acknowledgement

Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lock-

heed Martin Company, for the United States Department of Energy’s Na-

tional Nuclear Security Administration under Contract DE-AC04-94AL85000.

The work of the ﬁrst author (AKP) was supported in part by a grant from

the US Department of Energy under the award DE-FG07-02ID14336.

References

[JEM81] Morel, J.E.: Fokker-Planck Calculations Using Standard Discrete Ordi-

nates Codes. Nucl. Sci. Eng., 79, 340 (1981)

[CL83] Caro, M., Ligou, J.: Treatment of Scattering Anisotropy of Neutrons

Through the Boltzmann Fokker-Planck Equation. Nucl. Sci. Eng., 83,

242 (1983)

[GCP96] Pomraning, G.C.: Higher Order Fokker-Planck Operators. Nucl. Sci.

Eng., 124, 390 (1996)

[LL01] Leakeas, C.L., Larsen, E.W.: Generalized Fokker-Planck Approximations

of Particle Transport with Highly Forward-Peaked Scattering. Nucl. Sci.

Eng., 137, 236 (2001)

[PP01] Prinja, A.K., Pomraning, G.C.: A Generalized Fokker-Planck Model for

Transport of Collimated Beams. Nucl. Sci. Eng., 137, 227 (2001)

[PKH02] Prinja, A.K., Klein V.M., Hughes, H.G.: Moment Based Eﬀective Trans-

port Equations for Energy Straggling. Trans. Am. Nucl. Soc., 86, 204

(2002)

[FPKL1] Franke, B.C., Prinja, A.K., Kensek, R.P., Lorence, L.J.: Discrete

Scattering-Angle Model for Electron Pencil Beam Transport. Trans. Am.

Nucl. Soc., 86, 206 (2002)

[FPKL2] Franke, B.C., Prinja, A.K., Kensek, R.P., Lorence L.J.: Ray Eﬀect Mit-

igation for Electron Transport with Discrete Scattering-Angles. Trans.

Am. Nucl. Soc., 87, 133 (2002)

[RG41] Rossi, B., Greisen, K.: Cosmic Ray Theory. Rev. Modern Phys., 13, 240

(1941)

[PP98] Pomraning, G.C., Prinja, A.K.: The Pencil Beam Problem for Screened

Rutherford Scattering. Nucl. Sci. Eng., 130, 1 (1998)

[HMA81] Hogstrom, K.R., Mills, M.D., Almond, P.R.: Electron Beam Dose Calcu-

lations. Phys. Med. Biol., 26, 445 (1981)

Implicit Riemann Solvers for the P

Equations

Ryan McClarren

, James Paul Holloway

, Thomas Brunner

,and

Thomas Mehlhorn

Department of Nuclear Engineering and Radiological Sciences, College of

Engineering, University of Michigan, 2355 Bonisteel Boulevard, Ann Arbor, MI

48109 2104, USA

rmcclarr@umich.edu

Sandia National Laboratories

∗∗

, PO Box 5800, Albuquerque, NM 87185-1186

Summary. The spherical harmonics (P

) approximation to the transport equation

for time dependent problems has previously been treated using Riemann solvers and

explicit time integration. Here we present an implicit time integration method for

the P

equations using Riemann solvers. Both ﬁrst-order and high-resolution spatial

discretization schemes are detailed. One facet of the high-resolution scheme is that

a system of nonlinear equations must be solved at each time step. This nonlinearity

is the result of slope reconstruction techniques necessary to avoid the introduction

of artiﬁcal extrema in the numerical solution. Results are presented that show

auspicious agreement with analytical solutions using time steps well beyond the

CFL limit.

1 Introduction

Time dependent radiation transport is an important aspect in the problems

of radiation hydrodynamics and subcritical accelerator-driven nuclear sys-

tems. In both of these problems small time scales cause transport eﬀects

to be important in resolving macroscopic phenomena. In this article we

present the ﬁrst implicit, Roe-type Riemann solver for the spherical har-

monics (P

) approximation to the linear transport equation. In the past

explicit [BH2001a, BH2001b, Bru2000] and steady state [EPO2003] time in-

tegration methods were implemented and shown to be eﬀective in capturing

the propagation of sharp wavefronts. The most alluring characteristic of Rie-

mann solvers is the fact that they respect the characteristics of hyperbolic

equations by allowing information to travel only in the directions allowed by

the model. For ﬁrst-order spatial accuracy, the Riemann solver is linear for

the P

equations, however, higher order methods must use nonlinear slope

reconstruction. The implicit method we present uses the backward Euler

method and BDF2 for time integration. For the nonlinear equations of the

high resolution scheme a preconditioned inexact Newton method is used.

The computational results obtained show that it is not necessary to resolve

∗∗

Sandia is a multiprogram laboratory operated by Sandia Corporation, a

Lockheed Martin Company, for the United States Department of Energy’s National

Nuclear Security Administration under Contract DE-AC04-94AL85000.

458 R. McClarren et al.

the propagation of information across a spatial cell to produce numerical re-

sults tantamount to analytic solutions, that is, the CFL limit can be violated

and eﬀective results still ensue.

Equations

The one-dimensional spherical harmonic approximation of order n to the

linear Boltzmann equation is given by

∂

∂t

ψ(x, t)+A

∂

∂x

ψ(x, t)=−Sψ(x, t)+Q(x, t) . (1)

In the above, c is the particle speed, A is an n×n matrix, ψ is a vector of the

moments of the angular ﬂux, Q(x, t) is the inhomogeneous source term, and

S is a diagonal matrix that includes the collisional sources/sinks (namely

scattering and absorption). These equations are hyperbolic, meaning that

information propagates at ﬁnite speeds and that A has all real eigenvalues.

We seek to cast this equation in terms of cell averaged quantities. First we

prescribe a set of points, I =1, 2,...,i,...,N,whereeachi denotes a spatial

cell in our domain. Without loss of generality we assert that the spatial extent

of a cell is constant, ∆x, such that ∆x = X/N , where X is the length of our

spatial region of interest. Also, we shall designate the left and right edges of

a general spatial cell, i,byi − 1/2, and i +1/2, respectively.

Over a general cell, (1) is integrated and divided by ∆x (i.e. spatially

averaged). This yields

∂

∂t

ψψ

∆x



i+1/2

− F

i+1/2



= −Sψ

+ Q

(t) , (2)

wherewehavedeﬁned

(t)=

∆x



(i+1/2)∆x

(i−1/2)∆x

Q(x, t)dx , (3)

(t)=

∆x



(i+1/2)∆x

(i−1/2)∆x

ψ(x, t)dx , (4)

and

i±1/2

= Aψ



[i ± 1/2]∆x, t



. (5)

The F

i±1/2

terms are called the “ﬂuxes” across a cell interface. In the most

facile sense this is the ﬂow of information between cells. From (2) we see an

immediate issue, the quantities F

i±1/2

are not known in terms of cell-average

quantities. To express the F’s in these terms (and hence give us a numerical

method) we shall examine the solution to a model problem that exactly gives

the ﬂux between two cells–the Riemann problem.

Implicit Riemann Solvers for the P

Equations 459

3 Solving the Riemann Problem

The Riemann problem is to ﬁnd the ﬂux across a boundary, where the ﬂux

is given by a matrix, R, times a vector of state variables, u;thatistosay

F = Ru. The state variables are governed by the equation

∂u

∂t

+ R

∂u

∂x

=0, (6)

with the initial condition

u(x, 0) =



x<0

x>0

, (7)

where u

, u

are vectors of constants. Given the vectors l

, r

,thek

left

and right eigenvectors of R,andλ

,thek

eigenvalue of R, multiply (6) by

∂u

∂t

+ l

· R

∂u

∂x

∂a

∂t

+ λ

∂u

∂x

∂a

∂t

+ λ

∂a

∂x

=0, (8)

where a

= l

· u. Similarly the initial conditions are multiplied by l

(x, 0) =



ˆu

l,k

x<0

ˆu

r,k

x>0

, (9)

in the above ˆu

l,k

= l

· u

. The solution to (8) subject to IC (9) is

(x, t)=



ˆu

l,k

(x − λ

t) x −λ

t<0

ˆu

r,k

(x − λ

t) x − λ

t>0

. (10)

Next, the solution to (6) is reconstructed from the characteristic solutions

(10)

u =



(x, t)r



x−λ

t<0

ˆu

l,k

(x −λ

t)r



x−λ

t>0

ˆu

r,k

(x −λ

t)r

. (11)

At x =0

u(0,t)=



ˆu

l,k

(−λ

t)r



t<0

ˆu

r,k

(−λ

t)r

. (12)

The ﬂux at x =0is

F(0,t)=Ru(0,t)=



ˆu

l,k



t<0

ˆu

r,k



· u



· u

(13)