Multi-group unified nodal method with two-group coarse-mesh finite difference formulation

Annals of Nuclear Energy 35 (2008) 1975–1985

Contents lists available at ScienceDirect

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

Multi-group unified nodal method with two-group coarse-mesh finite difference formulation Tae Young Han a, Han Gyu Joo a,*, Hyun Chul Lee b, Chang Hyo Kim a a b

Department of Nuclear Engineering, Seoul National University, San 56-1, Sillim-dong, Seoul 151-744, Republic of Korea Korea Atomic Energy Research Institute, 150 Deokjin-dong, Daejeon 305-353, Republic of Korea

a r t i c l e

i n f o

Article history: Received 30 March 2008 Received in revised form 1 July 2008 Accepted 9 July 2008 Available online 22 August 2008

a b s t r a c t The one-node kernels of the unified nodal method (UNM) which were originally developed for two-group (2G) problems are extended to solve multi-group (MG) problems within the framework of the 2G coarsemesh finite difference (CMFD) formulation. The analytic nodal method (ANM) kernel of UNM is reformulated for the MG application by adopting the Padé approximation to avoid the similarity transform required to diagonalize the G  G buckling matrix. In addition, a one-node semi-analytic nodal method (SANM) kernel which is considered adequate for multi-group calculations is also integrated into the UNM formulation by expressing it in the form consistent with the other UNM kernels. As an efficient global solution framework, the 2G CMFD formulation with dynamic group condensation and prolongation is established and the performance of the various MG kernels is examined using various static and transient benchmark problems. It turns out that the SANM kernel is the best one for MG problems not only because it retains accuracy comparable to MGANM with a shorter computing time but also because its accuracy or its convergence does not depend on the eigenvalue range of the buckling matrix of the system. The 2G CMFD formulation with MG one-node UNM kernels turns out to be very effective in that it conveniently accelerates the MG source iteration. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction There is a long history in the development of the multi-group (MG) nodal methods (Lawrence, 1986). One of the earliest work was done by Shober (1979) 30 years ago for the solution of the time-dependent neutron diffusion problems by the two basic nodal methods – the analytical method and the polynomial method. One of the most recent work is the MG analytic coarse-mesh finite difference method (Analytic CMFD or ACMFD) (Aragones et al., 2007). The reason why there have been continued developments for the solved problem is that there is still the need for a faster and yet more accurate MG nodal code particularly for the transient core simulation involving many neutron groups. This work is another approach to establish an efficient and versatile MG nodal solution method. The unified nodal method (UNM) (Lee and Kim, 2001; Lee et al., 2004a) provides a unified formulation with which not only the transverse-integrated nodal methods (TINMs) such as the nodal expansion method (NEM) (Finnemann et al., 1977) and the analytic nodal method (ANM) (Smith, 1979) or the nodal integration method (NIM) (Fischer and Finnemann, 1981), but also the non-TINM such the analytic function expansion nodal (AFEN) method (Noh and Cho, 1994) can be realized in exactly the same way. It has a * Corresponding author. Tel.: +82 2 880 9241; fax: +82 2 889 2688. E-mail address: [email protected] (H.G. Joo). 0306-4549/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2008.07.004

couple of advantages over the conventional formulations. The nodal kernels can be easily switched simply by replacing the response matrices according to the local core condition. The singularity problem that may arise during the solution process of ANM can be removed easily and effectively (Lee and Kim, 2001). However, the original formulation of UNM has disadvantages in that it was limited to two-group (2G) problems and only the three nodal kernels mentioned above were considered. Lee proposed MG two-node UNM kernels (Lee et al., 2004b) and compared their performance with that of a two-node SANM (Esser and Smith, 1993) kernel. This work, however, was limited only to the local MG two-node problem, not considering the whole reactor problem. In this work, the original one-node 2G-UNM formulation for each kernel is extended to solve MG problems and the one-node SANM kernel is incorporated into UNM by reformulating it in the form consistent with the other UNM kernels. The underlying principle of UNM is to use a linear combination of 5 basis functions as the solution to the TI1D MG equation with different basis functions assigned for the different nodal option or kernel. For example, the transcendental functions consisting of trigonometric and/or exponential functions are used as the basis function in case of ANM. It was shown that the ANM option of UNM gives exactly the same solution as the original ANM solution (Lee and Kim, 2001). In principle, the 2G-UNM formulation can be directly extended to MG without any theoretical problem. In case of ANM, however, the direct extension of UNM involves the

1976

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

diagonalization of the buckling matrix, B2, which requires the solution of a G  G matrix eigenvalue problem with G being the number of groups. The diagonalization and associated ANM solution process are quite cumbersome if G is greater than 2 as will be shown in the next section. Although there are some codes such as SIMULATE-4 (Bahadir et al., 2005) which employs the direct diagonalization scheme to solve the TI1D MG equation by ANM, a simpler yet equally accurate way is desired to enhance the computational efficiency and reduce the complication, particularly when G is considerably large. As a practical method to overcome the diagonalization problem, the Padé approximation of the matrix exponential was introduced by Noh et al. (1997). Once the Padé approximation is adopted in MGUNM, the MG extension of the 2G-UNM can be done straightforwardly even for ANM although the inversion of a full G  G matrix is required because of the coupling between groups in the scattering and fission matrices. The first subject of this paper primarily deals with the MG extension of the ANM option of UNM. The Padé approximation can be a practical solution in case of MGANM, yet it may not be the best solution because the direct inversion of a G  G matrix becomes increasingly inefficient as G becomes larger. Note that the computational overhead required to invert a G  G matrix increases in the third order on G. There is also an inherent problem to the Padé approximation that it is limited to a certain range of problems depending on the types of the approximation. As a way to overcome the inversion problem and the inherent problem of the Padé approximation, one may consider group decoupling by introducing a polynomial approximation of the fission and scattering source. The analytic solution for a given polynomial form of the source will then be obtained easily as a combination of the exponential and polynomial functions, as in the ILLICO formalism (Rajic and Ougouag, 1989) and the semi-analytic nodal method (SANM) solution (Zimin and Ninokata, 1998; Fu and Cho, 2002; Yoon and Joo, 2008). The second subject of this paper is to integrate SANM into MGUNM for more efficient MG solution. Regarding the order of the polynomial approximation of the source distribution in SANM, a second order polynomial had been traditionally used (Rajic and Ougouag, 1989), but a fourth order scheme was later introduced (Kim et al., 1999). Yoon and Joo (2008) showed that the conventional second order scheme suffers nontrivial errors in some MOX loaded cores whereas the fourth order scheme does not. Since the addition of the two source terms is trivial in MGUNM, a fourth order source approximation will be used here. The original UNM was derived employing the one-node formulation in which the incoming currents are specified and the outgoing currents are updated as the result of the nodal solution. Repeated updates of the partial currents and node average fluxes would give the converged solution, but the convergence would be very slow. The coarse-mesh finite difference (CMFD) (Smith, 1983) formulation can be used to accelerate the one-node nodal solution. In principle, the CMFD formulation should be done with the MG structure in order to accelerate the MG nodal solution. However, it is possible to adopt a 2G CMFD in case of the one-node MG nodal formulation since the nodal MG spectrum updated from the MG one-node nodal solution can be used to generate dynamically condensed 2G nodal cross sections for the 2G CMFD calculation. This idea has been already demonstrated for the MG nodal calculation in the hexagonal geometry (Joo et al., 2001) as well as in the rectangular geometry (Han et al., 2006). Since a 2G CMFD solver is readily available in any 2G nodal code employing CMFD and can take the advantage of the Wielandt eigenvalue shift acceleration, the utilization of the 2G CMFD solver in the MG nodal calculation would be quite advantageous. The third subject of this paper is to incorporate the MGUNM one-node solution into the 2G CMFD formulation. Once this is accomplished, the MG calculation is performed only at the local one-node level.

The MGUNM with 2G CMFD can be applied to transient problems as well as the static problems. The UNCARDS (Unified Nodal Code for Advanced Reactor Design and Simulation) (Lee et al., 2004a) code, which had already the 2G transient calculation capability and a 2G CMFD solver, was extended to implement the MGUNM. The performance of the new method is examined with UNCARDS first for various MG static eigenvalue problems involving both thermal and fast reactor spectra. The MG transient calculation performance is then assessed for a MG control rod ejection benchmark problem. In the following section, the basic formulation of MGUNM is given first and the Padé approximation of the matrix exponential is introduced for the ANM kernel. The derivation of the SANM kernel is then followed in Section 3 and the 2G CMFD formulation with the one-node MGUNM kernel is presented in Section 4. The performance test in terms of solution accuracy and efficiency is given in Section 5. 2. Basic formulation of multi-group UNM Consider the neutron balance equation for a spatial node m: m m r  Jm g ðrÞ þ Rrg /g ðrÞ ¼

G vg X

keff

mRmfg0 /mg0 ðrÞ þ

G X

m m Rm sgg 0 /g 0 ðrÞ þ Sg ðrÞ;

g 0 ¼1 g 0 6¼g

g 0 ¼1

ð1Þ where m m Jm ðg ¼ 1; . . . ; GÞ: g ðrÞ ¼ Dg ðrÞr/g ðrÞ

The groupwise form of Eq. (1) can be put into the matrix form:

Dm r2 /m ðrÞ þ Am /m ðrÞ ¼ Sm ðrÞ

ð2Þ

where T m /m ðrÞ ¼ ½/m 1 ðrÞ; . . . ; /G ðrÞ ;

¼

T m Sm ðrÞ ¼ ½Sm 1 ðrÞ; . . . ; SG ðrÞ ;

Dm

m diag½Dm 1 ; . . . ; DG ;

and

2 6 Am ¼ 6 4

Rm r1

0 ..

0 2

0 6 6 m 6R 4 sg1

.

Rm rG Rm s1g ..

.

Rm Rm sG1 sGg

3

2

7 6 7 1 6 5 keff 4

Rm s1G

3

v1 mRmf1    v1 mRmfG .. .

..

.

.. .

vG mRmf1    vG vRmfG

3 7 7 5

7 7

7 Rm sgG 5 0

The superscript T denotes the transpose of a vector or a matrix. Eq. (1) can be integrated over the node to yield the nodal balance equation which involves two types of unknowns: the node average flux and the surface average currents. In order to solve the nodal balance equation, an additional coupling relation between the node average flux and the surface average currents is required. It can be obtained by the transverse-integrated nodal formulation given below. The node index m will be omitted from now on for brevity. Note that most of the derivation given in the next subsection follows the derivation of 2G-UNM given in Reference (Lee et al., 2004a) in terms of matrices and matrix functions. Slight changes in the notations are introduced for more clear presentation. 2.1. Generalized one-node solution After the transverse integration of Eq. (2) and multiplication by the inverse of the diffusion matrix, the following equation can be derived:

1977

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985 2



1 d / ðsÞ þ B2 /u ðsÞ ¼ Su ðsÞ  Lu ðsÞ a2u ds2 u

ð3Þ

where

s¼

u au

and B2 ¼ D1 A:

Here s is the normalized independent variable and au is the node width in the u-direction. Rigorously speaking, su instead of s should be used because au depends on the direction, but the subscript u of su is omitted for brevity. /u(s) is the transverse-integrated one-dimensional (TI1D) flux. Lu(s) and Su(s) denote the transverse leakage and the transverse-integrated source divided by diffusion coefficients, respectively. B2 is called the buckling matrix due to the similarity of the left hand side (LHS) of Eq. (3) to the scalar Helmholtz equation. Note that a minus () sign should be attached to B2 for literally correct definition of the buckling matrix. But we neglect the minus sign for simplicity in the later derivation. The transverse leakage term can be approximated either by a quadratic form or by an analytic form (Lee et al., 2004a) and is represented commonly by

Lu ðsÞ ¼ Lu þ

2 X

giu ðsÞLiu ¼

i¼1

2 X

giu ðsÞLiu

ð4Þ

i¼0

where Lu ¼ L0u ¼ ½L01u ; . . . ; L0Gu T and Liu ¼ ½li1u ; . . . ; liGu T ði ¼ 1; 2Þ. The definition of the basis function diagonal matrix giu(s) is given in the preceding UNM paper (Lee et al., 2004a). The transverse-integrated source term in Eq. (3) can be represented in exactly the same way as the transverse leakage. However, it should be noted that the transverse-integrated source appearing in Eq. (3) is a reduced source which contains the division by the diffusion coefficient. In order to distinguish the regular source from the reduced source, the subscript u is attached to the reduced source so that Su ¼ D1 S . The buckling matrix B2 is a full G  G matrix and thus Eq. (3) needs to be diagonalized if an analytic solution is sought. Assuming that the eigenvalues and eigenvectors of B2 are all known, one can define the similarity transform matrix R and the diagonal buckling b 2 ¼ R1 B2 R. Then the diagonalization leads to the followmatrix B ^ u ðsÞ ¼ R1 / ðsÞ: ing equation for the modal flux, / u 2



2 X 1 d ^ b 2/ ^ u ðsÞ ¼ ^iu ðsÞð b /u ðsÞ þ B S iu  b L iu Þ g 2 2 au ds i¼0

ð5Þ

^iu ðsÞ ¼ R1 gu ðsÞR, b where g S iu ¼ R1 Siu and b L iu ¼ R1 Liu . The modal ^ iu ðsÞ, is required to be a diagonal mabasis function for the source, g trix so that all modes can be independently solved. The analytic solution of the above equation consists of the homogeneous solution and the particular solution. The TI1D flux which is obtained as the linear combination of the modal flux as

^ u ðsÞ; /u ðsÞ ¼ R/

ð6Þ

has thus the homogeneous and particular solution components as well. In UNM, however, this distinction of the solution components is avoided and an expansion form of the solution is employed instead. Specifically, the one-dimensional (1D) flux is given in the similar form as the polynomial expansion of NEM:

þ /u ðsÞ ¼ /

4 X

f iu ðsÞCiu :

ð7Þ

i¼1

where the specific requirements of basis function matrices fiu(s) are also given in the preceding UNM paper (Lee et al., 2004a). Note that in case of the NEM option of UNM, f iu ðsÞ ¼ hi ðsÞI is a diagonal matrix with hi(s) being the ith NEM polynomial whereas fiu(s) becomes a full matrix for the ANM option although its modal counterpart ^f ðsÞ ¼ R1 f ðsÞR is a diagonal matrix. iu iu

The four flux coefficient vectors, Ciu ði ¼ 1; . . . ; 4Þ in Eq. (7) can be expressed in terms of the two surfaces values and the parameters derived from the two weighted residual conditions in exactly the same way as NEM. The resulting expressions for the coefficients were derived as follows in the preceding UNM paper (Lee et al., 2004a) in terms of the symbols defined in the paper:

/ur  /ul ; 2 / þ /ul   /; ¼ ur 2 ¼ M1 3u ðM1u C1u þ DL1u  DS1u Þ and

C1u ¼ C2u C3u

ð8Þ

C4u ¼ M1 4u ðM2u C2u þ DL2u  DS2u Þ: Once the coefficient vectors Ciu of Eq. (7) are determined, the outgoing and the incoming currents at the two nodal surfaces (s = l, r) are then obtained by using the following relation between the partial currents and the TI1D flux: 

jus ¼

 1 1 @/u ðsÞ /us  D 4 2au @ s s¼ss

ðs ¼ r; lÞ:

ð9Þ

Rearranging the outgoing partial currents in terms of the incoming partial currents leads the following response matrix relation: þ   C4u Þ  Q C3u  Q jþ þ Q j and jur ¼ Q 0u ð6/ 1u 2u ul 3u ur    C4u Þ þ Q C3u  Q j þ Q jþ jul ¼ Q 0u ð6/ 1u 2u ur 3u ul

ð10Þ

where the definition of the elements of the diagonal matrix Qiu is given in the first UNM paper (Lee and Kim, 2001). By incorporating Eq. (17) into the nodal balance equation, the following average nodal flux is obtained:

¼ /

X 1 A þ 12 Q a 0u u¼x;y;z u

!1 e S

ð11Þ

where

e S ¼Sþ

X 1  þ ð2Q 0u C4u þ ðI þ Q 2u  Q 3u Þðjur þ jul ÞÞ: a u¼x;y;z u

Eqs. (10) and (11) form the basic coupling relations of the MGUNM which determine the node average flux and the outgoing partial currents. They have exactly the same form as those of the conventional NEM. Therefore, the same solution scheme as the conventional NEM can be adopted to determine the nodal average MG fluxes. Namely, one can determine the outgoing partial currents from the incoming partial currents assumed initially or known from the previous iteration for each node. Sweeping through the whole nodes would result in updated partial currents, which are in turn used to update the nodal average group fluxes. It is possible though to treat the surface flux as unknown as derived in our previous work for SANM (Han et al., 2006). In the following subsections, the specific forms of the coefficients are derived for each nodal option. 2.2. NEM kernel The 2G NEM can be naturally extended to the MG NEM without any modification of the 2G formulation. With the NEM polynomials and corresponding moments, the response matrix components, Miu’s of Eq. (8) appear simply as

M1u ¼ M2u ¼ DB2 ¼ A; 6 1 A; and M3u ¼  2 D  au 10 10 1 A: M4u ¼  2 D  au 14

ð12Þ

1978

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

Therefore, Miu’s are easily determined by using the nodal group constants irrespective of the number of groups. Once Miu’s are determined, the expansion coefficients of the one-dimensional flux, Ciu, can be calculated from Eq. (8) and thus all the matrices appearing in Eq. (11) are also known. The node average flux of Eq. (11) can be determined by solving the G  G linear system using the Gaussian elimination method. This solution procedure for the NEM option of MGUNM is identical to that of 2G-UNM except for the Gaussian elimination. 2.3. ANM kernel With the ANM option, the transverse leakage can be approximated either by a quadratic polynomial or by an analytic function. The derivation of the MGUNM goes differently for the different transverse leakage treatment as shown in the following two subsections. 2.3.1. Quadratic transverse leakage With the quadratic transverse leakage option, the basis function for the transverse leakage is in the same way as the NEM option, namely:

giu ðsÞ ¼ hiu ðsÞI ði ¼ 1; 2Þ:

ð13Þ

With this choice, the particular solution of Eq. (3) becomes a quadratic polynomial as well so that some of the flux basis function should be polynomials. On the other hand, the homogeneous solution of Eq. (3) is obtained in terms of matrix functions, namely:

sÞ ¼ coshð Be sÞp þ sinhð Be sÞq

/hom ð u

ð14Þ

where

e ¼ au B B

ð15Þ

and p and q are coefficient vectors. Note that any matrix function is defined as

f ðAÞ ¼ R diagðf ðki ÞÞR1

ð16Þ

where diagðf ðki ÞÞ is a diagonal matrix consisting of the eigenvalues of A, and R is the matrix consisting of the eigenvectors. Noting this form of the homogeneous solution, we include e sÞ and sinhð B e sÞ in the basis functions. Specifically, the coshð B basis functions are given as

f iu ðsÞ ¼ hiu ðsÞI ði ¼ 1; 2Þ; e sÞ þ b3 s and f 3u ðsÞ ¼ a3 sinhð B e sÞ þ b4 s2 þ c4 f 4u ðsÞ ¼ a4 coshð B

ð17Þ

! !!1 e e e B B 1 B  sinh cosh and 2 2 2 2 ! ! !!1 e e e e B B B B sinh  cosh b3 ¼ sinh : 2 2 2 2

a3 ¼

Once all the coefficients are determined, the response matrix components are then obtained as

where

T4u

K 1 ðxÞ ¼ cosh x;

ð19Þ

K 2 ðxÞ ¼ sinh x

K 3 ðxÞ ¼ xK 1 ðxÞ  K 2 ðxÞ;

and

K 4 ðxÞ ¼ x2 K 2 ðxÞ  3K 3 ðxÞ: It should be noted that the matrix B is a full G  G matrix and the whole manipulation given above can be quite complicated as opposed to 2G-UNM in which the following similarity transform is easy:

" b 2 ¼ R1 B2 R ¼ B

j21

0

0

j22

# ð20Þ

;

where j2i is the eigenvalue of the buckling matrix B2. If j2i is negative, the exponential function involving ji becomes trigonometric functions. After the similarity transform, the response matrix component is obtained as

Miu ¼ DR c M iu R1

ði ¼ 1; 2; 3; 4Þ

ð21Þ

where

  0 b 1 H c b iu ¼ M i1u M iu ¼ G and ju 0 M i2u b ju ¼ R1 Gju R; H b iu ¼ R1 Hiu R G b ju and H b iu , are given in the paper The 2  2 diagonal matrices, G by Lee et al. (2004a). In MGUNM, however, the similarity transform itself becomes nontrivial because it requires the determination of all the eigenvalues of the G  G buckling matrix, B2. In order to avoid this complication, we employ an approximation directly to the matrix functions involved in Eq. (19). Specifically, Miu is approximated by a rational matrix which involves the powers of B2. Specifically, after defining

e¼B e 2 ¼ a2 D1 A; A u

ð22Þ

e in both numerator and we employ a second order polynomial of A denominator to obtain the following Padé approximation:

Miu ¼ DTiu  

1 DF1 Eiu a2u iu

ði ¼ 3; 4Þ

ð23Þ

where

where ai ; bi ; and ci are G  G matrices which can be determined by the constraint conditions given by Eq. (17) of the preceding UNM paper (Lee et al., 2004a). For example,

M1u ¼ M2u ¼ DB2 ¼ A; M3u ¼ DT3u and M4u ¼ DT4u :

   1 1 2 1e 1e B K2 B K3 B ; 2 2 2    1 1 1e 1e ¼ B2 K 3 ; B K4 B 2 2 2

T3u ¼

ð18Þ

e þ 1:97176814  104 A e 2; E3u ¼ 6:00051471I þ 1:33765862  101 A 3 e 6 e 2 F3u ¼ 1:00000000I þ 6:91360728  10 A þ 2:07584742  10 A ; e þ 1:32885391  104 A e 2 ; and E4u ¼ 10:0001490I þ 1:18937674  101 A 3 e 6 e 2 F4u ¼ 1:00000000I þ 4:97435033  10 A þ 1:29734201  10 A : ð24Þ

Note that the Padé approximation has an advantage of less truncation error than a Taylor expansion requiring similar amount of computation, particularly when the magnitude of the eigenvalues of the buckling matrix is considerably off from 0. The above Padé approximation was obtained by assuming 0:05 < k < 2:0 which is a sufficiently wide range of the bucking eigenvalue, k, applicable to light water reactor (LWR) cores. Note that the squared value can be negative for the imaginary numbers. Because Eq. (24) involves linear combinations of G  G matrices, e and A e 2 , Eq. (23) can be manipulated easily although it requires A the computation of an inverse matrix, F1 iu , which can be obtained using the Gauss-Jordan elimination method. In case of 2G, the ANM option of UNM produces exactly the same solution as that

1979

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

of the conventional ANM (Lee and Kim, 2001). On the contrary, the ANM option of MGUNM would involve some error due to the Padé approximation. The errors associated with this approximation, however, turn out to be very small as demonstrated in Section 5. 2.3.2. Analytic transverse leakage In the ANM solution for the transverse-integrated nodal method, the only approximation appears in the shape of the transverse leakage. As an effort to establish a better shape than the quadratic form, Lee introduced the analytic transverse leakage (ATL) approximation when deriving 2G-UNM (Lee and Kim, 2000). In the ATL option for 2G-UNM, the basis functions for the transverse leakage were determined from a 2D intra-nodal flux distribution assumed for the AFEN method (Noh and Cho, 1994). The resulting basis functions are

  1 e g1u ðsÞ ¼ a1 sinh pffiffiffi B s and 2   1 e g2u ðsÞ ¼ a2 cosh pffiffiffi B s þ b2 : 2

ð25Þ

The G  G coefficient matrices ai and bi are determined such that these basis functions satisfy the properties of giu(s) given by Eq. (13) of the preceding UNM paper (Lee et al., 2004a). Note that the factor of p1ffiffi2 which corresponds to sine of 45 degree comes from the corner point of a node (Noh and Cho, 1994). Once the basis functions of the ATL are determined, the particular solution can be represented with the same basis function. e sÞ Thus the basis functions for the flux solution contain sinhðp1ffiffi2 B e sÞ functions as their components. Consequently the and coshðp1ffiffi2 B four basis functions are now given by

   es ; e s þ b1 sinh p1ffiffiffi B f 1u ðsÞ ¼ a1 sinh B 2    1 e e f 2u ðsÞ ¼ a2 cosh B s þ b2 cosh pffiffiffi B s þ c2 ; 2   1 es e sÞ þ b3 sinh pffiffiffi B and f 3u ðsÞ ¼ a3 sinhð B 2    e s þ c4 : e s þ b4 cosh p1ffiffiffi B f 4u ðsÞ ¼ a4 cosh B 2

ð26Þ

and

e E2u ¼ 1:01125635I þ 8:45117290  103 A; e F2u ¼ 1:00000000I þ 5:28353735  103 A; e þ 2:76212731  103 A e 2; E3u ¼ 5:99971002I þ 1:90954942  101 A 2 e 4 e 2 F3u ¼ 1:00000000I þ 6:10005423  10 A þ 1:67088771  10 A ; e þ 1:49778246  103 A e 2 ; and E4u ¼ 9:99995218I þ 1:45991686  101 A 2 e 5 e 2 F4u ¼ 1:00000000I þ 3:81315883  10 A þ 8:15303625  10 A : ð30Þ

With these approximations of Miu ’s, the sole difference between the QTL and ATL options of MGUNM appears only in the coefficients of the Padé approximation. 3. SANM kernel for MGUNM The Padé approximation introduced for the MG-ANM formulation requires taking into account the radius of the convergence. Although the typical conditions of the current LWRs are considered to obtain the coefficients in Eqs. (24) and (30), these are not generally applicable, particularly, to reactor types other than LWR. Noting that all the complications of the analytic solution come from the group coupling, we can consider decoupling the group diffusion equation by moving all the source terms to the right hand side (RHS) and treating them as known from the previous iteration. In the semi-analytic nodal method (SANM) to be derived here, the source on the RHS is represented as a quartic polynomial through a projection. By moving the fission and scattering sources to the RHS, the MG TI1D equation appears as 2



1 d 1 1 / ðsÞ þ B2 /u ðsÞ ¼ D F/u ðsÞ þ D1 E/u ðsÞ a2u ds2 u keff

where the buckling matrix B is now diagonal, namely,

ð27Þ

M4u ¼ DT4u T2u ;

  Rr1 RrG B2 ¼ diag ;...; D1 DG and F and E are now the fission and scattering matrices. All the RHS terms are represented by a single quartic polynomial as

D1



where

0

T1u

T2u

!

!1

e e B B 1 ¼ @I  K 3 pffiffiffi K 2 pffiffiffi K2 2 2 2 2 2 0 ! !1 e e B B 1@ I  K 4 pffiffiffi K 3 pffiffiffi ¼ K3 2 2 2 2 2

!

e B K3 2 !

e B K4 2

ð31Þ

2

M2u ¼ AT2u ;

M3u ¼ DT3u T1u

e E1u ¼ 1:02954752I þ 1:90296243  102 A; e F1u ¼ 1:00000000I þ 1:15599348  102 A;

þ Su ðsÞ  Lu ðsÞ

With these basis functions, the response matrices, Miu, are derived the following form:

M1u ¼ AT1u ;

where

!1 11 e B A and 2 !1 11 e B A : 2 ð28Þ

¼

 1 F/u ðsÞ þ E/u ðsÞ þ Su ðsÞ  Lu ðsÞ  Pu ðsÞ keff

4 X

giu ðsÞPiu

ð32Þ

i¼0

where

g0u ðsÞ ¼ I;

giu ðsÞ ¼ hiu ðsÞI ði ¼ 1; 2; 3; 4Þ:

With this approximation of the RHS, the groupwise form of Eq. (31) appears as 2

with the definitions of Ki(x) functions given by Eq. (19). The Padé approximation yields the following approximate forms: 1 M1u ¼ AT1u  AF1 1u E1u ; M2u ¼ AT2u  AF2u E2u 1 M3u ¼ DT3u T1u   2 DF1 E3u and au 3u 1 E4u : M4u ¼ DT4u T2u   2 DF1 au 4u

ð29Þ



4 X 1 d /gu ðsÞ þ B2g /gu ðsÞ ¼ g giu ðsÞPgiu 2 2 au ds i¼0

ð33Þ

which can be solved independently for each group. The analytic solution of Eq. (33) consists of the cosh and sinh homogeneous solution and a fourth order polynomial particular solution so that there are seven independent functions. In contrast to the NEM or ANM option of MGUNM which requires five basis functions, the SANM option of MGUNM requires seven basis functions. In other words, the TI1D flux of the SANM is represented by the following vector form:

1980

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

þ /u ðsÞ ¼ /

6 X

f iu ðsÞCiu

ð34Þ

i¼1

where

f iu ðsÞ ¼ hiu ðsÞI ði ¼ 1; 2; 3; 4Þ; e sÞ þ b5 s and f 5u ðsÞ ¼ a5 sinhð B

ð35Þ

e sÞ þ b6 s2 þ c6 : f 6u ðsÞ ¼ a6 coshð B Here ai ; bi ; ci are G  G diagonal matrices satisfying the conditions of Eq. (17) of the preceding UNM paper (Lee et al., 2004a) with an extension of the conditions for the third and fourth basis functions to the additional fifth and sixth basis functions. The first four basis functions correspond to the particular solution while the fifth and sixth functions are for the homogeneous solution. The first two coefficients of Eq. (34) are determined by the same way as the other nodal options using the surface flux values. The other four expansion coefficients can be determined using the method of undetermined coefficients after inserting the groupwise form of Eq. (34) into Eq. (33). Then Ciu ’s ði ¼ 3; . . . ; 6Þ are expressed in terms of the coefficients of the effective source, Piu, as 1 C3u ¼ M1 3u DP3u ; C4u ¼ M4u DP4u ;   6 2 2 1 and C5u ¼ M1 5u C1u  Rr DP1u  2 Rr D P3u au   10 2 2 1 C6u ¼ M1 6u C2u  Rr DP2u  2 Rr D P4u : au

ð36Þ

where

M3u ¼ M4u ¼ DB2 Rr ! !1 e e B B 1 K3 M5u ¼ K 2 ; 2 2 2

M6u ¼

! !1 e e B B 1 K4 K3 2 2 2

Unlike the ANM option, the response matrix components, Miu, are easily obtained for SANM because the matrix functions are all diagonal matrices as follows,

K2

au ¼ diag sinh 2

sffiffiffiffiffiffiffi!

Rrg Dg

au ; . . . ; sinh 2

Dg

1u

¼ /

X 1 Rr þ 12 Q a 0u u¼x;y;z u

2u ur

3u ul

ð40Þ

where

X 1 1  e  þSþ S¼ F/ þ E/ ð2Q 0u ðC4u þ C6u Þ þ ðI þ Q 2u  Q 3u Þ keff a u¼x;y;z u 

12

ð42Þ

which produces

e 2u ¼ C2u þ T26u C6u ; e 1u ¼ C1u þ T15u C5u ; C C e 3u ¼ C3u þ T35u C5u and C e C 4u ¼ C4u þ T46u C6u

ð43Þ

where 0 ! !1 1   e e B B 1@ 1 4 2 1 2 e þ 28I  B e e þ 140I K 1 e Þ A; T15u ¼ 3B B I þ 60ð B K3 2 2 2 0 ! !1 1  1  1 e e 1@ B B 4 2 2 e e e A; I þ 420 B T26u ¼ B þ 36I  B þ 252I K 2 K4 6 4 2 2 0 ! !1 1   1 e e B B 1 4 2 2 e Þ e e þ 60I —5 B A and K2 B T35u ¼ 140@ð B K3 2 2 0 ! !1 1  e e B B 4 1 e 2 e @ A B þ 42I  7K 2 T46u ¼ 21 40ð B Þ K4 2 2

The source polynomial is then constructed using the approximated flux as 4 X

giu ðsÞPiu ¼ Pu þ

4 X

i¼0

i¼1 2 X

D1



 1 e iu F þ E giu ðsÞ C keff

giu ðsÞðSiu þ Liu Þ

ð45Þ

i¼1

Pu ¼ D1



 1  þ Su þ Lu : FþE / keff

This gives the source coefficients as follows:

8  > e iu þ Siu þ Liu < D1 k1 F þ E C eff  Piu ¼ > e iu : D1 1 F þ E C keff

!1 e S

h i ~ u ðsÞ ds ¼ 0 ði ¼ 1; 2; 3; 4Þ; giu ðsÞ /u ðsÞ  /

1 2

where

ð39Þ and

Z

ð38Þ

þ   C4u  C6u Þ  Q ðC3u þ C5u Þ  Q jþ þ Q j ; jur ¼ Q 0u ð6/ 1u 2u ul 3u ur    C4u  C6u Þ þ Q ðC3u þ C5u Þ  Q j þ Q jþ ; j ¼ Q ð6/ 0u

ð41Þ

e iu can be determined by the following The coefficients C weighted residual method equations,

þ :

Once all the flux expansion coefficients are determined, the outgoing partial currents and the node average flux are obtained in the same way as the other kernel of MGUNM. With the additional coefficients, C5u and C6u , however, the response relation and the node average RHS terms are modified to

ul

e iu giu ðsÞ C

i¼1

sffiffiffiffiffiffiffi!!

Rrg

4 X

ð44Þ

ð37Þ

!

~ u ðsÞ ¼ / þ /

and

e ¼ au B: B

e B 2

tion is a quartic polynomial, a projection of the semi-analytic function to the polynomial space should be performed to approximate the TI1D flux with a 4th order polynomial before constructing the source (Lee et al., 2004b). Let the approximate TI1D flux be

ði ¼ 1; 2Þ ð46Þ ði ¼ 3; 4Þ

By Eqs. (43) and (46) and, the source term can be updated using the TI1D flux of the previous iteration step under the framework of source iteration. For the SANM option in MGUNM, there is thus another level of iteration within the one-node solver. Since it is not necessary to achieve full convergence of the source shape in the one-node problem during the course of the incoming current iteration, the additional level of source iteration would not be burdensome. In fact, updates of the source shape and the incoming current can be performed in parallel, implying that no additional source iterations are in fact necessary.

þ

 ðjur þ jul ÞÞ: 4. Two-group CMFD formulation The semi-analytic solution represented by Eq. (34) contains the exponential function as well as the polynomial functions. Because the source which is to be constructed from this semi-analytic solu-

The MGUNM formulation presented above solves a one-node problem given the incoming partial currents or surface fluxes of

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

all surfaces. The result of the one-node calculation is the groupwise node average fluxes as well as the outgoing currents. In order to accelerate the MG one-node nodal calculation, the MG CMFD method can be applied in principle. But noting that the groupwise fluxes and partial currents updated by the one-node MGUNM solver can be used to condense the MG cross sections and nodal coupling coefficients into few-groups, we can apply the 2G CMFD directly to accelerate the MG one-node calculation. If 2G CMFD is used for the global neutron balance calculation, the MG calculation needs to be done only at the local one-node level. Therefore, extension of an existing 2G code to MG can be quite simple. The 2G CMFD formulation, however, needs prolongation of the 2G node average flux to the MG average flux and incoming current. The condensation and prolongation procedure is detailed below. Firstly, the multi-group cross sections are condensed into the two groups using the following relations in order to construct the 2G CMFD nodal balance equation.

RxG0 ¼

X

Rxg ug ðG0 ¼ 1; 2Þ

ð47Þ

from the CMFD, the surface flux which is needed to construct the partial current should be determined using an additional relation similar to Eq. (52):



 m r m l m r l /m G0 u ¼ xG0 /G0 u þ 1  xG0 /G0 u þ CG0 /G0 u þ /G0 u :

xmG0 ¼

blG0 u brG0 u

þ blG0 u

:

Similarly to Eq. (52), the first two terms are the approximation of the surface flux by the finite difference method and the third term is the corrective term. The surface flux correction coefficients should also be determined from the one-node MGUNM solution. Once the 2G surface net currents and surface average flux are determined, the two-group partial currents in the interface, m, is obtained by using the following relation. m

jG0 u ¼

1 m 1 / 0  J m0 4 Gu 2 Gu

ð54Þ

The prolongation to the MG structure then can be completed by using the stored flux and partial current spectra. Namely,

where

ug ¼ / g =/ G0

ð48Þ

and

X

ð53Þ

where

g2G0

 G0 ¼ /

1981

g ¼ /  G0  u / g

ð55Þ

and

g /

ð49Þ

g2G0

m

jg

m

¼ jG0  1m g

ð56Þ

0

0

and G is the group index in 2G, or the coarse group index. The summation in Eqs. (47) and (49) goes over the fine groups g’s belonging 0 to the coarse group G . In addition, the multi-group partial current spectrums can be defined as follows: m m 1m ¼ jg =jG0 g

ð50Þ

where m jG 0

¼

X

m

jg

ðG0 ¼ 1; 2Þ

ð51Þ

where G is the two-group index to which the fine group g belongs. After the MG incoming partial currents are updated using the 2G CMFD solution, the one-node MGUNM is performed again to determine new MG node average fluxes and outgoing currents which are to be used eventually to update the 2G cross sections b m0 and Cm0 . The alternative calculation and correction coefficients: D G G of 2G CMFD and one-node MGUNM continues until the solution converges. The calculation flow is depicted in Fig. 1.

g2G0

These flux and partial current spectra obtained after the onenode MGUNM calculation need to be stored for later use in prolongation. For the consistent 2G CMFD formulation, the 2G net current determined by condensing the MG net currents should be conserved in the 2G calculation. This can be achieved by employing the following CMFD relation between the node average fluxes and current:

e m r b m r l l Jm G0 u ¼  D G0 ð/G0 u  /G0 u Þ  D G0 ð/G0 u þ /G0 u Þ

ð52Þ

where r l e m0 ¼ 2bG0 u bG0 u D G brG0 u þ blG0 u

and bsG0 u ¼ DsG0 =asu ðs ¼ l; rÞ:

The superscripts, r and l, denote the left side and right side. The first term on the RHS represents the ordinary finite difference method (FDM) current while the second term is the correction term needed to conserve the nodal current. The current correction b m0 , is obtained from the 2G current of the MGUNM coefficients, D G solution. From the 2G constants and the current versus node average flux relation, the 2G CMFD linear system consisting of the nodal balance equations can be constructed and the solution of the 2G CMFD gives the 2G node average flux distribution. In order to go back to the MG one-node nodal calculation which needs the MG incoming currents as well as the MG sources, the MG outgoing currents and sources need to be constructed from the 2G CMFD result. This prolongation can be done by first finding the 2G partial currents. Since only the net current information is available

5. Performance examination and discussions As mentioned previously, there are four kernels integrated into MGUNM: NEM, MGANM/QTL, MGANM/ATL and SANM. In order to examine the performance of the four kernels and the 2G CMFD formulation, MGUNM was implemented into the UNCARDS code which already had a 2G-UNM and CMFD solution modules. The 2G CMFD solver employs the biconjugate gradient stabilized (BiCGSTAB) algorithm preconditioned with the incomplete lower and upper (ILU) preconditioner (Joo and Downar, 1996) for the inner iteration and the Wielandt shift method for the outer iteration method. Since the UNCARDS code has a 2G transient solution capability, the MG transient solution capability can be readily implemented within the 2G CMFD framework. Only the MG transient fixed source term needs to be specially treated in the one-node MGUNM kernel. At first, a group of 2G and MG eigenvalue benchmarks problems are analyzed. Some of the benchmark problems represent fast breeder reactors (FBR), but most problems are for light water reactors (LWRs). These are the IAEA3D (Argonne National Laboratory, 1977), L336 C5 (Lefebvre et al., 1991), KAIST4G (Cho and Noh, 1995), KOEBERG (Muller and Weiss, 1991), TAKEDA Model 2 (Takeda and Ikeda, 1991), MOX 7G (Lewis et al., 2001). Two transient benchmark problems, the NEACRP PWR core transient benchmark problems (Finnemann, 1992) with 2G cross sections and the NEA MOX transient benchmark problems (Kozlowski and Downar, 2006) with 4-group cross sections, are then analyzed. As the reference solution for the comparison, the result of the ANM/QTL of UNCARDS with 16 nodes per fuel assembly is used.

1982

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

Update XS

Prolong Flux,, Current

Setup CMFD Equation ti

Two-group CMFD

Multi-group

lc Calculation

Nodall Calculation lc Yes

Update Correction Factor? No No

Flux Error Fl rro < ε 1 ?

Flux Change <ε2 ? Min < N < Max ?

No

Yes Condense

Yes

Exit

Flux,, Current, XS Update Correction Factorr

Fig. 1. Calculation flow of 2G CMFD employing MG one-node kernel.

The 2G benchmark problems, IAEA3D and L336 C5, are solved to verify the newly implemented MG solver which can solve the 2G problems as well. Note that the ANM kernel of MGUNM is different from that of 2G-UNM because of the Padé approximation. Therefore, there can be slight differences in the result. Table 1 summarizes the accuracy of all the nodal kernels of UNM. It is worthy to note in Table 1 that the ANM kernel of MGUNM gives essentially the same solution as that of 2G-UNM. This assures that the Padé approximation with the coefficients given in Eqs. (24) and (30) are accurate. As presented in the previous work (Lee et al., 2004a; Lee and Kim, 2000) the ATL option is superior to the QTL option and the same trend is observed with MGUNM. Particularly for the L336C5 problem, the ATL option gives about 4 pcm better solution than QTL for the 4 nodes/FA case. However, the improvement is not considered remarkable in general since there is only negligible difference in the power distribution error. On the other hand, the SANM kernel shows accuracy comparable to ANM/QTL while the NEM error is quite large (more than 50 pcm with 4 nodes/FA) for the L336 C5 problem which involves a severely distorted intra-nodal flux distribution because of the MOX fuel loaded. Table 1 shows the results of the MG benchmark problems as well. In the KOEBERG problem and the TAKEDA Model 2 problem, the accuracy of the four nodal options are similar to each other. However, the NEM kernel is inferior again to MGANM and SANM in the KAIST problem and MGANM/ATL is consistently better than MGANM/QTL or SANM for the problems involving MOX fuels although not remarkable. It is also noted that SANM has almost the same accuracy as MGANM in the MG problems as well. The execution performance of the UNCARDS solvers is compared in Table 2 which shows the number of iterations and computing time. Because all the nodal options of UNM are derived from the same basic formulation, the number of outer iterations and the number of nodal updates, which is to update the current correction coefficients, are independent of the nodal options and

thus a consistent timing comparison between the nodal kernels is possible. However, for the SANM kernel which involves the additional iteration to update the source shape, the number of onenode sweeps can potentially be more than that for the other kernels in which the group coupling is resolved directly. One one-node sweep here means the spatial sweep over all nodes which is to update the incoming currents. In case of SANM, the one-node sweep involves a group sweep as well. The number of outer iterations here is the number of CMFD calculations to update the global 2G flux and fission source distributions, and also the eigenvalue. The fact that the number of outer iterations are within the range of 10–25 for all the problems signifies that the 2G CMFD acceleration is quite effective. The timing data here were obtained for the 4 nodes/FA runs on a PC having a Pentium 3.0GHZ CPU. It is noted in the result for the IAEA3D problem that MGANM requires almost two times more nodal computing times than ANM. It is because there are numerous matrix operations in the Padé approximation whereas the ANM operations for the 2G problems are trivial. SANM requires also more time than ANM because of the additional time required for updating the source terms. However, SANM always takes less time than MGANM regardless of the number of groups. The MGANM/ATL kernel is very good for the MG problems upto seven groups, but its stability and computing time deteriorate rapidly with more groups. The transient calculation results are summarized in Table 3 which shows the steady-state results as well for the 2G NEACRP A1 problem and the 4G MOX transient problem. Through the comparison of the critical soluble boron concentration, the peak power and the peak time for the NEACRP A1 problem, it is noted that the accuracy of MGANM and SANM is similar to that of ANM and the computational speed of SANM is slightly faster than that of MGANM. Since the difference due to different nodal options essentially disappears with 4 nodes/FA calculations for the NEACRP problem, only the results for 1 node/FA are compared in Table 3 and Fig. 2.

1983

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985 Table 1 Accuracy of UNM kernels for eigenvalue benchmark problems Name (Group) Ref. k-eff

Options

1 Node/FA keff

4 Nodes/FA

Dkeff (pcm)

Max. power error (%)

RMS power error (%)

keff

Dkeff (pcm)

Max. power error (%)

RMS power error (%)

IAEA3D (2G) 1.02907

NEM ANM/QTL ANM/ATL MGANM/QTL MGANM/ATL SANM

1.02906 1.02912 1.02909 1.02912 1.02910 1.02912

1.5 4.9 2.1 5.1 2.4 4.9

1.77 1.19 0.78 1.19 0.78 1.19

0.89 0.37 0.25 0.37 0.25 0.37

1.02909 1.02910 1.02909 1.02910 1.02909 1.02910

1.7 2.3 2.0 2.3 2.0 2.3

0.17 0.37 0.21 0.37 0.21 0.56

0.05 0.13 0.08 0.13 0.08 0.16

L336 C5 (2G) 0.93883

NEM ANM/QTL ANM/ATL MGANM/QTL MGANM/ATL SANM

0.93617 0.93906 0.93897 0.93907 0.93897 0.93904

265.8 22.9 13.7 23.4 13.4 21.2

3.11 0.38 0.25 0.40 0.25 0.38

2.03 0.20 0.13 0.21 0.13 0.20

0.93831 0.93888 0.93884 0.93888 0.93884 0.93888

52.1 4.7 0.6 4.7 0.6 4.8

0.55 0.13 0.04 0.13 0.04 0.13

0.36 0.08 0.03 0.08 0.03 0.08

KOEBERG (4G) 1.00796

NEM MGANM/QTL MGANM/ATL SANM

1.00840 1.00838 1.00838 1.00839

43.5 42.3 42.0 42.3

2.53 2.28 2.26 2.28

1.36 1.28 1.28 1.28

1.00800 1.00800 1.00800 1.00800

4.1 3.8 3.8 3.7

0.25 0.22 0.22 0.22

0.14 0.12 0.12 0.12

KAIST (4G) 1.06729

NEM MGANM/QTL MGANM/ATL SANM

1.06666 1.06756 1.06750 1.06756

62.4 27.7 21.3 27.8

1.54 0.27 0.28 0.27

0.90 0.13 0.15 0.13

1.06725 1.06737 1.06734 1.06737

3.3 8.4 5.9 8.3

0.28 0.09 0.04 0.10

0.14 0.05 0.03 0.06

TAKEDA Model 2 (4G 0.96885

NEM MGANM/QTL MGANM/ATL SANM

0.96886 0.96886 0.96886 0.96886

1.5 1.5 1.4 1.5

0.13 0.13 0.13 0.13

0.03 0.03 0.03 0.03

0.96885 0.96885 0.96885 0.96885

0.9 0.9 0.9 1.0

0.02 0.02 0.02 0.02

0.01 0.01 0.01 0.01

MOX 7G (7G) 1.21040

NEM MGANM/QTL MGANM/ATL SANM

1.21052 1.21061 1.21061 1.21067

12.7 20.9 21.5 27.0

5.04 1.62 1.48 1.82

2.22 0.65 0.66 0.74

1.21050 1.21044 1.21043 1.21045

9.9 4.8 3.8 5.1

0.93 0.50 0.39 0.49

0.44 0.18 0.15 0.18

Table 2 Iteration performance of UNM kernels for steady-state benchmark problems Problem

Options

Number of

CPU time (s) for

Nodal sweeps

Corr. factor updates

Outer iterations

IAEA3D (2G)

NEM ANM/QTL ANM/ATL MGANM/QTL MGANM/ATL SANM

34 34 34 34 34 37

10 10 10 10 10 9

12 12 12 12 12 11

0.95 1.03 1.02 2.02 2.08 1.86

1.19 1.17 1.15 1.20 1.19 1.06

2.80 2.80 2.80 4.08 4.17 3.88

0.028 0.030 0.030 0.059 0.061 0.050

L336 C5 (2G)

NEM ANM/QTL ANM/ATL MGANM/QTL MGANM/ATL SANM

28 32 28 32 28 27

8 8 8 8 8 8

10 10 10 10 10 10

0.063 0.062 0.047 0.063 0.062 0.047

0.015 0.016 0.015 0.015 0.016 0.015

0.094 0.110 0.141 0.157 0.141 0.140

0.002 0.002 0.002 0.002 0.002 0.002

KOEBERG (4G)

NEM MGANM/QTL MGANM/ATL SANM

29 29 29 25

8 8 8 8

10 10 10 10

0.219 0.205 0.218 0.202

0.016 0.032 0.016 0.016

0.375 0.390 0.406 0.343

0.008 0.007 0.008 0.008

KAIST (4G)

NEM MGANM/QTL MGANM/ATL SANM

27 29 31 24

7 7 7 7

10 10 10 10

0.062 0.063 0.078 0.048

0.016 0.015 0.015 0.015

0.110 0.109 0.125 0.109

0.002 0.002 0.003 0.002

TAKEDA Model 2 (4G)

NEM MGANM/QTL MGANM/ATL SANM

61 61 61 61

11 11 11 11

14 14 14 14

MOX 7G (7G)

NEM MGANM/QTL MGANM/ATL SANM

39 39 39 42

9 9 9 10

12 12 12 13

In the 4G MOX transient problem, the results are compared with the PARCS results which are given as the reference in

Nodal sweeps

11.6 12.4 12.8 10.3 0.59 0.63 0.69 0.51

CMFD

3.64 3.67 3.66 3.69 0.015 0.062 0.047 0.064

Total

17.4 18.3 18.7 16.4 0.750 0.813 0.843 0.734

A nodal sweep

0.19 0.20 0.21 0.17 0.015 0.016 0.018 0.012

the benchmark (Kozlowski and Downar, 2006). It is shown in Table 3 that all the nodal kernels of MGUNM except NEM produce

1984

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

Table 3 Results for transient benchmark problems Options NEACRP A1 1 Node/FA (2G)

MOX TR 1 Node/FA (4G)

MOX TR 4 Nodes/FA (4G)

Soluble boron (ppm)

Peak power (%)

CPU time (s)

561.0 566.8 562.6 562.7 562.5 562.7 562.6

833.6 807.4 818.3 821.1 820.9 821.1 818.7

0.540 0.735 0.615 0.610 0.610 0.610 0.620

132.8 69.7 98.7 99.8 99.3 99.9 99.2

–

Ref. (PARCS) NEM MGANM/QTL MGANM/ATL SANM

1337.0 1343.0 1344.8 1344.7 1344.8

68.2 67.9 68.1 68.1 68.1

0.334 0.352 0.344 0.346 0.348

152.3 153.1 162.2 158.2 157.2

–

Ref. (PARCS) NEM MGANM/QTL MGANM/ATL SANM

1337.0 1339.7 1340.2 1340.2 1340.2

68.2 68.2 68.2 68.2 68.2

0.334 0.334 0.332 0.332 0.332

152.3 171.9 174.3 173.9 174.6

– 10890.0 11119.5 11604.7 10713.1

Ref NEM ANMQTL ANMATL MGANMQTL MGANMATL SANM

120

100

27.18 27.78 28.44 45.88 46.81 43.32 2730.5 2788.6 2834.6 2384.6

essentially the same results. Compared with the PARCS results which are obtained with the 4 node/FA NEM solver, the UNCARDS results are very close in the peak time but the peak value is somewhat different as identified in Fig. 2. This difference is caused by the difference in the thermal hydraulic solution module of the two codes. Consequently, it is demonstrated from Table 3 and Fig. 3 that the multi-group UNM works properly for the transient problems as well.

140

Power (%)

Peak time (s)

Ref. (16 N/FA) NEM ANM/QTL ANM/ATL MGANM/QTL MGANM/ATL SANM

NEACRP A1 Problem

80

6. Summary and conclusions 60

40

20

0.0

0.5

1.0

1.5

2.0

2.5

Time (sec) Fig. 2. Transient core power for NEACRP A1 problem.

4G MOX Transient Problem 180

Ref. NEM MGANMQTL MGANMATL SANM

160 140 120

Power (%)

Rod worth (pcm)

100 80 60 40 20 0 0.0

0.2

0.4

0 .6

0.8

1 .0

Time (sec) Fig. 3. Transient core power for MOX transient 4-group problem.

In this work the original 2G-UNM formulations of the ANM/QTL, and ANM/ATL kernels were extended for MG applications and the one-node SANM kernel which is considered adequate for the MG calculation was also integrated into the UNM formulation by reformulating it into the form consistent with that of the other UNM kernels. Also an efficient global solution framework was established by employing the concept of dynamic group condensation and prolongation within the framework of 2G CMFD formulation. The performance of the various MGUNM kernels was examined for various static and transient benchmark problems. It was straightforward to extend the NEM formulation to MG applications because the NEM response matrices are defined by a linear function of the buckling matrix. In case of ANM, however, it was needed to introduce a different formulation for efficient MG applications because the direct MG extension of ANM requires the determination of all the eigenvectors of the buckling matrix which poses a considerable computational burden as the number of groups increases. The Padé approximation for the response matrices was thus introduced for MG-ANM with the approximation coefficients valid over the interval of 0.05 < k < 2.0, which represents a sufficiently wide range of the buckling matrix eigenvalues for LWR cores. In the SANM formulation, the neutron diffusion equation for each group was decoupled by introducing a quartic polynomial approximation for the fission source and the scattering source, which makes it natural and straightforward to apply the SANM to multi-group problems. The one-node SANM with a source approximation by a projection from the transcendental function space to the polynomial function space was reformulated into the form consistent with that of the other UNM kernels to integrate it into UNM. As an efficient global solution framework, a 2G CMFD formulation with dynamic group condensation and prolongation was established. In this framework, the MG nodal kernels are used for the intra-nodal MG solutions while the 2G CMFD solver is used

T.Y. Han et al. / Annals of Nuclear Energy 35 (2008) 1975–1985

for the global 2G solution. The 2G cross sections and nodal coupling coefficients needed for the 2G CMFD calculation are generated by dynamic condensation using the MG one-node solution. The MGUNM was implemented into the UNCARDS code and its performance were examined for various steady-state eigenvalue and transient benchmark problems representing fast breeder reactor (FBR) as well as LWR cores with 2G, 4G, 7G cross sections. The ANM kernels of MGUNM gave essentially the same solution as that of 2G-UNM for the 2G problems despite the approximation introduced by the Padé approximation in MGUNM. It was noted that the SANM kernel gives essentially the same accuracy as that of ANM kernels while the NEM kernel is inferior to the others in both 2G problems or MG problems. The MG kernels require much more time than 2G kernels in 2G problems. However, SANM takes less time than the MGANM for MG problems. From the comparisons above, it can be concluded that the SANM would be the best kernel for MG problems not only because it retains accuracy comparable to MGANM with a shorter computing time but also because its accuracy and convergence do not depend on the eigenvalue range of the buckling matrix of the system. Also the 2G CMFD formulation with MG one-node UNM kernels is very effective because it provides a good acceleration performance for the source iteration as identified with a small number of outer iterations. The 2G CMFD formulation is particularly beneficial for transient calculations since only 2G CMFD calculations can be performed at most of transient calculation time steps if a conditional nodal update scheme, which is to update conditionally the local MG spectrum and 2G current correction factors, is implemented. References Aragones, J.M., Ahnert, C., Garcia-Herranz, N., 2007. The analytic coarse-mesh finite difference method for multi-group and multidimensional diffusion calculations. Nucl. Sci. Eng. 157, 1. Argonne National Laboratory, 1977. Benchmark Problem Book. ANL-7416, Supplement 2, Argonne National University. Bahadir, T., Lindahl, S.-Ö., Palmtag, S.P., 2005. SIMULATE-4 Multigroup Nodal Code with Microscopic Depletion Model. In: Proc. Math. Comp. Supercom. Reac. Phys. Nucl. Bio. Appl. Palais des Papes, Avignon, France, September 12–15, on CDROM. Cho, N.Z., Noh, J.M., 1995. Hybrid of AFEN and PEN methods for multi-group diffusion nodal calculation. Trans. Am. Nucl. Soc. 73, 438. Esser, P.D., Smith, K.S., 1993. A semi-analytic two-group nodal model for SIMULATE3. Trans. Am. Nucl. Soc. 68, 220. Finnemann, H., Bennewitz, F., Wagner, M.R., 1977. Interface current techniques for multidimensional reactor calculations. Atomkernenergie 30, 123–128. Finnemann, H., 1992. NEACRP 3D LWR Core Transient Benchmark-Final Specifications. NEACRP-L-335 (Revision 1).

1985

Fischer, H.D., Finnemann, H., 1981. The nodal integration method – a diverse solver for neutron diffusion problems. Atomkernenergie 39, 229–236. Fu, X.D., Cho, N.Z., 2002. Nonlinear analytic and semi-analytic nodal methods for multi-group neutron diffusion calculations. Nucl. Sci. Tech. 39, 1015. Han, T. Y., Joo, H. G., Kim, C. H., 2006. Two-group CMFD accelerated multi-group calculation with a semi-analytic nodal kernel. In: Proc. PHYSOR 2006, ANS Topical Mtg. on Reactor Physics, Vancouver, BC, Canada, September 10–14, CDROM. Joo, H.G., Downar, T.J., 1996. An incomplete domain decomposition preconditioning method for nonlinear nodal kinetics calculations. Nucl. Sci. Eng. 123, 403. Joo, H.G., Cho, J.Y., Song, J.S., Zee, S.Q., 2001. Multi-group transient calculation within the framework of a two-group hexagonal CMFD formulation. M&C 2001, Salt Lake City. Kim, Y.I., Kim, Y.J., Kim, S.J., Kim, T.K., 1999. A semi-analytic multi-group nodal method. Ann. Nucl. Energ. 26, 699. Kozlowski, T., Downar, T., 2006. Pressurized WATER REACTOR MOX/UO2 Core Transient Benchmark, Final Report. NEA/NSC/DOC(2006)20, OECD Nuclear Energy Agency. Lawrence, R.D., 1986. Progress in nodal methods for the solution of the neutron diffusion and transport equations. Prog. Nucl. Energ. 17, 271–301. Lee, H.C., Kim, C.H., 2000. A new transverse leakage approximation for improved nodal reactor calculations. Trans. Am. Nucl. Soc. 83, 411–414. Lee, H.C., Kim, C.H., 2001. Unified formulation of nodal expansion method and analytic nodal method solutions to two-group diffusion equations. Nucl. Sci. Eng. 138, 192–204. Lee, H.C., Chung, K.Y., Kim, C.H., 2004a. Unified nodal method for solution to the space-time kinetics problems. Nucl. Sci. Eng. 147, 1–17. Lee, H. C., Joo, H. K., Kim, C. H., 2004b. A comparison between unified nodal method and semi-analytic nodal method for two-node multi-group neutron diffusion problem. In: Proc. Korean Nuclear Society 2004 Fall Mtg, Yongpyeong, Korea, October 28–29, 2004, CDROM (in Korean). Lefebvre, J.C., Mondot, J., West, J.P., 1991. Benchmark Calculations of Power Distribution within Assemblies. NEACRP-L-336, OECD Nuclear Energy Agency. Lewis, E.E. et al., 2001. Benchmark Specification for Deterministic 2-D/3-D MOX fuel Assembly Transport Calculations without Spatial Homogenization (C5G7MOX). OECD Nuclear Energy Agency. Muller, E.Z., Weiss, Z.J., 1991. Benchmarking with the multi-group diffusion highorder response matrix method. Ann. Nucl. Eng. 18, 535. Noh, J.M., Cho, N.Z., 1994. A new approach of analytic basis function expansion to neutron diffusion nodal calculation. Nucl. Sci. Eng. 116, 165. Noh, J.M., Pogosbekyan, L., Kim, Y.J., Joo, H.K., 1997. A general approach to multigroup extension of the analytic function expansion nodal method based on matrix function theory. In: Proc. Joint Intl. Conf. On Mathematical Methods and Supercomputing for Nuclear Applications, Saratoga Springs, 1, p. 144. Rajic, H.L., Ougouag, A.M., 1989. ILLICO: a nodal neutron diffusion method for modern computer architectures. Nucl. Sci. Eng. 103, 392–408. Shober, R.A., 1979. A nodal method for fast reactor analysis. Proceedings of the Computational Methods in Nuclear Engineering. ANS, Williamsburg. Smith, K.S., 1979. An analytic nodal method for solving the 2-Group, multidimensional, static and transient neutron diffusion equations. Nucl. Eng. Thesis, Department of Nucl. Eng., MIT. Smith, K.S., 1983. Nodal Method Storage Reduction by Nonlinear Iteration. Trans. Am. Nucl. Soc. 44, 265. Takeda, T., Ikeda, G., 1991. 3D neutron transport benchmarks. J. Nucl. Sci. Tech. 28, 656. Yoon, J.I., Joo, H.G., 2008. Two-level coarse mesh finite difference formulation with multigroup source expansion nodal kernels. J. Nucl. Sci. Technol. 45, 668–682. Zimin, V.G., Ninokata, H., 1998. Nodal neutron kinetics model based on nonlinear iteration procedure for LWR analysis. Ann. Nucl. Energ. 25, 507.