deterministic coarse mesh neutron transport method

Annals of Nuclear Energy 41 (2012) 1–11

Contents lists available at SciVerse ScienceDirect

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

An efficient hybrid stochastic/deterministic coarse mesh neutron transport method Dingkang Zhang, Farzad Rahnema ⇑ Nuclear & Radiological Engineering and Medical Physics Programs, Georgia Institute of Technology, 770 State Street NW, Atlanta, GA 30332-745, USA

a r t i c l e

i n f o

Article history: Received 14 June 2011 Accepted 27 September 2011 Available online 16 December 2011 Keywords: Coarse mesh neutron transport Incident flux response expansion Hybrid stochastic and deterministic Response function Whole core transport Whole core benchmark problem

a b s t r a c t A new incident flux response expansion method has been developed to significantly improve the accuracy of the hybrid stochastic/deterministic coarse mesh transport (COMET) method. Additionally, two acceleration techniques are introduced that significantly increase the computational efficiency of the method by several folds. The new expansion method removes singularities associated with the current method that degrade its accuracy and efficiency and ability to solve realistic problems with complexity and size that are inherent in operating commercial reactors. It also enables (paves the way for) the response method to be imbedded in low order transport methods (e.g., diffusion theory) for improving accuracy without degradation in efficiency. In general, the new expansion method also enables efficient and accurate coupling of different deterministic methods (e.g., characteristic to discrete ordinates and in general high order transport to high or low order transport). The new method improvements enable COMET to perform whole-core neutronics analysis in all light and heavy water operating reactors with Monte Carlo fidelity and efficiency that is several orders of magnitude faster than both direct Monte Carlo and fine mesh transport methods. A stylized CANDU-6 core benchmark problem with and without adjuster rods was used to test the accuracy and efficiency of the COMET method in whole (full) core configurations at two coolant states. The benchmark problem consisted of 4560 fuel bundles containing a total of 168,720 fuel pins and 21 adjuster rods. The COMET solutions were compared to direct Monte Carlo (MCNP) reference solutions. It was found that the core eigenvalue, bundle averaged and fuel pin power distributions predicated by COMET agree very well with the MCNP reference solution in all cases when the coarse mesh incident angular flux expansion in the two spatial and two angular (azimuthal and polar) variables is truncated at 4, 4, 2 and 2, respectively. These comparisons indicate that COMET can achieve accuracy comparable to that of the Monte Carlo method with a computational efficiency that is several orders of magnitude better. Ó 2011 Published by Elsevier Ltd.

1. Introduction The traditional steady-state neutronics methods for whole core criticality analysis consist of three steps. Initially, a series of detailed transport calculations are performed for each unique lattice cell (containing a single assembly/bundle) in the core with infinitemedium (specular reflective) boundary conditions. These calculations produce, in the pre-computation phase, a library of spatially homogenized and energy collapsed lattice cross sections by flux weighting the fine group or continuous energy cross sections in the lattice. This library which includes flux discontinuity factors (Smith, 1986) is finally used in a few-group coarse mesh diffusion theory based core calculations. The accuracy of this approach deteriorates as fuel assembly and core heterogeneity increases. For example, large flux gradients in the core as a result of strong heterogeneity break down the diffusion approximation and the ⇑ Corresponding author. E-mail address: [email protected] (F. Rahnema). 0306-4549/$ - see front matter Ó 2011 Published by Elsevier Ltd. doi:10.1016/j.anucene.2011.09.024

infinite-medium boundary condition approximation used in generating the homogenized cross sections and discontinuity factors. Also, because of the cross section homogenization, pin power reconstruction techniques are ad hoc and suffer from the same problems related to core environmental effects (e.g., assembly and core heterogeneity, infinite-medium boundary condition). To address these issues, an incident flux response expansion method for heterogeneous coarse mesh transport problems was developed by Mosher and Rahnema (2006). Similar to the traditional methods, they divide the core into a set of contiguous coarse meshes of fuel assembly size and then expand the angular flux at the mesh boundaries in terms of discrete Legendre polynomials. The discrete ordinates method is then used to generate a response function library for each unique coarse mesh. A two-level (inner and outer) iteration method is finally used to perform whole-core calculations. The benchmark calculations against the whole-core discrete ordinates method in 1-D slab geometry demonstrated that the method is highly accurate and efficient. To take advantage of the geometric flexibility of Monte Carlo methods, more recently

2

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

Forget and Rahnema (2006a, 2006b) expended the angular current in terms of tensor-products of (continuous) Legendre polynomials and used a stochastic method to generate the response function library. Although, the results (Forget and Rahnema, 2006a) for the small C5G7 MOX benchmark problem with control rods (Lewis et al., 2005) was encouraging, the maximum errors in fuel pin powers were very large (10%). Also, a highly simplified quarter core CANDU-6 benchmark problem with no adjuster rods (i.e., control material) and only four burnup points was considered by Forget and Rahnema (2006b) to evaluate the method. It was found that the results agreed well with direct Monte Carlo calculations. For example, the flux weighted average error was 0.7% with a maximum error of 3.5% for a 2nd order expansion. However, no appreciable improvement in the results was obtained when the expansion order was increased to 4; the computational speed, although appreciably better than MCNP (Briesmeister, 1997), was significantly increased from 4 h to about 22 h. Also, they did not consider any problems that are typical of a realistic operating reactor core, e.g., a CANDU-6 core with adjuster (absorber) rods. In this paper, the COMET method (Forget and Rahnema, 2006a,b,c) has been extended to resolves issues that affect its accuracy and computational efficiency. It also makes the method robust for imbedding into low order transport methods or for coupling different transport solution methods (e.g., characteristic to discrete ordinates). These extensions include an angular flux (instead of the usual current) expansion technique and deterministic iterative acceleration methods. For example, the angular flux expansion removes singularities issues that arise when the current-based method is coupled to or imbedded in low order methods such as those based on diffusion theory. More serious consequences of the current-based method are twofold. (1) In general, given the same truncation order a better accuracy is achieved with the flux expansion as opposed to the expansion of the incident current for any problem. This is because for grazing angles (i.e., cosine of the polar angle near zero), the current expansion method will require a much larger order to properly represent the angular distribution of the incident neutrons. (2) It is not possible for the current-based method to calculate the mesh interface scalar flux if needed for any purpose in the problem (e.g., into diffusion methods). The work in this paper results in a significant gain in computational efficiency and provides consistently high accuracy regardless of the core configuration (geometry and material distribution) and heterogeneity. The remainder of the paper is organized as follows. The hybrid stochastic/deterministic COMET method is presented in Section 2. Section 3 describes a stylized CANDU-6 core benchmark problem. The COMET benchmark results are presented in Section 4. Conclusions and future work are discussed in Section 5.

^  r þ rð~ H¼X r; EÞ 

The neutron angular flux distribution in a large heterogeneous system V with boundary o V is governed by the Boltzmann transport equation:

^ ; EÞ ¼ 1 Fuð~ ^ ; EÞ ~ Huð~ r; X r; X r2V k

ð1aÞ

with the boundary condition

^ ; EÞ ¼ Buð~ ^ 0 ; E0 Þ ~ ^ < 0 and n ^0 > 0 ^X ^X uð~r; X r0 ; X r;~ r 0 2 @V; n

ð1bÞ

^ ; EÞ represents the angular flux at space ~ ^ where uð~ r; X r, direction X and energy E, k is the eigenvalue (effective multiplication factor) ^ is the outward normal on the external boundof the large system, n ary. The symbol B represents the general boundary condition operator. The operators H and F are defined below.

dE

0

Z

^ 0 rs ð~ ^ 0 ; E0 ! X ^ ; EÞ dX r; X

ð2Þ

4p

F¼

vð~r; EÞ 4p

Z

dE

0

Z

^ 0 mrf ð~ dX r; E0 Þ

ð3Þ

4p

The spatial domain V is first divided into a number of nonoverlapping coarse meshes {Vi|i = 1, . . ., N}. Consistent with the incident flux response expansion method [Mosher and Rahnema, 2006], solving the whole-core transport problem (Eq. (1)) is then equivalent to solving a set of local fixed source problems given by

^ ; EÞ ¼ Hwi ð~ r; X

1 ^ ; EÞ ~ Fw ð~ r; X r 2 Vi k i

ð4aÞ

with the boundary condition

^ ; EÞ ¼ u ð~ ^ ^ <0 ~ ^þis  X wi ð~ r; X r 2 @V is and n is r; X; EÞ

ð4bÞ ^þ Vi, n is

where oVis represents surface s of coarse meshes denotes the ^ ; EÞ is the angular flux from the outward normal on oVis, and wis ð~ r; X adjacent coarse mesh which is next to surface s. 2.1. Flux expansion and orthogonal expansion function ^ ; EÞ; is not a Since the flux from the adjacent coarse mesh, wis ð~ r; X known priori, it is assumed that the angular flux at each bounding surface can be expanded in terms of orthogonal expansion functions in both outward and inward hemispheres.

^ ; EÞ ¼ wis ð~ r; X

X

~ ^ ~ J ;m r 2 @V is is Cm ðr; X; EÞ

^ >0 ^ is  X and n

ð5Þ

m

^ ; EÞ are orthogonal expansion coefficients which will where Cm ð~ r; X be described later, and J ;m represent the outgoing/incoming expanis sion coefficients. The choice of the symbol is motivated by the fact that 0th moment expansion coefficient is the partial current. Due to the linearity of the Boltzmann transport equation, the solution to the local problems can be constructed as:

^ ; EÞ ¼ wi ð~ r; X

X

m ~ ^ ~ J ;m r 2 Vi is Ris ðr; X; EÞ

ð6Þ

s;m

~ ^ where Rm is ðr; X; EÞ is the solution to the following local fixed-source problems:

^ ; E; kÞ ¼ HRis ð~ r; X

1 ^ ; E; kÞ ~ r; X r 2 Vi FRis ð~ k

ð7aÞ

with the boundary condition

^ ; E; kÞ ¼ Cm ð~ ^ ; EÞ ~ ^ <0 ^ þis  X Ris ð~ r; X r; X r 2 @V is and n

ð7bÞ

^ ; E; kÞ are the volumetric angular flux distribuClearly, X tions within coarse mesh i responding to a unit incoming flux with ^ ; EÞ imposed on surface s. As a rea phase space distribution Cm ð~ r; X m ^ ~ sult, Ris ðr; X; E; kÞ are called (surface-to-volume) response functions in this paper. It should be pointed out that Eq. (7) forms a local fixed-source problem with the fission source scaled by the global eigenvalue which is not known a priori. Theoretically, the choice of expansion functions is arbitrary. If the expansion function set is complete, one could use the incident flux response expansion method to solve the global transport problem (1) without approximation. However, for practical reasons, the expansion in Eq. (5) must be truncated at a finite expansion order. Ideally, both high accuracy and high efficiency can be achieved if one can choose an appropriate set of functions by which the actual flux distribution can be sufficiently represented with a low order expansion. Since the neutron distribution is isotropic on the average in the core, it is natural to require the 0th order expansion function to represent the isotropic flux. It is also very import to require that the magnitude of the partial current is always conserved no matter what ~ Rm is ðr;

2. Hybrid stochastic/deterministic COMET method

Z

3

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

expansion order is used in Eq. (5). As a result, the following orthogonality condition is required.

Z

dE

Z

d~ r

Z ^ >0 ^  X n is

@V is

^ ðn ^ ÞCm ð~ ^ ; EÞCm0 ð~ ^ ; EÞ ¼ Am dmm0 ^ is  X dX r; X r; X

ð8Þ

where Am is a constant and dmm0 is the Kronecker delta. It is found that the tensor-products of Legendre polynomial used in Mosher and Rahnema (2006) and Forget and Rahnema (2006a) do not satisfy these conditions. Specifically, in Mosher and Rahnema’s work the tensor-products of Legendre polynomial were used to expand the angular flux. However, this conserves only the scalar flux and not necessarily the particle balance (e.g., partial currents are not conserved). In Forget and Rahnema’s paper the tensor-products of Legendre polynomial were used to directly expand angular currents (rather than fluxes), this can guarantee that the partial currents are conserved, but its 0th expansion function does not represent an isotropic flux. This leads to a relatively high expansion order in the angular variables to achieve good accuracy degrading computational efficiency significantly. In this work, the expansion functions in Eq. (9) are used to expand angular fluxes on the mesh interface. The expansion functions are constructed as a product of Legendre polynomials Pn(x) and Chebyshev polynomials of the second kind Un(x) .

^ ; EÞ ¼ Pi ð~ Cijkg ð~ r; X rÞU j ðcos hÞPk ðcos uÞdðE  Eg Þ

ð9Þ

Here, polynomials Un(x) have the recurrence relation given below. The Dirac delta function in E is equivalent to the multigroup treatment of the energy variable.

It is obvious from Eqs. (7a), (7b), and (13) that the response functions depend on the coarse mesh geometry and material compositions (i.e., type) and consequently only the response functions for each unique coarse mesh type (t) need to be computed. The Monte Carlo method is chosen to solve the local problems given by Eqs. (7a) and (7b) because of its flexibility to handle arbitrarily complex geometry. In response function calculations, source particles are sampled based on the incoming flux phase space distribu^ ; EÞ. After the particles escape from a surface or enter tion Cm ð~ r; X into the region of interest, the expected value of response functions are tallied as the probability that initial neutrons and their progenies escape from the coarse mesh through a specific surface, scaled by a weight factor determined by the orders of the expansion and the location, energy and direction of the exiting neutrons. For example, the surface-to-surface response function for coarse mesh type t can be calculated as 0

Rmm tss0 ¼

1 X ^ n ; En Þwn Cm0 ð~ rn ; X N ~r 2@V

ð14Þ

ts0

n

^ n and wn are where N is the total number of source particles, ~ r n ; En ; X the position, energy, direction and weight of particle n when they escape through the coarse mesh surface. The response function uncertainties must be evaluated in addition to the mean value. Since the true sample mean is not known, an approximation of the true variance, also known as the sample variance, is evaluated according to the following equation: mm0

d2 ½Rtss0

2 32 X 1 X 1 ^ n ; En Þwn 2  4 ^ n ; En Þwn 5 ¼ ½Cm0 ð~ rn ; X Cm0 ð~ rn ; X N ~r 2@V N ~r 2@V n

ts0

n

ts0

U 0 ðlÞ ¼ 1

ð15Þ

U 1 ðlÞ ¼ 2l

ð10Þ

U nþ1 ðlÞ ¼ 2lU n ðlÞ  U n1 ðlÞ for n P 1 Clearly, the 0th order angular expansion is constant and consequently represents an isotropic flux in angle. From Eqs. (5) and (8), the flux expansion coefficients (moments) can be calculated as

J ;m ¼ is

Z

dE

Z

d~ r

Z ^ >0 ^  X n is

@V is

^ ðn ^ ÞCm ð~ ^ ; EÞw ð~ ^ ^ is  X dX r; X is r; X; EÞ:

ð11Þ

It is evident that the 0th expansion coefficient is identical to the total partial incoming/outgoing current crossing the mesh boundary when the new expansion functions are used. As a result, the total partial currents (or particles) based on these expansions are always conserved. 2.2. Response function generation The outgoing partial current and its high moments from surface s0 of coarse mesh i can be calculated by projecting Eq. (6) onto the expansion functions: 0

J þ;m ¼ is0

X

0

0

;m Rmm iss0 ðkÞJ is

ð12Þ

s;m 0

where Rmm iss0 is the surface-to-surface response function having the following relation with the surface-to-volume response function ~ ^ Rm is ðr; X; E; kÞ . 0

Rmm iss0 ðkÞ ¼

Z

dE

Z

@V is0

As evidenced from Eq. (7a) and (7b), response functions are only dependent on (1) the coarse mesh geometry and composition, (2) the phase space distribution of the incident flux and (3) the global (reactor core) eigenvalue k. The phase space distribution is known once the expansion functions Cm are chosen. The coarse mesh geometry and composition are also given (known) by the specification of the system volume (reactor core) and the choice of the coarse mesh grid. In theory, coarse meshes can be arbitrary in size. However, in order to take advantage of the repeated geometry of reactor cores, these meshes are chosen to be of the size of a fuel assembly/lattice as will be seen in the numerical results section. As it is well known, operating reactors use a limited number of fuel assemble/lattice types in the reactor core designs and therefore the response functions (RF) can be precomputed to form the library for solving the transport in equation in the core. Since k is not known a priori, the library must include response functions at a set of predefined values of k. The linear interpolation below is used to determine the RF for the calculated k in the outer iteration described in Section 2.3.

d~ r

Z

^ >0 ^ þ0 X n

^ ðn ^ ÞCm0 ð~ ^ ; EÞRm ð~ ^ ^ þis0  X dX r; X is r; X; E; kÞ

is

ð13Þ 0 Rmm iss0

It can be seen that represents the magnitude of the outgoing angular flux from surface s0 in moment m0 responding to a unit ^ ; EÞ impinging incoming flux with a phase space distribution Cm ð~ r; X on surface s.

Rmm tss0 ðkÞ ¼

1=k  1=k1 mm0 1=k2  1=k mm0 R 0 ðk1 Þ þ R 0 ðk2 Þ 1=k2  1=k1 tss 1=k2  1=k1 tss

ð16Þ

It is interesting to note that given the precomputed RF library, any combination of the unique coarse meshes making up the reactor core can be solved by the method very quickly. In this respect the COMET method is similar to the current industry methods which rely on a precomputed library of data (homogenized cross sections) that only depend on the lattice type (coarse mesh geometry and material composition). However, in contrast to industry methods in which an assumed boundary condition (e.g., specular reflection) is used to homogenize the cross sections, the COMET method retains the full heterogeneity of the lattice and the calculation of response functions do not require a pre-assumed boundary condition.

4

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

2.3. Deterministic iteration method Once the response function library is generated, an iterative deterministic method is used to perform the whole-core calculations. The algorithm consists of two level iterations: (1) outer iterations on the global eigenvalue, and (2) inner iterations on the incoming/outgoing partial currents crossing coarse meshes. We use indices (u, v) to denote the uth outer iteration and vth inner iteration. The numerical procedures without accelerations for solving whole core problems are described below. (1) Start an initial guess of the global eigenvalue k(1) and the expansion coefficients J ;m is ð1; 1Þ. A reasonable initial guess is k(1) = 1 since an operating reactor core is critical (k = 1). (2) Eq. (16) is used to update the response functions for the updated eigenvalue (3) Perform inner iterations to solve the eigenvalue of the numerical matrix given below.

X

0 ;m Rmm tðiÞss0 ðkÞJ is

¼

0 kJ þ;m is0

þ;m0

P

mm0

(a) Calculate the outgoing flux J is0 ðu; v þ 1Þ ¼ s;m RtðiÞss0 ðkðuÞÞJ ;m is ðu; v Þ (b) Use the interface or external boundary conditions to calculate the incoming flux J ;m is ðu; v þ 1Þ (c) Normalize the incoming flux J ;m is ðu; v þ 1Þ (d) Repeat steps (a)–(c) until the following convergence criterion is satisfied.

where eJ is the predefined convergence criterion on the scalar current. (4) Use the following neutron balance equation to update the global eigenvalue: kðu þ 1Þ R R R v Þ dE V i d~r 4p dX^ mrf RmtðiÞs ðr; X^ ; E;~kðuÞÞ R P m ;0 ^ v Þ dE V i d~r 4p dXra RtðiÞs ð~r; X^ ;E; kðuÞÞ þ @V is @V ðJþ;0 is ðu; v Þ  J is ðu; v ÞÞ P

;m i;s;m J is ðu;

R

;m i;s;m J is ðu;

R

ð19Þ (4) Repeat steps (2)–(4) until the following convergence criterion  is satisfied   kðu þ 1Þ  < ek   1   kðuÞ

L

ð21Þ (

ð1; 1Þ ¼ J ;m;H is

J ;m;L ðU; VÞ s ¼ 1; 2; . . . ; M L is 0

s ¼ ML þ 1; . . . ; M H

ð22Þ

where, the superscripts L and H denote a lower or high order expansion, ML and MH are the numbers of expansion functions, U and V represent the total number of outer and inner iterations in the low order calculations. Numerical experiments have shown that the improvement in the computational efficiency using this strategy is about a factor of 2–3.

ð18Þ

is

¼P

H

k ð1Þ ¼ k ðUÞ

ð17Þ

s;m

  þ;m0  J 0 ðu; v þ 1Þ   is  1 < eJ  þ;m0   J 0 ðu; v Þ

2.3.1. Low order acceleration (LOA) Clearly, if the initial guess of the global eigenvalue and partial currents are very close to the actual values, the iteration procedures converge very fast. Since the number of response function coefficients increases exponentially as the expansion order increases, a low order calculation is significantly faster than higher order calculations. For example, when the expansion order set {2, 2, 1, 1} is used, which denotes that the expansion orders in the two spatial and the two angular variables are 2, 2, 1 and 1, the whole-core computation time is about 2% of that when the expansion order set is {4, 4 , 2, 2}. Although the solution with a low order approximation may not be accurate, it would serve as a good estimate of the core effective multiplication factor and partial currents crossing coarse mesh interfaces. As a result, one can use a low order solution to accelerate the high order calculations. Since a low order expansion functions are a subset of high order expansion functions, the initial guess for the high order iterations can be easily written as:

ð20Þ

where ek is the predefined eigenvalue convergence criterion. Note in the above mentioned unaccelerated algorithm inner iterations (step 3) use the power method to converge on the largest eigenvalue k and its associated eigenvector (incoming/outgoing partial current moments). This largest eigenvalue is actually the discontinuity factor between the outgoing and incoming partial currents. When the unconverged effective multiplication factor k is less than the actual value, the discontinuity factor is greater than 1, i.e. the incoming partial current entering into a coarse mesh is less than the outgoing partial current exiting from its adjacent coarse mesh. When the exact effective multiplication factor is used, the discontinuity factor is equal to 1 and the incoming partial current matches the outgoing partial current from the adjacent region. It is well known that the convergence rate of the power iterations is very low when the dominant ratio is close to 1 for large reactors. This is especially true when the number of total coarse meshes and expansion order is very high. Two methods have been used to accelerate the iterative algorithm in this work.

2.3.2. Chebyshev acceleration method Since the convergence rate of the power iteration method might be unacceptably slow, the Chebyshev polynomial filtering (Saad, 1992) is often used to speed up the convergence of the standard method. The theoretical foundation is that a polynomial of matrix A, pq(A), has the same set of eigenvectors as the original matrix A, where pq is a polynomial of degree q. As a result, one can choose a good polynomial to converge on the desired eigenvectors faster. It has been shown that the Chebyshev polynomial filtering can efficiently converge on the eigenvector associated with the largest eigenvalue. To accelerate the inner iterations, we use the following equation to calculate the outgoing partial currents in the numerical step 3a): 0

J þ;m ðu; v þ 1Þ ¼ is0

X

0

;m pq ðRmm tðiÞss0 ðkðuÞÞÞJ is ðu; v Þ

ð23Þ

s;m

where pq and wi are defined as

pq ðRÞ ¼

q Y

ðR  wi IÞ;

ð24Þ

    2i  1 p þ 1  a; wi ¼ c cos 2q

ð25Þ

i¼1

c, q and a are the input parameters which can be determined empirically (Saad, 1992). Nevertheless, values of 0.985, 10 and 1.0 suggested in reference (Stamm’ler and Abbate, 1983) work very well in this study. 2.4. Uncertainty analysis Since the Monte Carlo method is used to generate the response function, their inherent uncertainty will propagate into core solutions when the response function library is used to perform

5

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

whole-core calculations. As a result, it is important to calculate the uncertainty of the whole core solution. In principle, the variance of the global eigenvalue k can be computed as:

X

@k

t;s;s0 ;m;m0

0 @Rmm tss0

2

d ½k ¼ þ

!2

X  @k 2 2 m d ½RAts  m t;s;m @RAts

0

d2 ½Rmm tss0  þ

X  @k 2 2 m d ½RF ts  @RF m ts t;s;m

ð26Þ

m where RAm ts and RF ts are the response absorption and fission rate m 2 functions for coarse mesh type t, d2 ½RAm ts  and d ½RF ts  are the corresponding variances, respectively. Note in the above equation the covariance between any two response function coefficients are ignored since we assume they are independent. The derivative terms can be approximately calculated as:

@k kðR þ DRÞ  kðRÞ  @R DR

ð27Þ

However, this is very time-consuming since it requires performing a whole-core calculation for a perturbation to each response coefficient. For example, if the expansion order is {4, 4, 2, 2} and the 2-group cross sections are used, the total number of response coefficients are about 3.5 million. As a result, it is impractical to use this method. In this work, we estimate the variances of the partial current moments and the global eigenvalue as: 0

d2 ½J þ;m ¼ is0

X mm0 2 2 ;m X 2 mm0 ;m 2 ðRtðiÞss0 Þ d ½J is  þ d ½RtðiÞss0 ðJ is Þ s;m

8 > < X>

RF m d ½k tðiÞs ¼ P 0 ;m0 2 > 0 0 RF m k 0J 0 i ;s ;m tðiÞs0 i;s;m > is : 2

þ

8 > < X> i;s

þ

> > :

þ

P

dm0 m0 ;m0 i0 ;s0 ;m0 J is0 RAtðiÞs0 þ

þP

1 0

0

i

;s0 ;m0

( X X

m0

J ;m RAtðiÞs0 þ is0

P V i0 s 0

J;m is P ;m0 m0 i2It i0 ;s0 ;m0 J i0 s0 RF t 0 s0

t;s;m

92 > > = ðJ þ;0  J ;0 Þ> ; is0 is0 >

>

P V i0 s0

92 > > = ðJ þ;0  J ;0 Þ> ; is0 is0 >

d2 ½J ;m is 

d2 ½J þ;0 is 

)2 d2 ½RF m ts 

8 > t;s;m > :

ð28Þ

s;m

J ;m is P m0 ;m0 i2It i0 ;s0 ;m0 J i0 s0 RAt0 s0 þ

P V i0 s0

92 > > = ðJ þ;0  J ;0 Þ> ; i0 s0 i0 s0 >

d2 ½RAm ts 

unrealistic boundary condition used in that reference which was intended to make the benchmark problem for benchmarking diffusion methods. (2) Pounders et al. (2011) described the problem in its half core configuration. We wish to expand the problem to full (whole) core configuration to verify the accuracy and robustness of the COMET method in calculating the detailed fuel pin fission density distribution. However, as discussed in that paper, generating the reference Monte Carlo pin fission density distribution even in the half core was proven to be challenging and therefore limited to a few selected channels. To avoid this issue, we chose the lower half of the core in the axial direction and then expanded the resulting configuration to full core with quarter core symmetry. This allowed modeling of only 1/4th of the core in MCNP in which case the entire fuel pin fission density distribution were tallied with acceptable statistical uncertainty. It is recognized that these simplifications are not realistic. However, it is clear that they do not affect the problem heterogeneity and therefore the benchmark problem serves as a valid evaluation of COMET’s accuracy and efficiency. The COMET solution was generated for the full core instead of quarter core as was the case in MCNP to determine its efficiency. A brief description of the problem is given below for the sake of completeness and the reader’s convenience. For a detailed description, the reader is referred to Pounders et al. (2011). The x–y (radial) and x–z (axial) planes of the core model are shown in Figs. 1 and 2. The problem consists of 380 (radial) 12 (axial) fuel bundles. The benchmark problem contains 4560 fuel bundles and 168,720 fuel pins. The problem is symmetric about the horizontal and axial mid-plane surfaces of the core. The peripheral fuel channels are surrounded by a reflector as shown in Fig. 1. Vacuum boundary condition is used on all external surfaces. We note that corrugated external boundary and arbitrarily general boundary conditions are not a limitation in COMET. The box type reflector geometry was intentional as discussed in change (1) above. The same burnup distribution as in Pounders et al. (2011) is used in this core. This reference used eight distinct bundle averaged burnup points to make up the core. These are 32.69, 78.38, 342.37, 818.87, 1638.73, 3608.15, 6381.44 and 8721.49 MWd/tU. For the convenience of the reader the detailed burnup distribution is provided in Appendix A. The CANDU-6 NU (natural uranium) bundle used in this benchmark problem is representative of the 37-fuel-element bundle as shown in Fig. 3. The dimension of the fuel bundles is

ð29Þ

where It is the set of all coarse meshes in the core that are of the same type t. It should be noted that the correlation among the response functions for the same coarse mesh type appearing at different locations in the core is taken into account. But the covariances among the incoming/outgoing partial current moments are ignored in the above calculations because its effect on the solution uncertainty is expected to be small and not worth the additional computational expense. 3. A CANDU benchmark problem A highly stylized 3D CANDU benchmark problem derived from Pounders et al. (2011) is used to evaluate the accuracy and computational efficiency of the COMET method. The radial geometric configuration and burnup distribution in this paper are the same as the benchmark problem in Pounders et al. (2011) except for the following changes. (1) The heavy water reflector surrounding the core was extended to a full box. This was done to avoid the

Blue: Fuel

Red: Adjuster Rod

White: Moderator/Reflector

Fig. 1. The x–y plane (radial slice) of the 3-D CANDU benchmark problem.

6

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11 Table 1 Horizontal locations of adjuster rods. Horizontal location

Distance to the mid-plane (cm)

1 2 3 4 5 6 7

171.45 114.3 57.15 0 57.15 114.3 171.45

Table 2 Adjuster rod dimensions. Fig. 2. The x–z plane (axial slice) of the 3-D CANDU benchmark problem.

Shim OR Steel tube IR Steel tube OR Guide tube IR Guide tube OR

Inner element [0.000–85.725 cm]

Outer element [85.725–171.450 cm]

0.650 3.607 3.725 4.519 4.572

0.710 3.607 3.690 4.519 4.572

is six-bundle long (171.45 cm). Each rod consists of the following two vertical segments with slightly different diameters: (1) the inner segment extending from the mid-plane to 85.725 cm (3-bundle lengths), and (2) the outer segment extending from the end of the inner segment to a point 171.45 cm from the mid-plane. The horizontal locations and the dimensions of the inner and outer adjuster rod segments are given in Tables 1 and 2. 4. Numerical results 4.1. MCNP reference calculation

Fig. 3. CANDU-6 lattice representation.

28.575 28.575 49.53 cm. The central fuel pin is surrounded by 3 rings of fuel pins. For each fuel element, natural uranium is contained in zirconium cladding and then surrounded by heavy water coolant. All fuel pins and surrounding heavy water coolant are contained in a pressure tube, beyond which is a calandria tube. The small air gap between the pressure and calandria tubes thermally isolates the pressurized heavy water inside the pressure tube from a low-pressure cool moderator surrounding the calandria tube. Two coolant states are considered in this paper; namely the cooled and fully voided states. The cooled state corresponds to the coolant and moderator at nominal hot operating condition, while the fully voided state corresponds to a postulated accident case in which the coolant density is reduced to 0.001 g/cm3 in all channels. The following two core configurations are modeled in this work. The first configuration corresponds to the case in which all adjuster rods are removed from the core. In the second configuration, 21 adjuster rods parallel with the y axis (vertical direction) are inserted as shown in Figs. 1 and 2. It is assumed that only one type of adjuster rods is present. This array of adjuster rods exists on three axial planes corresponding to z = 217.49 cm, 297.179 cm, and 377.18 cm, where z = 0 is the exterior boundary of axial plane 1. Each adjuster rod consists of a stainless steel cylinder and a steel shim. The adjuster rods are confined inside a zirconium guild tube. The gaps among the steel cylinder, shim and guide tube are filled with heavy water to provide adequate cooling to the rods. Their length

MCNP calculations were performed on a 22-CPU computer cluster. The 2-group cross sections provided in Pounders et al. (2011) were used in the MCNP reference calculations and the response

(a) z=19.059 cm

(b) z=49.53 cm

Fig. 4. Adjuster rod coarse mesh with varying position in (x–z) view.

7

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11 Table 3 Eigenvalue results including one-sigma uncertainty and running times. Core configuration

b c

Voided, without adjuster rods

Cooled, with adjuster rods

Voided, with adjuster rods

MCNP

Eigenvalue Timea (h)

1.01974 ± 0.00002 61.8

1.03613 ± 0.00002 62.1

1.00558 ± 0.00002 62.3

1.02218 ± 0.00002 62.7

COMET

Eigenvalue Timeb (h)

1.02065 ± 0.00028 1.13

1.03679 ± 0.00027 1.16

1.00653 ± 0.00028 1.15

1.02275 ± 0.00027 1.18

91 ± 28

66 ± 27

95 ± 28

57 ± 27

Difference (pcmc) a

Cooled, without adjuster rods

MCNP calculations for the core were performed on a 22-CPU computer cluster. COMET calculations for the whole core were performed on a single CPU. Per cent mille (one-thousandth of a percent).

Table 4 Summary of differences between MCNP and COMET predicted bundle averaged fission densitya.

a

Table 5 Summary of differences between MCNP and COMET predicted fuel pin fission densitya.

Core configuration

AVG (%)

MRE (%)

RMS (%)

MAX (%)

Core configuration

AVG (%)

MRE (%)

RMS (%)

MAX (%)

Cooled, without adjuster rods Voided, without adjuster rods Cooled, with adjuster rods Voided, with adjuster rods

0.24 0.20 0.37 0.48

0.21 0.18 0.33 0.44

0.32 0.24 0.43 0.58

1.7 1.3 1.4 1.5

Cooled, without adjuster rods Voided, without adjuster rods Cooled, with adjuster rods Voided, with adjuster rods

0.36 0.30 0.44 0.55

0.32 0.28 0.38 0.50

0.42 0.38 0.49 0.64

3.2 2.6 3.7 3.0

The total number of fuel bundles in the benchmark problems is 4560.

a

The total number of fuel pins in the benchmark problems is 168,720.

Fig. A.1. Core burnup distribution and bundle indexing.

8

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

function generation mentioned in Section 4.2. The fission source for each configuration was converged by running 2500 inactive cycles with 200,000 particles per cycle. For the first case, an initial source guess was created by specifying source particle locations within each fuel pin at each of the six axial levels using MCNP’s ‘‘ksrc’’ card. Since every case has the same geometry, a converged source from a previous run was used as the initial guess to accelerate the source convergence of the other cases. Using the fission source, 1000 inactive cycles and 4500 active cycles were followed with 200,000 particles per cycle (for a total of 900 million active particles) to generate the reference solution including the core eigenvalue, the bundle averaged and fuel pin fission density distributions.

4.2. COMET calculation The following 22 unique coarse meshes were used to model the cooled and voided core configurations in COMET: 8 coarse fuel meshes with the coolant at hot operating condition, 8 coarse fuel meshes with the coolant voided, 4 coarse meshes containing an adjustor rod and guide tube with different radii at different positions, and 2 coarse meshes of reflector/moderator with different sizes. All the fuel meshes used in COMET have the same size, namely, 23.775 cm 28.575 cm in the radial direction and 49.53 cm the

axial direction. The four fuel rings in each bundle were modeled exactly the same as in the whole core MCNP model. The adjuster rods were modeled in COMET using coarse meshes of 9.6 cm 28.575 cm 49.53 cm. There are two different rod dimensions (inner and outer segments) and two unique rod positions (z = 19.059 cm and z = 49.53). The adjuster rods at the other two locations (z = 30.471 cm and z = 0) are not unique because of the symmetry. This led to a total of four unique controlled meshes. Coarse meshes with two adjuster rod positions are shown in Fig. 4. It should be noted that the controlled meshes next to the top surface contain only half of an adjuster rod (i.e., the rod is cut in half), parallel to the direction of the guide tube as seen in Fig. 4b. For the COMET calculations, the method for estimating the incident flux expansion coefficients was implemented into the MCNP code. Using this modified version a library of response functions was generated that included all of the unique coarse meshes. For each response function (RF) calculation, 8 million histories were followed in about 6 min on a single CPU. The relative uncertainty of the 0th moment of response functions were less than 0.3%. It should be pointed out that the response functions for all unique coarse meshes were generated in the pre-computation phase with the maximum order of (4, 4, 2, 2). That is, the expansion orders in the two spatial and two angular variables were 4, 4, 2 and 2, respectively. A detailed description of the number of response function calculations and pre-computation time can be found in Appendix B. Note that the response function library generation is

Fig. A.2. Core burnup distribution and bundle indexing (continued).

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

a one-time pre-computation. Since each response function calculation is independent, it is very suitable for parallel computing since no communication are required among CPUs. This library generation process is similar to that of pre-computing the cross section library needed for the current whole core diffusion calculations. A number of techniques such as a perturbation method are currently under development to significantly improve the computational efficiency in response function generation. The discussion of these methods is outside of the scope of this paper. Using the response function library full (whole) core COMET calculations were performed to determine the eigenvalue, bundle averaged and fission density distributions in the benchmark problem in both cooled and voided states. The core eigenvalue results are compared to those calculated with MCNP in quarter core configuration with specular reflective boundary condition in Table 3. It is found that the core eigenvalues predicated by COMET agree very well with those of MCNP. The differences ranged from 57 to 95 pcm. The COMET whole-core calculations on a single CPU took slightly more than 1 h for each configuration, while the MCNP quarter-core calculations on a 22-CPU computer cluster took about 60 h each, in addition to 32 h to converge the fission source. Note that the pre-computation of the response function library took 130 h on a 22-CPU computer cluster. As mentioned in Section 2.4, the COMET code underestimates the uncertainty in its results. This indicates that the differences observed are likely within 2–3 standard deviation of the COMET uncertainty.

9

In addition to the core eigenvalue, the quantities of interest in the comparison between the COMET and MCNP solutions are the relative difference (AVG), root mean square (RMS), mean relative difference (MRE) and maximum (MAX) difference in the bundle and pin fission densities. These statistical quantities are defined in the following equations:

AVG ¼

X jFDCM  FDMC j i i FDMC i

i

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u 1 X FDCM  FDMC 2 i i RMS ¼ t N i FDMC i P MRE ¼

CM i jFDi

P

MAX ¼ max i

 FDMC i

MC i FDi CM jFDi  FDMC i j FDMC i

ð30Þ

ð31Þ

ð32Þ ð33Þ

In the above equations, N is the total number of fuel bundles or fuel pins in the problem, FDMC and FDCM are the bundle/pin fission i i densities computed by MCNP and COMET, respectively. The differences between MCNP and COMET calculated bundle averaged and fuel pin fission density values are summarized in Tables 4 and 5.

Fig. A.3. Core burnup distribution and bundle indexing.

10

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

It can be seen from Table 4 that the bundle fission density values predicated by the two codes are in excellent agreement, with the average relative difference of 0.24% and the maximum relative difference of 1.7% in all configurations, respectively. It can also be seen that when the adjuster rods are inserted, the discrepancies between the COMET and MCNP bundle fission density increase slightly. This is expected since the flux spatial gradients become stronger in the controlled cores. It is clear from Table 5 that the agreement in fuel pin fission density between COMET and MCNP is excellent. The average relative difference varies from 0.36% for the uncontrolled cooled case to 0.55% for the controlled voided case, while the maximum relative difference ranges from 2.6% to 3.7%. The uncertainty of the MCNP pin fission densities varies from 0.02% to 1.3%, while the COMET uncertainty is about 0.4%. The largest discrepancy between the COMET and MCNP pin fission densities occurs at a peripheral bundle next to the bottom boundary where fission density is low and calculation uncertainty is large. The difference is still within 3 standard deviations of the MCNP uncertainty as expected. These comparisons indicate that the expansion order {4, 4, 2, 2} is sufficient for representing the actual phase distribution of neutrons on each mesh boundary.

5. Conclusion A new incident flux response expansion method has been developed to significantly improve the accuracy of the hybrid stochastic/deterministic coarse mesh transport (COMET) method. Also, two acceleration techniques were introduced that increase the computational efficiency of the method significantly, by several folds (at least a factor of 20). The new expansion method addresses singularities associated with the original method in which the use of incident current expansion degrades accuracy and efficiency and the method’s ability to solve realistic problems with complexity and size that are inherent in operating commercial reactors. The new method also paves the way for COMET to be imbedded in low order transport methods (e.g., diffusion theory) for improving accuracy without degrading the efficiency of the combined method. In general, the new expansion method also enables efficient and accurate coupling of different deterministic methods (e.g., characteristic to discrete ordinates and in general high order transport to high or low order transport). The new method improvements have enabled COMET to perform whole-core neutronics analysis in all light and heavy water operating reactors with Monte Carlo fidelity and efficiency that is several orders of magnitude faster than both direct Monte Carlo and fine mesh transport methods. The computational efficiency and accuracy of the COMET method was assessed for neutronics calculations in a stylized CANDU-6 benchmark problem in full core configuration with and without adjuster rods. The core consisted of 4560 fuel bundles consisting 168,168 fuel pins. Two coolant sates were considered, namely the hot operating condition and a postulated accident condition. The model included a total of 21 adjuster rods. The COMET model used a very coarse-mesh grid. Comparisons to the MCNP reference solutions confirmed that the method has accuracy close to that of Monte Carlo, while achieving a significantly faster computational speed (about three orders of magnitude). Future work will include benchmarking COMET in full core configurations typical of current operating reactors including pressurized water reactors (PWR), boiling water reactors and additional heavy water reactors (CANDU-6) problems with realistic boundary condition and burnup distribution. Other interesting works are to extend the COMET method to hexagonal geometry typical of prismatic gas cooled and fast reactors and cylindrical geometry typical of pebble bed reactors. Since the method is based on the incident

flux response expansion, a natural extension would be to couple different solution methods (e.g., diffusion to transport) to improve the accuracy of low order methods while maintaining diffusion theory computational efficiency. Appendix A. Appendix A: Core burnup distribution The bundle-averaged burnup values in the present benchmark core have been rounded to one of eight discrete points. Figs. A.1– A.3 present the burnup distribution (taken from Pounders et al., 2011) in the core as indexed in Table A1. Appendix B. Appendix B: Response function generation All response functions were generated with 4th order expansion in space and 2nd order expansion in both polar and azimuthal angles. As mentioned in Section 2.2, the response functions for a coarse fuel mesh also depend on the core eigenvalue. As a result, the response functions for all unique fuel meshes were generated at three distinct core eigenvalues (0.95, 1.0 and 1.05). The total number of response function calculation and CPU times are shown

Table A1 Burnup legend. Index

Burnup (MWd/TU)

1 2 3 4 5 6 7 8

32.69 78.38 342.37 818.87 1638.73 3608.15 6381.44 8721.49

Table B1 Number and CPU time for response function (RF) calculations. Expansion order x y Number of spatial expansion functionsa Polar Azimuthal Number of angular expansion functionsb Energy (group) Total number of expansion functions Coarse fuel meshes Number of unique coarse meshes Number of unique incident Surfacesc Number of distinct core eigenvalues Number of RF calculations Adjuster rod coarse meshes Number of unique coarse meshes Number of unique incident surfaces Number of RF calculations Reflector/moderator coarse meshes Number of unique coarse meshes Number of unique incident surfacesc Number of RF calculations All coarse meshes Total number of RF calculations CPU time for each RF calculation (min) Total CPU time for all RF Calculations (h) Computational Time on 100 CPUs (h)

4 4 15 2 2 6 2 180 16 3 3 25920 4 6 4320 2 3 1080 31,320 6 3132 31

a All cross terms with a combined order higher than 4 are discarded since the numerical experiments have shown they do not improve the accuracy. b All cross terms with a combined order higher than 2 are discarded for the above-mentioned reason. c The symmetry of the coarse mesh is used.

D. Zhang, F. Rahnema / Annals of Nuclear Energy 41 (2012) 1–11

in Table B1. As it can be seen from Table B1, the total number of response function calculations performed is 31320. It took about 31 h on a 100-CPU cluster to generate the response function library. References Briesmeister, J.F., 1997. MCNP – A General Monte Carlo N-Particle Transport Code, Version 4B. Los Alamos National Laboratory, LA-12625-M. Forget, B., Rahnema, F., 2006a. COMET Solutions to the 3D C5G7 MOX Benchmark Problem. Progress in Nuclear Energy 48 (5), 467–475. Forget, B., Rahnema, F., 2006b. COMET Solutions to Whole Core CANDU-6 Benchmark Problem. PHYSOR-2006: Advances in Nuclear Analysis and Simulation, American Nuclear Society’s Topical Meeting on Reactor Physics, Vancouver, BC, Canada. Forget, B., Rahnema, F., 2006c. COMET solution in a highly heterogeneous boiling water reactor benchmark problem. Transactions of the American Nuclear Society 95, 709.

11

Lewis, E.E., Palmiotti, G., Taiwo, T.A., Blomquist, R.N., Smith, M.A., Tsoulfanidis, N., 2005. Benchmark specifications for deterministic MOX fuel assembly transport calculations without spatial homogenization (3-D extension C5G7 MOX). Nuclear Energy Agency. Mosher, S., Rahnema, F., 2006. The incident flux response expansion method for heterogeneous coarse mesh transport problems. Transport Theory and Statistical Physics 35 (1), 55–86. Pounders, J., Rahnema, F., Serghuita, D., Tholammakkil, J., 2011. A 3D stylized halfcore CANDU Benchmark problem. Annals of Nuclear Energy 38, 876–896. Saad, Y., 1992. Numerical Methods for Large Eigenvalue Problems. Halstead Press, New York. Smith, K.S., 1986. Assembly homogenization techniques for light water reactor analysis. Progress in Nuclear Energy 17, 303. Stamm’ler, R.J.J., Abbate, M.J., 1983. Methods of Steady-state Reactor Physics in Nuclear Design. Academic Press, London-New York.