- Email: [email protected]
Analytic function expansion nodal (AFEN) method for solving multigroup neutron simplified P3 (SP3) equations in hexagonal-z geometry
Annals of Nuclear Energy 98 (2016) 74–80
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Analytic function expansion nodal (AFEN) method for solving multigroup neutron simplified P3 (SP3) equations in hexagonal-z geometry Mohammad Hasan Jalili Bahabadi, Ali Pazirandeh ⇑, Mitra Athari Department of Nuclear Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 21 March 2015 Received in revised form 4 May 2016 Accepted 28 June 2016 Available online 3 August 2016 Keywords: Simplified P3 SP3 Nodal method Hexagonal AFEN MGHANSP3 code
a b s t r a c t In this paper, analytic function expansion nodal (AFEN) method was developed to solve multi-group and multi-dimensional neutron simplified P3 (SP3) equation in reactor cores with hexagonal fuel assembly. In this method, the intranodal fluxes are expanded into a set of analytic basis functions for each group and Legendre moment. Fourteen boundary conditions has been considered that constrain the intranodal flux distributions in the hexagonal-z node, which include twelve radial surface-averaged partial currents and two axial surface-averaged partial currents. The code takes few-groups cross sections produced by a lattice code and calculates the effective multiplication factor (keff), flux in multi-group energy, reactivity, and the relative power density at each fuel assembly. Finally, the solution accuracy is tested for the two and three dimensional IAEA benchmark problem. The 3D IAEA problem has been calculated for three different cases. The numerical results demonstrate that the SP3 AFEN method exhibits better accuracy compared to diffusion method especially when the control materials are inserted in the core. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The solution of linear transport problems is difficult because of the rich phase space; in general there are seven independent variables: three for space, two describing a direction on the unit sphere, one for speed or energy, and one for time (McClarren, 2011). Obtaining energy and time dependent solutions to 3D transport problems is still challenging, even on petascale computers (McClarren, 2011). The simplified PN (SPN) method was developed by Gelbard (1962) when computer resources were used to solve 3D diffusion problems. The SPN method is a simple extension of the one-dimensional PN equations to the multidimensional PN equations. Larsen et al. (1993, 1996) showed that the SPN approximation is an asymptotic correction to standard diffusion theory. This method greatly reduced complexity of the multidimensional spherical harmonics (PN) method. The simplified P3 (SP3) approximation to the multigroup neutron transport equation in arbitrary geometries was derived by (Brantley and Larsen, 2000) using a variational analysis. They showed that the SP3 approximation can eliminate much of the inaccuracy in the diffusion and SP2 calculations of mixed-oxide (MOX) fuel problems ⇑ Corresponding author. E-mail addresses: [email protected] (M.H. Jalili Bahabadi), [email protected], [email protected] (A. Pazirandeh), mitra_athari@ yahoo.com (M. Athari). http://dx.doi.org/10.1016/j.anucene.2016.06.034 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.
(Brantley and Larsen, 2000). Ryu and Joo (2013) solved SP3 equations by finite element method for general geometry applications. Their calculation results demonstrated that the SP3 formulation is particularly beneficial for fast neutron systems. Because the SP3 equations have higher accuracy than diffusion theory and also less computational expense than the solution of transport equation with the discrete ordinate (SN) or spherical harmonics (PN) approximations, they were utilized in several pin-by-pin whole core calculation codes such as PARCS (Downar et al., 2004, 2009), SCOPE2 (Sugimura et al., 2006), DYN3D (Beckert and Grundmann, 2008) and COCAGNE (Fliscounakis et al., 2012). For some calculations such as fuel burnup, fuel management or transient analysis, due to the need for many full core spatial solutions, direct few group finite difference solution is inefficient or even impractical (Stacey, 2007). For this reason the nodal solution methods for the SP3 equations were developed (Downar et al., 2004, 2009; Beckert and Grundmann, 2008; Kim et al., 2009; Cho and Cho, 2010, Jalili et al., 2015a). In development of the PARCS (Downar et al., 2004, 2009) code, the nodal expansion method (NEM) was utilized for the multigroup SP3 equation. In this method, the three-dimensional SP3 equations are changed to three type one-dimensional SP3 equations by transverse integration procedure. The one-dimensional flux moments were approximated in 4th-order polynomials and
75
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
the transverse leakage shape was approximated in a second-order polynomial for both 0th and 2nd flux moments. In development of the DYN3D (Beckert and Grundmann, 2008) code, the SP3 equations are reduced to two-dimensional equations in radial plane (x, y), and one-dimensional equations in axial zdirection by transverse integration over the node height and rectangular node area, respectively. Polynomial and exponential functions were used as intranodal flux expansions. In the work of (Kim et al., 2009), the conformal mapped nodal simplified P3 equations were derived and implemented for the two-dimensional neutronics analysis of fast reactor cores with hexagonal fuel assemblies. In the past twenty years, the analytic function expansion nodal (AFEN) method (Cho and Noh, 1994, 1995; Noh and Cho, 1996; Cho et al., 1997; Woo et al., 2001; Xia and Xie, 2006; Cho and Lee, 2006; Cho and Cho 2010; Jalili et al., 2015a,b) has been the most commonly used for analysis of neutronics core calculations. This method represents multidimensional intra nodal flux distribution in terms of analytic basis functions at any points in the node. In the work of Cho and Cho (2010), the AFEN method has been reformulated for SP3 equations and implemented in a version of the AFEN method code COREDAX for cubic geometry (Cho and Cho, 2010). In the work of Jalili et al. (2015a,b), the new AFEN method has been developed for SP3 equations in rectangular-z geometry. In this work, the nodes are coupled through the surface averaged partial currents at each nodal interface, in other words, six boundary conditions at each group and Legendre moments have been considered. In the present work, a multigroup AFEN code named MGHANSP3 is developed for solution of SP3 equations in hexagonal-z fuel lattice. The boundary conditions include twelve radial surface averaged partial currents across the half radial hexagonal surfaces and two axial surface averaged partial currents. Finally, the solution accuracy is examined for the two and three dimensional IAEA (Herbert, 2008) problem. The numerical results illustrated that the solution of SP3 equations by the AFEN method is an accurate method for predicting multiplication factor and power distribution in reactor core with hexagonal fuel assembly. 2. Description of the method
The relationship between partial currents, surface fluxes and net currents is as follows (Downar et al., 2009):
j1 ¼
1 1 3 U0s J1 /2s ; 4 2 16
r Þ þ ½A/ð r Þ r ½Dr/ð
r/2g þ 13 r/0g þ Rtrg /1g ¼ 0; r:/3g þ 25 r /1g þ Rtg /2g ¼ 0;
g0
v
X g0
ð1Þ 0
mRfg0 /g0
The zeroth and second Legendre moments of the angular flux in group g (/0g ; /2g ) are the scaler variables, and the first and third Legendre moments of the angular flux in group g (/1g ; /3g ) are the vector variables. The Eq. (1) in each group can be rewritten as follows (Downar et al., 2009):
"
D1 r2 þ Rr 23 Rr
2Rr D3 r2 þ Rrt
#
U0 /2
"
¼ S0
1 23
#
ð4Þ
½AG1 ½AG2
½AGG
and
½/g ¼ ½U0g ; /2g " g # D0 0 ½Dg ¼ 0 D2g " # 0 0 Rsg !g 2Rgs !g ½Agg0 ¼ 0 0 2 g !g R 43 Rgs !g " 3 gs # Rr 2Rrg ½Agg ¼ 23 Rrg 53 Rtg þ 43 Rrg h i g g ½F g T ¼ mRf 2mRf
ð5Þ
ð6Þ
1 ½A0 ¼ ½D1 ½A v ½FT keff
ð7Þ
For analytic solution of Eq. (1), the eigenvalues km and thus eigenvectors um of the matrix [A0 ] must be calculated. For general multigroup problems the eigenvalues (and thus eigenvectors) may be real or complex. In this project, the QR algorithm is used to find the eigenvalue evaluation. The QR decomposition is based on Householder transformations. If all eigenvalues km of the matrix [A0 ] are real, to reduce Eq. (6) to the following decoupled form:
r2 Wm km Wm ¼ 0;
m ¼ 1; 2; . . . ; 2 G
ð8Þ
we can use the following transformation
;
U0 ¼ /0 þ 2/2 ; Rrt ¼ 53 Rt þ 43 Rr ; J1 ¼ D1 rU0 ; J3 ¼ D3 r/2
Rr ; J1 ¼ D1 rU0 ; J3 ¼ D3 r/2
r Þ ¼ 0 ½F T /ð
We can rewrite Eq. (4) in the form of the following equation:
r:/1g þ Rrg /0g ¼ S0g ;
Rsg0 /0g þ keffg
keff
½FT ¼ ½½F 1 T ½F 2 T . . . ½F G T
where
X
v
r Þ ¼ f½/ ; ½/ ; . . . ; ½/ gT /ð 2 3G 2 1 ½D1 0 6 . .. .. 7 7 ½D ¼ 6 . . 5 4 .. 0 ½DG 3 2 ½A11 ½A12 ½A1G 6 .. 7 6 ½A . . . . 7 7 6 21 ½A ¼ 6 7 .. .. 7 6 .. 4 . . . 5
The governing equations of the multidimensional SP3 equation is as follows (Brantley and Larsen, 2000):
S0g ¼
ð3Þ
where
rÞ þ ½A0 /ð r Þ ¼ 0 r2 /ð
r/2g þ Rtg /3g ¼ 0;
7 1 1 / J U0s 16 0s 2 3 16
Eq. (2) can be rewritten in the form of matrix equation as follows (Cho and Cho, 2011):
2.1. Intranodal flux expansion
2 3 3 5 3 7
j3 ¼
ð2Þ
U ¼ ½u1 ; . . . ; um ; . . . ; u2G ; um ¼ ½u1m ; . . . ; u2G m
WðrÞ ¼ U 1 UðrÞ where WðrÞ denote auxiliary flux.
T
ð9Þ
76
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
The elements of um vector can be obtained as: 2 3 a1ð2GÞ a11 6 7 . .. If ½A00 ¼ ½A km I ¼ 4 ... 5 (I is the identity . .. að2GÞ1 að2GÞð2GÞ matrix) then
2
.. .
a11 6. 6. 4. aðð2GÞ1Þ1
32
3
2
3
a1ð2GÞ u1m 76 . 7 6 7 .. 76 . 7 ¼ 6 7; . 54 . 5 4 5 2G1 aðð2GÞ1Þðð2GÞ1Þ aðð2GÞ1Þð2GÞ um a1ðð2GÞ1Þ .. .
u2G m ¼ 1; m ¼ 1; . . . ; 2 G
where uim is the element of the matrix U that is defined in Eq. (10). If some eigenvalues km are complex numbers, the corresponding eigenvectors will be complex and cannot be handled effectively. For this case, we used the following transformation
U ¼ ½Reðu1 Þ; Imðu1 Þ; . . . ; um ; . . . ; uG ; m ¼ j þ 1:
and the following auxiliary fluxes that are presented by Cho et al., 1997:
W1 ðrÞ ¼
Wm ðrÞ ¼
1 h X pffiffiffiffi pffiffiffiffi i km el r km el r A1 þ B1 ml e ml e
l¼1
el ¼ cos½ð2l 3Þhex þ sin½ð2l 3Þhey ; h ¼ tan1
pffiffi
ð12Þ
e7 ¼ ez
l¼1
l¼1
SN ¼ CS ¼
l¼1
ð13Þ
km < 0
The zero eigenvalue of matrix [A] in a node has not been considered in our calculations. It is noteworthy that the eigenvalues are rarely equal to zero. By using Eqs. (9) and (13) the intranodal flux distribution Ung ðrÞ can be obtained:
( ) 2G 7 7 X X X g i Un ðrÞ ¼ um Aml SNðkm el rÞ þ Bml CSðkm el rÞ ; m¼1
l¼1
U
g n ðrÞ
( ) j 2G 7 7 X X X X i i ¼ um Wm ðrÞ þ um Aml SNðkm el rÞ þ Bml CSðkm el rÞ ; m¼J
m¼1
l¼1
l¼1
ð17Þ
To calculate the volumetric averaged flux, the volumetric averaged auxiliary flux at the volume of node must be found.
ð18Þ
3x2H 3 3 3
where H and h are shown in Fig. 1. When km is the real number, the Eq. (13) is inserted in Eq. (18) and when km is the complex number, the Eq. (16) is inserted in Eq. (18). Therefore the volumetric averaged flux can be obtained as: ¼g
Un ¼
j 2G X X ¼ ¼ ui uim fWm g; m fWm g þ m¼1
m¼j
j ¼ Number of complex eigen v alues; n ¼ 0; 2ðLegendre momentsÞ; g ¼ group number and i ¼ n þ g; ð19Þ 2.2. Nodal boundary conditions
l¼1
i ¼ n þ g And n ¼ 0; 2;
l¼1
In this case, the intermodal flux distribution Ung ðrÞ can be written as follows:
0
sin; km < 0 cosh; km > 0 cos;
l¼1
Z h Z 0 Z pffi3xþ2Hpffi3 3 3 1 Wm ¼ pffiffiffi 2 pffi pffi Wm ðx; y; zÞdydxdz 2 3H h h H 33x2H3 3 ! Z H Z pffi3xþ2Hpffi3 3 3 Wm ðx; y; zÞdydxdz þ pffi pffi
7 7 X X Aml SNðkm el rÞ þ Bml CSðkm el rÞ
pffiffiffiffiffiffiffiffi jkm j; sinh; km > 0
l¼1
7 7 X X B2l sinhðpel rÞcosðqel rÞ A2l coshðpel rÞsinðqel rÞ þ
¼
So that Wm ðrÞ can be written in the following form:
km ¼
ð16Þ
i ¼ n þ g And n ¼ 0; 2;
3 6
l ¼ 1; 2; . . . ; 6;
l¼1
7 7 X X W2 ðrÞ ¼ A1l sinhðpel rÞsinðqel rÞ þ B1l coshðpel rÞcosðqel rÞ
ð11Þ
where el is an arbitrary unit vector. In practical calculations, the number of terms in Eq. (11) mainly depends on the number of the available nodal boundary conditions. Increasing the number of the boundary conditions will result in higher accuracy, but the computing time will also increase. Considering the coordinate system in Fig. 1, the following seven unit vectors are chosen for Eq. (11) (Jalili et al., 2015b).
l¼1
7 7 X X B2l coshðpel rÞsinðqel rÞ þ A2l sinhðpel rÞcosðqel rÞ þ
l¼1
Wm ðrÞ ¼
7 7 X X A1l coshðpel rÞcosðqel rÞ þ B1l sinhðpel rÞsinðqel rÞ l¼1
ð10Þ The analytic solution to Eq. (8) can be found in the following form:
ð15Þ
ð14Þ In this project, we consider 14 boundary conditions in each group and Legendre moment that constrain the intranodal flux distributions in the node shown in Fig. 1, which include twelve radial surface averaged partial currents and two axial surface averaged partial currents. g
g
The radial surface averaged partial currents j1;Sr and j3;Sr at given k
k
surface Srk are defined as:
1g 1 g 3 g 7 g 1 1 g g g / r J g r U J r / j1;Sr ¼ U r ¼ r r; j r k 4 0;Sk 2 1;Sk 16 2;Sk 3;Sk 16 2;Sk 2 3;Sk 16 0;Sk
ð20Þ
g r and / g r represent the scalar neutron fluxes and J g r and where U 0;S 2;S 1;S k
k
k
g
j3;Sr represent the net current averaged over the surface Srk , respeck g rþ and J g rþ ðrÞ can be found tively. For example at surfacesSrþ , U 2
Fig. 1. Coordinate system for one node.
in the following forms:
n;S2
n;S2
77
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
Fig. 2. IAEA benchmark configuration.
g rþ U n;S
1 ¼ 2Hh
2
J g rþ ðrÞ n;S2
R h R H2
/ng ðx; y; zÞjy¼p3ffi x2Hp3ffi dxdz
h
0
m¼1
R h R H2 @/ng ðx;y;zÞ
3 3 ( ) j G X X R R H2 R h R H i h um h 0 Wm ðx; y; zÞðFrom Eq: ð17ÞÞjy¼p3ffi x2Hp3ffi dxdz þ uim h 0 2 Wm ðx; y; zÞðFrom Eq: ð13ÞÞjy¼p3ffi x2Hp3ffi dxdz ;
¼
1 2Hh
¼
1 Dng 2Hh
3
h
( 1 ¼ Dng 2Hh
@x
0
þ
g
@/n ðx;y;zÞ @y
p
y¼
ffi
3
G X
uim
ffi dxdz
j X R h R H2 @ Wm ðx;y;zÞðFrom Eq: ð16ÞÞ @Wm ðx;y;zÞðFrom Eq: ð17ÞÞ
ui þ
m h 0 @x @y
R h R H2 @Wm ðx;y;zÞðFrom Eq: ð13ÞÞ h
@x
0
3
p 3x2H 3 3 3
m¼1
þ
3
m¼jþ1
Eq: ð13ÞÞ
þ @ Wm ðx;y;zÞðFrom
@y
p
y¼
m¼jþ1
ffi3
3
ffi
ffi dxdz
p p 3 x2H3 3 3
y¼
p
x2H3
)
ffi3 dxdz ; i ¼ n þ g: ð21Þ
It is obvious that all the partial currents are functions of the intranodal flux expansion coefficients Aml and Bml
where j is number of complex eigenvalues. g
In the axial direction, the surface averaged partial currents j1;S z
1
g
z and j3;S z at the surfaces S1 are defined as:
2.3. Iterative solution of the SP3 equation
1
1g 1 3 g U z J g z / z ; 4 0;S1 2 1;S1 16 2;S1 7 1 1 g / z J g z ¼ U z 16 2;S1 2 3;S1 16 0;S1
g
j1;S z ¼ 1
g
j3;S z
1
ð22Þ
Two sets of surface averaged incoming and outgoing partial currents correlations (21) and (22) in all groups and all surfaces of a node give the following two matrix equations:
J ¼ R X;
ð24:1Þ
g z and J g z ðrÞ can be found in the following forms: where U n;S n;S 1
g z U n;S
1
1
¼ 2p1ffiffi3H2 ¼ 2p1ffiffi3H2 þ
pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 g R H R y¼ 33xþ2H3 3 g pffi pffi / ðx; y; zÞdydxj pffi pffi þ / ðx; y; zÞdydxj z¼h z¼h n n 3 2H 3 3 2H 3 H 0
( j X
G X
1
3
x
y¼
3
3
x
3
pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 R H R y¼ 33xþ2H3 3 pffi pffi W ðx; y; zÞðFrom Eq: ð16ÞÞdxdyj pffi pffi ui þ W ðx; y; zÞðFrom Eq: ð16ÞÞdxdyj m m z¼h z¼h m 3 2H 3 3 2H 3 H 0 y¼
m¼1
3
x
y¼
3
3
x
3
) pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 R H R y¼ 33xþ2H3 3 i ffi ffi ffi ffi p p p p um H Wm ðx; y; zÞðFrom Eq: ð13ÞÞdxdyjz¼h þ 0 Wm ðx; y; zÞðFrom Eq: ð13ÞÞdxdyjz¼h 3 2H 3 3 2H 3 y¼
m¼jþ1
J g z n;S
y¼
3
x
y¼
3
3
x
3
pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 @/ng ðx;y;zÞ R H R y¼ 33xþ2H3 3 @/ng ðx;y;zÞ pffi pffi pffi pffi ¼ 2p1ffiffi3H2 H dydxj þ dydxj z¼h z¼h 0 @z @z y¼ 33x2H3 3 y¼ 33x2H3 3 ( pffi pffi pffi pffi j X R 0 R y¼ 33xþ2H3 3 @ Wm ðx;y;zÞðFrom Eq: ð16ÞÞ R H R y¼ 33xþ2H3 3 @ Wm ðx;y;zÞðFrom Eq: ð16ÞÞ pffi pffi pffi pffi ¼ p1ffiffi 2 ui dxdyj þ dxdyj z¼h z¼h m 3 2H 3 3 2H 3 H 0 @z @z 2 3H
y¼
m¼1
3
x
pffi pffi G X R 0 R y¼ 33xþ2H3 3 pffi pffi uim H þ 3 2H 3 m¼jþ1
y¼
3
x
3
y¼
3
3
y¼
3
ffi3
x
pffi R H R y¼ 3 xþ2H3 3 @ Wm ðx;y;zÞðFrom Eq: ð13ÞÞ pffi pffi dxdyj þ z¼h 3 2H 3 0 @z p
x
3
3
@ Wm ðx;y;zÞðFrom Eq:ð13ÞÞ dxdyjz¼h @z
) ; i ¼ n þ g:
ð23Þ
78
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
Table 1 Two-group cross sections for the IAEA problem. Region
Group
D (cm)
Ra (cm1)
tRf (cm1)
Rs;1!g (cm1)
Rs;2!g (cm1)
Fuel 1
1 2
1.5 0.4
0.01 0.08
0.0 0.135
0.1922222 0.02
– 0.7533333
Fuel 2
1 2
1.5 0.4
0.01 0.085
0.0 0.135
0.1922222 0.02
– 0.7483333
Fuel 2 + Rod
1 2
1.5 0.4
0.01 0.13
0.0 0.135
0.1922222 0.02
– 0.7033333
Reflector
1 2
1.5 0.4
0.0 0.01
0.0 0.0
0.1822222 0.04
– 0.8233333
Reflector + Fuel
1 2
1.5 0.4
0.0 0.055
0.0 0.0
0.1822222 0.04
– 0.7783333
Fig. 3. Comparison of assembly powers for the IAEA problem.
are obtained by all nodal sweeps from the initially given or previously updated incoming currents by using Eqs. (24.2) and (25). After that, the nodal average fluxes are calculated from Eq. (19). Finally, at the end of the inner iteration the effective multiplication factor (keff) is calculated. The above iteration is continued to satisfy the keff convergence criterion. 3. Numerical results
Fig. 4. Comparison of assembly powers for the IAEA problem.
J þ ¼ Rþ X;
ð24:2Þ
1 1 1 where J is the partial currents and J ¼ ½ðj1;1 , . . ., j1;14 ; j3;1 , . . ., 1
G
G
G
G
j3;14 ),. . ., ðj1;1 , . . ., j1;14 ; j3;1 , . . ., j3;14 )] and X is the flux expansion coefficients and X = [(A11, . . ., A17, B11, . . ., B17), . . ., (A2G1, . . ., A2G7, B2G1, . . ., B2G7)]T. Elements of both matrices R+ and R are constant numbers that depend on keff and the group constants of the node. In addition to Eq. (24), we require another correlation as follows:
J g; ¼ RJ g;þ s s
The analytic function expansion nodal (AFEN) method developed above has been implemented into a computer code written in C#, which solves the multi-groups neutron simplified P3 equation in hexagonal-z geometry. The accuracy and capability of this method have been tested against 2D and 3D IAEA (Herbert, 2008) benchmark problem. We will compare SP3 solutions obtained with the new AFEN method (MGHANSP3 code) to results that were obtained by MCNP4C code.
ð25Þ
where R is the core external boundary conditions. Given an initial value for multiplication factor, the matrix eigenvalues problem of the matrix A0 is first solved for eigenvalues km and corresponding eigenvectors um. X is obtained from Eq. (24.1) by Gauss elimination method. Then the outgoing partial currents
3.1. Benchmark problem (A): Two dimensional IAEA benchmark problem The 2D benchmark of IAEA (Herbert, 2008) models a core with hexagonal fuel assembly in steady state. The core geometry is shown in Fig. 2. In this figure, the radial fuel assembly lattice pitch is 20 cm. Also the cross sections are given in Table 1. In MGHANSP30 s calculations, each fuel assembly was considered as a node and a convergence criterion of 107 was used for keff. In Fig. 3 our AFEN SP3 (MGHANSP3) relative assembly core power density distribution is compared to McCARD code values that presented by Ryu and Joo, 2010. Only three relative power density errors are higher than 1% and the maximum error is 1.77%.
79
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80 Table 2 Numerical results of the IAEA problem. McCARDa
MCNP4C
DcCARTa
SP3 FDMa
SP3 FEMa
1.00622
1.00623 (+1)
1.00675 (+53)
1.00624 (+2)
Linear Quadratic
1.00633(+11) 1.00625(+3)
Run timeb (s)
Diffusion AFEN
SP3
1.00536(86)
1.00610(12)
AFEN
Diffusion AFEN 0.836
SP3AFEN 9.540
Deviation from McCARD code in PCM. a Ryu and Joo (2010). b Sony vaioVPCCW21FX.
Fig. 5. Axial shapes of three cases of the three dimensional IAEA problem.
Table 3 Results of 10 statistical check for the estimated answer for the tally fluctuation chart bin of tally 4. Tfc bin
Mean
Relative error
Behavior
Behavior
Value
Decrease
Decrease rate
Variance of the variance Value
Decrease
Decrease rate
Figure of merit Value
Behavior
Pdf Slope
Desired Observed Passed?
Random Random Yes
<0.10 0.00 Yes
Yes Yes Yes
1/sqrt (nps) Yes Yes
<0.10 0.00 Yes
Yes Yes Yes
1/nps Yes Yes
Constant Constant Yes
Random Random Yes
>3.00 10.00 Yes
Table 4 Numerical results of the three dimensional IAEA benchmark problems. MCNP4C (Reference)
Case 1 Case 2 Case 3
1.04067 1.00253 1.01774
ea (eMax )
keff (ek )
Run timea (s)
SP3
Diffusion
SP3
Diffusion
SP3
Diffusion
1.04064 (3) 1.00244 (9) 1.01769 (5)
1.04052 (15) 1.00167 (86) 1.01719 (55)
0.13 (0.43) 0.33 (0.81) 0.40 (1.24)
0.29 (0.67) 0.45 (1.67) 0.91 (2.57)
563.980 555.991 604.971
31.214 31.897 29.301
ek is the keff relative error in PCM. ea and eMax are average and maximum percentage relative errors of assembly powers. a
Sony vaioVPCCW21FX.
In Fig. 4, the diffusion and SP3 relative assembly core power density distribution are compared to MCNP4C code values. From this figure, we calculated the average percentage relative errors for SP3 and diffusion approximations. The average percentage relative errors for SP3 and diffusion equations are 0.51 and 0.61 respectively. Therefore, it is obvious that the SP3 power density results is more accurate than the diffusion theory. Table 2 gives the keff of several methods for the IAEA problem. It is seen that the MGANMSP3 keff deviation from McCARD code is only -12 PCM. 3.2. Benchmark problem (B): Three dimensional IAEA benchmark problem Second benchmark problem is based on the same core for problem (A) in which the original two dimensional core is generalized to a three dimensional core. The total core height is 380.0 cm including 20.0 cm thick axial reflectors. Vacuum boundary conditions are applied to the entire reflector external boundaries. In this problem, the following three cases are considered:
Case 1: All control rods are withdrawn. Case 2: All control rods are fully inserted. Case 3: Seven control rods are fully inserted, and 100.0 cm of the six control rods are inserted in the core. The axial shapes of the above cases are shown in Fig. 5. The reference results were obtained by MCNP4C code (MonteCarlo method). In all cases, the results of 10 statistical check for the estimated answer for the tally fluctuation chart bin of tally 4 are like Table 3. This table shows that the MCNP4C results are trustable. MGHANSP3 has been ran for the three-dimensional IAEA problem using coarse nodes i.e. one node per a fuel assembly in radial shape and axial sections with 20.0 cm height. Also the keff convergence criterion is 107 . The numerical results presented in Table 4 demonstrate that the SP3 method exhibits better accuracy compared to diffusion method especially in cases 2 and 3. Radial power distributions obtained from the SP3 and diffusion methods are compared with those from MCNP4C code in Figs. 6–8. These figures demonstrate that the SP3 solution is more accurate
80
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
Acknowledgment This work was supported by the Islamic Azad University Science and Research Branch. References
Fig. 6. Comparison of assembly powers for the 3D IAEA problem in case 1.
Fig. 7. Comparison of assembly powers for the 3D IAEA problem in case 2.
Fig. 8. Comparison of assembly powers for the 3D IAEA problem in case 3.
than the diffusion solution in the cases 1, 2 and 3. This more precision of SP3 solution is tangible in cases 2 and 3. 4. Conclusions A new SP3 AFEN code (MGHASP3) has been successfully developed for reactor cores with the hexagonal fuel assembly. This method represents multidimensional intra nodal flux distribution in terms of analytic basis functions at any point in a node. The surface averaged partial currents in half of the surfaces of the hexagon adopted at the nodal boundary coupling conditions. The code was verified by two benchmark problems. The numerical results demonstrate that the SP3 AFEN method exhibits better accuracy compared to diffusion method especially when the control materials are inserted in the core. These results also confirm the validity of MGHANSP3 code.
Beckert, C., Grundmann, U., 2008. Development and verification of a nodal approach for solving the multigroup SP3 equations. Ann. Nucl. Energy 35, 75–86. Brantley, P.S., Larsen, E.W., 2000. The simplified P3 approximation. Nucl. Sci. Eng. 134, 1–21. Cho, B.H. Cho, N.Z., 2010. Extension of Analytic Function Expansion Nodal (AFEN) Method for Multigroup Simplified P3 (SP3) Equations, Spring Meeting of the KNS; Pyongchang (Korea, Republic of); 27–28 May 2010. Available from KNS, Daejeon (KR). Cho, B.H., Cho, N.Z., 2011. Analytic Function Expansion Nodal (AFEN) Method for Multigroup Simplified P3 (SP3) Equations (Master’s thesis). Department of Nuclear and Quantum Engineering KAIST. Cho, N.Z., Lee, J., 2006. Analytic function expansion nodal (AFEN) method in hexagonal-Z three-dimensional geometry for neutron diffusion calculation. J. Nucl. Sci. Technol. 43 (11), 1320–1326. Cho, N.Z., Noh, J.M., 1994. The AFEN method for hexagonal nodal calculation and reconstruction. Trans. Am. Nucl. Soc. 71, 466–468. Cho, N.Z., Noh, J.M., 1995. Analytic function expansion nodal method for hexagonal geometry. Nucl. Sci. Eng. 121, 245–253. Cho, N.Z., Kim, Y.H., Park, K.W., 1997. Extension of analytic function expansion nodal method to multigroup problems in hexagonal-z geometry. Nucl. Sci. Eng. 126, 35–47. Downar, T., Lee, D., Xu, Y., Kozlowski, 2004. User Manual for the PARCS v2.6 U.S. NRC Core Neutronics Simulator. School of Nuclear Engineering Purdue University, USA. Downar, T., Xu, Y., Seker, V., 2009. Theory Manual for the PARCS Kinetics Core Simulator Module. Department of Nuclear Engineering and Radiological Sciences University of Michigan, USA. Fliscounakis, M., Girardi, E., Courau, T., Couyras, D., 2012. Potential of pin-by-pin SPN calculations as industrial reference. In: Proceedings of ICAPP 2012, Chicago, USA, June 24–28, 2012, on CD-ROM. Gelbard, E.M., 1962. Application of the Simplified Spherical Harmonics Equation in Spherical Geometry, WAPD-TM-294. Bettis Atomic Power Laboratory. Herbert, A., 2008. A Raviart-Thomas-Schneider solution of diffusion equation in hexagonal geometry. Ann. Nucl. Energy 35, 363–376. Jalili, M.H., Pazirandeh, A., Athari, M., 2015a. New analytic function expansion nodal (AFEN) method for solving multigroup neutron simplified P3 equation. Ann. Nucl. Eng. 77, 148–160. Jalili, M.H., Pazirandeh, A., Athari, M., 2015b. Development of a new analytic function expansion nodal code, HexDANM, for solving the neutron diffusion equation in hexagonal-Z geometry. Kerntechnik 80, 529–536. Kim, Y.I. et al., 2009. A conformal mapped nodal SP3 method for hexagonal core analysis. Ann. Nucl. Energy 36, 498–504. Larsen, E.W., Morel, J.E., McGhee, J., 1993, Asymptotic derivation of the simplified Pn equations. In Proceedings of Joint International Conference on Mathematical Methods and Supercomputing in Nuclear Applications, Portland, Oregon. Larsen, E., Morel, J., McGhee, J., 1996. Asymptotic derivation of the multigroup P1 and simplified PN equations with anisotropic scattering. Nucl. Sci. Eng. 123 (3), 328–342. McClarren, R.G., 2011. Theoretical aspects of the simplified Pn equations. Transp. Theory Stat. Phys. 39, 73–109. Noh, J.M., Cho, N.Z., 1996. A multigroup diffusion nodal scheme in rectangular and hexagonal geometry: hybrid of AFEN and PEN methods. In: Proceedings of the International Conference on the Physics of Reactor (PHYSOR 96), vol. 1, Mito, Ibaraki, Japan, 16–20 September 1996. Atomic Energy Society of Japan. pp. A50– A59. Ryu, E.H., Joo, H.G., 2010. Development of a 2D simplified P3 FEM solver for arbitrary geometry application. Transaction of the Korean Nuclear Society Autumn Meeting Jeju, Korea, October 21–22, 2010. Ryu, E.H., Joo, H.G., 2013. Finite element method solution of the simplified P3 equations for general geometry applications. Ann. Nucl. Energy 56, 194–207. Stacey, W.M., 2007. Nuclear Reactor Physics. John Wiley & Sons Inc., New York, USA. Sugimura, N., Ushio, T., Yamamoto, A., Tatsumi, M., 2006. Calculation models of AEGIS/SCOPE2, a core calculation system of next generation. In: Proceedings of PHYSOR-2006, ANS Topical Meeting on Reactor Physics, Vancouver, BC, Canada, and September 10–14, 2006. Woo, S.W., Cho, N.Z., Noh, J.M., 2001. The analytic function expansion nodal method refined with transverse gradient basis functions and interface flux moments. Nucl. Sci. Eng. 139, 156. Xia, B., Xie, Z., 2006. Flux expansion nodal method for solving multigroup neutron diffusion equations in hexagonal-z geometry. Ann. Nucl. Energy 33, 370–376.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Analytic function expansion nodal (AFEN) method for solving multigroup neutron simplified P3 (SP3) equations in hexagonal-z geometry Mohammad Hasan Jalili Bahabadi, Ali Pazirandeh ⇑, Mitra Athari Department of Nuclear Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 21 March 2015 Received in revised form 4 May 2016 Accepted 28 June 2016 Available online 3 August 2016 Keywords: Simplified P3 SP3 Nodal method Hexagonal AFEN MGHANSP3 code
a b s t r a c t In this paper, analytic function expansion nodal (AFEN) method was developed to solve multi-group and multi-dimensional neutron simplified P3 (SP3) equation in reactor cores with hexagonal fuel assembly. In this method, the intranodal fluxes are expanded into a set of analytic basis functions for each group and Legendre moment. Fourteen boundary conditions has been considered that constrain the intranodal flux distributions in the hexagonal-z node, which include twelve radial surface-averaged partial currents and two axial surface-averaged partial currents. The code takes few-groups cross sections produced by a lattice code and calculates the effective multiplication factor (keff), flux in multi-group energy, reactivity, and the relative power density at each fuel assembly. Finally, the solution accuracy is tested for the two and three dimensional IAEA benchmark problem. The 3D IAEA problem has been calculated for three different cases. The numerical results demonstrate that the SP3 AFEN method exhibits better accuracy compared to diffusion method especially when the control materials are inserted in the core. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The solution of linear transport problems is difficult because of the rich phase space; in general there are seven independent variables: three for space, two describing a direction on the unit sphere, one for speed or energy, and one for time (McClarren, 2011). Obtaining energy and time dependent solutions to 3D transport problems is still challenging, even on petascale computers (McClarren, 2011). The simplified PN (SPN) method was developed by Gelbard (1962) when computer resources were used to solve 3D diffusion problems. The SPN method is a simple extension of the one-dimensional PN equations to the multidimensional PN equations. Larsen et al. (1993, 1996) showed that the SPN approximation is an asymptotic correction to standard diffusion theory. This method greatly reduced complexity of the multidimensional spherical harmonics (PN) method. The simplified P3 (SP3) approximation to the multigroup neutron transport equation in arbitrary geometries was derived by (Brantley and Larsen, 2000) using a variational analysis. They showed that the SP3 approximation can eliminate much of the inaccuracy in the diffusion and SP2 calculations of mixed-oxide (MOX) fuel problems ⇑ Corresponding author. E-mail addresses: [email protected] (M.H. Jalili Bahabadi), [email protected], [email protected] (A. Pazirandeh), mitra_athari@ yahoo.com (M. Athari). http://dx.doi.org/10.1016/j.anucene.2016.06.034 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.
(Brantley and Larsen, 2000). Ryu and Joo (2013) solved SP3 equations by finite element method for general geometry applications. Their calculation results demonstrated that the SP3 formulation is particularly beneficial for fast neutron systems. Because the SP3 equations have higher accuracy than diffusion theory and also less computational expense than the solution of transport equation with the discrete ordinate (SN) or spherical harmonics (PN) approximations, they were utilized in several pin-by-pin whole core calculation codes such as PARCS (Downar et al., 2004, 2009), SCOPE2 (Sugimura et al., 2006), DYN3D (Beckert and Grundmann, 2008) and COCAGNE (Fliscounakis et al., 2012). For some calculations such as fuel burnup, fuel management or transient analysis, due to the need for many full core spatial solutions, direct few group finite difference solution is inefficient or even impractical (Stacey, 2007). For this reason the nodal solution methods for the SP3 equations were developed (Downar et al., 2004, 2009; Beckert and Grundmann, 2008; Kim et al., 2009; Cho and Cho, 2010, Jalili et al., 2015a). In development of the PARCS (Downar et al., 2004, 2009) code, the nodal expansion method (NEM) was utilized for the multigroup SP3 equation. In this method, the three-dimensional SP3 equations are changed to three type one-dimensional SP3 equations by transverse integration procedure. The one-dimensional flux moments were approximated in 4th-order polynomials and
75
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
the transverse leakage shape was approximated in a second-order polynomial for both 0th and 2nd flux moments. In development of the DYN3D (Beckert and Grundmann, 2008) code, the SP3 equations are reduced to two-dimensional equations in radial plane (x, y), and one-dimensional equations in axial zdirection by transverse integration over the node height and rectangular node area, respectively. Polynomial and exponential functions were used as intranodal flux expansions. In the work of (Kim et al., 2009), the conformal mapped nodal simplified P3 equations were derived and implemented for the two-dimensional neutronics analysis of fast reactor cores with hexagonal fuel assemblies. In the past twenty years, the analytic function expansion nodal (AFEN) method (Cho and Noh, 1994, 1995; Noh and Cho, 1996; Cho et al., 1997; Woo et al., 2001; Xia and Xie, 2006; Cho and Lee, 2006; Cho and Cho 2010; Jalili et al., 2015a,b) has been the most commonly used for analysis of neutronics core calculations. This method represents multidimensional intra nodal flux distribution in terms of analytic basis functions at any points in the node. In the work of Cho and Cho (2010), the AFEN method has been reformulated for SP3 equations and implemented in a version of the AFEN method code COREDAX for cubic geometry (Cho and Cho, 2010). In the work of Jalili et al. (2015a,b), the new AFEN method has been developed for SP3 equations in rectangular-z geometry. In this work, the nodes are coupled through the surface averaged partial currents at each nodal interface, in other words, six boundary conditions at each group and Legendre moments have been considered. In the present work, a multigroup AFEN code named MGHANSP3 is developed for solution of SP3 equations in hexagonal-z fuel lattice. The boundary conditions include twelve radial surface averaged partial currents across the half radial hexagonal surfaces and two axial surface averaged partial currents. Finally, the solution accuracy is examined for the two and three dimensional IAEA (Herbert, 2008) problem. The numerical results illustrated that the solution of SP3 equations by the AFEN method is an accurate method for predicting multiplication factor and power distribution in reactor core with hexagonal fuel assembly. 2. Description of the method
The relationship between partial currents, surface fluxes and net currents is as follows (Downar et al., 2009):
j1 ¼
1 1 3 U0s J1 /2s ; 4 2 16
r Þ þ ½A/ð r Þ r ½Dr/ð
r/2g þ 13 r/0g þ Rtrg /1g ¼ 0; r:/3g þ 25 r /1g þ Rtg /2g ¼ 0;
g0
v
X g0
ð1Þ 0
mRfg0 /g0
The zeroth and second Legendre moments of the angular flux in group g (/0g ; /2g ) are the scaler variables, and the first and third Legendre moments of the angular flux in group g (/1g ; /3g ) are the vector variables. The Eq. (1) in each group can be rewritten as follows (Downar et al., 2009):
"
D1 r2 þ Rr 23 Rr
2Rr D3 r2 þ Rrt
#
U0 /2
"
¼ S0
1 23
#
ð4Þ
½AG1 ½AG2
½AGG
and
½/g ¼ ½U0g ; /2g " g # D0 0 ½Dg ¼ 0 D2g " # 0 0 Rsg !g 2Rgs !g ½Agg0 ¼ 0 0 2 g !g R 43 Rgs !g " 3 gs # Rr 2Rrg ½Agg ¼ 23 Rrg 53 Rtg þ 43 Rrg h i g g ½F g T ¼ mRf 2mRf
ð5Þ
ð6Þ
1 ½A0 ¼ ½D1 ½A v ½FT keff
ð7Þ
For analytic solution of Eq. (1), the eigenvalues km and thus eigenvectors um of the matrix [A0 ] must be calculated. For general multigroup problems the eigenvalues (and thus eigenvectors) may be real or complex. In this project, the QR algorithm is used to find the eigenvalue evaluation. The QR decomposition is based on Householder transformations. If all eigenvalues km of the matrix [A0 ] are real, to reduce Eq. (6) to the following decoupled form:
r2 Wm km Wm ¼ 0;
m ¼ 1; 2; . . . ; 2 G
ð8Þ
we can use the following transformation
;
U0 ¼ /0 þ 2/2 ; Rrt ¼ 53 Rt þ 43 Rr ; J1 ¼ D1 rU0 ; J3 ¼ D3 r/2
Rr ; J1 ¼ D1 rU0 ; J3 ¼ D3 r/2
r Þ ¼ 0 ½F T /ð
We can rewrite Eq. (4) in the form of the following equation:
r:/1g þ Rrg /0g ¼ S0g ;
Rsg0 /0g þ keffg
keff
½FT ¼ ½½F 1 T ½F 2 T . . . ½F G T
where
X
v
r Þ ¼ f½/ ; ½/ ; . . . ; ½/ gT /ð 2 3G 2 1 ½D1 0 6 . .. .. 7 7 ½D ¼ 6 . . 5 4 .. 0 ½DG 3 2 ½A11 ½A12 ½A1G 6 .. 7 6 ½A . . . . 7 7 6 21 ½A ¼ 6 7 .. .. 7 6 .. 4 . . . 5
The governing equations of the multidimensional SP3 equation is as follows (Brantley and Larsen, 2000):
S0g ¼
ð3Þ
where
rÞ þ ½A0 /ð r Þ ¼ 0 r2 /ð
r/2g þ Rtg /3g ¼ 0;
7 1 1 / J U0s 16 0s 2 3 16
Eq. (2) can be rewritten in the form of matrix equation as follows (Cho and Cho, 2011):
2.1. Intranodal flux expansion
2 3 3 5 3 7
j3 ¼
ð2Þ
U ¼ ½u1 ; . . . ; um ; . . . ; u2G ; um ¼ ½u1m ; . . . ; u2G m
WðrÞ ¼ U 1 UðrÞ where WðrÞ denote auxiliary flux.
T
ð9Þ
76
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
The elements of um vector can be obtained as: 2 3 a1ð2GÞ a11 6 7 . .. If ½A00 ¼ ½A km I ¼ 4 ... 5 (I is the identity . .. að2GÞ1 að2GÞð2GÞ matrix) then
2
.. .
a11 6. 6. 4. aðð2GÞ1Þ1
32
3
2
3
a1ð2GÞ u1m 76 . 7 6 7 .. 76 . 7 ¼ 6 7; . 54 . 5 4 5 2G1 aðð2GÞ1Þðð2GÞ1Þ aðð2GÞ1Þð2GÞ um a1ðð2GÞ1Þ .. .
u2G m ¼ 1; m ¼ 1; . . . ; 2 G
where uim is the element of the matrix U that is defined in Eq. (10). If some eigenvalues km are complex numbers, the corresponding eigenvectors will be complex and cannot be handled effectively. For this case, we used the following transformation
U ¼ ½Reðu1 Þ; Imðu1 Þ; . . . ; um ; . . . ; uG ; m ¼ j þ 1:
and the following auxiliary fluxes that are presented by Cho et al., 1997:
W1 ðrÞ ¼
Wm ðrÞ ¼
1 h X pffiffiffiffi pffiffiffiffi i km el r km el r A1 þ B1 ml e ml e
l¼1
el ¼ cos½ð2l 3Þhex þ sin½ð2l 3Þhey ; h ¼ tan1
pffiffi
ð12Þ
e7 ¼ ez
l¼1
l¼1
SN ¼ CS ¼
l¼1
ð13Þ
km < 0
The zero eigenvalue of matrix [A] in a node has not been considered in our calculations. It is noteworthy that the eigenvalues are rarely equal to zero. By using Eqs. (9) and (13) the intranodal flux distribution Ung ðrÞ can be obtained:
( ) 2G 7 7 X X X g i Un ðrÞ ¼ um Aml SNðkm el rÞ þ Bml CSðkm el rÞ ; m¼1
l¼1
U
g n ðrÞ
( ) j 2G 7 7 X X X X i i ¼ um Wm ðrÞ þ um Aml SNðkm el rÞ þ Bml CSðkm el rÞ ; m¼J
m¼1
l¼1
l¼1
ð17Þ
To calculate the volumetric averaged flux, the volumetric averaged auxiliary flux at the volume of node must be found.
ð18Þ
3x2H 3 3 3
where H and h are shown in Fig. 1. When km is the real number, the Eq. (13) is inserted in Eq. (18) and when km is the complex number, the Eq. (16) is inserted in Eq. (18). Therefore the volumetric averaged flux can be obtained as: ¼g
Un ¼
j 2G X X ¼ ¼ ui uim fWm g; m fWm g þ m¼1
m¼j
j ¼ Number of complex eigen v alues; n ¼ 0; 2ðLegendre momentsÞ; g ¼ group number and i ¼ n þ g; ð19Þ 2.2. Nodal boundary conditions
l¼1
i ¼ n þ g And n ¼ 0; 2;
l¼1
In this case, the intermodal flux distribution Ung ðrÞ can be written as follows:
0
sin; km < 0 cosh; km > 0 cos;
l¼1
Z h Z 0 Z pffi3xþ2Hpffi3 3 3 1 Wm ¼ pffiffiffi 2 pffi pffi Wm ðx; y; zÞdydxdz 2 3H h h H 33x2H3 3 ! Z H Z pffi3xþ2Hpffi3 3 3 Wm ðx; y; zÞdydxdz þ pffi pffi
7 7 X X Aml SNðkm el rÞ þ Bml CSðkm el rÞ
pffiffiffiffiffiffiffiffi jkm j; sinh; km > 0
l¼1
7 7 X X B2l sinhðpel rÞcosðqel rÞ A2l coshðpel rÞsinðqel rÞ þ
¼
So that Wm ðrÞ can be written in the following form:
km ¼
ð16Þ
i ¼ n þ g And n ¼ 0; 2;
3 6
l ¼ 1; 2; . . . ; 6;
l¼1
7 7 X X W2 ðrÞ ¼ A1l sinhðpel rÞsinðqel rÞ þ B1l coshðpel rÞcosðqel rÞ
ð11Þ
where el is an arbitrary unit vector. In practical calculations, the number of terms in Eq. (11) mainly depends on the number of the available nodal boundary conditions. Increasing the number of the boundary conditions will result in higher accuracy, but the computing time will also increase. Considering the coordinate system in Fig. 1, the following seven unit vectors are chosen for Eq. (11) (Jalili et al., 2015b).
l¼1
7 7 X X B2l coshðpel rÞsinðqel rÞ þ A2l sinhðpel rÞcosðqel rÞ þ
l¼1
Wm ðrÞ ¼
7 7 X X A1l coshðpel rÞcosðqel rÞ þ B1l sinhðpel rÞsinðqel rÞ l¼1
ð10Þ The analytic solution to Eq. (8) can be found in the following form:
ð15Þ
ð14Þ In this project, we consider 14 boundary conditions in each group and Legendre moment that constrain the intranodal flux distributions in the node shown in Fig. 1, which include twelve radial surface averaged partial currents and two axial surface averaged partial currents. g
g
The radial surface averaged partial currents j1;Sr and j3;Sr at given k
k
surface Srk are defined as:
1g 1 g 3 g 7 g 1 1 g g g / r J g r U J r / j1;Sr ¼ U r ¼ r r; j r k 4 0;Sk 2 1;Sk 16 2;Sk 3;Sk 16 2;Sk 2 3;Sk 16 0;Sk
ð20Þ
g r and / g r represent the scalar neutron fluxes and J g r and where U 0;S 2;S 1;S k
k
k
g
j3;Sr represent the net current averaged over the surface Srk , respeck g rþ and J g rþ ðrÞ can be found tively. For example at surfacesSrþ , U 2
Fig. 1. Coordinate system for one node.
in the following forms:
n;S2
n;S2
77
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
Fig. 2. IAEA benchmark configuration.
g rþ U n;S
1 ¼ 2Hh
2
J g rþ ðrÞ n;S2
R h R H2
/ng ðx; y; zÞjy¼p3ffi x2Hp3ffi dxdz
h
0
m¼1
R h R H2 @/ng ðx;y;zÞ
3 3 ( ) j G X X R R H2 R h R H i h um h 0 Wm ðx; y; zÞðFrom Eq: ð17ÞÞjy¼p3ffi x2Hp3ffi dxdz þ uim h 0 2 Wm ðx; y; zÞðFrom Eq: ð13ÞÞjy¼p3ffi x2Hp3ffi dxdz ;
¼
1 2Hh
¼
1 Dng 2Hh
3
h
( 1 ¼ Dng 2Hh
@x
0
þ
g
@/n ðx;y;zÞ @y
p
y¼
ffi
3
G X
uim
ffi dxdz
j X R h R H2 @ Wm ðx;y;zÞðFrom Eq: ð16ÞÞ @Wm ðx;y;zÞðFrom Eq: ð17ÞÞ
ui þ
m h 0 @x @y
R h R H2 @Wm ðx;y;zÞðFrom Eq: ð13ÞÞ h
@x
0
3
p 3x2H 3 3 3
m¼1
þ
3
m¼jþ1
Eq: ð13ÞÞ
þ @ Wm ðx;y;zÞðFrom
@y
p
y¼
m¼jþ1
ffi3
3
ffi
ffi dxdz
p p 3 x2H3 3 3
y¼
p
x2H3
)
ffi3 dxdz ; i ¼ n þ g: ð21Þ
It is obvious that all the partial currents are functions of the intranodal flux expansion coefficients Aml and Bml
where j is number of complex eigenvalues. g
In the axial direction, the surface averaged partial currents j1;S z
1
g
z and j3;S z at the surfaces S1 are defined as:
2.3. Iterative solution of the SP3 equation
1
1g 1 3 g U z J g z / z ; 4 0;S1 2 1;S1 16 2;S1 7 1 1 g / z J g z ¼ U z 16 2;S1 2 3;S1 16 0;S1
g
j1;S z ¼ 1
g
j3;S z
1
ð22Þ
Two sets of surface averaged incoming and outgoing partial currents correlations (21) and (22) in all groups and all surfaces of a node give the following two matrix equations:
J ¼ R X;
ð24:1Þ
g z and J g z ðrÞ can be found in the following forms: where U n;S n;S 1
g z U n;S
1
1
¼ 2p1ffiffi3H2 ¼ 2p1ffiffi3H2 þ
pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 g R H R y¼ 33xþ2H3 3 g pffi pffi / ðx; y; zÞdydxj pffi pffi þ / ðx; y; zÞdydxj z¼h z¼h n n 3 2H 3 3 2H 3 H 0
( j X
G X
1
3
x
y¼
3
3
x
3
pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 R H R y¼ 33xþ2H3 3 pffi pffi W ðx; y; zÞðFrom Eq: ð16ÞÞdxdyj pffi pffi ui þ W ðx; y; zÞðFrom Eq: ð16ÞÞdxdyj m m z¼h z¼h m 3 2H 3 3 2H 3 H 0 y¼
m¼1
3
x
y¼
3
3
x
3
) pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 R H R y¼ 33xþ2H3 3 i ffi ffi ffi ffi p p p p um H Wm ðx; y; zÞðFrom Eq: ð13ÞÞdxdyjz¼h þ 0 Wm ðx; y; zÞðFrom Eq: ð13ÞÞdxdyjz¼h 3 2H 3 3 2H 3 y¼
m¼jþ1
J g z n;S
y¼
3
x
y¼
3
3
x
3
pffi pffi pffi pffi R 0 R y¼ 33xþ2H3 3 @/ng ðx;y;zÞ R H R y¼ 33xþ2H3 3 @/ng ðx;y;zÞ pffi pffi pffi pffi ¼ 2p1ffiffi3H2 H dydxj þ dydxj z¼h z¼h 0 @z @z y¼ 33x2H3 3 y¼ 33x2H3 3 ( pffi pffi pffi pffi j X R 0 R y¼ 33xþ2H3 3 @ Wm ðx;y;zÞðFrom Eq: ð16ÞÞ R H R y¼ 33xþ2H3 3 @ Wm ðx;y;zÞðFrom Eq: ð16ÞÞ pffi pffi pffi pffi ¼ p1ffiffi 2 ui dxdyj þ dxdyj z¼h z¼h m 3 2H 3 3 2H 3 H 0 @z @z 2 3H
y¼
m¼1
3
x
pffi pffi G X R 0 R y¼ 33xþ2H3 3 pffi pffi uim H þ 3 2H 3 m¼jþ1
y¼
3
x
3
y¼
3
3
y¼
3
ffi3
x
pffi R H R y¼ 3 xþ2H3 3 @ Wm ðx;y;zÞðFrom Eq: ð13ÞÞ pffi pffi dxdyj þ z¼h 3 2H 3 0 @z p
x
3
3
@ Wm ðx;y;zÞðFrom Eq:ð13ÞÞ dxdyjz¼h @z
) ; i ¼ n þ g:
ð23Þ
78
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
Table 1 Two-group cross sections for the IAEA problem. Region
Group
D (cm)
Ra (cm1)
tRf (cm1)
Rs;1!g (cm1)
Rs;2!g (cm1)
Fuel 1
1 2
1.5 0.4
0.01 0.08
0.0 0.135
0.1922222 0.02
– 0.7533333
Fuel 2
1 2
1.5 0.4
0.01 0.085
0.0 0.135
0.1922222 0.02
– 0.7483333
Fuel 2 + Rod
1 2
1.5 0.4
0.01 0.13
0.0 0.135
0.1922222 0.02
– 0.7033333
Reflector
1 2
1.5 0.4
0.0 0.01
0.0 0.0
0.1822222 0.04
– 0.8233333
Reflector + Fuel
1 2
1.5 0.4
0.0 0.055
0.0 0.0
0.1822222 0.04
– 0.7783333
Fig. 3. Comparison of assembly powers for the IAEA problem.
are obtained by all nodal sweeps from the initially given or previously updated incoming currents by using Eqs. (24.2) and (25). After that, the nodal average fluxes are calculated from Eq. (19). Finally, at the end of the inner iteration the effective multiplication factor (keff) is calculated. The above iteration is continued to satisfy the keff convergence criterion. 3. Numerical results
Fig. 4. Comparison of assembly powers for the IAEA problem.
J þ ¼ Rþ X;
ð24:2Þ
1 1 1 where J is the partial currents and J ¼ ½ðj1;1 , . . ., j1;14 ; j3;1 , . . ., 1
G
G
G
G
j3;14 ),. . ., ðj1;1 , . . ., j1;14 ; j3;1 , . . ., j3;14 )] and X is the flux expansion coefficients and X = [(A11, . . ., A17, B11, . . ., B17), . . ., (A2G1, . . ., A2G7, B2G1, . . ., B2G7)]T. Elements of both matrices R+ and R are constant numbers that depend on keff and the group constants of the node. In addition to Eq. (24), we require another correlation as follows:
J g; ¼ RJ g;þ s s
The analytic function expansion nodal (AFEN) method developed above has been implemented into a computer code written in C#, which solves the multi-groups neutron simplified P3 equation in hexagonal-z geometry. The accuracy and capability of this method have been tested against 2D and 3D IAEA (Herbert, 2008) benchmark problem. We will compare SP3 solutions obtained with the new AFEN method (MGHANSP3 code) to results that were obtained by MCNP4C code.
ð25Þ
where R is the core external boundary conditions. Given an initial value for multiplication factor, the matrix eigenvalues problem of the matrix A0 is first solved for eigenvalues km and corresponding eigenvectors um. X is obtained from Eq. (24.1) by Gauss elimination method. Then the outgoing partial currents
3.1. Benchmark problem (A): Two dimensional IAEA benchmark problem The 2D benchmark of IAEA (Herbert, 2008) models a core with hexagonal fuel assembly in steady state. The core geometry is shown in Fig. 2. In this figure, the radial fuel assembly lattice pitch is 20 cm. Also the cross sections are given in Table 1. In MGHANSP30 s calculations, each fuel assembly was considered as a node and a convergence criterion of 107 was used for keff. In Fig. 3 our AFEN SP3 (MGHANSP3) relative assembly core power density distribution is compared to McCARD code values that presented by Ryu and Joo, 2010. Only three relative power density errors are higher than 1% and the maximum error is 1.77%.
79
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80 Table 2 Numerical results of the IAEA problem. McCARDa
MCNP4C
DcCARTa
SP3 FDMa
SP3 FEMa
1.00622
1.00623 (+1)
1.00675 (+53)
1.00624 (+2)
Linear Quadratic
1.00633(+11) 1.00625(+3)
Run timeb (s)
Diffusion AFEN
SP3
1.00536(86)
1.00610(12)
AFEN
Diffusion AFEN 0.836
SP3AFEN 9.540
Deviation from McCARD code in PCM. a Ryu and Joo (2010). b Sony vaioVPCCW21FX.
Fig. 5. Axial shapes of three cases of the three dimensional IAEA problem.
Table 3 Results of 10 statistical check for the estimated answer for the tally fluctuation chart bin of tally 4. Tfc bin
Mean
Relative error
Behavior
Behavior
Value
Decrease
Decrease rate
Variance of the variance Value
Decrease
Decrease rate
Figure of merit Value
Behavior
Pdf Slope
Desired Observed Passed?
Random Random Yes
<0.10 0.00 Yes
Yes Yes Yes
1/sqrt (nps) Yes Yes
<0.10 0.00 Yes
Yes Yes Yes
1/nps Yes Yes
Constant Constant Yes
Random Random Yes
>3.00 10.00 Yes
Table 4 Numerical results of the three dimensional IAEA benchmark problems. MCNP4C (Reference)
Case 1 Case 2 Case 3
1.04067 1.00253 1.01774
ea (eMax )
keff (ek )
Run timea (s)
SP3
Diffusion
SP3
Diffusion
SP3
Diffusion
1.04064 (3) 1.00244 (9) 1.01769 (5)
1.04052 (15) 1.00167 (86) 1.01719 (55)
0.13 (0.43) 0.33 (0.81) 0.40 (1.24)
0.29 (0.67) 0.45 (1.67) 0.91 (2.57)
563.980 555.991 604.971
31.214 31.897 29.301
ek is the keff relative error in PCM. ea and eMax are average and maximum percentage relative errors of assembly powers. a
Sony vaioVPCCW21FX.
In Fig. 4, the diffusion and SP3 relative assembly core power density distribution are compared to MCNP4C code values. From this figure, we calculated the average percentage relative errors for SP3 and diffusion approximations. The average percentage relative errors for SP3 and diffusion equations are 0.51 and 0.61 respectively. Therefore, it is obvious that the SP3 power density results is more accurate than the diffusion theory. Table 2 gives the keff of several methods for the IAEA problem. It is seen that the MGANMSP3 keff deviation from McCARD code is only -12 PCM. 3.2. Benchmark problem (B): Three dimensional IAEA benchmark problem Second benchmark problem is based on the same core for problem (A) in which the original two dimensional core is generalized to a three dimensional core. The total core height is 380.0 cm including 20.0 cm thick axial reflectors. Vacuum boundary conditions are applied to the entire reflector external boundaries. In this problem, the following three cases are considered:
Case 1: All control rods are withdrawn. Case 2: All control rods are fully inserted. Case 3: Seven control rods are fully inserted, and 100.0 cm of the six control rods are inserted in the core. The axial shapes of the above cases are shown in Fig. 5. The reference results were obtained by MCNP4C code (MonteCarlo method). In all cases, the results of 10 statistical check for the estimated answer for the tally fluctuation chart bin of tally 4 are like Table 3. This table shows that the MCNP4C results are trustable. MGHANSP3 has been ran for the three-dimensional IAEA problem using coarse nodes i.e. one node per a fuel assembly in radial shape and axial sections with 20.0 cm height. Also the keff convergence criterion is 107 . The numerical results presented in Table 4 demonstrate that the SP3 method exhibits better accuracy compared to diffusion method especially in cases 2 and 3. Radial power distributions obtained from the SP3 and diffusion methods are compared with those from MCNP4C code in Figs. 6–8. These figures demonstrate that the SP3 solution is more accurate
80
M.H. Jalili Bahabadi et al. / Annals of Nuclear Energy 98 (2016) 74–80
Acknowledgment This work was supported by the Islamic Azad University Science and Research Branch. References
Fig. 6. Comparison of assembly powers for the 3D IAEA problem in case 1.
Fig. 7. Comparison of assembly powers for the 3D IAEA problem in case 2.
Fig. 8. Comparison of assembly powers for the 3D IAEA problem in case 3.
than the diffusion solution in the cases 1, 2 and 3. This more precision of SP3 solution is tangible in cases 2 and 3. 4. Conclusions A new SP3 AFEN code (MGHASP3) has been successfully developed for reactor cores with the hexagonal fuel assembly. This method represents multidimensional intra nodal flux distribution in terms of analytic basis functions at any point in a node. The surface averaged partial currents in half of the surfaces of the hexagon adopted at the nodal boundary coupling conditions. The code was verified by two benchmark problems. The numerical results demonstrate that the SP3 AFEN method exhibits better accuracy compared to diffusion method especially when the control materials are inserted in the core. These results also confirm the validity of MGHANSP3 code.
Beckert, C., Grundmann, U., 2008. Development and verification of a nodal approach for solving the multigroup SP3 equations. Ann. Nucl. Energy 35, 75–86. Brantley, P.S., Larsen, E.W., 2000. The simplified P3 approximation. Nucl. Sci. Eng. 134, 1–21. Cho, B.H. Cho, N.Z., 2010. Extension of Analytic Function Expansion Nodal (AFEN) Method for Multigroup Simplified P3 (SP3) Equations, Spring Meeting of the KNS; Pyongchang (Korea, Republic of); 27–28 May 2010. Available from KNS, Daejeon (KR). Cho, B.H., Cho, N.Z., 2011. Analytic Function Expansion Nodal (AFEN) Method for Multigroup Simplified P3 (SP3) Equations (Master’s thesis). Department of Nuclear and Quantum Engineering KAIST. Cho, N.Z., Lee, J., 2006. Analytic function expansion nodal (AFEN) method in hexagonal-Z three-dimensional geometry for neutron diffusion calculation. J. Nucl. Sci. Technol. 43 (11), 1320–1326. Cho, N.Z., Noh, J.M., 1994. The AFEN method for hexagonal nodal calculation and reconstruction. Trans. Am. Nucl. Soc. 71, 466–468. Cho, N.Z., Noh, J.M., 1995. Analytic function expansion nodal method for hexagonal geometry. Nucl. Sci. Eng. 121, 245–253. Cho, N.Z., Kim, Y.H., Park, K.W., 1997. Extension of analytic function expansion nodal method to multigroup problems in hexagonal-z geometry. Nucl. Sci. Eng. 126, 35–47. Downar, T., Lee, D., Xu, Y., Kozlowski, 2004. User Manual for the PARCS v2.6 U.S. NRC Core Neutronics Simulator. School of Nuclear Engineering Purdue University, USA. Downar, T., Xu, Y., Seker, V., 2009. Theory Manual for the PARCS Kinetics Core Simulator Module. Department of Nuclear Engineering and Radiological Sciences University of Michigan, USA. Fliscounakis, M., Girardi, E., Courau, T., Couyras, D., 2012. Potential of pin-by-pin SPN calculations as industrial reference. In: Proceedings of ICAPP 2012, Chicago, USA, June 24–28, 2012, on CD-ROM. Gelbard, E.M., 1962. Application of the Simplified Spherical Harmonics Equation in Spherical Geometry, WAPD-TM-294. Bettis Atomic Power Laboratory. Herbert, A., 2008. A Raviart-Thomas-Schneider solution of diffusion equation in hexagonal geometry. Ann. Nucl. Energy 35, 363–376. Jalili, M.H., Pazirandeh, A., Athari, M., 2015a. New analytic function expansion nodal (AFEN) method for solving multigroup neutron simplified P3 equation. Ann. Nucl. Eng. 77, 148–160. Jalili, M.H., Pazirandeh, A., Athari, M., 2015b. Development of a new analytic function expansion nodal code, HexDANM, for solving the neutron diffusion equation in hexagonal-Z geometry. Kerntechnik 80, 529–536. Kim, Y.I. et al., 2009. A conformal mapped nodal SP3 method for hexagonal core analysis. Ann. Nucl. Energy 36, 498–504. Larsen, E.W., Morel, J.E., McGhee, J., 1993, Asymptotic derivation of the simplified Pn equations. In Proceedings of Joint International Conference on Mathematical Methods and Supercomputing in Nuclear Applications, Portland, Oregon. Larsen, E., Morel, J., McGhee, J., 1996. Asymptotic derivation of the multigroup P1 and simplified PN equations with anisotropic scattering. Nucl. Sci. Eng. 123 (3), 328–342. McClarren, R.G., 2011. Theoretical aspects of the simplified Pn equations. Transp. Theory Stat. Phys. 39, 73–109. Noh, J.M., Cho, N.Z., 1996. A multigroup diffusion nodal scheme in rectangular and hexagonal geometry: hybrid of AFEN and PEN methods. In: Proceedings of the International Conference on the Physics of Reactor (PHYSOR 96), vol. 1, Mito, Ibaraki, Japan, 16–20 September 1996. Atomic Energy Society of Japan. pp. A50– A59. Ryu, E.H., Joo, H.G., 2010. Development of a 2D simplified P3 FEM solver for arbitrary geometry application. Transaction of the Korean Nuclear Society Autumn Meeting Jeju, Korea, October 21–22, 2010. Ryu, E.H., Joo, H.G., 2013. Finite element method solution of the simplified P3 equations for general geometry applications. Ann. Nucl. Energy 56, 194–207. Stacey, W.M., 2007. Nuclear Reactor Physics. John Wiley & Sons Inc., New York, USA. Sugimura, N., Ushio, T., Yamamoto, A., Tatsumi, M., 2006. Calculation models of AEGIS/SCOPE2, a core calculation system of next generation. In: Proceedings of PHYSOR-2006, ANS Topical Meeting on Reactor Physics, Vancouver, BC, Canada, and September 10–14, 2006. Woo, S.W., Cho, N.Z., Noh, J.M., 2001. The analytic function expansion nodal method refined with transverse gradient basis functions and interface flux moments. Nucl. Sci. Eng. 139, 156. Xia, B., Xie, Z., 2006. Flux expansion nodal method for solving multigroup neutron diffusion equations in hexagonal-z geometry. Ann. Nucl. Energy 33, 370–376.