Discontinuous finite element formulations for neutron transport in spherical geometry

Annals of Nuclear Energy 64 (2014) 244–255

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Discontinuous finite element formulations for neutron transport in spherical geometry Mehmet Mercimek ⇑, H. Atilla Özgener Energy Institute, Istanbul Technical University, Ayazag˘a Campus, 34469 Maslak, Istanbul, Turkey

a r t i c l e

i n f o

Article history: Received 20 June 2013 Received in revised form 9 October 2013 Accepted 11 October 2013 Available online 7 November 2013 Keywords: Discrete ordinates method Discontinuous finite element method Diamond differencing Implicit method

a b s t r a c t We have developed the linear and quadratic Galerkin discontinuous finite element methods for the solution of both time-independent and time-dependent spherical geometry neutron transport problems. Discrete ordinates method is used for the angular discretization while the implicit method is utilized for temporal discretization in time-dependent problems. In order to assess the relative performance of the newly developed linear and quadratic discontinuous finite element spatial differencing methods relative to the previously developed linear discontinuous finite element and diamond difference discretizations, a computer code is developed and numerical solutions of the neutron transport equation for some benchmark problems are obtained. These numerical applications reveal that the newly developed quadratic discontinuous finite element method produces the most accurate results while the newly developed linear discontinuous finite element method follows as the second best discontinuous finite element method. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The discontinuous finite element method (DFEM) is among the most flexible numerical methods in discrete ordinates (SN) formulations of neutral particle transport. The method was first introduced for the solution of the time-independent neutron transport equation (Reed and Hill, 1973). The developments in general DFEM have been recently reviewed (Cockburn et al., 2000; Cockburn, 2001). DFEM uses piecewise polynomial spaces similar to the (continuous) finite element method (FEM) but with relaxed continuity conditions at interelement boundaries. Trial functions can be chosen so that the field variable, its derivative or both are discontinuous across interelement boundaries. The method includes as its subsets both the FEM and finite difference methods (FDM). It has been reported that DFEM has the advantages of both the FDM and FEM (Li, 2006). Several computer programs have been developed for neutron transport using the linear discontinuous finite element method (LDFEM) (Reed et al., 1973, 1977; Seed et al., 1977, 1978; Wareing et al., 1996). ONETRAN and TIMEX are examples of these codes for the solution of time-independent and time-dependent transport, respectively (Hill, 1975, 1977; Hill et al., 1976). Recently, the LDFEM formulation, similar to the formulation used in TIMEX, has been employed in spherical geometry (Hong et al., 2010). In another study, piecewise LDFEM is applied to the two-dimensional ⇑ Corresponding author. Tel.: +90 538 641 90 36. E-mail addresses: [email protected] (M. Mercimek), [email protected] (H. Atilla Özgener). 0306-4549/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2013.10.012

cylindrical geometry and accurate results are reported (Bailey et al., 2009). Higher-order DFEM has been extensively investigated in other disciplines (Hesthaven and Warburton, 2008). But, in neutron transport, limited research has been carried out using elements of higher degree. For example, a higher order DFEM has been developed for the time-independent transport equation by using hierarchical basis functions (Wang, 2009). In another study, higher order DFEM was used in both angular and spatial differencing for the solution of time-independent, pure absorber problems, in spherical geometry (Machorro, 2007). A comparison of linear and quadratic DFEM for triangular lattice discrete ordinates calculations was also presented (Chang and Warsa, 2007). The first objective of this study is to reformulate and implement LDFEM for time-independent and time-dependent neutron transport problems in spherical geometry. The time dependent formulation will be based on the prompt neutron transport equation and the delayed neutron effects will not be treated. Our LDFEM formulation will be henceforth called as LD2 while the spherical geometry LDFEM formulation of TIMEX will be called as LD1 (Hill, 1975, 1977; Hill et al., 1976; Hong et al., 2010). The second objective of this study is to develop a higher-order DFEM formulation for the solution of neutron transport problems in spherical geometry both for the time-independent and time dependent variety. As a higher-order DFEM, we employ the quadratic discontinuous finite element method (QDFEM). In this paper, we will first present the derivations of linear systems of equations for LD1, LD2 and QDFEM for the solution of time-independent and time-dependent prompt neutron transport

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

equation in spherical geometry. We will use the classical Galerkin method in which the weight functions are taken to be the same as the shape functions which are Lagrange type polynomials with compact support in LD2 and QDFEM. We have developed a discrete ordinates computer code, SPDOT, in which the spatial discretization can be carried out using diamond difference (DD), LDFEM (LD2) or QDFEM approaches. Also, we have developed another version of the SPDOT for LD1 calculations of this study. Since DD formulation is given elsewhere (Lewis and Miller, 1984), it would not be repeated here. Implicit method is used for the temporal differencing in timedependent problems. Both codes are run for several time-independent and time-dependent benchmark problems to assess the performance of various spatial discretization methods in the solutions of neutron transport problems by the discrete ordinates method. The remainder of this paper is organized as follows. In Section 2, LDFEM and QDFEM formulations are described. Section 2.1 describes the discrete ordinate spherical geometry transport equation and the weighted residual forms of this equation for spatial discretization methods LD1, LD2 and QDFEM. LD2 and QDFEM approximations are given in Section 2.2 while LD1 approximation is given in Section 2.3. Finite element formulations of LDFEM and QDFEM are given in Sections 2.4 and 2.5 respectively. In Section 2.6, time differencing scheme is applied to the derived equations in previous sections. Section 2.7 describes the boundary condition applied at the centre of the sphere. In Section 3, numerical results of DD, LDFEM and QDFEM are compared. Section 3.1 is devoted to time-independent criticality problems while Section 3.2 is for the analysis of time-dependent behavior of the methods. In the final section, some concluding remarks are given. 2. Discontinuous finite element formulation First, we describe the linear and quadratic discontinuous finite element methods for the spatial discretization of the time-independent within group neutron transport equation in spherical geometry using the discrete ordinates form. Details of the sweep of the space-angle mesh in the spherical geometry discrete ordinates neutron transport are given elsewhere (Lewis and Miller, 1984).

lm

i @ðr 2 wm Þ 2r h þ amþ12 wmþ12 ðrÞ  am12 wm12 ðrÞ þ r2 rðrÞwm ðrÞ @r wm

¼ r 2 qm ðrÞ

wmþ1 ðrÞ ffi 2wm ðrÞ  wm1=2 ðrÞ 2

rðnÞ ¼ r i þ

Dr i n 2

ð4Þ

where r i ¼ ðr i1=2 þ riþ1=2 Þ=2 and Dri = ri+1/2  ri1/2 represent the midpoint and interval size of the ith mesh interval which is assumed to be homogeneous. Fig. 1 illustrates the arrangement of angular flux node points in (i, m)th cell of the spherical domain. Using Eqs. (3) and (4), Eq. (2) is written with respect to local variable as:

  i 2l m d 2 2rðnÞ h ½r ðnÞwm ðnÞ þ 2amþ1 wm ðnÞ  amþ1 þ am1 wm1 ðnÞ 2 2 2 2 wm Dr i dn þ ri r 2 ðnÞwm ðnÞ ¼ r 2 ðnÞqm ðnÞ; 1 6 n 6 1 ð5Þ 2.1.1. Weighted residual form for the LD2 and QDFEM formulations LD2 and QDFEM formulations are based directly on Eq. (5). Since wm(n) and wm1/2(n) are approximate solutions then they cannot be expected to satisfy Eq. (5) at all points in the mesh interval. We require that these approximate solutions satisfy Eq. (5) only in an integral sense by integrating Eq. (5) after multiplication with a weight function, w(n), and arrive at the weighted residual (Galerkin) form:

2lm Dr i

Z

1

d 2 ½r ðnÞwm ðnÞdn dn   Z 1 Z 1 2 þ 2amþ1 wðnÞrðnÞwm ðnÞdn  am wðnÞrðnÞwm1 ðnÞdn 2 2 wm 1 1 Z 1 wðnÞr 2 ðnÞwm ðnÞdn þ ri wðnÞ

1

Z

The spherical geometry time-independent within group neutron transport equation in conservation form is written as:

where

We omit group indices for the sake of clarity. In Eq. (1), angular variable l denotes the radial component of particle direction, spatial variable r is the distance from the origin, w is the angular neutron flux, r is the macroscopic total cross section and q is the emission density which includes the inhomogeneous (or fixed) source, scattering source and fission contribution to the group source. In order to develop the discrete ordinates equations, the angular domain l 2 [1, +1] is discretized into M (m = 1, 2, . . . , M) quadrature points lm with weights wm. We introduce the spherical geometry discrete ordinates form of transport equation by introducing the angular differencing coefficients, am+1/2 and am1/2 (Lewis and Miller, 1984) into Eq. (1):

ð3Þ

In order to make the derivation easier, the local spatial variable n 2 [1, +1], is introduced as:

¼

ð1Þ

ð2Þ

Diamond differencing (DD) in angle is used to relate the edge (wm+1/2(r), wm1/2(r)), and cell-centered (wm(r)) angular fluxes:

2.1. The discrete ordinates spherical geometry transport equation and the weighted residual forms

  @ ð1  l2 Þw @ðr 2 wÞ l þr þ r 2 rðrÞwðr; lÞ ¼ r2 qðr; lÞ @l @r

245

1 1

wðnÞr2 ðnÞqm ðnÞdn

1

ð6Þ

am ¼ amþ1=2 þ am1=2

ð7Þ

Applying integration by parts to the first term on the left hand side of Eq. (6), we obtain:

Fig. 1. A unit cell for DFEM formulations.

246

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255



 Z 1 2lm dw 2 wð1Þr2iþ1 wm ð1Þ  wð1Þr 2i1 wm ð1Þ  r wm ðnÞdn 2 2 Dr i 1 dn   Z 1 Z 1 2 2amþ1 wðnÞrðnÞwm ðnÞdn  am wðnÞrðnÞwm1 ðnÞdn þ 2 2 wm 1 1 Z 1 wðnÞr 2 ðnÞwm ðnÞdn þ ri 1

Z

¼

1

wðnÞr2 ðnÞqm ðnÞdn

Using the transformation to local coordinates, Eq. (4), multiplying with a weight function and integrating by parts we arrive at the weighted residual form for the special direction:



  Z 1 2 dwðnÞ wð1Þw1 ð1Þ  wð1Þw1 ð1Þ  w1 ðnÞ dn þ ri 2 2 2 dn Dr i 1 Z 1  wðnÞw1 ðnÞdn

1

ð8Þ Eq. (8) constitutes the final weighted residual form from which LD2 and QDFEM formulations follow. 2.1.2. Weighted residual form for the LD1 formulation For the LD1 formulation Eq. (5) is further simplified by the additional approximations:

wm ðnÞ ffi wi;m ðnÞ ffi

2

1 1

wðnÞq1 ðnÞdn

ð15Þ

2

1

Eq. (15) is used in all (LD1, LD2, QDFEM) formulations. 2.2. The linear (LD2) and quadratic (QDFEM) discontinuous finite element approximations The finite element approximation is introduced by assuming

1 ðw þ wi1=2;m Þ 2 iþ1=2;m

wm1=2 ðnÞ ffi wi;m1=2

¼

Z

ð9Þ ð10Þ

When Eqs. (9) and (10) are utilized in the second term on the lefthand side of Eq. (5), we obtain:

wm ðnÞ ¼

Nþ1 X N hj ðnÞwjm

ð16Þ

j¼1

wm1=2 ðnÞ ¼

Nþ1 X N hj ðnÞwjm1=2

ð17Þ

j¼1

  2l m d 2 2rðnÞ h ½r ðnÞwm ðnÞ þ amþ12 wiþ12;m þ wi12;m wm Dri dn   i  amþ1 þ am1 wi;m1 þ ri r 2 ðnÞwm ðnÞ ¼ r 2 ðnÞqm ðnÞ 2

2

2

qm ðnÞ ¼

2l m Dr i

ð11Þ

1

d wðnÞ ½r 2 ðnÞwm ðnÞdn dn   Z 1 2 þ amþ12 wiþ12;m þ wi12;m wðnÞrðnÞdn wm 1  Z 1 Z 1 wðnÞrðnÞdn þ ri wðnÞr 2 ðnÞwm ðnÞdn am wi;m1 Z

¼

1

1

wðnÞr2 ðnÞqm ðnÞdn

ð12Þ

1 1

wðnÞr 2 ðnÞqm ðnÞdn

Eq. (13) constitutes the final weighted residual form from which LD1 formulation follows. For starting the calculations in spherical geometry discrete ordinates method, we need to derive the weighted residual form also for the starting direction that is l1/2 = 1. Transport equation for the radially inward direction is:

dr

þ ri w1 ðrÞ ¼ q1 ðrÞ 2

That

is

w1m ¼ wi1=2;m

and

w2m ¼ wi;m only for quadratic elements. The finite element trial functions have compact support, that is:

i; j ¼ 1; . . . ; N þ 1

ð19Þ

2ði  1Þ N

ð20Þ

2

The finite element trial functions which are linear and quadratic polynomials over the element are given in Table 1. For LD2 and QDFEM, as weight functions in Eq. (8), we choose the finite element trial functions. That is we take

wðnÞ ¼ hk ðnÞ;

k ¼ 1; . . . ; N þ 1

ð21Þ

Inserting Eqs. (16)–(18) and (21) into Eq. (8) and using Eq. (19), we obtain the following equation for k = 1, . . . , N + 1:

ð13Þ

2

coordinates.

1

1

dw1

local

wNþ1 ¼ wiþ1=2;m for both linear and quadratic elements while m

ni ¼ 1 þ

  Z 1 2lm dw 2 wð1Þr 2iþ1 wm ð1Þ  wð1Þr 2i1 wm ð1Þ  r wm ðnÞdn 2 2 Dr i 1 dn   Z 1  Z 1 2 amþ12 wiþ12;m þ wi12;m wðnÞrðnÞdn  am wi;m1 wðnÞrðnÞdn þ 2 wm 1 1 Z 1 2 wðnÞr ðnÞwm ðnÞdn þ ri



to

where

Applying integration by parts to the first term on the left hand side of Eq. (12), we obtain:

¼

global

N

1

Z

ð18Þ

where N = 1 corresponds to linear and N = 2 corresponds to quadratic discontinuous finite element approximations. Here we introduced a change of notation as we changed the formulation from

hj ðni Þ ¼ dij

1

2

N

hj ðnÞqjm

j¼1

LD1 formulation is based on the approximate equation, Eq. (11). For the LD1 formulation, we require wm(n) to satisfy Eq. (11) only in an integral sense by integrating of Eq. (11) after multiplication with a weight function w(n) and arrive at the weighted residual form:

Z

Nþ1 X

ð14Þ

" # Nþ1 X 2lm 1 j 2 dkðNþ1Þ r 2iþ1 wNþ1  d r w  k w k1 i1 m k;j m m 2 2 Dr i j¼1 " # N þ1 N þ1 X X 2 j j 2amþ1 lk;j wm  am lk;j wm1 þ 2 wm 2 j¼1 j¼1 þ ri

Nþ1 Nþ1 X X mk;j wjm ¼ mk;j qjm j¼1

ð22Þ

j¼1

Table 1 Finite element trial functions. N

h1(n)

h2(n)

h3(n)

1

1n 2 n2 n 2

1þn 2

–

2

1  n2

n2 þn 2

247

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

The matrices K, L, M for LD2 and QDFEM are defined as:

Z

kk;j ¼

1 1

Z

lk;j ¼

1

Explicit form of these matrices for the linear finite elements is presented below:

dhk 2 r ðnÞhj ðnÞdn dn

ð23Þ

"



K¼

hk ðnÞrðnÞhj ðnÞdn

1 2

ð24Þ

"

1

mk;j ¼

Z



1 2

hk ðnÞr ðnÞhj ðnÞdn

ð25Þ

 12  12

2

  r2 3 r 2 Dr 2 Dr 2    2i  ri D6 ri  24i  2i þ ri D6ri  24i 6 7 K ¼ 4 2  r2  5 Dr 2i Dr 2i ri ri Dr i r i Dr i i  þ þ þ 2 6 24 2 6 24

M¼

2 6 L¼4

2r i 3

 D6ri



symmetric

2 6 M¼4

2r 2i 3

i

3



2r i 3



 ri D3ri þ

Dr 2i 15

þ

Dr i 6



ð27Þ

r2 i

3

2r2

symmetric

þ

Dr 2i 60



3



þ ri D3 ri þ

i

3

ð28Þ

i

15

  2r2  r2  3 r2 3Dr 2 Dr 2 Dr 2 r i Dr i i Dr i  Dr i i  3i  4r15  40i  30i  ri30  40i  2i þ 1130 6 7 2r2     2r2 7 Dr 2i Dr 2i 7 2r i Dr i 8r i Dri 2r i Dri i i       7 3 5 10 15 3 5 10 7  r2      5 2 2 2 2 2   D r 2 3 D r 3 D r r r    Dr i 4r i Dr i 11r i Dr i i i i i  6i þ ri30   40i þ  þ  3 15 30 2 30 40

2

2 6 6 L¼6 6 4

4r i 15

 D10ri





2r i 15

 D15ri  



3 ri  15  7 7 2r i Dr i 7 þ 15 15 7  5 4r i Dr i þ 15 10

symmetric

6 6 6 M¼6 6 4

4r 2i 15



 ri D5ri þ

3Dr2i 70

6 K ¼4

 12  23 2 3

0

 16

2 3

1 6

7  23 5

ð37Þ

1 2

4 15

6 M¼4

3

symmetric

2 15 16 15

1  15 2 15 4 15

3 7 5

ð38Þ

In LD1 formulation, weight functions are different for lm > 0 and lm < 0. That is:

w1 ðnÞ ¼ 1; w2 ðnÞ ¼ h2 ðnÞ w1 ðnÞ ¼ 1; w2 ðnÞ ¼ h1 ðnÞ

lm > 0 lm < 0

ð39Þ ð40Þ

For this reason it proves to be convenient to define a variable n such that:



1; 2;

lm > 0 lm < 0

ð29Þ

n¼

ð30Þ

Inserting Eqs. (16), (18), (39) and (40) into Eq. (13), we obtain the following equation:

ð41Þ



16r i 15

2

ð36Þ

2.3. The linear discontinuous finite element approximation (LD1)

7 5 Dr 2

and for quadratic finite elements (QDFEM), these matrices assume the form: 6 6 6 K ¼6 6 4

2



7 5

#

For quadratic elements, these matrices assume the form:

2

3

r

1 3 2 3

symmetric



ð26Þ

ð35Þ

1 2 2 3

1

Explicit form of these matrices for linear finite elements (LD2) is presented below:

#

 2r2 i

15



i Dr i  2r15 þ

16r2 i

15

þ

4Dr 2i 105

symmetric

Dr2i 70





 r2 Dr2   15i  140i

3

7 7 7 þ þ 7 15 7 4r2  5 2 3 D r  ri Dr i i i þ 5 þ 70 15 2r2 i

2r i Dr i 15

Dr2i 70

" # Nþ1 X 2lm 2 1 j 2 2 ð1  d2k d2n Þr iþ1 wm  ð1  d2k d1n Þr i1 wm  kk;j wm 2 2 Dr i j¼1 " # 2 2 X X 2 amþ12 lk wjm  am lk w m1=2 þ ri mk;j wjm þ wm j¼1 j¼1 ¼

2 X mk;j qjm

ð42Þ

j¼1

ð31Þ

where For the treatment of the special direction (m = 1/2) for LD2 and QDFEM, we insert Eqs. (16), (17) and (21) into Eq. (15) and use Eq. (19) to obtain the equations for the radially inward direction. That is for k = 1, . . . , N + 1, we have:

"

N þ1 X

~ k;j qj1 m

~ k;j ¼ m

1 Z 1

dwk 2 r ðnÞhj ðnÞdn dn

ð44Þ

wk ðnÞr 2 ðnÞhj ðnÞdn

ð45Þ

1

ð33Þ

The square matrix L defined for LD2 and QDFEM becomes a column matrix of dimension two in the LD1 formulation and is defined as:

ð34Þ

lk ¼

where

~ ¼ k k;j

mk;j ¼

1

2

j¼1

Z

Z

kk;j ¼

ð32Þ

ð43Þ

The matrices K and M are defined for LD1 as:

#

Nþ1 Nþ1 X X 2 ~ wj þ r ~ k;j wj1 k m  dkðNþ1Þ wNþ1  dk;1 w11  1 i k;j 1 2 2 Dr i 2 2 j¼1 j¼1

¼

 . 2  m1=2 ¼ w1 w 2 m1=2 þ wm1=2

1

1 Z 1 1

dhk hj ðnÞdn dn hk ðnÞhj ðnÞdn

Z

1

1

wk ðnÞrðnÞdn

ð46Þ

248

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

Explicit form of these matrices is presented below:

" K ¼ ð1Þ

r2

n1

i

2

" L¼

0 

 ri D6ri þ

2r i r i þ ð1Þn1 D6ri

 r2

24

2

i

0 

þ ri D6 ri þ

Dr 2i



ð47Þ

24

#

2 M¼4

Dr 2i

ð48Þ

r 2i  ri D3 ri þ r 2

ð1 þ dn2 Þ 3i þ ð1 þ 3dn2 Þ

Dr 2i 12

Dr2i 60

4amþ1 2lm 2 k1;j þ l1;j þ ri m1;j ; j ¼ 1; 2 Dr i wm   4amþ1 2lm 2 2 r 1 d2j k2;j þ l2;j þ ri m2;j ; j ¼ 1; 2 a2;j ¼ Dr i iþ2 wm 2 2 X 2lm 2 2am X r 1 w dk1 þ lk;j wjm1 þ mk;j qjm ; k ¼ 1; 2 fk ¼ Dri i2 b wm j¼1 2 j¼1

a1;j ¼ 

#

r 2i þ ri D3 ri þ

þ dn2 ð1Þn1 ri D3 ri 

r 2

ð1 þ dn1 Þ 3i þ ð1 þ 3dn1 Þ

The starting direction (l1/2 = 1) equations are obtained similarly after we insert Eqs. (16), (18) and (39), (40) into Eq. (15):

" # 2 2 2 X X X 2 j 2 1 ~ ~ k;j wj1 ¼ ~ k;j qj1 kk;j w1 þ ri m  dk1 w1  w1  m 2 2 Dr i 2 2 2 j¼1 j¼1 j¼1 Z 1 dwk ~ ¼ hj ðnÞdn k k;j 1 dn Z 1 ~ k;j ¼ wk ðnÞhj ðnÞdn m

ð50Þ ð51Þ

þ dn1 ð1Þn1 ri D3ri 

5

ð49Þ



a2;2

a2;1 a1;1

"

w2m w1m

#

 ¼

f2

 ð59Þ

f1

ð52Þ 4amþ1 2lm 2 k2;j þ l2;j þ ri m2;j ; j ¼ 1; 2 Dr i wm  4amþ1 2l  2 l1;j þ ri m1;j ; j ¼ 1; 2 a1;j ¼  m r 2i1 d1j þ k1;j þ 2 Dr i wm 2 2 X 2l 2am X lk;j wjm1 þ mk;j qjm ; k ¼ 1; 2 fk ¼  m r2iþ1 wb dk2 þ 2 2 Dr i wm j¼1 j¼1

a2;j ¼  After taking the integrals, the following explicit form of these matrices are obtained for LD1:

e ¼ K

0

" f M¼

0

# ð53Þ

 12  12 1 1 2 3

ð54Þ

1 3

The only difference between LD1 and LDFEM in Hill’s study is that a different special direction method (step-start) is used in the TIMEX code to calculate the starting angular fluxes than the above mentioned method which has been used in all spatial discretization methods in this study. Experience has shown that the step-start method is less accurate than using Eq. (45) but the loss in accuracy is quite small (O’Dell and Alcouffe, 1987). For the LD1 formulation, it is worth to note that if the weight functions are chosen according to Eq. (21) as in the LD2 formulation instead of using Eqs. (39) and (40), the calculated effective multiplication factor and the group fluxes would be identical since the spanned vector space is the same. This assertion has been confirmed by numerical experimentation in this study. 2.4. The linear discontinuous finite element equations 2.4.1. LD2 Formulation For lm > 0 in LD2 method, and k = 1 in Eq. (22), we take the angular flux w1m ¼ wb in the second term with the Kronecker delta where wb represents the angular flux at the second node of the previous ((i  1)st) element in the marching scheme. Thus the discontinuity is introduced. For lm > 0 and k = 2, we leave Eq. (22) unchanged. The resulting linear system is:

a1;1 a2;1

a1;2 a2;2

ð60Þ ð61Þ ð62Þ

# For the special direction, the analogous equations are:





ð58Þ

For lm < 0 and k = 2 in Eq. (22), in the first term with the Kronecker delta, we take the angular flux w2m ¼ wb where wb represents the angular flux at the first node of the previous ((i + 1)st) element in the marching scheme. For lm < 0 and k = 1, we leave Eq. (22) unchanged. The resulting linear system is:

a1;2

1

"

ð57Þ

3

Dr 2i 12

Dr 2i 60

ð56Þ

"

w1m w2m

#

 ¼

f1 f2

 ð55Þ

~2;2 a ~1;2 a

2 3 " #  w21 ~ ~2;1 a 4 2 5 ¼ f2 1 ~f 1 ~1;1 a w1

ð63Þ

2

2 ~ ~ 2;j ; j ¼ 1; 2 k2;j þ ri m Dr i 2 ~1j Þ þ ri m ~ 1;j ; j ¼ 1; 2 ¼ ðd1j þ k Dr i

~2;j ¼ a

ð64Þ

~1;j a

ð65Þ

fk ¼

2 X 2 ~ k;j qj1 ; wb dk2 þ m Dr i 2 j¼1

k ¼ 1; 2

ð66Þ

2.4.2. LD1 Formulation Similarly for lm > 0 in LD1 method, and k = 1 in Eq. (42), we take the angular flux w1m ¼ wb in the second term with the Kronecker deltas. Thus the discontinuity is introduced. For lm > 0 and k = 2, we leave Eq. (42) unchanged. The resulting linear system is:

 2amþ1 2lm  2 2 riþ1 d2j k1;j þ l1 þ ri m1;j ; 2 Dr i wm  2amþ1 2lm  2 2 ¼ riþ1 d2j k2;j þ l2 þ ri m2;j ; 2 Dr i wm

a1;j ¼

j ¼ 1; 2

ð67Þ

a2;j

j ¼ 1; 2

ð68Þ

fk ¼

2 X 2lm 2 2am  r i1 wb dk1 þ lk wm1=2 þ mk;j qjm ; 2 Dr i wm j¼1

k ¼ 1; 2

ð69Þ

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

For lm < 0 and k = 1 in Eq. (22), in the first term with the Kronecker delta, we take the angular flux w2m ¼ wb . For lm < 0 and k = 2, we leave Eq. (22) unchanged. The resulting linear system is:

 2amþ1 2l m  2 2 r 1 d1j þ k1;j þ l1 þ ri m1;j ; j ¼ 1; 2 Dr i i2 wm  2amþ1 2l  2 l2 þ ri m2;j ; j ¼ 1; 2 a2;j ¼  m r2i1 d1j þ k2;j þ 2 Dr i wm 2 X 2l 2am  lk wm1=2 þ mk;j qjm ; k ¼ 1; 2 fk ¼  m r 2iþ1 wb dk1 þ 2 Dr i wm j¼1

a1;j ¼ 

ð70Þ

2 ~ Þþrm ðd1j þ k j ¼ 1; 2 1j i ~ 1;j ; Dr i 2 ~2j Þ þ ri m ~2;j ¼ ~ 2;j ; j ¼ 1; 2 a ðd1j þ k Dr i 2 X 2 ~ k;j qj1 ; k ¼ 1; 2 wb dk1 þ fk ¼ m 2 Dr i j¼1

For the special direction, the analogous equations are:

2

~3;3 a 6~ a 4 2;3

~3;2 a ~2;2 a

~1;1 a

~1;2 a

ð71Þ

2 ~ ~ 3;j ; j ¼ 1; 2; 3 k3;j þ ri m Dr i 2 ~ ~2;j ¼ ~ 2;j ; j ¼ 1; 2; 3 a k2;j þ ri m Dr i 2 ~ Þþrm ~1;j ¼ a ðd1j þ k j ¼ 1; 2; 3 1;j i ~ 1;j ; Dr i 3 ~f ¼ 2 w d þ m ~ k;j qj1 ; k ¼ 1; 2; 3 m k k3 j¼1 Dr i b 2 ~3;j ¼ a

ð72Þ

ð73Þ

ð75Þ

For lm > 0 and k = 1 in Eq. (22), in the second term with the Kronecker delta, we take the angular flux w1m ¼ wb where wb represents the angular flux at the third node of the previous ((i  1)st) element in the marching scheme. Thus the discontinuity is introduced. For lm > 0 and for k = 2 or 3, we leave Eq. (22) unchanged. The resulting linear system is:

a1;2 a2;2

32 1 3 2 3 wm a1;3 f1 76 2 7 6 7 a2;3 54 wm 5 ¼ 4 f2 5

a3;1

a3;2

a3;3

4amþ1 2l 2 a1;j ¼  m k1;j þ l1;j þ ri m1;j ; j ¼ 1; 2; 3 Dr i wm 4amþ1 2l 2 a2;j ¼  m k2;j þ l2;j þ ri m2;j ; j ¼ 1; 2; 3 Dr i wm  4amþ1 2l m  2 2 r 1 d3j  k3;j þ l3;j þ ri m3;j ; j ¼ 1; 2; 3 a3;j ¼ Dri iþ2 wm 3 3 X 2lm 2 2am X r i1 wb dk1 þ lk;j wjm1 þ mk;j qjm ; k ¼ 1; 2; 3 fk ¼ 2 2 Dr i wm j¼1 j¼1

1

r2

¼ q2 qðr; l; tÞ

6 4 a2;3 a1;3

a3;2

a3;1

32

w3m

3

2

f3

76 7 6 7 a2;1 54 w2m 5 ¼ 4 f2 5

a1;2

a1;1

w1m

1

v

r 2 ðwjþ1  wj Þ   Z tjþ1  @ ð1  l2 Þw @ðr 2 wÞ þ l þr þ r 2 rðrÞwðr; l; tÞ dt @l @r tj Z tjþ1 r 2 qðr; l; tÞdt ¼

ð92Þ

ð77Þ ð78Þ

We use the implicit method for time integration. That is the integrands are evaluated at t = tj+1, thus:

ð79Þ

v

1

r 2 ðwjþ1  wj Þ þ

ð80Þ

l

! @ðr2 wjþ1 Þ @½ð1  l2 Þwjþ1  þr þ r 2 rðrÞwðr; l; t jþ1 Þ Dt j @r @l

ð81Þ

f1

ð93Þ

Collecting the angular flux terms at t = tj+1 on the left side of the equation and dividing by Dtj, we obtain:

l

@ðr2 wjþ1 Þ @½ð1  l2 Þwjþ1  ~ ðrÞwjþ1 ðr; lÞ þr þ r2 r @r @l ~jþ1 ðr; lÞ ¼ r2 q

4amþ1 2l m 2 k3;j þ l3;j þ ri m3;j ; j ¼ 1; 2; 3 ð82Þ Dr i wm 4amþ1 2l 2 l2;j þ ri m2;j ; j ¼ 1; 2; 3 ð83Þ a2;j ¼  m k2;j þ Dr i wm  4amþ1 2l  2 l1;j þ ri m1;j ; j ¼ 1; 2; 3 ð84Þ a1;j ¼  m r2i1 d1j þ k1;j þ 2 Dr i wm 3 3 X 2l 2am X lk;j wjm1 þ mk;j qjm ; k ¼ 1; 2; 3 ð85Þ fk ¼  m r 2iþ1 wb dk3 þ 2 Dr i wm j¼1 2 j¼1

a3;j ¼ 

ð91Þ

tj

3

a2;2

ð90Þ

where ‘‘v’’ is the average speed of group neutrons. We define time steps as tj+1 = tj + Dtj and integrate over the jth time step:

ð76Þ

For lm < 0 and k = 3 in Eq. (22), in the first term with the Kronecker delta, we take the angular flux w3m ¼ wb where wb represents the angular flux at the first node of the previous ((i + 1)st) element in the marching scheme. For lm < 0 and for k = 1 or 2, we leave Eq. (22) unchanged. The resulting linear system is:

a3;3

ð89Þ

@w @ðr 2 wÞ @½ð1  l2 Þw þl þr þ r2 rwðr; l; tÞ @t @r @l

¼ r 2 qðr; l; tjþ1 ÞDtj

2

ð88Þ

The spherical geometry within group prompt neutron transport equation for the time dependent case can be written as:

f3

w3m

ð87Þ

2.6. Time differencing

v

a1;1 6 4 a2;1

ð86Þ

ð74Þ

2.5. The quadratic discontinuous finite element (QDFEM) equations

2

2 3 3 3 2~ 3 ~3;1 6 w12 7 a f3 6 7 7 27 ~ ~2;1 56 w1 7 ¼ 6 a 4f2 5 4 25 ~1;3 ~f 1 a w11 2

For the special direction, the following equations are obtained:

~1;j ¼ a

249

ð94Þ

where

r~ ¼ r þ

1

ð95Þ

v Dt j

and

~jþ1 ¼ qjþ1 þ q

wj V Dt j

ð96Þ

Since Eq. (94) is formally identical to Eq. (1), all previous formulation regarding spatial discretization is valid. The accuracy of the implicit method can be shown to be of order Dt (Lewis and Miller, 1984).

250

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

2.7. Centre condition

Table 2 Cross-sections (cm1) for the U-D2O system.

Since the direction cosines and weights are ordered so that

1 < l1 <    < lM < 0 < lMþ1 <    < lM < 1 2

ð97Þ

2

m

rf

rs

rt

1.70

0.054628

0.464338

0.54628

and

ð98Þ

wm ¼ wMmþ1 ;

ð99Þ

for the Gauss–Legendre quadrature we use the following centre condition:

wm ðr ¼ 0Þ ¼ wMmþ1 ðr ¼ 0Þ;

DD

k-eff

M 2 M m ¼ 1; . . . ; 2

lm ¼ lMmþ1 ; m ¼ 1; . . . ;

LD2 QDFEM LD1

M m ¼ 1; . . . ; 2

ð100Þ

3. Numerical applications Number of points

Our FORTRAN 90 program SPDOT which has been developed in a personal computer platform is capable of solving both time-independent (fixed source and criticality) and time-dependent problems in spherical geometry multigroup neutron transport. The Gauss–Legendre quadrature up to M = 256 is allowed. The main formulations (LD2 and QDFEM) which have been developed in this work are implemented as spatial discretization options. The TIMEX type LDFEM, that is LD1, and DD are also implemented in another version of the code for comparison purposes. When N is the number of spherical shells into which the spherical system is divided, the number of points for which the scalar flux is evaluated is equal to N for DD discretization. On the other hand, the linear discontinuous finite element discretization requires scalar flux evaluation at N + 1 points. If the quadratic discontinuous finite element method is employed and each spherical shell is taken as a quadratic finite element, the number of points for which scalar fluxes are evaluated equals 2N + 1. Linear (LD1, LD2) and quadratic formulations require the solution of 2  2 and 3  3 linear systems, respectively. Due to the low dimensionality of these linear systems, Cramer’s rule is employed for the solution. In this study comparisons are made between different spatial discretization methods holding the number of points for which scalar flux evaluations are made possibly constant.

Fig. 2. The variation of k-eff with respect to the number of points for the one-group problem with isotropic scattering.

3.1.1. One-group model of a critical bare homogeneous U-D2O system with isotropic scattering The one-group macroscopic cross sections for the uranium-D2O system are given in Table 2. The critical dimension (Rc) of the system is calculated analytically to be 22.017156 cm. In Fig. 2, the variation of the calculated effective multiplication factors (k-eff) is given as a function of the number of points for which scalar fluxes are evaluated in various spatial discretization methods. The convergence criterion for the inner and outer iterations is chosen to be 106. For outer iterations it is on the absolute fractional difference between k-eff estimates in consecutive iterations. In the inner iterations it is on the maximum absolute fractional difference in pointwise flux estimates in consecutive inner iterations. The quadratic discontinuous finite element method appears to yield more accurate k-eff values than all of its alternatives when the number of points is held constant. The newly developed linear discontinuous finite element method (LD2) seems to be more accurate than its linear alternative, LD1. The performance of the DD discretization is in between the two LDFEM discretizations. Next, we compare the spatial discretization methods with respect to the number of fission source (or outer) iterations to convergence and the computation times to assess their convergence behavior and computational efficiencies. Table 3 shows the number of outer iterations (NIT), the number of elements (NEM), the number of points (for which scalar fluxes are evaluated) (NP) and relative CPU times (t). t is defined as the ratio of the CPU time to the shortest CPU time observed in Table 3 (DD with NEM = 2).

3.1. Time independent (criticality) problems In order to test the developed methods in the time-independent case, we use the benchmark problems defined in the literature (Sood et al., 2003). Some of the benchmark problems are defined for spherical geometry in this set. In all cases, vacuum boundary condition is applied at the boundary, r = R. S32 calculation is performed using the Gauss–Legendre quadrature.

Table 3 Number of outer iterations and relative CPU times for different spatial discretization methods. NEM

2 4 5 7 10 15 20 30 50

DD

LD1

LD2

QDFEM

NIT

NP

t

NIT

NP

t

NIT

NP

t

NIT

NP

t

17 30 33 34 35 35 35 36 36

2 4 5 7 10 15 20 30 50

1.0 2.0 2.7 3.7 5.1 7.6 10.9 15.1 25.4

22 33 34 35 36 36 36 36 36

3 5 6 8 11 16 21 31 51

1.9 3.7 5.0 6.6 9.3 13.7 18.3 28.1 45.7

28 34 35 35 36 36 36 36 36

3 5 6 8 11 16 21 31 51

1.7 3.4 4.4 6.9 8.6 13.7 18.0 26.6 43.6

36 36 36 36 36 36 36 36 36

5 9 11 15 21 31 41 61 101

2.7 6.6 7.6 10.1 14.0 21.4 28.4 41.7 69.3

251

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255 Table 4 Number of elements, corresponding number of points, number of outer iterations and CPU times for different spatial discretization methods. Method

k-eff  0.999932

DD LD1 LD2 QDFEM

k-eff  0.999990

k-eff  0.999999

NEM

NP

NIT

t

NEM

NP

NIT

t

NEM

NP

NIT

t

10 21 5 2

10 22 6 5

35 36 35 36

5.1 19.9 4.4 2.7

28 72 15 5

28 73 16 11

36 36 36 36

13.9 65.9 13.7 7.6

59 176 55 21

59 177 56 43

36 36 36 36

29.7 162.3 45.1 28.9

k-eff

Table 7 Cross-sections (cm1) for U-D2O Reactor.

QDFEM LD2 DD LD1

Fig. 3. The variation of k-eff with respect to the number of points for the two-group problem.

m

rf

rs,gg

rs,g+1,g

rt

v

2.50 2.50

0.0010484 0.050632

0.62568 2.44383

0.029227 –

0.65696 2.52025

1.0 0.0

rs0

rs1

rt

0.054628

0.464338

0.056312624

0.54628

3.1.2. Two-group model of a critical bare, homogeneous system made of 93% enriched uranium The scattering is also assumed to be isotropic in this problem. The two-group cross sections for the system are presented in Table 5. The critical radius of the system is calculated analytically to be 16.049836 cm. Fig. 3 shows again the variation of the calculated k-eff values with respect to number of points for various spatial discretization methods. For a fixed number of points QDFEM yields more accurate k-eff values than the other methods. Our observations regarding the performance of various methods in k-eff calculation are qualitatively identical to those in the previous problem. QDFEM shows the best performance with LD2 following it closely as the second best. Table 6 shows the relative CPU times (t) for each method to obtain approximately given k-eff values. t is here defined as the ratio of the CPU time to the shortest CPU time observed (QDFEM with NEM = 3). As before, QDFEM results in shortest t compared to the other methods.

Table 5 Two-group cross-sections (cm1) for the 93% enriched uranium system.

1 2

rf

1.808381

codes are run to get approximately the given k-eff values. t is again defined as the ratio of the CPU time to the shortest CPU time observed (DD with NEM = 2). It is obvious that a coarser mesh is usually sufficient for QDFEM to get an equivalent accuracy to LDFEM and DD. If QDFEM is used, t is usually significantly reduced. Thus it appears to be the computationally most efficient spatial discretization method among the four.

Number of points

Group, g

m

Although DD requires the least number of outer iterations for convergence in coarse meshes the number of iterations becomes identical for all methods when spatial mesh is refined. When the number of elements is held constant, t is the longest for QDFEM since the number of points is approximately twice that of LD and DD methods. DD gives the shortest t when the number of elements is held constant as expected. While it is necessary to calculate the cell-edge fluxes and the cell-centre flux by solving a 3  3 linear system in QDFEM and only cell-edge fluxes by solving a 2  2 linear system in LDFEM, merely the cell-centre flux needs to be calculated in the DD method. When the number of points is held constant, QDFEM gives the shortest execution time among all methods. Indeed QDFEM execution times are at least 20% shorter than those of the LD methods for fixed number of points. The number of linear systems to be solved is reduced by a factor of two in the quadratic formulation relative to the linear one. Obviously this effect overrides the computational load incurred by the necessity of solving 3  3 systems instead of 2  2 ones. Table 4 shows the relative CPU times, the number of outer iterations to convergence and the number of points required when the

3.1.3. One-group model of a critical bare homogeneous U-D2O system with anisotropic scattering The cross sections which are given in Table 7 are identical to the first problem we have considered with the exceptions of scattering anisotropy introduced by the nonzero first Legendre component of the scattering cross section and the increase in the average number of neutrons emitted per fission. The modifications in the cross section data result in a decrease in the analytically calculated critical radius Rc, which is 18.30563081 cm in this case.

Table 6 Number of points, outer iterations and relative CPU times for different spatial discretization methods. Method

DD LD1 LD2 QDFEM

k-eff  0.999918

k-eff  0.999990

k-eff  0.999999

NEM

NP

NIT

t

NEM

NP

NIT

t

NEM

NP

NIT

t

19 44 10 3

19 45 11 7

16 16 16 16

2.3 9.4 2.0 1.0

44 119 33 12

44 120 34 25

16 16 16 16

5.2 25.4 6.6 3.8

68 195 63 23

68 196 64 47

16 16 16 16

8.2 42.2 12.6 7.4

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

k-eff

252

U-235 core

QDFEM LD2 DD LD1

H2 O reflector Fig. 5. Water reflected U-235 fuelled critical system.

Number of points Fig. 4. The variation of k-eff with respect to the number of points for the one group problem with anisotropic scattering.

DD QDFEM LD2 LD1

k-eff

The variation of k-eff with respect to the number of points for various methods is given in Fig. 4. A study of the figure reveals the fact that QDFEM yields the most accurate results with LD2 following it closely. Table 8 gives the relative CPU times (t) for comparison of the computational performance of various methods. t is here defined as the ratio of the CPU time to the shortest CPU time observed (QDFEM with NEM = 2). As before, QDFEM has the smallest relative CPU times to obtain the given k-eff values.

Number of points Fig. 6. The variation of k-eff with respect to the number of points in the core region for the reflected reactor problem.

3.1.4. One-group model of a critical, water reflected U-235 fuelled system Fig. 5 shows a two-region reactor which consists of the core and the reflector regions. Scattering is again assumed to be isotropic. The cross sections of the core and the reflector are presented in Table 9. The outer radii of the core and the reflector regions are 6.12745 cm and 15.318626 cm, respectively. With these dimensions the core and the reflector are 2 and 3 mean free paths thick, respectively. In the numerical solutions, the number of spherical shells in the reflector region is taken to be one and a half times that of the core region. The variation of k-eff with respect to the number of points in the core region for various methods is given in Fig. 6. The performance of various discretization methods follows the same order as in the previous problems. The only notable difference is the convergence of the DD method from above while the other techniques continue to converge from below to the analytical k-eff. Fig. 7 shows the percent error in k-eff with respect to CPU times to achieve this error values for different methods. Here, S256 quadrature is used. Also convergence criterion in k-eff is chosen as 1010. In this two-region problem, QDFEM appears to be by far the most computationally efficient method compared to the DD and LD methods.

Table 9 Cross sections (cm1) for the water reflected U-235 fuelled system.

U-235 H2O

m

rf

rs

rt

2.679198 0.0

0.065280 0.0

0.248064 0.293760

0.32640 0.32640

integral transport equation analytically (Olson and Henderson, 2004). One of these is the spherical shell source problem which is shown in Fig. 8. There is a bare spherical system of radius 10 cm. The outermost spherical shell part of the system of thickness 2 cm contains an isotropic neutron source of strength, 1 neutron/(cm3 s) which is activated at time t = 0. Initially there are no neutrons in the system. The total and scattering cross sections are 1 and 0.9 cm1 for the whole sphere. The neutrons are assumed to be monoenergetic with a speed of v = 1 cm/s. The exact scalar flux values are given at discrete points which are Dr = 0.5 cm apart for t = 1, 2.5, 5, 7.5, 10 s (early times) and 20, 50, 70, 100 s (late times) (Olson, 1999). Mesh cell sizes (Dr) are chosen identical in both source and non-source regions. Our calculations are conducted using the S8 approximation with the Gauss Legendre quadrature. The number of spherical shells used is 20 for both the DD and LDFEM (LD1 and LD2) methods. For

3.2. Time dependent problems Some benchmark solutions for time dependent neutron transport in spherical geometry have been generated by solving the

Table 8 Number of points, outer iterations and relative CPU times for different spatial discretization methods. Method

DD LD1 LD2 QDFEM

k-eff  0.999925

k-eff  0.999990

k-eff  0.999999

NEM

NP

NIT

t

NEM

NP

NIT

t

NEM

NP

NIT

t

12 29 5 2

12 30 6 5

28 28 28 28

1.9 8.0 1.5 1.0

24 69 14 5

24 70 15 11

28 28 28 28

3.6 19.2 3.9 2.3

33 104 27 10

33 105 28 21

28 28 28 28

4.8 29.5 7.0 4.4

253

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

DD

Error in k-eff (%)

LD2

L2 error norm

QDFEM LD1

S4 S8 S16 S32

Time (s)

CPU time (s) Fig. 7. The variation of percent error in k-eff with respect to the CPU time necessary to achieve this accuracy.

Fig. 10. The variation of the L2 error norms with time for various angular approximations in the QDFEM.

Δt=10 s Δt=1 s

L2 error norm

8 cm

2 cm

Δt=0.1 s Δt=0.01 s Δt=0.001 s

Fig. 8. Spherical shell source problem.

Time (s) DD

Fig. 11. Variation of the L2 error norms with time for different time step sizes, Dt.

LD1

L2 error norm

LD2 QDFEM

L2 error norm

10 el.

Time (s)

20 el. 40 el. 80 el.

Fig. 9. The variation of L2 error norms with time for various discretization methods in the S8 approximation.

the quadratic case (QDFEM) the number of spherical shells (or quadratic elements) is 10. This means 20 points for DD and 21 points for LDFEM and QDFEM. For time integration, we have chosen Dt = 0.1 s. To assess the relative merits of various methods we define the L2 norm as:

k/exact

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XN  2  /k ¼ £exact  £n r 2n Dr n n n¼1

ð101Þ

Fig. 9 shows the evolution of the L2-error norm with respect to time for various spatial discretization methods. A study of the results reveals that the best method is QDFEM with LD2 following it as the second best which is consistent with our observations in time independent problems. The major difference between time dependent and time independent cases is the superiority of LD1 over the DD in time dependent problems in contrast to the situation met in time independent problems. Since QDFEM appears to be the most accurate approximation among those considered, only its convergence properties with respect to SN order, time step size and number of elements will be investigated in the rest of this subsection.

Time (s) Fig. 12. L2 error norms of QDFEM for different number of elements.

First of all, the convergence behavior of QDFEM with respect to SN order will be studied. A time step size of Dt = 0.001 s and 500 quadratic elements are used. The variation of the L2 error norm with time is presented graphically in Fig. 10 for various SN orders. From the figure, the reduction in the L2 error norm with respect to increasing SN order is clear, indicating convergence with respect to refinement in angular approximation. Secondly, the variation of L2 error norms with respect to time for different time step sizes, Dt, will be considered. In order to show the convergence of QDFEM with smaller time steps, we run SPDOT with S8 quadrature order and 20 elements. The relevant results are given in Fig. 11. L2 error norms for Dt = 0.01 s and Dt = 0.001 s are almost indistinguishable at all times. This shows the convergence of the method as time step size is decreased. Also we note that the error norms for all time step cases become nearly identical after t = 100 s since convergence to the asymptotic, timeindependent solution is reached.

254

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

time ind.

Scalar flux (n 0/cm2s)

t=100 s t=70 s t=50 s t=20 s t=10 s t=5 s t=1 s

r (cm) Fig. 13. Comparison of time-dependent and time-independent flux profiles.

Table 10 Temporal–spatial–angular mesh cases considered. Case

Dt (s)

Quad. order

Number of elements

1 2 3 4

0.5 0.1 0.01 0.001

S8 S16 S32 S64

10 20 50 100

travel from the source region into the inner regions. More neutrons will be able to arrive at the inner regions of the sphere with increasing time. Finally, the flux profile will reach a time-independent distribution. The spatial variation of the flux at various elapsed times, t is given in Fig. 13. These results are obtained by running code with Dt = 0.5 s, S8 quadrature order and 10 elements. At t = 50, 70 and 100 s, flux profiles almost overlap with time-independent scalar flux profile. Thus the convergence of time-dependent solutions to the time-independent distribution is obvious. Finally, we will consider the radial variation of the difference between the time dependent flux at a very late time (t = 100 s) and the time-independent flux for various temporal–spatial–angular meshes which are defined in Table 10. Mesh cases in Table 10 vary from the coarsest case (Case 1) to the finest case (Case 4). The results are presented graphically in Figs. 14 and 15 which show absolute and relative error variation with respect to radial distance respectively. The maximum difference between timedependent and time-independent scalar fluxes occurs at the centre of the sphere in all cases. The refinement in the temporal–spatial– angular mesh results in convergence as expected. Time dependent fluxes are closer to the time-independent fluxes in the source region than the central non-source region.

ind-

dep

(n0/cm2s)

4. Conclusion

case 1 case 2 case 3 case 4

r (cm)

100x(Φ Φ ind-Φ dep)/Φ ind

Fig. 14. Difference between time-independent and time-dependent scalar fluxes as a function of radial distance (r) for different cases.

r (cm) Fig. 15. Variation of relative errors with respect to radial distance for different cases.

Thirdly, the variation of L2 error norms with respect to time for different number of elements will be investigated. The relevant results are given in Fig. 12. S8 and 0.1 s are chosen as the quadrature order and time step size respectively. The convergence is clear with increasing number of elements. Physically, there can be no neutrons in the interior regions of the sphere at early times since it takes some time for neutrons to

In this study, we have developed linear and quadratic discontinuous finite element methods for the solution of the spherical geometry neutron transport equation. Both time-independent and time-dependent cases are investigated. The developed formulations (LD2 and QDFEM) are compared with the results of formulations previously developed (LD1 and DD). QDFEM turns out to be the spatial discretization method which can generate the most accurate k-eff values and proves itself to be the fastest method for time-independent problems. LD2 gives the second best k-eff results for a given number of points in time-independent problems. When we try to assess computational efficiency in terms of CPU times required for a fixed accuracy in k-eff, LD2 turns out to be the second best when the required accuracy in k-eff is not too great. But DD becomes the second best if the required accuracy in k-eff is enhanced when the k-eff convergence criterion is 106. But in water reflected U-235 problem LD2 is the second best with 1010 convergence criterion and S256 quadrature. LD1 becomes the worst method in all time-independent problems. For time-dependent problems, QDFEM and LD2 come out as the best and the second best methods respectively when the accuracy is measured in terms of the L2 error norm in the scalar flux. In contrast to time-independent problems, LD1 outperforms DD when accuracy is measured in terms of the L2 error norm in scalar flux. In both time-dependent and time-independent problems LD2 formulation outperforms LD1 in all cases studied. This stems from the fact that the neutron spatial dependence is largely ignored in the angular redistribution term in the LD1 method while the finite element expansion is utilized in the LD2 method for this term. By numerical experimentation we have found out that if LD1 is modified by using the finite element expansion for the angular fluxes in the angular redistribution term, it gives identical results to LD2. Our study of the convergence properties of QDFEM reveals that this method possesses the desired convergence properties under the refinement of temporal, spatial and angular meshes. In short, numerical experiments indicate that QDFEM is a well suited method for spherical geometry neutron transport problems when high accuracy and computational efficiency is required.

M. Mercimek, H. Atilla Özgener / Annals of Nuclear Energy 64 (2014) 244–255

References Bailey, T.S., Adams, M.L., Chang, J.H., 2009. A piecewise linear discontinuous finite element spatial discretization of the transport equation in 2D cylindrical geometry. In: International Conference on Advances in Mathematics, Computational Methods, and Reactor Physics, 3–7 May, 2009. Saratoga springs, NY, United States. Chang, J.H., Warsa, J.S., 2007. Comparison of linear and quadratic discontinuous spatial finite element methods for parallel SN transport on triangles. In: Joint International Topical Meeting on Mathematics and Computation and Supercomputing in Nuclear Applications, April 15–19, Monterey, California. Cockburn, B., 2001. Devising discontinuous Galerkin methods for non-linear hyperbolic conservation laws. J. Comput. Appl. Math. 128, 187–204. Cockburn, B., Karniadakis, G., Shu, CW., 2000. The development of discontinuous Galerkin methods. In: Cockburn, B., Karniadakis, G., Shu, CW. (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Lecture Notes in Computational Science and Engineering, vol. 11. SpringerVerlag, New York, pp. 3–50. Hesthaven, J.S., Warburton, T., 2008. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Springer Texts in Applied Mathematics, vol. 54. Springer-Verlag, New York. Hill, T.R., 1975. ONETRAN, A Discrete Ordinates Finite Element Code for the Solution of the One-dimensional Multigroup Transport Equation. Los Alamos Scientific Laboratory Report LA-5990-MS (June). Hill, T.R., Reed, W.H., 1976. TIMEX: A Time-dependent Explicit Discrete Ordinates Program for the Solution of Multigroup Transport Equations with Delayed Neutrons. Technical Report LA-6201-MS. Hill, T.R., 1977. Application of the ONETRAN and TIMEX codes to shielding problems. In: Fifth International Conference on Reactor Shielding, April 18– 22, Knoxville, Tennesse. Hong, Z., Yuan, G., Fu, X., 2010. Oscillation of numerical solution for time-dependent particle transport equations. Prog. Nucl. Energy 52, 315–320. Lewis, E.E., Miller, W.F., 1984. Computational Methods of Neutron Transport. John Wiley, New York. Li, B.Q., 2006. Discontinuous finite elements in fluid dynamics and heat transfer. Computational Fluid and Solid Mechanics Series, ISSN 1860-482X, Germany.

255

Machorro, E., 2007. Discontinuous Galerkin finite element method applied to the 1D spherical neutron transport equation. J. Comput. Phys. 223 (1), 67–81. O’Dell, R.D., Alcouffe, R.E., 1987. Transport Calculations for Nuclear Analysis: Theory and Guidelines for Effective Use of Transport Codes. Report LA-10983-MS, Los Alamos National Laboratory. Olson, K.R., Henderson, D.L., 2004. benchmark solutions for time-dependent neutral particle transport in one-dimensional homogeneous media using integral transport. Ann. Nucl. Energy 31, 1495–1537. Olson K.R., 1999. Neutral Particle Integral Transport in Inertial Confinement of Fusion Systems using Time Dependent Integral Transport Methods. Doctoral Dissertation, University of Wisconsin-Madison. Reed, W., Hill, T., 1973. Triangular Mesh Methods for the Neutron Transport Equation. Technical Report la-ur-73-479, Los Alamos Scientific Laboratory. Reed, W.H., Hill, T.R., Brinkley, F.W., Lathrop, K.D., 1973. TRIPLET: A Two Dimensional, Multigroup, Triangular Mesh, Planar Geometry, Explicit Transport Code. Technical Report LA-5428-MS, Los Alamos National Laboratory. Reed, W.H., Hill, T.R., Brinkley, F.W., Lathrop K.D., 1977. TRIDENT: A Two Dimensional Multigroup, Triangular Mesh, Explicit Neutron Transport Code. Technical Report LA-6735-MS, Los Alamos National Laboratory. Seed, T.J., Miller, W.F., Brinkley, F.W., 1977. TRIDENT: A Two-dimensional, Multigroup Triangular Discrete Ordinates, Explicit Neutron Transport Code. Technical Report LA-5428-M, Los Alamos Laboratory. Seed, J., Miller, W.F., Bosler, G.E., 1978. TRIDENT: a new triangular mesh discrete ordinates code. In: Topical Meeting on Advances in Reactor Physics, Gatlinburg, TN (April 9–12). Sood, A., Forster, R.A., Parsons, D.K., 2003. Analytical benchmark test set for criticality code verification. Prog. Nucl. Energy 42 (1), 55–106. Wang, Y., 2009. Adaptive Mesh Refinement Solution Techniques for the Multigroup SN Transport Equation using a Higher-order Discontinuous Finite Element Method. Dissertation, Texas A&M University (May). Wareing, T.A., McGhee, J.J., Morel, J.E., 1996. ATTILA: A three-dimensional unstructured tetrahedral mesh discrete ordinates transport code. In: American Nuclear Society Annual Winter Meeting, vol. 75, ANS, Washington, DC, p. 146.