A comparative study of analytical solutions for some one-dimensional transport equation approximations

Progress in Nuclear Energy, Vol. 33, No. 3, pp. 289-300, 1998

Pergamon

© 1998 Elsevier Science Ltd All rights reserved. Printed in Great Britain 0149-t970/98 $19.00 + 0.00

PII: S0149-1970(96)00023--6

A COMPARATIVE

S T U D Y OF A N A L Y T I C A L

ONE-DIMENSIONAL

TRANSPORT

SOLUTIONS FOR SOME

EQUATION

APPROXIMATIONS

Augusto V. Cardona 1 and Marco Tt~llio de Vilhena 2

1 IMA - PUCRS Av. Ipiranga, 6681 - pr~dio 15 90619-900, Porto Alegre - Brazil E-Mail avcardona~music.pucrs.br

2PROMEC - UFRGS Av. Osvaldo Aranha, 99 - 4o andar - DENUC 90046-900, Porto Alegre - Brazil Abstract

The aim of this work is to study the numerical performance of the WN, ChN, AM and LDN methods to solve one-dimensional transport problems. Numerical comparisons for small and large thickness slab are presented.

© ]998 Elsevier Science Ltd.

Introduction

The key feature of deterministic methods to solve one-dimensional transport equations relies on the fact that they transform this equation into a set of linear differential equations (for instance SN, P., WN, ChN, AM and LDN) [1] or integral equation as in the FN method [2]. For the first type of approximations, the Laplace transform was well succeeded to generate analytical solutions for some approximations, well known as LTSN [3], LTPN [4], LTWN [5], LTChN [6], LTAN [7] and LTLDN [8] solutions. The procedure to attain these solutions consists in the application of the Laplace transform to the first order differential matrix equation, solution of the resulting algebraic linear system, which hinges on the complex parameter s, through the use of matrix inversion definition and reconstruction of the flux by the Heaviside expansion technique. It is important to notice 289

290

A.V. Cardona and M. T. de Vilhena

that besides the analytical character, these solutions posses a continuous, piecewise continuous and discrete dependency on the angular variable. In this work, our objective is to study the performance of these methods to solve transport problems in a slab, with small and large thickness, considering strong and weak absorption, constant and piecewise continuous incident angular flux at boundary. Furthermore, it is shown that for small thickness problems, all the methods here considered are equivalent regarding the accuracy of the results. On the other hand, for large thickness problems the LTCh. is the one which furnishes more accurate results for low order approximation. This conclusion is bolstered by the fact that the dominant eigenvalue of the LTChN matrix is the one which has better coincidence with the Case's discrete eigenvalue [9]. The paper is outlined as follows: in section 1, is shown a generic solution for the ChN, WN, PN and LDN approximations, in section 2, are depicted numerical simulations and comparisons and in section 3, is reported a discussion of the results encountered.

1. Analysis Let us consider the following transport problem in a homogeneous slab with anisotropic linear scattering:

t~

(X,l~)+9(x,p.) = 2 ~-lq~(x'l~')dP'+

I~J'!ll~9(x'l~')dl~' +S(x't~)'

(1)

where the coefficients fK are the parameters of the anisotropic linear model scattering; e(x,}~) and S(x,~) are respectively defined as the angular flux and the source term, subject to the boundary conditions: m(0, i~)= f(l~), if p. > 0

(la)

and e(a,-I~) = 0, if p. > 0.

(1 b)

In what follows, is briefly described how the mentioned approximations are obtained. Indeed, The W. [5] and ChN [6] approximations are achieved expanding the angular flux, in the angular variable, respectively in a truncated series of Walsh functions and Chebyshev polynomials, replacing these ansatz in equation (1) and taking moments. The AN approximation [7] is attained transforming the range of the angular variable of the

One-dimensional transport equation approximations

291

angular flux (p(x,l~) from the interval [-1,1] to [0,1] and evaluating the integral term of equation (1) by Gaussian quadrature scheme. The LDN equations [8] are obtained considering a linear expansion of the angular flux in the angular variable inside each node of the angular grid, replacing these expansions in equation (1) and taking moments. It turns out that all these approximations might be recasted in a matrix form as:

A;

(2)

~_(x)+ B:_~(x) : S(x),

where the matrices A and B and the vectors (l)(x) and S(x) depend strictly on the approximation considered, whose definition will come after. Henceforth, for the WN approximation, the elements of vectors (I)(x) and S(x) are respectively the coefficients of expansion of the angular flux and source in terms of Walsh functions in the angular variable. The entries of the matrices A and B ((2N+2) x (2N+2)) are respectively defined as [5]: I Dn_l,m_N_2, if n _ N + 1 an, m =

Dn_N_2,m_l, if n > N + 1 and m < N + 1

(3)

L0, at another case and

8n,0, if n_ N + 1 and m > N + 1

bn m = ~Sn m - 3 f l '

'

4

0

'

(3a)

'

L0, at another case where

(~n,m stands for the delta

of Kronecker and:

1/2, i f n = m gl3,m =

-2 -(k+2), if (n + m) mod2 = 2 k, k integer positive, 0, at another case

with (n + m) mod 2 denoting the mod 2 sum of the binary digits n and m [10].

(3b)

A . V . Cardona and M. T. dc Vilhena

292

For the ChN approximation, the elements of the vectors ¢(x) and _S(x) are respectively the coefficients of the expansion of the angular flux and source in Chebyshev polynomials in the angular variable. The components an, m and bn, m of the matrices A and =B ((N+I) x (N+I)) are respectively defined as [6]:

al.- ,1

(4)

an'm - 2 ( 2 - an+rn,3) and

I 8nrn _T..:-~, foam1 2-8m,1

bn,m= 8nm

IIUI

if n odd

~ )

3 f18rn2

' + ' , if n even 2 2 (n2 - 2 n + 3)

(4a)

For the AN approximation, the elements of the vectors ¢(x) and _S(x) are respectively the transformed angular flux and the source term at discrete directions. The matrix A is the identity matrix of order 2N, and the components bn,m of the matrix B (2N x 2N) are written as [7]: 1 l~n 1 I.Lm bn, m :

3flC0nI.Ln i f n = m _ N 2 f0mn , i f n - N = m 2 ~n

f0COn , if m _ N and n - N ~ m

(5)

2 I.Ln_N 3 fl0)rn-N tJ'm-N1 i f n < N , m > N a n d m - N . n 2 0, at another case The parameters

t~n and en denote the abscissas and weights of the Gaussian

quadrature scheme for the interval [0,1]. For the LDN approximations, the elements of the vectors ~(x) and S(x) are respectively the coefficients of the linear approximation of the angular flux and moments of the source term at each angular node. The components an, m and bn, m of the matrices A and B (2N x 2N) are respectively expressed as [8]:

One-dimensional transport equation approximations

293

~n, if n= m_ N aN, m --

-~,if n-N= m

(6)

3~,if m-N= n [0, at another case where P-n denotes direction at the center of the node, and 3 fl I~n I~m if n < N a n d m < N N ' fl , - 5 if n_N

&n,m bn, m =

(6a) 6nm '

f0 N

2fl i f n < N a n d m < N 3N 5 ' -

2 fl ILm , i f n > N a n d m < _ N

N4

Returning to the initial point of solving equation (2), notice that this equation has the well known solution [11 ]: ~_(x)= exp[-A-1Bxl~_(0) + foeXp[-A-1B(x-~) 1S(~)d~,

(7)

once the components of the column vector @(0) are known. Having established

an

analytical formulation for the exponential appearing in equation (7), the components of the vector ¢(0) for the boundary problem (1) can be readily obtained applying the boundary conditions (la) and (lb) and solving the resulting linear system [5-8]. To derive an analytical formulation for the exponential of the matrix A-1B, appearing in equation (7), let us solve the homogeneous version of equation (2), namely:

A ~x ~_(x) +B__¢(x) : O.

(8)

Applying the Laplace transform in equation (8), it turns out an algebric linear system, which has the solution:

A. V. Cardona and M. T. de Vilhena

294

(9)

= (sA +B )-1,(0).

Here the bar stands for the transformed component of the angular flux. The matrix (s A + B )-1 is readily obtained using the Trzaska's method [12]:

M

1

( s A + B ) -1 = ~ =

=

k=l S-

___kA

(9a)

Sk

Notice that, M = 2N + 2; N + 1; 2N; 2N respectively for the WN, ChN, AN and LDN approximations. The coefficients sk denote the eJgenvalues of the matrix A-1B and the matrices A are the ones resulting from the application of Trzaska's method. The =k

inversion is performed using the Heaviside expansion technique.

Following this

procedure, we get the insuing analytical expression for the exponential of the A-1B matrix: M

exp[ -A-1B x] = ~ eSkX=kA

(10)

k=1

Finally, replacing (10) into (7), it turns out the generic analytical solution: M

M

~_(x)= ~ eSk×_A__k~_(0)+ ~ SO eSk(X-~) S(~) d~ &___k k=l

(1I)

k=l

for all approximations here considered. 2. Numerical Comparisons

To check the performance of these methods, some one-dimensional standard transport problems were selected. (i) Initially let us consider an isotropic problem in homogeneous slab without source with constant and sectionally continuous incident flux at x = 0. The slab parameters are: a = 10 mfp and fo = 0.99; 0.5. The numerical results attained by the

One-dimensional transport equation approximations

295

LTWs, LTChl~, LTAo and LTLD6 methods to the scalar flux at x = 0, 5 and 10 mfp, defined as: 1

,l,(x) : j'_.~e(x, 1.0 d~.,

(12)

are presented in table 1.

Table 1. Numerical simulation to the scalar flux by LTW., LTChN, LTA. and LTLD. formulations for a = 10 mfp.

fo = 0.99 LTWs

Constant incident flux f(~,) = 1 at x = 0. x=0 x=5

x=10

1,810261

0.6547454

5.712323x10 -2

LTChll

1,812462

0.6555920

5.670869x10 -2

LTA8

1,810210

0.6554202

5.730714x10 -2

LTLDe

1,810212

0.6554537

fo = 0.5 LTWs

x=0 1,171573

x=5

5.730133x10 -2 x=10

4.555377x10 -3

2.603669x10 -5

LTCh11

1,184449

4.618824x10 -3

2.752705x10 -5

LTA6

1.171573

4.615734x10 -3

2.887815x10 -5

LTLD6

1.171573

LTLDe

re =0.5

0.7314091 x=0

0.1192809 x=5

1.042388x10-2 x=10

4.616003x10 -3 2.750903x 10-5 Incident flux f(F) = 1, if 0 < ~ < 0.5, or f(F) = 0, if 0.5 < p. < 1, at x = 0. x=0 x=5 x=10 fo = 0.99 LTWs 0.7311334 0.1185494 1.033988x10-2 LTChll 0.7456863 0.1204080 1.041136x10-2 LTAe 0.7280950 0.1151696 1.006703x10-2

LTWs

0.5600353

3.561281x10-4

1.970637x10-6

LTChll

0.5754633

3.699549x10-4

2.097712x10-6

LTA6

0.5595543

3.445009x10-4

2.291806x10-6

LTLD6

0.5600022

3.689565x10-4

2.084138x10-6

(ii) Let us now consider an isotropic problem in a homogeneous slab with thickness of 40 mfp and the following parameters: fo = 0.999; 0.95; 0.9; 0.8; f(I-0 = 1 and without source. The numerical results encountered by LTWs, LTCh11, LTA6 and LTLD6 methods, for the scalar flux at x = 20 mfp, are compared with the LTS12o solutions [13] and reported in table 2.

296

A . V . Cardona and M. T. de Vilhena

Table 2. Numerical comparisons of the LTVVN, LTChN, LTAN, LTLDN and LTSN results. 0.999

0.950

0.900

0.800

LTWs

0.5822254

7.523920x10 4

3.538306x10-5

6.866019x10-7

LTCh~

0.5832499

7.652124x10 4

3.636029x10 4

7.244695x10-7

LTAs

0.5837271

7.648984x10 4

3.632051x10 4

7.237517x10 q

LTLD6

0.5832071

7.646748x10 4

3.631408x10-5

7.216179x10-7

LTS~2o

0.5832087

7.647777x10 4

3.633341x10-5

7.238371x10-7

(iii) Finally, let us consider an isotropic and linearly anisotropic problem in homogeneous slab with the parameters: isotropic case, a = 100 and 1000 mfp, fo = 0.8, f(l~)--1

and a unitary source at the slab region ( a - l ) _ < x_< a. For the anisotropic

problem, a = 100 mfp, fo = 0.99, fl = 0.8 and f(l~) -- 1. The numerical results for the scalar flux, obtained by the LTW5, LTCh11, LTA6 and LTLD6 methods, at x = 0, a/2, a are shown in table 3. The reason for considering the anisotropic problem relies on the fact that the LTS12o solution for this problem, considered exact, is available.

Table 3. Numerical comparisons of the L'FVVN,LTChN, LTAN, LTLDN and LTSN results. a = 100 mfp LTWs

x=0 1.381966

Isotropic Problems x = 50 1.176230x10-16

x = 100 1.826786

LTChll

1.390848

1.367732x10-16

1.791950

LTA6

1.381965

1.377307x10-16

1.826792

1.361093x10-16 x = 500 0.386264x10-1~

1.826529 x = 1000 1.826786 1.791950

LTLD6

1.381966

a : 1000 mfp LTW5

x=0 1.381966

LTChll

1.390848

LTA6

1.381965

LT LDe a = 100 mfp LTWs LTChll

1.988780x10 -155 1.995946x 10-155

1.826529

1.381966 1.730356x 10-155 Linearly Anisotropic Problem x=50 x=0 1.648013 3.174995 x 10-2 1.649783 3.292111 x 10-2

1.826786

2.257939 x 10-4 2.421066 x 10.4

x = 100

LTA6

1.645990

3.298146 x 10~

0.724299 x 104

LTLD6

1.656435

2.711847 x 102

1.591500 x 104

LTSI~o

1.645988

3.290443 x 10-2

2.444444 x 104

One-dimensional transport equation approximations

297

For one side, examining the results for the isotropic case appearing in table 3, we promptly realize a discrepancy between the LTChl~ solutions and the remaining ones. On the other hand, looking at the results for the linearly anisotropic case, we notice a better coincidence of the L T C h , solution with the LTS12o results. These facts, suggest that for large thickness slabs the LTChl~, among the methods herein studied,

is the one which

generates more accurate results This thing is bolstered by the interesting result that the dominant eigenvalue of the LTChll matrix is the one which better approximated the discrete eigenvalue of the Case's solution [9].

Table 4. Dominant eigenvalue comparison and conditioning number. fo -- 0.99

Eigenvalue

Error (%)

Conditioning Number

LTWs

0.1729

4.5 x 10-5

45.83

LTChll

0.1725

1.0 x 104

38.69

LTA6

0.1725

1.0 x l 0 -6

164.33

LTLD6 fo = 0.8

0.1725 Eigenvalue

1.0 x 104 Error (%)

75.40 Conditioning Number

LTW5

0.7137

1.8 x 10-3

11.41

LTChl~

0.7104

6.0 x l 0 -6

9.66

LTA8

0.7104

6.0 x l 0 -6

40.24

LTLD6 fo = 0.5

0.7107 Eigenvalue

1.4 x l 0 -4 Error (%)

18.63 Conditioning Number

LTWs

0.9830

2.2 x 10-2

8.64

LTChll

0.9570

4.4 x l 0 -4

7.48

LTA8

0.9589

1.2 x 10-3

30.22

LTLD6 fo = 0.1

0.9671 Eigenvalue

8.3 x l 0 -3 Error (%)

14.06 Conditioning Number

LTWs

1.1081

9.8 x l 0 -2

8.08

LTChll

1.0067

6.7 x l 0 -3

7.51

LTAe

1.0297

2.9 x l 0 -2

28.63

LTLD6

1.0655

6.2 x 10-2

13.21

Our previous comments as well the conditioning number for the matrices of all methods considered are illustrated in table 4, which also increase our confidence in the perform of the LTChN formulation.

A. V. Cardona and M. T. de Vilhena

298

3. Conclusions The foregoing results point to the following conclusions. First of all, it was possible to develop an analytical generic method to solve the WN, ChN, AN and LDN approximations of the one-dimensional transport equation. It is also worth to point out that this solution shall be valid for all approximations of first type of the one-dimensional transport equation, for instance, the SN and PN approximations [3-4]. The analytical character of this solution allow us to change basis of the space solution, in such manner that it will become appropriated to handle problems with large thickness without discretisation of the domain. Analyzing the results appearing in previous tables, we realize that for small thickness slab, roughly speaking due to the low order of approximation considered, the methods are equivalent regarding the accuracy of the results, meanwhile, for large thickness slab, the LTChll was the one which generates best results. This fact is reinforced by the best coincidence of its dominant eigenvalue with the discrete Case's eigenvalue [9]. Concluding, we would like to emphasize that our attention is now focused to infer a conclusive answer to the question which method is the best approximation of onedimensional transport problem in a large thickness slab. To reach this goal, we will employ the recursive method [13-14] to invert the LTWN, LTChN, LTAN and LTLDN matrices, which will allow us to solve problems with increasing N (N = 120) and we will include in the analysis the SN and PN approximations. We hope that following this procedure, it will be possible to definitively confirm the efficiency of the LTChN formulation to solve one-dimensional transport problems for large thickness slab. Our attention is now focused in this direction.

Acknowledgments The first and second authors are gratefully indebted to CNEN (Comiss~o Nacional de Energia Nuclear) and CNPq (Conselho Nacional de Desenvolvimento Tecnol6gico e Cientifico) for having partially supported this work.

One-dimensional transport equation approximations

299

References

[1] CARDONA, A. V.: A Generic Method of Analytical Solution for Some Approximations of the Linear Transport Equation. PhD thesis, Mechanical Enginnering, Federal University of Rio Grande do Sul, Porto Alegre, Brazil (1996). [2]

GARCIA, R. D. M.: A Review of the Facile (FN) Method in Particle Transport Theory. Transport Theory and Statistical Physics 14(1985) 391-435.

[3]

BARICHELLO, L. B.; VILHENA, M. T.: A General Approach to One Group One Dimensional Transport Equation. Kemtechnik 58(1993) 182-184.

[4] VILHENA, M. T.; STRECK, E. E.: An Approximated Analytical Solution of the OneGroup Neutron Transport Equation. Kemtechnik 57(1992) 196-198. [5] CARDONA, A. V. & VILHENA, M. T.: A Solution of Linear Transport Equation using Walsh Function and Laplace Transform. Annals of Nuclear Energy 21(1994) 495505. [6]

CARDONA, A. V. & VILHENA, M. T.: A Solution of Linear Transport Equation using Chebyshev Polynomials and Laplace Transform. Kerntechnik 59(1994) 278-281.

[7]

CARDONA, A. V.; VlLHENA, M. T.: Analytical Solution for the AM Approximation. Progress in Nuclear Energy. In press.

[8]

DE BARROS, R. C.; CARDONA, A. V. & VILHENA, M. T.: Analytical Numerical Methods Applied to Linear Discontinuous Angular Approximations of the Transport Equation in Slab Geometry. Kerntechnik. In press.

[9] CASE, K. M.: Elementary Solutions of the Transport Equation and Their Applications. Annals of Physics 9(1960) 1-23.

[10] CORRINGTON, M. S.: Solution of Differential and Integral Equations with Walsh Functions. IEEE Transactions of Circuit Theory 20(1973) 470-476. [11] SCHILLING, R. J. & LEE, H.: Engineering Analysis - A Vector Space Approach.

John Wiley and Sons, New York (1988). [12] TRZASKA, Z.: An Efficient Algorithm for Partial Fraction Expansion of the Linear Matrix Pencil Inverse. Joumal of the Franklin Institute 324(1987) 465-477.

300

A.V. Cardona and M. T. de Vilhena

[13] BRANCHER, J. D.; CARDONA, A. V.; VILHENA, M. T.: A Recursive Method to Invert the LTSN Matrix. Progress in Nuclear Energy. Submitted for publication. [14] BRANCHER, J. D.; SEGATTO, C. F.; VILHENA, M. T.: The One-Dimensional LTSN Formulation for High Order Anisotropy. In preparation.