Solutions of diffusion equations by Fourier expansions

Annals of Nuclear Science and Engineering, Vol. 1, pp. 547 to 553. Pergamon Press 1974, Printed in Northern Ireland

SOLUTIONS OF DIFFUSION EQUATIONS BY FOURIER EXPANSIONS NOBUO OHTANI, JUNGCHUNG JUNG, KEISUKE KOBAYASHIand HIROSHI NISHIHARA Department of Nuclear Engineering, Kyoto University, Kyoto, Japan

(Received 18 March 1974; is revisedfrom 19 May 1974) Abstract--One- and two-dimensional diffusion equations in slab geometry are solved by a method of Fourier expansion. In this method, at first, equations for the fluxes on the boundaries and their normal derivatives are derived. Applying boundary conditions, these equations are solved and all boundary values are determined. Then using these boundary values, the Fourier coefficients of the flux in the region are calculated. Different from the eigenfunction expansion method, the function series used for the expansion is independent of the boundary conditions. Therefore multi-regional problems are also solved by this method. The results of the numerical calculations are given and compared with the results by the usual finite difference method. 1. INTRODUCTION One- and two-dimensional diffusion equations in slab geometry are solved by a method of Fourier expansion. Diffusion equations, especially in two-dimensional geometry, are solved in the following ways; the finite difference method (e.g. Fowler, Tobias and Vondy, 1967), the finite element method (Kaper, Leaf and Lindeman, 1972), the flux synthesis method (Kaplan, 1962; Kobayashi, 1968), the Green function method (Kobayashi and Nishihara, 1967), the eigenfunction expansion method (Garabedian and Thomas, 1962), and the other function expansion method (Cohen, 1956). In this paper solutions of the diffusion equations are given in the Fourier expansion series. An outline of the method is as follows; first we multiply the non-homogeneous diffusion equation by the general solutions of the homogeneous diffusion equation, and integrate the resulting equations over the homogeneous region. Then the volume terms of the integral are cancelled out and equations which include the fluxes on the boundaries and their normal derivatives are obtained. Applying the given boundary conditions to these equations, all fluxes on the boundaries and their normal derivatives are determined. In one-dimensional case, boundary values are given as numerical values and two independent solutions of the homogeneous equation are given explicitly. In two-dimensional case, boundary values are one-dimensional functions, which are expanded in one-dimensional Fourier series. As for solutions of the homogeneous equation, we use the function series which are used in the solutions of a Dirichlet

problem for the homogeneous equation (Churchill, 1963). Secondly, expecting that the flux inside the region can be expressed in terms of a Fourier series, we multiply the diffusion equation by the general term of the Fourier series, and integrate it over the region. Using the boundary values already obtained, the Fourier coefficients of the flux are determined. In two-dimensional case, we use a double Fourier series. This method is different from the usual eigenfunction expansion method as follows. In the eigenfunction expansion method, the eigenfunctions are determined with the boundary conditions; the expansion cannot be defined without specifying all the boundary conditions at the regional and outer boundaries. In muhi-regional problems, boundary conditions at material interfaces are continuity conditions and are not given in explicit values, then we can not determine the eigenfunction series for each homogeneous region separately. Therefore, for example, Garabedian and Thomas (Garabedian and Thomas, 1962) expanded the flux over the whole system in one eigenfunction series, which is determined only from the outer boundary conditions. The effects of the heterogeneity of the system are introduced by the transient terms. In the present method, each of the functions used for expanding the flux at interior points of a homogeneous region is chosen independently of the flux boundary conditions. Therefore, it is always possible to expand the flux in regionwise Fourier series. It should be noticed that each term of the resulting Fourier series does not generally satisfy the boundary conditions, and that the series does not uniformly converge at boundaries; the boundary values of the Fourier series do not generally coincide with the flux 547

548

NOBUO OHTANI, JUNGCHUNG JUNG, KEISUKE KOBAYASm and H m o s m NISHn-IARA

values at the boundaries. This fact is not an obstacle for this method because we do not apply boundary conditions to the boundary values of the Fourier series of the flux inside the region. However the convergence of the series is rather slow and we cannot differentiate the series term by term. In the following section, general equations in oneand two-dimensional problems are derived. In Section 3, numerical results of a simple twodimensional problem are given. 2. BASIC E Q U A T I O N S

by parts over 0 < x < a. Using equation (4) ~v(a)9'(a) - V,(0)tp'(0) - ~p'(a)cp(a) + ~p'(0)9(0)

¢u

where 9' (a) = I d q~(x) l = = a etc.

The two independent solutions of equation (4) are e ~= and e-~=. Substituting ~o(x) = e ~= into equation (5), we get e~aqJ(a) - 9'(0) - Ke-~ag(a) + Jcq~(O) =

A one-group steady-state diffusion equation in the homogeneous region is

-

- DV~9(r) + Y'a~(r) = koo~q~(r) + S(r),

Q(x):

®dx.

(6.1)

(1) Similarly substituting ~o(x) = e- ~ ,

where

e-~aqJ(a) - q/(0) + Ke-~ag(a) - Kqo(O) =

~o(r) : the neutron flux, D Za k~o S(r)

(5)

= - IiaQ(x)vJ(x) dx,

: the diffusion coefficient, : the macroscopic absorption cross section, : the infinite medium multiplication factor, : the external neutron source.

Equation (1) can be rewritten as V2~(r) - ~0(r) = - O(r),

(2)

where

--

Q(x)e -~= dx.

(6.2)

Among four unknowns 9(0), 9(a), 9'(0) and 9'(a) in equations (6.1) and (6.2), two of them, or two relations are given as boundary conditions. Thus all boundary values are obtained. Next, we multiply equation (2) by exp [ - - i(2=mla) x] and perform integration by parts over 0 < x < a.

1 K2 = ~ ( 1

a(r) =

-k=),

S(r) D

x

9'(a) -- q0'(0) + i a - {9(a) -- ~p(0)} + aQra

,

(7) 1.

One-dimensional case

Now we study a one-dimensional problem in the slab geometry. We consider equation (2) in the following region

where 9m and Qm are the Fourier coefficients of the flux and the external source respectively, i.e. +co

9(x) =

O~x<=a.

Q(x) =

Equation (2) becomes d~

q~(x) -- K~qJ(x) = -- Q(x).

(3)

the region; x = 0 and x = a. Let ~p(x) a general solution of the homogeneous equation of equation (3). (4)

We multiply equation (3) by ~(x) and integrate it

9m exp [i(2=m]a) x], Qm exp [i (2=m/a) x].

~ m=

d First, we calculate 9(x) and ~x ~o(x)on boundaries of

da dx ~ ~o(x) -- ,:~o(x) = 0.

~

m ~ --aO +00 --00

With 9(0), 9(a), ~0'(0) and 9'(a) already obtained, we get q'm from equation (7), and the Fourier expansion of the flux is determined. 2.

Two-dimensional case We set the region V in two-dimensional slab geometry as

V : 0 ~ x =
9(x, y ) - , : 9 ( x , y )

Q(x, y).

(8)

Moreover we use the general solution ~0(x,y) of the

Solutions of diffusion equations by Fourier expansions homogeneous equation, +

549

where

~v(x, y) - ,J~v(x, y) = O.

(9)

We multiply equation (8) by ~(x, y) and integrate it over the region V. Applying Gauss' theorem, we get

~m

2~rm c--~

2

(m = O, ~1, :52 . . . . ) g m (x') =

~ (x', y)

, y=0

sn. {~o(x, y).v~(x, y)} dS

0

g,2) (y') = [~_xCP(X, y )lz=a, - js,," {~(x, Y)VV(x, y)}dS = - j;O.(x, y)~(x, y ) d r , (lO)

g(S'(X') ~

'p (X ,y

u=b"

where

0

gm(Y') = f ~-x~(x'" V');*=0'

S: the surface of the region V, n : the normal unit vector at N. The general solution of equation (9) can be expanded in the following functions, exp [

-i(2~m/a) x] exp [ ~ / ~

+ (2~rn]a)~y], (m ~ 0 , zk 1, 4-2 . . . . . ) (11.1)

f m (x') = ~ (x', 0),

f(2) (y,) = ~ (a,y'),

fc~) (x') = v (x', b),

f~4~(y,) = ~ (0,y').

When we use equation (11.2), we get

- foadX ' exp ( ~ fl~,Ux')gm (x') + exp (~f~.)a)

and exp [ ± V ~ + (2~n/b) ~x] exp [ - i (2~n/8) y], (n ~ 0 , 4- I, 4-2 . . . . . ) (11.2) Therefore, we use equations (11.1) and (11.2) as the function ~o(x, y) in equation (10). Whenwe use equation (11.1) as ~o(x, y) in equation (10), we get -

+

~0~dx

t

e x p ( ~ / / .111 x ) gr

£

+ exp ( ± ~ ' b )

f~

13

)(x')

-fobdY ' exp

dx exp( - , % . ,a, x ), g (n ( x,) +

dy' exp (--tfl. . cz~y ,)g (2) (y),

"<

y exp(:k%y)g

(y)

. (11 v dx t exp ( -~g,~ x )3"(~)( x t)

~ i~

II

~ Jo

r

dx exp ( ~

q: f12' exp (4-¢2'a)

(--i/¢.')")ov~)(y')

~2'x')fm (x') f]dy' exp (-ip~."v')f '~' (y')

d0

"~) Jo(~"ax ' e x p i"4r e 0~m '.~(3~.tx),. +~P, , x~f ~=2'

dx' exp ( - ~ %. x )fI~)(x ) ' ~

(b

4-;L"' .

.

(a) , . , - ( 2 1

~'exp(±%

(2)

t

b)

•

{1) ¢

dx e x p ( - t % x ) f

(3)

dx'

£

dy' exp ( -i~.'~ y') f

~4~(y')

dy'exp(±/~, x ) e x p ( - ~ / ~ , y )

t

Q ( xt, y )¢,

(x)

0

--lea' f O dy' exp(

(n = 0, :LI, ~2 . . . . . )

4- o~,. t~) y t .)j~ ( 4 ) (y')

(12.2)

where

= -

£;o dx

q), exp (-~% " ~" x )' e x p ( = '% y~)

'

0

(m -- 0, :~ 1, ± 2 . . . . . )

× Q (x',y'), (12.1)

+ (2bn)z - - ,/~°,2, ~ - 2~n ~-, (n = 0 , ~1, i 2 . . . . . . )

550

NOBUO OHTANI, JUNGCI-I'UNOJUNO, K~str~ KOBAYASHIand HIROSm NISHIHARA

We expand the boundary valuesg (a) ~ ~.(4),f(1) ...f(4) andthe external neutron source Q(x,y) in the Fourier series,

where (z~l) = 9',..

ft~ ~exp

• (1)x ~"). exp (-r-E.(1) x t )" d x ¢, (ta,.

+oo

g(t) (x) =

g(~"exp [i (2zrfl/a) x] (l = 1, 3),

~ ~9= --(x)

" (~) (x) =

~

( ::E:2) = ~,..

f ' f exp [i (2~rp/a) x] (l

=

. f a ~exp

(3) y t~). exp (tfl'.~'y) • ( ± ~,.. dy'.

3),

1,

~=--oo

g(t) (y) =

g~, exp [i (2zrq/b) y] (l = 2,4),

~ q = --oo +cO

f(t) (y) =

f~" exp [i (2=q/b) y] (l = 2, 4),

~ q=--oo +co

Q (x, y) =

+~o

~

Q~q exp [i (2rrfl/a) x]

~

~=--o~

q=--oo

× exp [i(2,q/b)y]

(13)

Then we substitute equation (13) into equation (12.1) and (12.2). In equation (13), we restrict the summation indices to those such as - M _~ p _-
Applying the boundary conditions given (e.g. f(1) = f ( 2 ) = f ( 3 ) = f ( 4 ) = 0 ) ' we solve equation (14.1) and (14.2), then Fourier coefficients of all boundary values fix) ~f(a), g(1) .~g(a) are obtained. Next, we calculate the Fourier coefficients of the flux in the region. We multiply equation (8) by exp [-- i (2zrm/a) x] exp [ - - i(2zrn/b) y] and integrate the equation over the region V, then we get the Fourier coefficient ~0,nn of the flux in the region, i.e. f~

=

;F dx

dy~o(x,y)exp [-i(2~rm/a)x] × exp [ - i(2rrn/b)y] 1

-

37

(1)

-ag,,, +

(-;-2)

(2)

~. 9',,~ g~ q = - - ~V

+ exp (+%,b)ag,, (~' (~, (3)

:-(1)

iV

• (1)

q=-~

~,.~('~2'g~("

27rm ...(2,

( ::~Z2) ~t-(2)

X i(f(2' - f 2 ' ) + 7

q=--N

Z~2 (Z(2}

z

1. (3, (1,. 1 (g(z, ., 2rrn ~(g,, - - g ~ ) +a--g" ) +--~

X

(2)~\

t(y,

_ft,,) + Q,,,, }, (15)

,,-(3)

,, exp t ± %, o) af,~

M •

(I)

--,cz,,

N

( - I - 2 ) z.(a)

2 7,,~ J. q = --N

N

=--a

(-I-3)

Z ~q q = --N

Qmq,

~p(x,y) =

~r

~. m=--M

q~mnexp [i (2rrm/a) x]

~ n=--N

X

(14.1) (-M

_~ m _--
-

(-]-1)

~,

(1)

g.

(I)

+ e x p ( - 4 - f l . a)bg.

+

M

~,

12}

~ (.: E l )g~(s) - b g (4)

~o = - - M

M •-(~) ~ 7(4-x~f(I) p=--M q: ~ . ) e x p t-± p -.( I , PC

(1)

[4)

M

Numerical calculations are carried out for the simple two-dimensional homogeneous system shown in Fig. 1. We impose the boundary conditions that the fluxes are zero on boundaries, that is f(1)(x ) =ft3)(x ) = 0,f(2)(y) =ft4)(y) = 0.

(17)

Further the external neutron source is choosen as ..,.(2,

a)oj. • (21

+fl,, bf,, = - - b

(16)

3. NUMERICAL RESULTS*

Equation (12.2) becomes M

exp [i (2rrn/b) y]

3/

a(x, y) = (z~zl)tq3)

(~zl)

~_, ~,,,, Q•n, p=--M

(--N
~N),

(14.2)

1 :

a--c 2

<=x,y

a+c 2

<= - -

0 : otherwise• * These computations were performed by using FACOM 230-75 digital computer of the Data Processing Center of Kyoto University.

Solutions of diffusion equations by Fourier expansions

551

The right-hand side of equation (14.3) is calculated directly, without expanding Q(x, y) into Fourier series. Hence 2

a+c

(i)

2

(I)

ff'm 17/

-- ~

cos ---f- sin

~m C

2

m

× {exp Lr--~%: (2~(a + c)] Source

~ C~~ ~

- exp L--2%

~m C

--

O-C

So

2

{expFl - - g1% (~)(a + c)] i ,2, (a -- c)] - exp L--z%

I ~.-c 2

I o+c 2

(m :~ O)

(m = 0).

For the flux in the region, we get from equation (16),

(11

1

Fig. 1. The calculated system. --

+ K2

From the property of symmetry of the system, g(l)(s) = -g(Z)(s) = -g(3)(s) =g(4)(s).

(18)

Substituting equations (17) and (18) to equations (12.1) and (12.2), two equations become identical with each other.

×

- - (ag o r e + g o . ) + Qc~,,,~ ,

(16.1)

where ctVeemn = qPmn -I Cpm_n + q)-mn -t- Ct) m_n,

- [ i + e x p ( - %t2~aa. I

dscosc~(.:I

-

2

~

'a

sgre(s)

Qccmn : Qmn + a m - n + Q-ran + Q - m - n

ds exp ( - %(~)s)g (i)

4 = - ~ Qman,

(s)

0

=-

dx'

dy'cos% x expt-%y) × Q(x',y')

Defining

(1)

(i)

"nT/'/

(12.3)

+g_m,

(m = 0).

M

+2

q=0

(m = 0,1,2 . . . . .

M)

(14.3)

where

=:

%emo COS +

(2) \

(I}

__

12~

exp(-%.a)]

2rrm x + %c0m COS T y

n=l

a

The dimensions and the constants of the calculated system are

ds • exp t - : % s) • cos % s (Xm

a

M M 2ran 2~rn ]~ ~ %cm,~COS XCOS--~--y. ~=I

= ~ [ 1

( m :~ 0 )

~(x, y) = ~qP 1 ccO0

expt-%a)~g~,.

I

a

The flux at each point in the region is calculated by the following equation,

we get from equations (14.1) or (14.2),

f:

sin

(I)

go~ = g ~

-(1 +

1)" . .-.-. .

(~m ~

(q¢0)

a =8,

c =4,

x°- = 0 - 0 4

The flux along the line x = 3.0 are shown in Fig. 2. l f o a ds exp

(-%.~2,s)

1 = ~1 ~,

X [1 - e x p ( - % (~ a ) ]

(q = 0)

*The calculation is performed by the program EXTERMINATOR-2 developed by Fowler et al. (1967).

NOBUO OHTANI, JONGCHUNG JUNG, KEISUKE KOBAYASm and HIROSHI NISHIHARA

552

calculating the total neutron balance rather than spatial distributions. Further, the integral of the flux can be obtained readily as follows: foadx fO ~dy ~o(x, y) = T aS" ¢~oo. "k~, i la_

4

5

6

7

8

Y

Fig. 2. Fluxes along the line x = 3. Present method, 25 x 25 terms. Present method, 4 x 4 terms. Present method, 2 x 2 terms. • EXTERMINATOR-2, 33 × 33 meshes.

Results of the calculations by the present method are compared with those by the finite difference method.* The calculation by the finite difference method is performed in a quarter of the region because of the symmetrical property of the system. For this system, the series converges sufficiently at each point with 7 x 7 terms, and the results by two different methods show an excellent agreement. In Table 1 we show the integral of the flux and the computational time by each method. By the present method, the integral of the flux converges very rapidly with increased number of terms and we can attain the same accuracy of the integral in shorter computational time than by the finite difference method. This indicates that the present method will provide us a powerful means for treating a criticality problem where the problem is essentially that of Table 1. The integral of the flux and the computational time by the Present Method No. of Terms 2x 4x 7x 13 x 25 x

2 4 7 13 25

by EXTERMINATOR-2

Integral T i m e No. of Integral Time of flux (msec) meshes of flux (msec) 574.4 572.3 572.2 572.2 572.2

64 77 115 236 677

5x 9x 17 x 33 x

5 9 17 33

545.4 566.4 570.8 572.1

452 1066 4809 26959

4. DISCUSSION AND CONCLUSION The diffusion equations in slab geometry have been solved by the Fourier expansion. Formulations were restricted to one-regional problems, but multi-regional problems can also be solved in the same way. In multi-regional problems of two dimensions, equations for fluxes and their normal derivatives at the outer and interface boundaries are derived by applying the same procedures to each of the homogeneous regions. Using given boundary conditions to these equations, we obtain all the boundary values. Next, the flux in each homogeneous region is expanded in a double Fourier series. With boundary values already obtained for each region, the Fourier coefficients of the flux in the region are determined. In the sample problem shown in this paper, the flux took periodic boundary values, and consequently, the Fourier series converged uniformly even at the boundaries. However, in multi-regional problems, the flux in each homogeneous region does not generally assume periodic boundary values, and the Fourier series of the flux does not converge uniformly at boundaries; the boundary values of the Fourier series of the flux in each homogeneous region do not take acculate values. In such problems the series converges rather slowly at each point in the region and cannot be differentiated term by term. In this method, the boundary conditions are not directly applied to the boundary values of the Fourier series. Each of the coefficients of the Fourier expansion of the flux at internal points is determined separately with the relation (7) or (16), which include the boundary values of flux and its normal derivatives. Thus we can solve such multi-regional problems. The algebraic equations to be solved contain only boundary values as the unknowns. Fourier coefficients of the flux in each region are determined by equation (7) or (16), without solving any algebraic equations. Therefore, the number of the unknowns of the algebraic equations in this method is far less than that in the finite difference method. Moreover, the present method is feasible to be applied for multiregional and multi-group problems in two dimensions.

Solutions of diffusion equations by Fourier expansions

553

Multigroup Neutron Diffusion Equations In Two Dimensions, ORNL-4078. Garabedian H. L. and Thomas D. H. (1962) Nucl. Sci. Engng 14, 266. REFERENCES Kaper H. G., Leaf G. K. and Lindeman A. J. (1972) Churchill R. V. (1963) Fourier Series and Boundary Value NucL ScL Engnff 49, 27. Kaplan S. (1962) Nucl. ScL Engng 13, 22. Problems, 2nd edn. McGraw-Hill, New York. Kobayashi K. and Nishihara H. (1967) Nucl. Sci. Engng Cohen E. R. (1956) NucL Sci. Engng 1,268. Fowler T. B., Tobias M. L. and Vondy D. R. (1967) 28, 93. EXTERMINATOR-2: A Fortran IV Code for Solving Kobayashi K. (1968) Nucl. Sci. Enffng 31, 91.

Acknowledgement--The authors wish to express their thanks to Mr. M. Nakagawa of JAERI for making available the computer program EXTERMINATOR-2.