Finite-moments approaches to the time-dependent boltzmann equation

Progress in Nuclear Energy, Vol. 25, No. 2-3, pp. 127-157,1991.

0149-1970/91 $0-00+ .50 Copyright© 1991 PergamonPress plc

Printed in Great Britain. All rights reserved.

FINITE-MOMENTS APPROACHES TO THE TIME-DEPENDENT BOLTZMANN EQUATION* AHMED BADRUZZAMAN S a n d i a N a t i o n a l L a b o r a t o r i e s , R a d i a t i o n E f f e c t s D e p a r t m e n t , A l b u q u e r q u e , N M 8 7 1 8 5 - 5 8 0 0 , U.S.A.

ABSTRACT Two f i n i t e - m o m e n t s t e c h n i q u e s to solve the t i m e - d e p e n d e n t d i s c r e t e o r d i n a t e s t r a n s p o r t e q u a t i o n are discussed. with

a

finite-moments

expansion

finite-moments expansions algorithms

is

compared

o r d i n a t e s methods.

One of the t e c h n i q u e s is i m p l i c i t in time of

the

flux

in space,

in both time and space. to

that

of

conventional

and

certain

can

other

uses

The p e r f o r m a n c e of these time-dependent

discrete

The a c c u r a c y of the basic i m p l i c i t t i m e d i f f e r e n c i n g and

the e f f e c t of a n g u l a r d i s c r e t i z a t i o n are b r i e f l y discussed. techniques

the

offer

a significant

time-dependent

transport

advantage

problems,

over

Finite-moments

conventional

especially

methods

in

in m u l t i d i m e n s i o n a l

geometry. i. I N T R O D U C T I O N

Time-dependent with

pulsed

sources,

transport problems nuclear

bursts

are e n c o u n t e r e d

in

the

in o i l - w e l l

atmosphere,

reactor

logging

kinetics,

r a d i a t i v e transfer, etc. 1"4 As in s t e a d y - s t a t e problems, b o t h M o n t e Carlo and d e t e r m i n i s t i c t r a n s p o r t m e t h o d s have been u s e d to solve such problems. Carlo m e t h o d s can h a n d l e the time d e p e n d e n c e naturally. errors m a y r e m a i n u n a c c e p t a b l e ,

especially

Monte

However, s t a t i s t i c a l

in i n d i v i d u a l

t i m e bins.

Since

the e a r l y 1970s, two o n e - d i m e n s i o n a l d i s c r e t e - o r d i n a t e s codes h a v e e x i s t e d to solve t i m e - d e p e n d e n t t r a n s p o r t problems. 5°6 One uses a p o s i t i v e - f l u x w e i g h t e d finite-difference

algorithm

in both space and time. 5 This is s i m i l a r to the

scheme u s e d in the t w o - d i m e n s i o n a l the o t h e r t i m e - d e p e n d e n t

s t e a d y - s t a t e code,

discrete-ordinates

code,

DOT

(reference 7). In

the t e m p o r a l

*Work p e r f o r m e d u n d e r U S D O E C o n t r a c t No. D E - A C 0 4 - 7 6 D P 0 0 7 8 9 .

127

variable

is

128

A. B ADRUZZAMAN

t r e a t e d i m p l i c i t l y w i t h e i t h e r a linear d i a m o n d d i f f e r e n c e scheme to t r e a t the spatial variable. 6 Here, the

steady-state

effective

transport

source

cros s s e c t i o n interval.

term.

and 1/vAt,

The

equation

The

with

effective

an e f f e c t i v e

cross

section

w h e r e v is the p a r t i c l e

effective

source

is the

to the flux at the initial p o i n t

method,

a linear

with used

Pearce

treatment

can

in large

result

follow

time

the

steps

steps.

This

yielded

developed

and the The

discrete-ordinates

In this

dependence

expanded

in

terms

is of

stability

properties in a p r o b l e m

of

an

scheme

inability

a term

was

problem. 8

with that

a finitethe m e t h o d

of the m e t h o d

appropriately

accurate

to

constructed

results.

However,

for each p a r t i c l e

the

direction,

< i, w h e r e Ax is the s p a t i a l - c e l l inordinately

were

and

The implicit

derivatives,

transfer

concluded

that

more

performed

of the

similar

the

to

large

on

dependences

for the s t e a d y - s t a t e

treated

of

number

a nodal

of time

scheme

flux w e r e

where

expanded

finite-moments

transport

in

algorithm

equation, ~I and

it

time-dependent

implicitly,

and

the

The

second

i0.

algorithms

as

a short

assumption also

approaches.

polynomials.

a scattering

problem

will

finite-moments

reference

the

with

implicit We

two

Legendre

of

of the

the

is

examine

algorithm

analytically.

showed

source

m o r e a c c u r a t e r e s u l t s than c o n v e n t i o n a l

properties

accuracy

and

an

and At is the time

spatial

implicit

due to the

to

tests

method

we

temporal

edge

the

and

methods.

paper

finite-moments

of the

vAt/Ax

lead

spatial

in the m i d - 1 9 8 0 s

significantly

on at one

can

velocity

radiative

are restricted,

condition,

preliminary

m o m e n t s . I0

yielded

They

section

in the time step.

variable

primarily

schemes

restriction

Recently,

examined

spatial

considerably

stability

the t e m p o r a l

finite

errors

in e x p l i c i t

by the C o u r a n t

both

of the

characteristic.

schemes

width.

differencing

a two-dimensional

and M i t c h e l l 9 have

difference

explicit

diamond

in s o l v i n g

cross

is a sum of the actual

sum of the a c t u a l

proportional

recently

or a d i s c o n t i n u o u s

the r e s u l t a n t e q u a t i o n r e s e m b l e s

briefly

We and

medium pulse.

with

the

discuss

spatial

will then

is

spatial the

their

source then

of

discretization.

2. C O N V E N T I O N A L

The t i m e - d e p e n d e n t

TIME-DEPENDENT

Boltzmann

DISCRETE-ORDINATES

transport

equation

METHODS

is g i v e n

is

the error

is turned

discuss

variable

accuracy

the

space-time examine

examine

will

one,

variation the

briefly

a f t e r the We

In

by Iz

the

treated angular

Finite-moments approach

129

A

(z)

v 8t where

the variables

consists

have

of a scattering

their usual meanings

and the source,

term,

a n d a n e x t r a n e o u s s o u r c e term. A d i r e c t i o n , ~ , is d i s c r e t i z e d ,

In t h e d i s c r e t e - o r d i n a t e s such that the transport

a fission

scheme equation

term,

the particle is s o l v e d

for each member

S,

in g e n e r a l ,

of a discrete

set

{~m} "

2a.

Implicit

temporal

In a time such that

differencinq

interval

equation

A t = tk+1--tk, t h e

--

v

where

S" = S k÷1 f o r a f u l l y

explicitly.

temporal

derivative

is d i f f e r e n c e d

1 becomes

Equation

÷

At

~

implicit

• V~k+z+a~

scheme

2 can be written

A

k+1

=

(2)

S*

o r S* = S k if t h e s o u r c e

as a s t e a d y - s t a t e

• V ~ k*1 + a . ~ : ~ 1ffk÷l = S ~ : ~

is t r e a t e d

equation,

,

(3)

where 1 ~3 e f f

=

0

+

-----7 ChlV

(4 )

and Se~~ = S" +

i

@k

vAt

Note the

that time

reference

equation step

6, t h e

conventional discontinuous slab

geometry

3 is a s t e a d y - s t a t e

with

an

effective

spatial

linear

derivative

diamond

finite-element as

we

cross

consider

equation section

for

and

in e q u a t i o n

finite-difference algorithm. the

(5) "

former

the

an

flux

2 is t r e a t e d approximation

We will

confine

approach

at t h e

effective

here.

our

end

of

source.

In

by

either

the

or

the

linear

discussions

Then

to

equation

becomes

r

a__~[~_ (x,~) Ox

+ a~ze~k÷1(x,p)

= S~z z

(6) .

3

130

A. BADRUZZAMAN

Spatially

differencing

equation

6 as

k*l

(7)

~'eff

k+l

Ax

where

Ax = x]+ I - xj is t h e

cell

width,

=

where

the

quantity

distribution

at

a = 1/2, w h i l e of t h e

temporal

a is t h e

the

end

a~j+1

each

time

spatial

:

2b.

la.

Implicit

Space-time

temporal

9 resembles

differenced "cells" and

by

the

or " n o d e s " .

spatial

width, 1

vAt

Thus,

~

approach.

[~k+1 _ ~k]

+ ~

j+l

with

-~x

+ O~

for

the

flux

scheme,

the c o n n e c t i o n s

x

spatial

finite

differencing.

equation

equation

in s l a b g e o m e t r y ,

: S .

transport Here,

in a s p a c e - t i m e

Ax = Xj+l-Xj,

solve

linear-diamond

la d i s p l a y s

transport

a two-dimensional usual

relation,

alqorithms

the time-dependent

l a Ot ~+

Equation

can

the

Figure

differencing

v

auxiliary

points.

weiqhted-difference

L e t us c o n s i d e r

we

In

j j+} Fig.

the

,

parameter, step.

a = 1 .

mesh

assuming

(l-a)

+

weighting

of

in t h e s t e p s c h e m e and

and

we

(9)

equation

divide

and

space

n o d e of t i m e - i n t e r v a l , 9 takes

[<~9+i > _ <~j~]

the

hence

and

can

time

be

into

At = tk+1-tk,

form

+ 0<7>

= <~>

(I0)

Finite-moments approach

where

bars

averaging. boundary

denote

spatial

There

are

averaging

three

conditions).

and

unknowns

To s o l v e

<~

the

131

symbol

in e q u a t i o n

for these,

= a<~j÷1~

i0

we assume

(l-a)~j~

+

<

>

denotes

(given

the

temporal

initial

the auxiliary

and

relations

,

(ii)

and <~ in t h e

spatial

and temporal

a = d = 1 yields linear

diamond

points

in t h e

the

"cell"

spatial

if t h e

cell width

weighting

positive.

This

are shown

diamond than

coefficients results

analogous

are

chosen

such

following

above is t h e

connectivity

can

free

f r e e path, that

space-time of the m e s h

yield

time,

the

12 for

negative

[ =

i/ov,

"edge" and

=

"edge"

fluxes

-~-E~ken + _Ax _ 2 <0/j>(a + --~-{) for ~ < 0.

8 s = 1 and

(14)

13 a n d 8n = 0.

O

k ÷ ~1 k

O

!

j+l lb.

Space-time

finite

are

(13)

,

In e q u a t i o n s 5,

5,

for ~ > 0

VAt

In r e f e r e n c e

the

In r e f e r e n c e

05+ n-~--~x<~~ On + i ~k

k+l

Figure

equations,

in r e f e r e n c e

i = i/a.

relations

=

expressions

parameters.

The

scheme

a mean

than a mean

in t h e

In t h e

lb. As d i s c u s s e d

larger

1-a

adjustable

in F i g u r e

linear

is l a r g e r

(12)

a n d a = d = 1/2

scheme.

is

l-d

with

scheme

the

step

,

respectively.

Crank-Nicholson

case,

time

channels,

implicit-step

or t h e

steady-state

fluxes

the

the

= d~k+1 + ( l _ d ) ~ k

+1

X

differencing.

14, 8 s a n d

8 n are

132

A. B ADRUZZAMAN

3. F I N I T E - M O M E N T S

Two c a t e g o r i e s is

treated

of t h e s e a l g o r i t h m s

implicitly

and

the

arise.

spatial

expansion.

In

treated

equivalently

u s i n g finite m o m e n t s

Implicit

Let A M

other,

finite-moments

us r e w r i t e

the

the

In one, the t e m p o r a l

variation

moments

3a.

the

ALGORITHMS

is

temporal

treated

and

spatial

expansions.

by

variable a

finite-

variables

We d i s c u s s

are

both next.

scheme

implicit

difference

equation

6 for

a cell

of w i d t h

Xj+l--Xj,

=

2[/ Ax

so that -i ~ x ~ i. equation

15

can

over x to y i e l d

@~k÷z + aeff~ k÷: = Self ,

As d i s c u s s e d

be m u l t i p l i e d the m o m e n t s

2 • r , k÷l (--1) AxL,gj+=.

(15)

ax

in r e f e r e n c e

by

Legendre

balance

ii for the s t e a d y - s t a t e

polynomials

in x

and

case,

integrated

equation

nx÷l 1,i];+l--{(2/]X--1),b "k÷l +

}]

S elf

11rk+l

+ 0 effTnx

:

(16)

n . .

where 1

~

:

f ~ (x) P~xdx

(17)

-1

are

the

equation factor,

spatial 15

moments

exp(xa,f~/~).

polynomials

of

is i n t e g r a t e d

in space

the

over

The

flux.

To

x after

integrands

derive

being

the

auxiliary

multiplied

are then

expanded

relations,

by the

integrating

in t e r m s

of L e g e n d r e

to o b t a i n k+l : ~ ~j÷z

k÷l

exp(-OeffAx/~)

+ ~F

n xS nelf x ,

(18)

nx

where

1 Fn x _

Ax2~ 2nx+lexp(-aeffAx/2~)2

In m u l t i d i m e n s i o n a l othe r

spatial

-i f l~nx(x) e x p ( ~ x ) d x

geometry,

derivatives,

and

additional one

s c h e m e 11 to r e d u c e the m u l t i d i m e n s i o n a l one in each

spatial

channel.

thos e g i v e n

in r e f e r e n c e

would

t e r m s will a r i s e r e f l e c t i n g the use

equation

The d e t a i l s

(z9)

"

the

transverse

to s i n g l e - c h a n n e l

of the p r o c e d u r e

ii for the s t e a d y - s t a t e

transport

integration equations,

are identical equation,

to

except

Finite-moments approach

133

that h e r e we h a v e an e f f e c t i v e cross s e c t i o n g i v e n by e q u a t i o n 4 and h e n c e an e f f e c t i v e m e a n free path,

3b.

Space-time

leff =

finite-moments

i/aef f.

scheme

As in S e c t i o n 2b, here we d i v i d e space and time into "cells" or "nodes". Thus,

in a s p a c e - t i m e node of t i m e - i n t e r v a l ,

At, and spatial width,

Ax,

we

r e w r i t e e q u a t i o n 9, as

w h e r e -i ~ x,t ~ i. differential spatial

2

a~

vAt

3t

+

2~

o_~ + o ~

Ax

: s,

(20)

3x

To solve e q u a t i o n 20, we r e d u c e it to two s i n g l e - c h a n n e l

equations,

variable.

one for the t e m p o r a l v a r i a b l e and the other for the

For

example,

multiplying

equation

20

by

Legendre

p o l y n o m i a l s in x and i n t e g r a t i n g over the cell-width, we o b t a i n the t-channel differential exp(tavAt)

equation

which

is t h e n m u l t i p l i e d

by the

integrating

factor,

and i n t e g r a t e d from -i to t to o b t a i n

~n (t) = @ nkx e x p [ - o v ~ t ( l + t )

/2] +

VAt exp(-tovAt/2) 2 (21)

f exp(tlov~t/2) [Snx(t/)

- T x(r/,n x)] dt I ,

1

where 1

Tx(t,nx)

=

A2}~ X

f

~-~xPn x (x/) dx/ "

(22)

-1

We then e x p a n d the s o u r c e and the t r a n s v e r s e terms in e q u a t i o n 21 in L e g e n d r e polynomials

in t to o b t a i n ~nkx+I =

~:xeXp

(-(Jv A t) + E She [Sn .... - Tx(n:' nx) ] ' nt

(23)

where 1

Fn t _

v A2t

2n~+l •-ovAt, 2 exp~----~--) f P n

(t) exp(

t) dt

(24)

-i

and

Tx(nt,n x) = 2~[,i,J+i /%x[7n ~ + (-I) nx+l ~ nj~ - (2nx-l)~nc,nx-1----'] -

A s i m i l a r p r o c e d u r e y i e l d s the f o l l o w i n g e q u a t i o n ~j÷1 = ~Jtexp(_oAx/~) nt

+ Efnx[Sn~ nx

,nx

(25)

for the x - c h a n n e l

_ T~(nt,nx)]

,

(26)

134

A. BADRUZZAMAN

w h e r e Tt(nt,nx) channel

is the t - c h a n n e l

counterparts

We t h e n o b t a i n a m o m e n t s Legendre t hes e

polynomials

variables

equation

in

counterpart

Fnt in e q u a t i o n

of

of e q u a t i o n

balance equation

by m u l t i p l y i n g

in t and x and i n t e g r a t i n g the

space-time

node.

The

general

example,

solution

algorithms

by s e t t i n g n t = n x = 0 in e q u a t i o n

n x = 0,i

in e q u a t i o n s

moments,

we recover

23-27

the

and

4.

of

the

over

balance

time-dependent emerged

4a.

Implicit

The

transport

the

only

the

linear

of the s t a b i l i t y

Here

we

review

finite-difference

implicit

scheme

and

zeroth-order

characteristics

the

in e q u a t i o n

8, w i t h a = 1 and a = 1/2,

linear

diamond

algorithms,

conclusions

the w e i g h t i n g

parameter,

a =

analyze

the

a,

is g i v e n

1 l_exp(_oe:~Ax/~)

stability

properties

7 and 8 can be r e w r i t t e n ~ kj+i/2-~'j÷i/2 +l ,I,k

7 and

the

represent

of

that

auxiliary

the implicit

respectively.

e q u a t i o n s w i l l also r e p r e s e n t the i m p l i c i t z e r o t h - o r d e r

equations

only the zeroth-

By s e t t i n g n t = 0,i and

algorithms

g i v e n by e q u a t i o n

and

scheme.

For

to date.

implicit

relation

from the a b o v e equations.

PROPERTIES

analysis

algorithms.

(27)

•

scheme.

STABILITY

We h a v e b e g u n a s y s t e m a t i c

have

retaining

linear-moments

S n ....

=

27 and r e t a i n i n g

orde r m o m e n t s we r e c o v e r the C r a n k - N i c h o l s o n

To

form

20 by

equation

is

We can d e r i v e v a r i o u s

where

equation

the r e s u l t i n g

Tt(nt,n x) + Tx(nt,n x) + o ~ ....

step

25, and Fn. is the x-

24.

These

two

finite m o m e n t s scheme

by -

~

(28)

o,ffAx

of

these

algorithms

we

note

that

as

+ P [ ~ j +k+l I - ~k+l]

+ Hu/S÷I/2N'I'k~I = pC~)9÷i/2~ --k~1

(29)

and ~;jk+i/2 = a~gk+1 + ( l - a ) ~ j

where

k

,

(30)

Finite-moments approach

135

p = ~vAt/Ax

(31)

= ovAt

(32)

and c is the s c a t t e r i n g ratio, as/a.

= At/T

In e q u a t i o n 32, r is the m e a n free time

g i v e n by (33)

= I/av.

In the a b o v e e q u a t i o n s we u s e d ~ i / 2 f o r

~k.

The t e r m on the r i g h t - h a n d side

of e q u a t i o n 29 r e p r e s e n t s the s c a t t e r i n g source t e r m w i t h the scalar flux M

~.w.

where

w m are

the

discrete

ordinates

,

(34)

quadrature

weights.

The

stability

a n a l y s i s is p e r f o r m e d by u s i n g the F o u r i e r t e c h n i q u e 13 w h e r e we let the errors in ~ and ~ be e~ = A ( ~ m , ~ ) ~ k e i~je~

(35)

e~ = ~ keilJe~

(36)

81 = ~A____~x , Lx

(37)

and

with

where

L x is the

width

of

the

slab,

and

note

that

the

errors

also

satisfy

e q u a t i o n s 29 and 30, we o b t a i n the f o l l o w i n g e x p r e s s i o n for the error growth factor,

i~i~

=

(3sa)

i

(I+~-~C) 2 + 2 k ~ ( l - C O S Q 8 I)

where vAt A(~ m

,~) w m .

(38b)

136

A. B ADRUZZAMAN

Note

that

the maximum

value

of

I~l ~ o c c u r s

J~IL

We note

the

i.

following

properties

2.

I~I~ decreases

3.

I~I2max a p p r o a c h e s

We should

are

We

also

moments

scattering

sources.

Space-time

quantity:

as c ~

in the

diamond

factor

29 and 30

be

nature

evaluated.

scheme.

source

We

Numerical

show

that

the

stable.

the

implicit for

the iterative

is b e i n g

iteration

scattering

evaluated

stable

ignored

sources

of t h e

30 w i l l

technique

growth

radius

29 a n d

the algorithm

Fourier

done above of s u c h

iterative

be

by e q u a t i o n s

i.

an

retained to

represented

increases

spectral

linear

Analyzing

error

~

numerically

are

appears

(Bg)

stable.

Inclusion

with

in e q u a t i o n s

algorithm

the

the

problems

have

spatial

by t h e

of t h i s

unity

source.

examining

on

algorithms

4b.

1

note that the analysis

scattering

also

tests

as

= 1 so t h a t

(1+~_~c) 2

l~12~x ~ 1 so t h a t t h e a l g o r i t h m s are unconditionally

of t h e

-

at cos£81

stability

when

finite-moments

both

iterative

higher-order

algorithm.

and

The

non-iterative

scheme

given

discussed

by e q u a t i o n s above,

we

10-12,

find

for a = 1/2,

that

the maximum

d = 1/2, value

of

is

[ 1 - ~ (1-c)] 2 l~l~

--

~

(40)

[ I + ~ 2 (l-c) ] 2 Note

that

l~12maxr e a c h e s

increases problems 4c.

for of

Space-time

equations lowest

a

minimum

Also,

interest,

For t h e

the

~ > 2.

finite

the

value

moments

constant-flux

23-27, order

of

zero

it a p p r o a c h e s

the Fourier

of

~ can

at

~ =

unity

2

when

as c ~ i.

be m u c h

larger

c

Note than

=

0

that

and

then

in m a n y

2.

alqorithms

assumption analysis

(~ = 0) h a r m o n i c

in b o t h shows

is g i v e n

time

that by

and

space

the error

in the m o m e n t s

growth

factor

for

Finite-moments approach

137

I~l : p e " + c [ l - e ~ - p e " ] + c[1-ePWe note

the

following

i.

I~l ~

2.

I{I

=

3.

I{I

~

We

are

examining moments

the

moments

e-P 1

for

as c "

for

1.

of

and

the

stability

harmonics.

depends

moments.

time)

with

no

In a d d i t i o n ,

we are

space-time

finite

higher-order

properties

in g e n e r a l ,

higher-order

(in s p a c e the

the

factor,

and

the higher

properties

Here,

growth

zerothscheme

expression

c = 0.

evaluating

stability

algorithms.

expression,

of t h i s

1 in g e n e r a l .

currently

the

since the error both

properties

(41)

~]

become

more

on t h e e r r o r

For

example,

sources

we

for

obtain

complicated

amplitudes the

the

of

linear-

following

~ = 0 harmonic,

--~[

sn Vx ~] + (Bnl2vr vx

~+Ana

= 1 -

~An) U~ d

de

(42)

Bn Vt An d

where 1 l - e x p (-ovA t)

d =

1 avA t

(43)

1 - 2yt 1 - y t/d

(44)

(Jv A t 2 Y r/(1-Y J d )

(45)

Ut -

with Vr =

Y t =

_ _3

20

~

_L

~ oAx

[i - expl - o ~ 1 1 . _ L . 2 ~ ~& ]] v A t

a~ = 2 d The

quantities

respectively. order

moments.

a

and

A n and Figure

At/T for d i f f e r e n t

Vx

B n are

are the

2 displays

values

the error the

of t h e ratio,

,

(46)

1 .

(47)

x-channel

counterpart

amplitudes behavior B./A,.

of of

the

I~I

Numerical

of

zeroth-

d

and

and

as a f u n c t i o n tests with

Vt,

firstof

~ =

iterative

138

A. BADRUZZAMAN

sources

show

satisfies

that

the

constant

"edge"

exhibits

the

Courant flux

superior

linear-moments

condition, assumption,

stability

algorithm

r = vAt/Ax with

a

will

i.

We

be also

higher-order

stable found

interior

if

one

that

the

expansion,

properties.

3,2

i

~//B 2,8

<

i

i

n / A n = 1.50

t

Bn : LINEAR M O M E N T ERROR ~'J' 2.4

nO I-

AMPLITUDE

I

An : Z E R O T H M O M E N T ERROR

AMPLITUDE

I

O

2.0

|

< IJ. "I-

~ ,

~

B n / A n = 1.00

1.6

0 ¢~

1.2

%% I

_ B n / A n = 0.00

0 nw

0.8

0.4

--

'

% s

. ...........

"%s, S %.

~ /

I 0.4

0.0 0.0

-

I 0.8

I 1.2

1.6

At/~

Fig.

Error growth algorithm.

2.

factor

5.

Here

we

test

slab 20.5 mean free time

distribution scattering quadrature assumed. spatial

scheme

right

(mfps)

is

of

used

=

analytical

as t h e

finite-moments

vacuum.

an interval

and

a

Os/a

=

is t u r n e d

pointwise

and

a

flux

for the

ratio We

r

then

other

We

then

of 2 0 . 5 mft.

0.2,

the space-time

m e t h o d . 13014

reference

the

A source

algorithms

0.5

and

observe Three

0.9,

convergence

criteria scheme,

v~t/Ax

use

this

algorithms.

=

1/2,

with

fine-mesh

which the

values

were

linear moments =

in

a

o n f o r a 0.5 m e a n

0.5 m f t w i d e a t t h e l e f t b o u n d a r y

medium,

moments

EXAMPLE

is

the

I/4)

space-time

in a c e l l

We first compare (Ax

and

wide.

linear

boundary

in t h e s l a b a f t e r ratio

mesh

collision

interval

the

source-free

A NUMERICAL

conventional

free paths

(mft)

reflective;

the

for the

is

flux of t h e

used.

An

S8

of

10 .5 w a s

using a fine a

multiple-

linear

moments

Finite-moments approach

In linear with The

the

flux

latter

source a

Figure

3 the

moments

of

distribution

using

distribution

uses

of

width

from

turned

nearly

14%

Consequently,

the

the

as/a

the

Since

from

r =

fluxes

calculated

in t h e f l u x s h a p e w a s s e e n f o r as/G = 0.2 a n d 0.9.

fine-mesh

the

conventional

space-time

moments

and

results

finite

as t h e

shape

is

moments

m e t h o d . 14

scheme

used

there by

in t h e s o u r c e

agreement

of

compared

interval,

The

accuracy

flux

space-time is

analytical

finite

normalized

in t h e

the

1/2),

the moments

but

pointwise

fluxes were

agreement

0.5

i/4,

multiple-collision

s o u r c e . 15

on for a small in

=

(Ax =

problem.

the

excellent

for

a fine-mesh

a delta-function

finite

difference

methods.

flux

method,

139

the

a

was two

c e l l of t h e

apparent.

Similar

We next discuss

algorithms

using

the

reference.

l@S O - - - O A N A L Y T I C A L METHOD H SPACE-TIME LINEAR MOMENTS METHOD (z~x = ~./4, r = 1/2)

10-6

10-7

X

..J 14. 10-6

10-9

10-1o 1

3

5

7

9

11

13 15

DISTANCE Fig.

3.

Flux distribution analytical method as/o = 0.5.

5a.

ImPlicit

Figure

linear

4

linear diamond spatial that

mesh

for

Ax

significantly with different

diamond

displays scheme,

widths, =

I/4

the

and

flux

various r

(x/k)

the space-time linear moments 20.5-mfp-wide slab at time 20.5

and an mft with

alqorithm

discussed and

underestimates values

by in a

17 19 21

= the

distribution in S e c t i o n values

1/2,

the

flux.

of r s h o w t h a t

2a,

of

obtained

ratio,

implicit

The

the

f o r as/a = 0.5 u s i n g

the

next

a t r = 1/2,

with

r = vAt/Ax.

linear

three

diamond

curves

the flux

for

implicit various We

see

scheme A x = I/2

is u n d e r e s t i m a t e d ,

140

and

A. BADRUZZAMAN

then

diamond

as

r

scheme

is

increased

grossly

midpoint

of

presents

the error distribution.

scheme

the

the

performs

For example, from nearly

system

and

even worse

22% t o 550%,

F o r as/a = 0.9,

still

very

shape

changes

while

it

beyond

F o r as/a = 0.2,

r = i,

from the

Figure

that

that

more

Table

1

linear diamond

flux distribution.

with

Os/a = 0.5 r a n g e

is s o m e w h a t

implicit

point.

the implicit

errors

the

less than near the

5 for the

absolute

those with

the method

such

the flux at distances

underestimates

as s e e n

f o r A x = I/2,

81%.

flux

overestimates

as/a = 0.2

from nearly

accurate

range 10% t o

but the errors

are

large. 10-5

~" . ~ < . ~ Z L x = )../2, r = 2 ~'~-.,k ~ , x - -

= )J2, r = 1

10-6

r=1/2/ N~ Zix = k / 2 . r = 1/2

10-7

~

X

::3 -J LL 10-8

o',/o" = 0.5 TIME = 20.5 mft 10 -g

IMPLICIT LINEAR DIAMOND

~

e.--~l, SPACE-TIME LINEAR MOMENTS METHOD = ~./4, r = 1/2 (REFERENCE) 10-10

1.5

5.5

g.5

13.5

DISTANCE

4.

5b.

Space-time

Figure

weighted

flux

shape

somewhat

positive-flux

6 displays

flux

21.5

(X/Z)

shifts

as

weiqhted

difference

the flux distribution

difference

scheme

the

ratio

discussed r

is

are larger.

as/a = 0.2,

distribution.

For example,

range

from

in a 2 0 . 5 -

f o r as/a = 0.5 f r o m t h e p o s i t i v e in

Section

2b.

However,

linear diamond

Figure

scheme

scheme

increased.

less than that with the implicit

2 for the error

errors with

17.5

Flux distribution by the implicit linear diamond m f p - w i d e s l a b a t t i m e 2 0 . 5 m f t a n d as/a = 0.5.

Fig.

Table

'~

7 shows

We

see the

scheme,

that

that effect

42%

to

215%,

while

f o r as/G = 0.2,

those

is

as s e e n f r o m

f o r A x = I/4 a n d r = I, t h e a b s o l u t e

nearly

the

with

the

errors,

Gs/a = 0.5

Finite-moments approach

range much

from less

Table

nearly for i.

o,/o

16% =

to approximately

42%.

Error in scalar flux by (Os/O = 0 . 5 ; T i m e = 2 0 . 5

Ax r =

errors

(not shown

here)

were

0.9. the implicit mft). ERROR

X/l

The

141

1/2

=

linear

AX

1

r =

scheme

(%)

I/2

r =

diamond

2

r =

1/2

=

r =

1

r =

2

0.25

-14.44

49.30

274.90

49.27

274.92

1339.93

1.5

-14.94

47.52

266.73

47.20

265.63

1378.60

3.5

-16.97

39.90

226.01

40.34

228.12

971.90

5.5

-20.75

26.84

166.71

27.60

168.86

609.10

7.5

-25.86

i0.ii

97.04

11.23

98.43

259.97

9.5

-32.46

-9.85

32.03

-10.48

27.19

52.18

11.5

-40.65

-29.79

-18.45

-34.87

-30.82

-44.02

13.5

-49.31

-48.81

-52.61

-58.57

-72.69

-77.57

15.5

-60.36

-65.58

-72.53

-85.70

-90.58

-91.44

17.5

-69.24

-74.99

-80.39

*

-90.40

*

19.5

-80.91

-81.21

-82.56

*

-86.47

*

* FLU:

WAS

NEGATIVE

10 .7

Ax=)J4 r = 1

10 -s

~x:AJ4 r = 1 / 2 " * ' ~ A , %

10"9

x

10-1o (~,/G = 0.2 TIME = 20.5 mft

I M P L I C f f LINEAR DIAMOND METHOD

10-~1 O-O

8 P A C E - T I M E LINEAR MOMENT8 METHOD ~ x = L ' 4 r = 1/2 (REFERENCE)

10-12

1.5

5.5

g.5

13.5

17.5

21.5

DISTANCE (x/k) Fig.

5.

Flux distribution mfp-wide slab at

by the implicit linear diamond t i m e 2 0 . 5 m f t a n d as/a = 0 . 2 .

scheme

in a 2 0 . 5 -

142

A. BADRUZZAMAN

10-5

2 .,r=l

10-6

= ~./4, r = 1/2

10-7

10-8 o',/o" = 0.5 TIME = 20.5 mft

10-9

[D---O POSITIVE-FLUX WEIGHTED DIFFERENCE METHOD

L~

e--.-e SPACE-TIME LINEAR MOMENTS METHOD Ax = ~/4, r = 1/2 (REFERENCE) 10-10 1.5

5.5

9.5

DISTANCE

Fig.

6.

Table

Flux distribution in a 20.5-mfp-wide

2.

17.5

Error in scalar flux by the positive-flux scheme (Gs/a = 0 . 5 ; T i m e = 2 0 . 5 m f t ) .

~x r =

1/2

21.5

by the positive-flux weighted difference s l a b a t t i m e 2 0 . 5 m f t a n d as/a = 0 . 5 .

ERROR

x/l

13.5

(x/~,)

=

i/4

r =

1

weighted-difference

(% ) ~x

r =

2

scheme

r =

1/2

=

r =

i/2 1

r = 2

0.25

19.06

39.42

62.50

17.44

25

23

4.11

1.5

19.44

41.83

70.96

16.94

25

65

5.56

3.5

18.97

41.74

74.45

17.08

28

23

11.05

5.5

18.73

40.66

75.64

18.97

33

25

20.74

7.5

19.15

37.53

72.48

23.12

39

54

33.70

9.5

20.78

34.85

60.54

30.59

45

40

47.91

11.5

24.88

29.60

40.07

43.93

49

65

51.10

13.5

33.22

22.77

12.94

67.53

51

27

26.13

15.5

48.47

12.82

-16.07

111.18

50

26

-16.69

17 • 5

99.93

13.40

-31.26

241.83

70

81

-41.41

19.5

207.15

15.63

-38.61

585.76

115

84

-54.36

Finite-momentsapproach

143

10 ~ Ax=~J4 r = 1

. Ax=;~J4 r = 1 / 2

10-9

10-1o

T I M E = 20.5 rnft

POSITIVE- F L U X W E I G H T E D DIFFERENCE METHOD C--O

SPACE-TIME LINEAR MOM ENT8 METHOD Ax=;d4 r = 1/2 (REFERENCE)

10-11

1,5

5.5

9.5

13.5

17.5

21.5

DISTANCE (x/k)

F l u x d i s t r i b u t i o n by the p o s i t i v e - f l u x w e i g h t e d d i f f e r e n c e in a 2 0 . 5 - m f p - w i d e slab at time 20.5 m f t and os/a = 0.2.

Fig.

7.

5c.

Space-time

Table Section

3 displays

2b,

exhibited when

linear

for

diamond

the

varying

negative

error

fluxes.

While

and as one m o v e s

Ax

flux

l,

the

in

distribution

widths

the

and

the

the errors

made coarser =

scheme

mesh

Ax = i/4 and r = 1/2,

scheme

for this

r =

1/2.

agreement increase

two

points

r ~

is e x c e l l e n t

near

right

in

the

method

at all

points

as the mesh

(reflecting) the

discussed i,

significantly

away from the left

last

method,

For

boundary.

boundary

is For

becomes

negative.

5d.

Implicit-in-time

This T able the

left

implicit are

algorithm

4 displays

error

Also,

this m e t h o d

by s e t t i n g

distribution the

errors

scheme,

while

near

the

flux

shape

n x = 0,i

are

as

in e q u a t i o n s

method. similar

the r i g h t

changes

u n l i k e the i m p l i c i t

for Ax = i are p o s i t i v e

alqorithm

for this

boundary,

diamond

However,

spatial-moments

is o b t a i n e d

the

(reflecting) linear

smaller.

changed.

linear

the

We to

at all points.

those

(vacuum) ratio,

linear d i a m o n d scheme,

16 and

see that from

boundary

r

=

vAt/Ax,

18. near the they is

the fluxes from

144

A. BADRUZZAMAN

Table

3.

Error in scalar flux by (as/o = 0.5; Time = 20.5

the space-time mft).

ERROR Ax x/A

r

0.25

= =

~/4

Ax

1/2

-4.12

x

1 0 "I 1 0 "I

=

-1.18

diamond

scheme

(%)

=

r

linear

~/2

Ax

1/2 x

r

=

= 1/2

i00

-5.39

x

i00 i0 °

1.5

-3.32

x

-1.13

x

i00

-5.73

x

3.5

-2.53

x

1 0 -I

-9.32

x

1 0 -I

-5.03

x

I0 °

5.5

-1.35

x

1 0 -I

-6.03

x

10 I

-3.81

x

i00

7.5

-1.05

x

1 0 .2

-2.44

x

1 0 -I

-2.55

x

i0 °

9.5

-1.25

x

1 0 "I

-4.72

x

1 0 .2

-1.67

x

i00

11.5

-2.30

x

10 I

-7.90

x

1 0 .2

-2.99

x

i0 °

13.5

-1.58

x

1 0 -I

-2.33

x

i0 °

-7.57

x

i0 °

15.5

-5.27

x

1 0 -I

-6.02

x

i0 °

-8.50

x

i0 °

17.5

-4.68

x

1 0 -I

-1.06

x

101

19.5

-8.88

x

10 I

-6.02

x

101

* FLUX

Table

WAS

NEGATIVE

4.

Error (as/a

in scalar flux by the implicit = 0.5; Time = 20.5 mft).

ERROR AX X/l 0.25

r

=

1/2

=

~/2

r

=

1

r

=

2

linear

moments

algorithm

(%)

r

=

1/2

AX

=

1

r

=

1

r

=

2

-14.29

49.57

275.72

49.62

275.89

1344.08

1.5

-14.82

47.76

267.47

47.75

267.35

1283.43

3.5

-17.01

39.76

225.47

39.79

225.66

964.56

26.46

165.36

26.47

165.38

586.62

5.5

-20.93

7.5

-26.14

9.59

96.13

9.61

96.13

274.85

9.5

-32.74

-9.89

33.03

-9.91

32.80

86.08

11.5

-40.13

-28.69

-15.12

-28.87

-15.57

-9.85

13.5

-47.62

-45.82

-47.30

-46.24

-47.93

-54.26

15.5

-55.99

-59.96

-66.19

-60.78

-66.91

-73.89

17.5

-59.33

-65.48

-72.31

-66.96

-73.22

-79.31

19.5

-61.36

-65.60

-70.78

-68.26

-72.11

-77.18

Finite-moments approach

5e.

Space-time

This

ratio

constant-edqe-constant-interior

scheme

5 displays

the

r =

is o b t a i n e d error

At/Ax to

order

magnitude

diamond

the

and

positive

Table

and

than

stability

see

also

that

for

moments for

r

the

the

<

1 and fn

for

r ~

from in

the

general,

method

yields

thus

confirming

C-C

moments

algorithm

(%)

=

~/2

r =

1

r =

2

10 "I

4.01

x

i00

2.10

x

1.5

-5.61

x

10 I

3.52

x

i0 °

2

ii

x

101

3.5

1.59

x

10 .2

4.08

x

i0 °

2 00

x

101

5.5

1.17

x

i0 °

5.18

x

i0 °

2 96

x

101

7.5

2.80

x

i00

6.69

x

I0 °

1 16

x

101

9.5

4.99

x

i0 °

8.65

x

i0 °

3 02

x

101

11.5

7.99

x

i0 °

1.12

x

101

2 36

x

101

13.5

1.13

x

101

1.35

x

101

2 41

x

101

15.5

1.44

x

101

1.54

x

101

1.88

x

101

17.5

2.32

x

101

1.82

x

101

3.60

x

i0 °

19.5

1.72

x

i0 °

8.53

x

101

-3.62

x

101

0,i

This

moments. the

However,

also

interior

is

since

obtained

source

and

6 displays

here

when

note I.

the

Table

errors

better

constant-edqe-linear-interior

scheme in

are

the it

that

With

r

represents

of

the

=

the

errors 2

slab.

(not A

equations

to

the not

shown similar

from

was

interior grow in

~C-L)

and

error

increased, flux

behavior

nt = ~x

the

shape

table)

setting

the was

r

C-L

is

errors seen

in

= I/2.

We

see

that

performed

From

Table

increased increased Ax

and

"edge"

algorithm.

scheme

for

0,i

the

C-C

better. as

nt =

0

space-time

appreciably

the

by nx = for

the

101

alqorithm

23-27

distribution

those

ratio

do

moments,

error

similar

scattering

the

from interior

the

4.

x

nx =

an

linear

-7.09

Space-time

the left

implicit the

i,

Table as

0.25

5f.

to

are,

Also,

the space-time mft).

1/2

moves

errors

Section

Ax r =

one

23-27.

increases

corresponding

ERROR

x/l

alqorithm

error

as

schemes.

analysis

E r r o r in s c a l a r f l u x b y (os/a = 0 . 5 ; T i m e = 2 0 . 5

the

increases

However,

those

linear

fC-C)

n t = n x = 0 in equations

We

it

solutions

the

setting

boundary.

less

stable of

5.

right

implicit

and

conclusions

by

distribution.

increases,

boundary of

145

=

6 we

from

1/2

in

the

i/4.

146

A. BADRUZZAMAN

Table

6.

Error (as/a

in scalar flux by the space-time = 0.5; Time = 20.5 mft).

ERROR Ax x/l

r =

=

C-L

moments

(%)

i/2

Ax

1/2

r =

algorithm

1

r =

=

1/2

r =

1

0 25

-9.22

x

10 "I

-1.40

x

i0 °

-3.92

x

i00

-5

80

x

i00

1 5

-8.31

x

i0 "I

-1.25

x

i0 °

-3.29

x

i00

-4

79

x

i00

3 5

-4.81

x

10 -I

-6.86

x

10 I

-1.92

x

i00

-2

85

x

i00

5 5

4.83

x

10 -I

3.91

x

10 -I

7.33

x

10 -I

1 13

x

I00

7 5

1.19

x

i0 °

1.83

x

i00

4.60

x

i00

6 88

x

i00

9 5

2.44

x

i00

3.61

x

i0 °

9.87

x

i0 °

1 45

x

101

ii

5

4.11

x

i0 °

5.78

x

i0 °

1.70

x

101

2 42

x

101

13

5

5.46

x

i00

7.17

x

i00

2.51

x

101

3 25

x

101

15

5

6.46

x

I00

7.81

x

i00

3.11

x

101

2 61

x

101

17

5

i.i0

x

101

5.80

x

i00

4.13

x

101

6 40

x

101

19

5

-1.46

x

101

6.32

x

i00

-3.79

x

101

2 89

x

101

(L-L)

alqorithm

5g.

Space-time

This

linear-edqe-linear-interior

scheme

follows

in t h e

edge

and

7 that

this

algorithm

more the

accurate results

from

the

from

more

Nicholson) the

in

significantly time

steps

method C-L

and than

with

Section

especially

other

accurate

linear

(of

a

diamond 2b).

the

interior accurate

still

of

23-27

the

results

which

one

C-C

algorithms

and

to

from

mesh, the

are For

are

those

two

and

a

positive-flux errors

of

the

slab.

the

other

the

Courant

grew

We

of

Ax

finer

weighted

more

an

condition,

can r <

1/2, those

order

use i.

than

difference

when the

of

(Crank-

accurate

finite

since one

r =

than

diamond

considerably

However,

Table

magnitude

= l,

mesh,

linear

n x = 0,i

from

accurate

orders

schemes,

and

see

orders with

more

space-time

the

satisfying

several

orders

several

n t = 0,i

node.

example,

with

the

Here,

than

by setting

space-time

algorithms.

finer

more

while

equations

moments

yields

the

this

scheme

implicit

scheme

than

space-time

magnitude

from

interior

r

=

method far

i, is

fewer

147

Finite-moments approach

Table

E r r o r in s c a l a r f l u x b y t h e s p a c e - t i m e (Us/a = 0.5; T i m e = 2 0 . 5 m f t ) .

7.

ERROR Ax

xlX

5h.

r =

10 .3

4.54

x

I0 "z

x 10 .2

3.5

2.39

x

10 .3

-2.39

5.5

1.01

x 10 .3

-1.65

x

10 .2

7.5

-2.60

x

10 .3

-7.80

x

10 .2

9.5

-7.80

x 10 .3

-1.61

x

10 I

11.5

-1.70

x

10 .2

-2.73

x

10 "I

13.5

-I.00

x

10 .2

-4.27

x

10 I

15.5

-1.12

x

10 -I

-7.43

x

10 I

x 10 "I

8.91

x

10 I

x

9.22

x

i0 °

17.5

6.75

19.5

-9.31

Cost

schemes

processing this

The

However,

the

increased here,

with

and

added

119 s e c t o y i e l d 20.5

With

where

the

A x = i/2,

cost

of

absolute

error

be

the

cited

are

for

For

us

fluxes linearon a

unaccelerated

scheme

scheme

weighting

in t h e Let

respectively,

diamond

weighted

savings

space-time

sec,

than

coarser.

as/a = 0.5.

342

linear

a net

can

with

and

times

implicit

from was took

the

higher-order

obtained.

A x = A,

r = 1/2,

yield mesh

computations

took

took 432

coefficients

ii0

sec

sec.

The

after

each

step.

space-time

ranging

441

CPU

the

more

a n d r = 1/2,

took

All

they the

above

positive-flux

accuracy

took

but

Ax = i/4

schemes

evaluates

time

time

errors

the

require

because

discussed

(Note:

the

mft.

mesh,

time

problem

mft,

scheme

in e a c h

discussed

a given

conventional

mesh

positive-flux

the

20.5

computer.

same

iteration

the

algorithms

(CPU)

constant-linear

calculations.) the

for

unit

for

at time

and

10 "I

Considerations

higher-order-moments

11/780

right

10 .2

x

computed

with

x

1.96

illustrate

spatial

1/2

7.22

1.5

central

by

=l

r =

1/2

0.25

conventional

with

Ax

x 10 .3

The

VAX

algorithm

(% )

1/2

=

3.77

Computinq

linear

L-L

For

linear-linear spatial

r = 1/2, 0.072%

nearly

errors the

to 9%.)

30 s e c w i t h

method

example, scheme

0.9% The

with

ranging

algorithm

the Ax

from took

(except implicit

absolute

is m o r e

in

in

offset being

= i/2,

0.004% 35 s e c

the

last

linear

spatial

than

problem

1/2,

to 0.68% for

at

absolute

cell

diamond

errors

r =

on

the

scheme

ranging

from

148

A. B ADRUZZAMAN

14.4%

to

errors

80.9%

while

ranging

6.

In the large

arise

only

These

to

the

scheme

from

section,

implicit

the

the

difference we

52

sec

for

SOLUTIONS

implicit

implicit

scheme, briefly

assumption

needed

r = 2.

IN T I M E - D E P E N D E N T

arise

this

weighted A x = l/2,

5, w e s a w t h a t

spatial

In

with

ISSUES

errors

the

etc.

due

54.4%

in S e c t i o n

variable,

approximation,

positive-flux ~o

ACCURACY

example

errors.

temporal

the

from 4.1%

schemes

differencing

the

discrete

examine

and

produced

the

the

of

the

ordinates

errors

discrete

that

ordinates

approximation.

6a.

Implicit

Let with

us

temporal

consider

differencing

equation

6,

the

transport

equation

in

slab

geometry,

S = 0,

a°Jk+1(x'~)ax + °eee°Jk+~(x'~)= ~i ~ k(x,~) This equation the

Integrating right

can be solved

implicit

assumption equation

analytically

irrespective

in s p a c e t o e x a m i n e of

48 f r o m x = x L o n t h e

the

(48)

the accuracy

of

spatial

difference

scheme.

left boundary

to a point

x on the

we obtain $k÷I (x, ~) = ~k+1 (xL, ~) exp[- (x-x L) o efflbt] +

1

~tvA t

(49)

x

f

x L

F o r x L = 0, w r i t i n g

down

equation

49 s u c c e s s i v e l y

f o r k = 0,

i,

2,

... etc.,

we obtain ~k+1(x,~)

= exp(-XOeaJ~)

(s0) 1

+

(t,vA

o

0

Ik+1(x) ], t) k+~

where X

Ik*l(x)

Note that boundary

=

X I

f dxl 0

in e q u a t i o n condition

X (k)

f d x 2 ... o 50 t h e

while

f o

d21(k.l)~lo(X(k+l),~)exT~(x(k.l)Geff/~)

first term

in t h e b r a c k e t s

the second term arises

is d e t e r m i n e d

due to the initial

.

(51)

by the

condition.

Hmte-momen~ ~pmach

The m u l t i p l e condition.

EXAMPLE

integrals We

in e q u a t i o n

illustrate

149

51 can be e v a l u a t e d

for a g i v e n

initial

this w i t h two examples.

1

Let us c o n s i d e r such that the

a linearly

solution

increasing

to the t r a n s p o r t

function

initial

in a v a c u u m

equation,

! a ~ + ~a_~ v 8t ax

w i t h the

of p o s i t i o n

= o

(52)

0 < x < 3 ,

(53)

condition

1 = ~ x

~(x,~,0)

Z~9(

for ~ = 1 is (x,t)

-

l

(x-

(54)

vt)

2Ax

The

left b o u n d a r y

condition

is ~(0,t)

W h e n we d i f f e r e n c e

equation

52 i m p l i c i t l y

(as in e q u a t i o n

o eff = I / v A t

since

a

=

domain.

0.

In

The

condition

time,

implicit

in e q u a t i o n

the

spatial

solution

in

1

Let us e v a l u a t e Table

8 displays

shape

time intervals, the

implicit

large.

These

f

o

o

flux d i s t r i b u t i o n

the e x a c t At = Ax/v,

and the

given

implicit

2 A x / v and 3Ax/v.

scheme

overestimated

errors

decrease

but t h e y are still

propagates 51

across

the

spatial

using

the

initial

becomes,

x (k~

f dx 1 ...

(vAt) k+1

the

(56)

53,

exp(-x/vAt)

2Ax

6),

,

equation

x

%k+1(x ) _

(55)

: 0 .

the

by e q u a t i o n analytical

57 at a time solutions

We a s s u m e Ax = 0.5.

results

considerably

substantial.

dx(k+'>X(k÷l)exp[x
and

the

overall

as the t e m p o r a l

mesh

. (57)

6Ax/v.

w i t h three We see that errors

are

is refined,

150

A. BADRUZZAMAN

Table

8.

Exact and implicit analytical s e c t i o n (t = 6 A x / v ) .

solutions

IMPLICIT

for example

1 in t h i s

ANALYTICAL

x

EXACT

0.0

0

0.0

0.0

0.0

1.0

0

0.00296

0.02334

0.05367

2.0

0

0.09772

0. 2 1 8 0 2

0.31799

3.0

0

0.48187

0. 6 7 2 1 3

0.81201

3.5

0.5

0.79665

0.98693

1.13032

4.0

1.0

1. 1 7 5 1 8

1.34800

1.48638

&t

=

&x/v

&t

= 2AX/V

&t

= 3&x/v

4.5

1.5

1.59974

1.74579

1.87340

5.0

2.0

2.05499

2. 1 7 1 8 2

2.28539

5.5

2.5

2.52948

2. 6 1 9 0 3

2.71727

6.0

3.0

3.01542

3. 0 8 1 8 0

3.16484

EXAMPLE

2

In this

problem,

the

initial

~(x,~,0)

and the

boundary

condition

condition

= 1 + e -°x/" ,

the

solution

(58)

~ > 0

is @'(0,1 ~,t)

so t h a t

is

to the

=

source-free

1

+

(59)

e -°vt

transport

equation,

la~_ + a¢ +o,~=o vat ~ ax

(60)

(x, ~ , t)

(61)

is

Using

the

solution

given

initial

in e q u a t i o n

and

= e -°x/" + e -"vt

boundary

49 b e c o m e s ,

after

conditions,

the

performing

implicit

some

analytical

algebra,

k

~ k . l ( x ,Is ) = e_OX/. + e _ . V t e - X ° . z z / . ~ m=o

1 ( ~ v & t) m

x = m! (62)

+

1 - e -x, .~I~

( v A C a , ~ ) k+1

Note that the first term of

the

exact

modification transient

solution,

term

~o

in e q u a t i o n equation

of t h e t r a n s i e n t

is u n d e r e s t i m a t e d .

62 c o r r e s p o n d s

61.

part

/m !

The

of t h e

second exact

The additional

to the steady-state term

in

solution. term

equation As

62

a result,

in e q u a t i o n

part is

a

the

62 a r i s e s

Finite-momentsapproach

purely

due to the implicit

equations

61 a n d 62 v a n i s h

k ~ ~).

Thus,

However,

in

transients

both

are

Table

3/v,

significantly the

time

9.

a =

overestimates

interval

Table

the exact

far

the

and the

i. W e

states

results

the results,

that

as one moves arises

but

Exact and implicit analytical s e c t i o n (t = 3/v, ~ = i).

errors

solutions

IMPLICIT

as

when

may

yield

considered

next.

for At = i/2v and implicit

into the

from the

large

(such

the

scheme

interior

last term

of t h e s e c o n d

(since

as t ~ ~.

assumption

solutions

in b o t h

62

solution

example

see

terms

in e q u a t i o n

implicit

implicit

~ =

the underestimation

improves

steady

by the numerical

1 and

the

term

steady-state

from

The overestimation

compensating

the

in a r e a c t o r ) ,

assuming

left boundary.

62 o v e r l y

are

is i l l u s t r a t e d

9 displays

at t =

approach

which

introduced This

A s t ~ ~, t h e t r a n s i e n t

as d o e s t h e a d d i t i o n a l

equations

problems

large errors.

i/v

assumption.

151

term.

from

in e q u a t i o n Refining

remain.

for

example

2 in t h i s

ANALYTICAL

At = i/2v

At = i/v

x

EXACT SOLUTION

SOLUTION

0.0

1.04979

1.04979

0.0

1.04979

0.0

0.5

0.65632

0.63710

-2.93

0.64633

-1.52

1.0

0.41767

0.39326

-5.84

0.42514

1.80

2.0

0.18512

0.18929

2.25

0.23513

27.01

3.0

0.09957

0.12853

29.08

0.16809

68.81

4.0

0.06810

0.10450

53.44

0.14181

108.23

5.0

0.05653

0.09431

66.84

0.13143

132.52

6.0

0.05227

0.09024

72.66

0.12742

143.79

A

similar

analytical

conclusion

solution

ERROR

is

reached

in e q u a t i o n

the

SOLUTION

(%)

for

the

scalar

ERROR

flux

which

(%)

for

the

61 is

i

¢(X,t)

= f

~(x,~,t)

d~ = E2(xu)

(63)

+ e -"vt

o

Integrating

the

implicit

solution,

equation

62,

the

scalar

flux

is g i v e n

by

k

4~k+l(X)

= E2(XO)

+

e-art E ~ m=-o

"-~. I,,,(xOe~f)

(o .~x) -i.(xo .±z)] ( V ~ t O eft"

where

k+l

m-o

In!

(64)

152

A. BADRUZZAMAN

1

Im(xg~tf) = f dN e x p ( - x s e : : / M ) / M m

(65)

.

0

From

Table

i0 w e

see

that

those

shown

in T a b l e

Table

i0.

that

the

errors

in t h e

9 for t h e

scalar

angular

Angle-integrated scalar (a = I, t = 3/v).

flux

flux,

flux

are

somewhat

example

2,

section

x

EXACT

IMPLICIT At = i / 2 v

0.i

0.77233

0.72663

-5.92

0.5

0.37643

0.35457

-5.81

1.0

0.19830

0.20594

3.85

2.0

0.08732

0.12025

37.71

3.0

0.06043

0.09783

61.88

4.0

0.05299

0.09093

71.61

5.0

0.05078

0.08878

74.83

6.0

0.05011

0.08811

75.85

that and

The

examples

can

result

the

sign

differencing with

6b.

a grossly of

accurate. discrete equations

errors

can

of

which

the Let

the

63 a n d

discrete

can

the

The

large

When

further

assumption

(%)

actual

errors

magnitude

the

In a d d i t i o n ,

increase

implicit

spatial

for p r o b l e m s

leading

to

the

in time.

time-dependent in o n e d i m e n s i o n .

angular

Also,

not

be

this

assumption,

16, even

spatial

flux.

may

illustrate

64

can i n c r e a s e .

the

effects

scalar

ordinates

alone.

problem-dependent.

errors

in t h e i r

quadrature us

be

reference ray

results

erroneous

demonstrate

6a,

assumption

in

exhibit

subsection

the e r r o r s

5 using

ERROR

approximation

would

these

ordinates

discussed

in t h i s

implicit

sources,

propagation

velocities

order

the

in S e c t i o n

Discrete

As

the

is i n c l u d e d ,

seen

algorithms the

of

iterative

results

discussed

from

larger

at ~ = i.

the

0

with

separation depending

sufficient

using

Example

steady-state

is a p p r o x i m a t e d

/ e-a~d~

modes

discrete This

different

in time. on

for

the

the

2 of part

ordinates arises

of

effective

T h i s can y i e l d

problem, scalar

Section the

from

a given

flux

6a.

to

With

be the

solutions

in

by

" E n=l

exp(-xa/~n)wn "

(66)

Finite-moments approach

153 N/2

where

the

angular

quadrature

weights

are

normalized

as

E

For

wn = 1

the

nll

implicit

solution,

the

integral

in equation

65 i s

approximated

as

N/2

I~(xo ,rz) " SN/2,m (x° ,rr) = ~

(67)

(-xo ,z:/l*,) wn/ (~n) m

exp

rill

Table t = 3/v,

ii d i s p l a y s

with

analytical

At

and the

boundary.

The

these

errors

large

distances

scalar The

=

implicit

from

Table

ii.

in t h e d i s c r e t e - o r d i n a t e s

5)

in

cases,

scheme

the

approximates

primarily

=

left its

implicit

errors

are

somewhat

of t h e a n g u l a r

boundary,

flux

from

the

implicit

the

counterpart exhibit

flux,

example

N=

8

9.47

4.46

0.5

4.74

-7.63

x 10 "I

2.75

1.0

-2.61

-5.31

x 10 I

2.71

3.00

2.92

8.15 x 10 -I

-2.31

2.0

-1.85

1.46

x 10 "I

32.30

34.56

3.0

2.15

x 10 "I

-1.31

x 10 .3

60.76

60.71

4.0

2.82

x 10 "I

-7.99

x 10 .3

71.71

71.40

-8.55

x 10 .4

74.91

74.80

3.52

x 10 .4

75.87

75.85

1.88

x

* Relative

to the

1 0 .2

angle-integrated

7.

In this the

paper

time-dependent

finite-moments equivalently The

we

and expands

implicit-moments

properties growth

the

scheme

in

general,

AND

two

equation.

in s p a c e .

of the space-time

factor,

examined

transport

expansion

exact

DISCUSSION

have

flux is

2,

N = 8

4

0.i

6.0

large

IMPLICIT - SN At = i/2v

- SN N=

x 10 .2

well.

(% )

N=

9.91

At

analytical

x

5.0

left

expected,

scheme

(S,) s c a l a r

ERROR*

4

As

at the

is r e f i n e d .

discrete-ordinates analytical

both

near the

larger.

quadrature

implicit

fluxes

In

assumption.

E r r o r s in d i s c r e t e o r d i n a t e s s e c t i o n 6a (a = i, t = 3/v).

ANALYTICAL

scalar

solution.

are considerable

angle-integrated

scalar

due to the

the

the errors

as the order

discrete-ordinates

errors

(k

implicit

decrease

flux

the errors

i/2v

scalar

flux

finite-moments One

in t e r m s

is

of t e m p o r a l

schemes on

the

approaches

implicit

unconditionally

depends

i0.

CONCLUSIONS

The other procedure

moments

in T a b l e

in t i m e

treats and

relative

and

solve uses

a

time and space

spatial

moments.

The

stability

stable.

are more

to

complex, values

and the error of

the

error

154

A. B ADRUZZAMAN

a m p l i t u d e s of the v a r i o u s moments. r e s u l t s t h a t the C-C s p a c e - t i m e numerical

evaluation

with

a

We saw from b o t h a n a l y t i c a l and n u m e r i c a l

algorithm

scattering

is u n c o n d i t i o n a l l y

stable.

source,

that

we

found

a l g o r i t h m was s t a b l e for r ~ 1 and the L-L a l g o r i t h m was s t a b l e Stability

properties

with

higher

harmonics

are

being

From

the

C-L

for r < I.

investigated

analytically.

N u m e r i c a l t e s t s on a p r o b l e m w i t h an i t e r a t i v e s c a t t e r i n g s o u r c e showed that b o t h the c o n v e n t i o n a l

i m p l i c i t linear d i a m o n d scheme and the p o s i t i v e -

flux w e i g h t e d a l g o r i t h m s p e r f o r m poorly, of the two,

especially

for

large

a l t h o u g h the latter was the better

scattering

ratios.

The

implicit

finite-

m o m e n t s s c h e m e w i t h a linear e x p a n s i o n d i s p l a y e d s i m i l a r a c c u r a c i e s near the r e f l e c t i n g b o u n d a r y but was c o n s i d e r a b l y Also,

using

expansion,

the

one

can

time-dependent solutions,

finite-moments prevent

solutions

especially

negative more

to

with

an

fluxes.

severely

adequate

Negative

than

they

away from it. order

fluxes

degrade

in

can

the

degrade

steady-state

in n o n l i n e a r p r o b l e m s such as r a d i a t i v e t r a n s f e r with

t e m p e r a t u r e e v o l u t i o n in time. needed

b e t t e r as one m o v e s

scheme

accelerate

the

In addition,

convergence

a c c e l e r a t i o n schemes, w h i c h are

in i t e r a t i v e

scattering,

emission

or

f i s s i o n s o u r c e problems, e s p e c i a l l y if the m e d i u m is o p t i c a l l y thick, may not f u n c t i o n w e l l w i t h n e g a t i v e fluxes.

The C r a n k - N i c h o l s o n scheme scheme)

provided

negative path.

fluxes

The

more arise

unless

and the C-L

C-C

both the c o n v e n t i o n a l time, linear m o m e n t s test problem.

(which is also the s p a c e - t i m e linear diamond

accurate

solutions

the

mesh

than

spacing

algorithms

the is

performed

and the f i n i t e - m o m e n t s

implicit

less

than

one

significantly

implicit

scheme,

but

mean

free

better

than

schemes.

The space-

(L-L) scheme y i e l d e d e x t r e m e l y a c c u r a t e s o l u t i o n s to the

It was one or m o r e orders of m a g n i t u d e m o r e a c c u r a t e than the

o t h e r a l g o r i t h m s tested.

Thus,

in spite of the r e q u i r e m e n t that this m e t h o d

s a t i s f y the C o u r a n t condition, one can use c o n s i d e r a b l y fewer time steps with it

than

greater

with in

other

algorithms.

multidimensional

The

advantage

problems.

Also,

of

the

algorithm

we

saw

earlier

will that

be the

s t a b i l i t y of a g i v e n f i n i t e - m o m e n t s a l g o r i t h m d e p e n d s on the the o r d e r of the e x p a n s i o n and the r e l a t i v e a m p l i t u d e s of the errors in the moments.

Thus, it

may

obtain

accurate

s o l u t i o n s by r e t a i n i n g a p p r o p r i a t e terms in the m o m e n t s expansions.

This is

be

possible

currently being equivalently,

to

relax

the

investigated.

stability

criteria

and

yet

We should note t h a t t r e a t i n g

time and space

one i n c r e a s e s the d i m e n s i o n a l i t y of the p r o b l e m r e l a t i v e to the

i m p l i c i t scheme, w h i c h is e q u i v a l e n t to an e f f e c t i v e s t e a d y - s t a t e algorithm. Thus,

the

space-time

finite-moments

intensive than impicit-in-time

scheme

will

be

f i n i t e - m o m e n t s schemes,

computationally

more

e s p e c i a l l y in three-

Finite-momentsapproach

d i m e n s i o n a l geometry.

155

However, the i n c r e a s e d c o m p u t i n g costs are e x p e c t e d to

be m o r e t h a n o f f s e t by the gain in accuracy. We b r i e f l y e x a m i n e d the a c c u r a c y of the i m p l i c i t a s s u m p t i o n in time with the s p a t i a l v a r i a b l e t r e a t e d a n a l y t i c a l l y . i m p l i c i t s c h e m e e x h i b i t e d large errors, difference

scheme

and s o u r c e

iteration.

In the e x a m p l e s c i t e d here, the

w h i c h may be m a g n i f i e d by a spatial However,

the

implicit

approach

is

a t t r a c t i v e due to its f a v o r a b l e s t a b i l i t y properties. The m e t h o d s d e v e l o p e d h e r e can be a p p l i e d to a v a r i e t y of t i m e - d e p e n d e n t problems.

We

have

developed

a

finite-moments

algorithm

to

solve

the

r a d i a t i o n t r a n s p o r t e q u a t i o n c o u p l e d to the p r e c u r s o r e v o l u t i o n e q u a t i o n in r e a c t o r kinetics.

We are also t e s t i n g the i m p l i c i t f i n i t e - m o m e n t s scheme on

p r o b l e m s of r a d i a t i v e t r a n s f e r c o u p l e d to t e m p e r a t u r e e v o l u t i o n in a medium. Our p r e l i m i n a r y t e s t r e s u l t s on m o d e l r a d i a t i v e t r a n s f e r p r o b l e m s 17, w h i c h are extremely

thick

optically,

show

that

the

implicit

p e r f o r m s s i g n i f i c a n t l y b e t t e r t h a n the c o n v e n t i o n a l d i f f e r e n c e scheme.

For example,

finite-moments

scheme

i m p l i c i t linear diamond

in one e s p e c i a l l y d i f f i c u l t m o d e l problem,

the f i n i t e - m o m e n t s s c h e m e w i t h a q u a d r a t i c spatial expansion, y i e l d e d results c o m p a r a b l e to t h o s e from the I m p l i c i t M o n t e C a r l o m e t h o d 17, w h i l e the linear diamond

and

temperatures.

the

linear

Detailed

discontinuous results

of

schemes

these

failed

tests

will

to be

yield

positive

reported

in

the

future. ACKNOWLEDGMENT The a u t h o r g r a t e f u l l y

acknowledges

B. D.

University

Ganapol

Scrivner,

of

and W.

the C.

Fan

of

the m a n y v a l u a b l e

of A r i z o n a

Sandia

National

and

Drs.

comments J.

H.

Laboratories.

of P r o f e s s o r Renken,

He

also

G.

J.

thanks

P r o f e s s o r A. J. H. G o d d a r d of the I m p e r i a l C o l l e g e and the e d i t o r of P r o q r e s s in N u c l e a r E n e r a v for t h e i r p a t i e n c e d u r i n g p r e p a r a t i o n of the manuscript.

REFERENCES

i.

Wahl J. S., et al.

(1970) "The T h e r m a l D e c a y T i m e Log," Soc. Pet. Eng. J,

p. 365. 2. S a n d m e i e r H. A., D u p r e e S. A. and H a n s e n G. E. (1972) " E l e c t r o m a g n e t i c Pulse and T i m e - D e p e n d e n t Escape of N e u t r o n s and G a m m a Rays from a N u c l e a r Explosion," Nucl. Sci. Eng., 48, 343.

3.

L e w i n s J.

Oxford.

(1978)

Nuclear Reactor Kinetics

and Control,

P e r g a m o n Press,

156

4.

A. B ADRUZZAMAN

P o m r a n i n g G. C.

Press,

5.

(1973) The E q u a t i o n of R a d i a t i o n H y d r o d y n a m i c s ,

Dupree

S.

Calculations

A.,

et

Using

Reed

W.

H.

7.

(May

1971)

"Time-Dependent

of D i s c r e t e

Ordinates,"

(also cited as O R N L - 4 6 6 2

(July

ordinates Program LA-4800,

al.

the M e t h o d

Scientific Laboratory

6.

Pergamon

Oxford.

1972)

"TIMEX:

A

Neutron

and

LA-4557,

Photon

Los Alamos

and SAI-70-125).

Time-Dependent

Explicit

for the S o l u t i o n of the M u l t i g r o u p T r a n s p o r t

Discrete-

Equations,"

Los A l a m o s S c i e n t i f i c Laboratory.

R h o a d e s W. A.

and Childs R.

L.

(April

1982)

"DOT IV, V e r s i o n

4.3:

An

U p d a t e d V e r s i o n of the DOT4 One- and T w o - D i m e n s i o n a l N e u t r o n / P h o t o n T r a n s p o r t Code," O R N L / C C C - 4 2 9 ,

8.

A l c o u f f e R. E.

Oak R i d g e N a t i o n a l Laboratory.

(April 1988)

"A D i f f u s i o n A c c e l e r a t e d S, T r a n s p o r t Method

for R a d i a t i o n T r a n s p o r t on a G e n e r a l Q u a d r i l a t e r a l M e s h , " Proc.

ANS Topical

Mtg. on A d v a n c e s in N u c l e a r E n g i n e e r i n g C o m p u t a t i o n and R a d i a t i o n Shielding, Santa Fe, NM, Vol 2, P a p e r 48.

9.

P e a r c e R.P.

and M i t c h e l l A. R.

(1962),

"On F i n i t e D i f f e r e n c e M e t h o d s of

S o l u t i o n of the T r a n s p o r t E q u a t i o n , " Math.

i0.

Badruzzaman

Equations,"

ii.

B a d r u z z a m a n A.

Nuclear

Science

P l e n u m Press,

12.

A.

(1988)

"A

Trans. Am. Nucl.

(1990)

and

Method

for

16, 155.

Time-Dependent

Transport

56, 303.

"Nodal M e t h o d s

Technology,

N e w York.

Nodal

Soc.,

Comput.

Vol

21,

in T r a n s p o r t T h e o r y , " A d v a n c e s J.

Lewins

and M.

Becker

in

Editors,

(This is a r e v i e w chapter.)

Bell G. I. and G l a s s t o n e S.

(1970) N u c l e a r R e a c t o r Theory, Van N o s t r a n d

R e i n h o l d Co., N e w York.

13. C l a r k Jr., M. and H a n s e n K. F. Analysis, A c a d e m i c Press, N e w York.

14.

G a n a p o l B. D.

(1986)

(1964)

Numerical

Methods

of

Reactor

" S o l u t i o n of the O n e - G r o u p T i m e - D e p e n d e n t Neutron

T r a n s p o r t E q u a t i o n in an I n f i n i t e M e d i u m by P o l y n o m i a l R e c o n s t r u c t i o n , " Nucl. Sci.

Eng.,

92, 272.

15.

G a n a p o l B. D.

(Sept 1987), U n i v e r s i t y of Arizona, p r i v a t e communication.

Finite-moments approach

16.

M a r t i n W. R., et al.

(1981)

157

"Phase-space Finite Element Methods Applied

to the F i r s t - o r d e r F o r m of the T r a n s p o r t E q u a t i o n , " A n n a l s of Nucl. 8, 17.

Energy,

633. Fleck

Scheme

for

Transport,"

Jr.,

J.

A.

Calculating J. Comp.

and

Cummings

Time

Phy., 8,

and 313.

J.

D.

Frequency

(1971)

"An I m p l i c i t

Dependent

Monte

Nonlinear

Carlo

Radiation