- Email: [email protected]
Three-dimensional high order nodal code, ACNECH, for the neutronic modeling of hexagonal-z geometry
Annals of Nuclear Energy 68 (2014) 172–182
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Three-dimensional high order nodal code, ACNECH, for the neutronic modeling of hexagonal-z geometry N. Poursalehi ⇑, A. Zolfaghari, A. Minuchehr Nuclear Engineering Department, Shahid Beheshti University, GC, P.O. Box 1983963113, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 10 July 2013 Received in revised form 21 November 2013 Accepted 8 January 2014 Available online 5 February 2014 Keywords: Diffusion equation Average current nodal expansion method Hexagonal geometry High order calculation
a b s t r a c t In this work, we developed a new high order nodal code for the neutronic analysis of hexagonal-z geometry using the first order accuracy of Average Current Nodal Expansion Method (ACNEM) for radial direction and the second order solution of ACNEM for axial direction. For this purpose, we prepared Average Current Nodal Expansion Code for three-dimensional Hexagonal geometries (ACNECH) in which it calculates with coarse meshes i.e. one node per hexagon assembly. In this code, we performed an adopted iterative approach for the solution of neutron diffusion equation in order to overbear some divergence situations. The results of ACNECH are validated against some well-known benchmark problems. Results show the accuracy of high order ACNEM calculation consuming relatively suitable computational time. In addition, for obtaining high order expansion coefficients, the computations are done using Galerkin and Moments weightings for the radial direction in which results confirm the higher accuracy of Moments weighting. At last, we can conclude that ACNECH can be used in the neutronic analysis of problems where the run time is important such as loading pattern optimization. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The fast and accurate determination of the instantaneous, threedimensional, spatial core power distribution and the total core reactivity can contribute significantly to the safety and control of a large commercial nuclear power plant. The more accurately the power distribution is obtained, the more optimally a plant can be operated. Also, the faster the power distribution and the reactivity are determined, the easier the power plant can be controlled. Conventionally, the power distribution is calculated by the diffusion theory approximation to the rigorous Boltzmann transport theory, which is impractical for use on a routine basis. Before the maturation of the nodal diffusion models, finite-difference methods for solving the diffusion equations were used. Although they can generate very accurate results, the requirement of very fine meshes (on the order of 1 cm) renders three-dimensional computations very expensive. Today, modern nodal methods are applied in the nuclear industry where use of relatively large meshes or nodes (about the size of a fuel assembly) is common and the computing speed is about two orders of magnitude faster than for the finite difference methods. This improvement in computing efficiency is due to a substantial reduction in the size of the system of equations to be solved. The basic strategy of nodal methods is ⇑ Corresponding author. Tel.: +98 21 22431595; fax: +98 21 29902546. E-mail address: [email protected] (N. Poursalehi). http://dx.doi.org/10.1016/j.anucene.2014.01.007 0306-4549/Ó 2014 Elsevier Ltd. All rights reserved.
to derive nodal balance equations and nodal coupling equations, and then solve these two sets of equations iteratively. The nodal balance equations involve node-averaged fluxes and surface-averaged currents. The nodal coupling equations provide the relationships between the node-averaged fluxes and the surface-averaged currents. In solving the coupling equations, the spatial shapes of the transverse leakages are not known, and various assumptions can be made. Procedures used to solve these coupling equations lead to the classification of the modern nodal method into two separate branches: analytic nodal methods, Smith (1979), and the polynomial nodal methods, Smith (1985) and Zerkle (1992). In the analytic nodal method, a quadratic shape of the transverse leakage is generally used and the coupling equations are solved analytically. On the other hand, in the modern polynomial nodal method, a special polynomial expansion technique is used to solve the coupling equations, Kuo (1993). Some advanced nodal methods have been developed for the reactor core analysis such as Nodal Expansion Method (NEM), Finnemann (1975), Nodal Green’s Function Method (NGFM), Lawrence (1979) and Analytic Function Expansion Method (AFEN), Noh and Cho (1993, 1994) in which AFEN was refined with transverse gradient basis functions and interface flux moments, Woo et al. (2001). Extension of AFEN for hexagonal-z geometry was also employed by Cho et al. (1997) and improved features and verification of this in hexagonal-z 3-D geometry presented by Cho et al. (2005). Furthermore, Chao and Shatilla (1995) and Chao and Tsolfanidis (1995) have been
173
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
provided a mathematical basis for an extension of nodal methods to hexagonal geometry using conformal mapping of the hexagon into a rectangle. In this work, a transverse integrated procedure is applied in the new coordinate system in which the methods developed for Cartesian geometry can be used, Zimin and Baturin (2002). Nodal Expansion Method (NEM) plays an important role in nodal diffusion calculations, Finnemann et al. (1977) and Wagner et al. (1977). NEM, which is the one of polynomial nodal methods, originated from the Nodal Synthesis Method by elimination of the finemesh calculations and improvement of the calculation of inter-node coupling coefficients. Furthermore, the NEM presents very suitable features, such as simplicity, computational efficiency and acceptable accuracy for most realistic nuclear reactor simulations. Recently, Poursalehi et al. (2012) exploited the zeroth order nodal expansion methods for three-dimensional Cartesian geometries in order to solve the neutron balance equation and then, compared their results together as respect to accuracy and convergence treatment. Moreover, they provided a high order Average Current Nodal Expansion Code (ACNEC) for 3D rectangular geometries in which this code can use the zeroth to second orders of ACNEM for performing the neutronic calculation of reactor core with fuel assemblies by rectangular shape, Poursalehi et al. (2013a). But in this research, we developed a high order ACNEM for hexagonal-z geometries. This method exploits the first order ACNEM in the radial direction using Galerkin and Moments weightings and the second order ACNEM in the axial direction using Moments weighting. For this purpose, we provided a 3D nodal code i.e. Average Current Nodal Expansion Code for Hexagonal-z geometries (ACNECH). In this code, an iterative scheme is also used for the diffusion equation solution. Numerical results of ACNECH demonstrate the efficiency of method because of gaining the relatively good accuracy consuming the suitable run time. The organization of the remainder of this paper is as follows. Section 2 presents the summary of high order average current nodal expansion method in 3D hexagonal geometry, in Section 3, the proposed iterative solution algorithm is given, Section 4 presents the results of ACNECH against some popular benchmarks and at last, conclusions are given in Section 5.
must be sufficient to provide the necessary additional equations to determine the values. Thus, the fundamental problem in deriving a hexagonal geometry nodal expansion method is to select a nodal expansion which satisfies the above requirements, and which leads to a system of nodal equations which are consistent with the inherent properties of a hexagon. The nodal equations for a 3D hexagonal geometry average current method, which has first order accuracy in the radial direction and second order accuracy in the axial direction, are developed below. 2.1. Nodal coordinate systems A typical node in 3D hexagonal geometry is shown in Fig. 1. The fundamental point to note is that, despite its name, the node only possesses hexagonal geometry characteristics in the radial directions; in the axial direction its characteristics are of an essentially rectangular geometry nature. Thus, the derivation of a nodal expansion method for this geometry requires the combination of a 2D hexagonal geometry analysis in the radial directions with a 1D rectangular geometry analysis in the axial direction. Evidently, this necessitates the use of a nodal expansion which, in the radial directions, reflects the inherent symmetry of a hexagon, but in the axial direction models the properties of a 1D slab. Irrespective of the geometry, in order to expand the neutron flux in each node it is necessary to introduce a set, or sets, of axes in its dependent variables. For our treatment of hexagonal geometry in radial direction, we shall introduce into each node both the local ‘‘rectangular Cartesian axes’’ (x, y) and local ‘‘hexagonal axes’’ (x, u, v) shown in Fig. 2. Thus, for 3D hexagonal geometry, we have the local rectangular-z Cartesian axes (x, y, z) and the local hexagonal-z axes (x, u, v, z). On each axis we define the local dimensionless variables:
nw ¼
w ; H
nz ¼
z m hz
w ¼ x; u; v ; y ð1Þ m
2. Average current nodal expansion method for hexagonal-z geometry In this section, we present the high order ACNEM for threedimensional hexagonal geometries, using the fundamentals of ACNEM in rectangular geometry, Putney (1984), Poursalehi et al. (2012, 2013a). However, we can point out that ACNEM can be interpreted as the solution of the (continuous) neutron diffusion equation via the incomplete determination of the coefficients of a polynomial flux expansion in each node, using a combination of weighted residual techniques and continuity conditions. Indeed, sufficient conditions to define the coefficients were obtained through the imposition of the neutron balance and continuity conditions then led to the necessary additional nodal equations to determine the nodal values. It is clearly unreasonable to expect all characteristics of the rectangular geometry nodal expansion methods to carry over to other geometries. Nevertheless, it is apparent that the basic requirements for a method to be classed as a nodal expansion method are as follows: (a) In its zeroth order variant, when the nodal expansions are fitted to the nodal values, the neutron balance and continuity conditions must be sufficient to determine them. (b) In the higher order variants, the higher order polynomials must be chosen to retain the fit of the zeroth order expansion and, after the local weighted diffusion equations have been applied, the neutron balance and continuity conditions
where H is the across flats width of a hexagon node and hz is the thickness of node Pm in the z-direction. Also, rectangular and hexagonal variables in radial directions are related by:
u¼
pffiffiffi 3y x ; 2
v¼
pffiffiffi 3y x 2
ð2Þ
The local hexagonal-z axes form the basis of our notation: m Cm ws , left (s = l)/right (s = r) w-surface of node P , w = x, u, v, z, m Um g , average flux for group g in P , m Wm gws , average flux for group g at Cws ,
Fig. 1. A three-dimensional hexagonal geometry.
174
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
1 1 n2w ; h5 ðnw Þ ¼ nw n2w 4 12
w ¼ x; u; v
ð6Þ
The various order polynomials gi(nz) for the axial direction also are obtained as:
g 1 ðnz Þ ¼ nz g 2 ðnz Þ ¼ n2z
1 12
1 g 3 ðnz Þ ¼ nz n2z 4 1 1 n2z n2z 4 20
g 4 ðnz Þ ¼ Fig. 2. Cartesian and hexagonal axes systems for the radial direction of a hexagon node.
wm gz , one dimensional spatially averaged flux in the z-direction of node Pm,
þm
m
jgws ; jgws , average outgoing (+) and incoming () partial currents for group g at Cm ws , km gws ,
value of boundary condition at C
Am ws ,
area of Cm ws ,
m ws ,
em ws , unit vector in the direction of the outward normal to C
m ws ,
m Pmws, node adjacent to surface Cm ws (of node P ).
2.2. A high order ACNEM formulation for three-dimensional hexagonal geometry The starting point for the derivation is the 3D multi group neutron continuity equation integrated over each hexagon node Pm:
n o 1 Xn o X 2 þm m þm m m jgws jgws þ m jgzs jgzs þ Rm rg Ug 3H w ¼ x; u; v hz s¼l;r ¼
s ¼ l; r G X
m Rm gg 0 Ug 0 þ
0
g ¼1 g 0 –g m ¼ 1; 2; . . . ; M;
Pm
/g dPm ¼ Um g
pffiffiffi Z 3 m /g dCm m ws ¼ Wgws ; Hhz Cmws 2 pffiffiffi 2 3H
Z Cm zs
w ¼ x; u; v
m /g dCm zs ¼ Wgzs ;
; s ¼ 1; r
s ¼ l; r
ð8Þ
For determining the higher order expansion coefficients of axial direction in Eq. (4) i.e. dgz1 and dgz2, the nodal expansions are forced to satisfy the weighted residual equations: 8 > > > > > <
Z P
m
9 > > > > > =
G G X v X W k ðzÞ div Dg r/g þ Rrg /g Rgg0 /g0 g tRfg0 /g0 dPm ¼ 0; > > k eff g 0 ¼1 > > 0 > > g ¼1 > > > > : ; g 0 –g
g ¼ 1; 2; . . . ; G;
k ¼ 1; 2
ð9Þ
where using the Moments weighting scheme, we have:
W 1 ¼ g1
ð10Þ
W 2 ¼ g2
g 0 ¼1
g ¼ 1; 2; . . . ; G
ð3Þ
In order to solve Eq. (3), we define the neutron flux distribution in each node and energy group g as the polynomial expansion i.e. Eq. (4) using the first order accuracy in the radial directions and the second order accuracy in the axial direction:
/g ðnx ; nu ; nv ; nz Þ ¼ Ag h0 þ agx h1 ðnx Þ þ bgx h2 ðnx Þ þ agu h1 ðnu Þ
2
Wgz ¼ pffiffiffi ð4Þ
where hi(nw) and gi(nz) are polynomials of degree i in nw(w = x, u, v) and nz, respectively. Using orthogonal characteristics of hi(nw) polynomials, they are obtained as:
h1 ðnw Þ ¼ nw ;
h2 ðnw Þ ¼
1 12
ð11Þ m
where if (P )2 denotes a horizontal plane through P
þ dgu h5 ðnu Þ þ dgv h5 ðnv Þ þ agz g 1 ðnz Þ
n2w
9 8 > > > > > > > > > > Z 2 = < G G X d Wm vg X gz m m m m m m m 0W 0 W ðzÞ D þ R W þ L R W t R dz ¼ 0; 0 0 rg k gz gz gz g gz gg fg 2 hm > > z k eff g 0 ¼1 2 d z > > > > g0 ¼ 1 > > > > ; : g 0 –g hm z 2
m
þ cg h1 ðnx Þh1 ðnu Þh1 ðnv Þ þ dgx h5 ðnx Þ þ bgz g 2 ðnz Þ þ dgz1 g 3 ðnz Þ þ dgz2 g 4 ðnz Þ
Due to the orthogonality properties of the expansion functions, Eq. (9) takes on the form of the one-dimensional weighted diffusion equation
m ¼ 1;2; .. .;M; g ¼ 1; 2;.. .; G; k ¼ 1;2
þ bgu h2 ðnu Þ þ agv h1 ðnv Þ þ bgv h2 ðnv Þ
h0 ¼ 1;
Z
1
m ¼ 1; 2; . . . ; M;
G vg X tRmfg0 Umg0 ;
keff
In order to determine the zeroth expansion coefficients, the polynomial of Eq. (4) is fitted to the average flux in the node and the average fluxes on its surfaces. That is, we impose the conditions:
vm
m Lm gz , z-dependent transverse leakage in node P ,
ð7Þ
ð5Þ
But in Eq. (4), dgw (w = x, u, v) are the first order expansion coefficients and the suitable polynomials h5(nw) are considered as:
Z
3H 2
ðPm Þ2
/g dðPm Þ2
ð12Þ
is the one-dimensional spatially averaged flux in the z-direction, and
2 Lgz ¼ pffiffiffi 2 3H
Z ðPm Þ 2
" Dg
2
d /g 2
d x
2
þ Dg
d /g 2
d y
# dðPm 2Þ
ð13Þ
is the z-dependent transverse leakage. The one dimensional spatially averaged flux for second order solution in the z-direction is also defined as:
175
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Wgz ¼ Um g þ agz g 1 ðnz Þ þ bgz g 2 ðnz Þ þ dgz1 g 3 ðnz Þ þ dgz2 g 4 ðnz Þ
ð14Þ
For a nodal expansion of up to 4th degree in the z-direction, the z-dependent transverse leakage variation will be of second degree, so that a consistent approximation is:
Lgz ðnz Þ ¼ Lgz0 þ Lgz1 g 1 ðnz Þ þ Lgz2 g 2 ðnz Þ
ð15Þ
In Eq. (15), Lm gz0 can be identified with the average z-dependent transverse leakage in the node, Lm gz , with the following relation:
" # Z hmz Z 2 2 2 d /g d /g m 1 2 ffiffiffi p D þ D d P2 dz g g m m 2 2 hz 3H2 h2z ðPm2 Þ d x d y o 2 X n þm m þm m j jgwr þ jgwl jgwl ¼ 3H w¼x;u;v gwr
Lgz ¼
ð16Þ
The most convenient means of determining the first and second m transverse leakage moments Lm gz1 and Lgz2 , is to fit the expansion (15) to the values of the average z-dependent; thereby enabling them to be calculated from the interface currents in these nodes. The resulting expressions for the transverse leakage moments take the form:
"
Lm gz1
#
Lm gz2
br 1 ¼ al br ar bl ar
bl al
"
m Lmzl gz Lgz
1 1 gx ¼ agx agu agv a 2 2
Dm g
m2 hz
)
m
þ Ak R m rg dgzk ¼
G X
n
8 > > > > > <
Pm
9 > > > > > =
G G X v X h5 ðnw Þ Dg r2 /g þ Rrg /g Rgg0 /g0 g tRfg0 /g0 dPm ¼ 0; > > k eff g 0 ¼1 > > 0 > > g ¼ 1 > > > > : ; 0 g –g
m ¼ 1; 2; . . . ; M;
g ¼ 1; 2; . . . ; G;
w ¼ x; u; v
ð19Þ
After evaluating the weighted integrals of Eq. (19) for the x-direction, we have:
61 13 13 d d d gx gu g v 9 36 36 H2
233 1 179 49 49 m gx cg þ dgx þ dgu þ dgv þ Rrg a 36 8 4752 19; 008 19; 008
G X 233 1 179 49 49 g 0 x cg 0 þ ¼ Rm dg 0 x þ dg0 u þ dg0 v a gg 0 36 8 4752 19; 008 19; 008 g0 ¼ 1 0 g –g
G vg X 233 1 179 49 49 g 0 x c g 0 þ þ tRmfg0 dg0 x þ dg 0 u þ dg0 v ; a 36 8 4752 19;008 19; 008 keff g0 ¼1
Dm g
m ¼ 1; 2;. . . ; M; g ¼ 1; 2; . . . ; G
ð24Þ
ð18Þ
m where Ak, Bk are constant coefficients. Note that em gz1 and egz2 are the m zeroth order nodal expansion coefficients Am and B , respectively, gz gz Poursalehi et al. (2013a). Here we need to calculate the first order expansion coefficients of radial directions in Eq. (4) i.e. dgx, dgu and dgv. Using Galerkin weighting, these coefficients are the solutions of the weighted diffusion equations
Z
ð23Þ
( ) ( ) m m 61 Dg 179 m 13 Dg 49 m þ R dgx þ þ R dgu 9 H2 4752 rg 36 H2 19; 008 rg ( ) m 13 Dg 49 þ þ Rm dgv ¼ Q m gx 36 H2 19; 008 rg
m ¼ 1; 2; .. .; M; g ¼ 1;2;. .. ;
G; k ¼ 1; 2
will be known if the groups are treated in the order of decreasing energy (‘‘-’’ denotes the value from the previous iteration). Rearranging gives:
g0 ¼ 1 g 0 –g G n o vg X þ tRm0 Ak dmg0 zk Bk emg0 zk þ Bk Rmrg emgzk keff g0 ¼1 fg þ Bk Lm gzk ;
ð22Þ
233 ~ 1 g0 x þ ~cg a 36 8
G X 233 ~ 1 179 ~ 49 ~ 49 ~ g0 x ~cg0 þ þ Rm dg0 x þ dg 0 u þ dg0 v a gg 0 36 8 4752 19; 008 19; 008 g0 ¼ 1 g 0 –g
G vg X 233 ~ 1 179 ~ 49 ~ 49 ~ g0 x ~cg0 þ þ dg0 x þ dg 0 u þ dg0 v tRm0 a 36 8 4752 19; 008 19; 008 keff g0 ¼1 fg
m Qm gx ¼ Rrg
ð17Þ
o
g ¼ 1; 2; . . . ; G
in which
mzs
m m Rm gg 0 Ak dg 0 zk Bk eg 0 zk
61 13 13 d d d gx gu gv 9 36 36 H2
179 49 49 ¼ Qm d þ d þ d þ Rm gx gu g v gx ; rg 4752 19; 008 19; 008
Dm g
m ¼ 1; 2; . . . ; M;
which as and bs are the function of hz and hz . Using Moments weighting in Eq. (11), the resulting system for the 1st and 2st order expansion coefficients in the z-direction takes the form:
(
ð21Þ
The first order nodal equations can be solved using a nodal fission source iteration scheme analogous to that adopted for the rectangular geometry methods. In the fission source iteration scheme effectively the fission source is frozen at each outer iteration and similarly for Eq. (20), the weighted integral of fission source can be frozen and computed using values of expansion coefficients from the previous iteration. Thus, Eq. (20) can be written:
#
m Lmzr gz Lgz m
where
Making the same assumptions for the weighted residual equations using the weighting functions h5(nu) and h5(nv), leads to a system of the form:
2
32
2
3 ð25Þ
where m 61 Dg 179 m þ R 9 H2 4752 rg m 13 Dg 49 bm þ Rm g ¼ 36 H2 19; 008 rg
amg ¼
ð26Þ
m Qm gu and Q g v follow by symmetry from Eq. (23). Solving the system of Eq. (25), we obtain:
dgx ¼
ð20Þ
3
Qm amg bmg bmg dgx gx 6 6 bm am bm 76 7 6 m 7 7 4 g g g 54 dgu 5 ¼ 4 Q gu 5 m m m m dgv bg bg ag Q gv
amg þ bmg Q mgx bmg Q mgu bmg Q mgv
amg bmg
amg þ 2bmg
m m m m m m bm g Q gx þ ag þ bg Q gu bg Q g v dgu ¼ amg bmg ðamg þ 2bmg Þ m m m m m m bm g Q gx bg Q gu þ ag þ bg Q g v dgv ¼ amg bmg amg þ 2bmg
ð27Þ
176
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Moments weighting also can be used in Eq. (19). In this work, we select the weighting function n3w as a Moments weighting. Therefore, the first order x-directed weighted diffusion equations arising from these Moments weighting scheme is given below, this for u- and v-directions follow by symmetry.
283 251 251 dgx dgu dgv 84 336 336 H G G X v X m 0m m 0m Rgg0 Ug0 þ g tRm0 U0m0 ; þ Rrg Ug ¼ keff g0 ¼1 fg g 0 g ¼1 g 0 –g m ¼ 1; 2; . . . ; M; g ¼ 1; 2; . . . ; G Dm g 2
Fig. 3. The flowchart of implemented adopted iterative scheme in ACNECH.
ð28Þ
177
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Similarly at surface Cm xl :
where
7 3
gx U0m a g ¼
1 1 1 1 cg þ dgx dgu dgv 32 112 896 896
ð29Þ
Having determined the high order expansion coefficients for the radial and axial directions, we can now proceed to calculate the interface currents in the integral neutron balance equation (3) for a hexagon node analogous to process in rectangular geometry treatment, Putney (1984) and Poursalehi et al. (2013a). Considering surface Cm xr we have, from Fick’s law: pffiffiffi Z pffiffiffi Z hmz
Z pHffi 2 @/g
2 3 3 3 þm m m jgxr jgxr ¼ m Dg r /g e m d C ¼ Dg dy m dz m xr xr h Hffi @x x¼H Hhz Cmxr Hhz 2z p 2 2 3 7 29 Dm g 47 þm m þm m þm m ¼ j þ jgxr þ jgxl þ jgxl þ j þj 5 10 gur gur H 5 gxr 11 þm m 29 þm m 11 þm m 36 m 1 jgul þ jgul þ jg v r þ jgv r jg v l þ jgv l Ug þ dgx 10 10 10 5 12 ð30Þ
þm
m
jgxl jgxl ¼
7 29 Dm g 47 þm m þm m þm m jgxl þ jgxl þ jgxr þ jgxr þ jgul þ jgul 5 10 H 5
29 11 þm m þm m jgur þ jgur þ jg v l þ jg v l 10 10 36 11 þm 1 m jg v r þ jg v r Um dgx g 10 5 12
Alike expressions for the other six surfaces (in the radial and axial directions) of node Pm follow by symmetry. The complete system, after separating the incoming and outgoing currents, takes the form of Eq. (32) as the interface current equations:
2
2
þm
3
2 m cg1 6 þm 7 6 jgxl 7 6 cm 7 6 g2 6 6 þm 7 6 m 6 jgur 7 6 cg4 7 6 6 6 þm 7 6 cm 6 jgul 7 6 g3 7 6 6 6 jþm 7 ¼ 6 cm 6 gv r 7 6 g4 6 þm 7 6 cm 6 j 7 6 g3 6 gv l 7 6 6 þm 7 6 0 6j 7 4 4 gzr 5 0 þm jgzl jgxr
cm g2
cm g4
cm g3
cm g4
cm g3
cm g1 cm g3 cm g4 cm g3 cm g4
cm g3 cm g1 cm g2 cm g4 cm g3
cm g4 cm g2 cm g1 cm g3 cm g4
cm g3 cm g4 cm g3 cm g1 cm g2
cm g4 cm g3 cm g4 cm g2 cm g1
0
cm g5
cm g6
cm g7
cm g7
0
0
0
cm g7
0
0
0
cm g7 cm g6 cm g6 cm g7 cm g7
0
0
cm g6 cm g7 cm g7 cm g7 cm g7
cm g7
0
cm g5 cm g5 cm g5 cm g5 cm g5 Am gz Am gz
cm g7
0
cm g6 cm g6
0 0
0 0
0 0
0
0
0
0
0
Cm gz
Bm gz
0
0
0
Dm gz
0
0
0
0
0
Bm gz
Cm gz
0
0
0
Dm gz
0
ð31Þ
m
jgxr
3
6 m 7 6 jgxl 7 7 6 6 m 7 6 jgur 7 36 m 7 7 0 6 6 jgul 7 6 m 7 7 0 76 j vr 7 7 76 7 6 gm 0 7 7 76 j 6 gv l 7 7 0 76 m 7 6 76 jgzl 7 7 0 7 7 76 7 6 jm 7 gzr 7 6 0 76 m 7 7 U 7 76 6 g 7 Em gz 56 7 6 dgx 7 Em 7 gz 6 6 dgu 7 7 6 6 dgv 7 7 6 7 6 4 dgz1 5 dgz2
Fig. 4. Fuel assembly relative power distribution of LWR-2D problem including reference solution and ACNECH results.
ð32Þ
178
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
(
Table 1 Cross-sections of LWR-2D problem. Region
Group, g
Dg (cm)
Rag (cm1)
tRag (cm1)
R1?2 (cm1)
1
1 2 1 2 1 2 1 2
1.230991 0.43165 1.23996 0.426453 1.24486 0.421773 1.24567 0.254215
0.0108167 0.0793297 0.011486 0.0100011 0.011943 0.11373 0.000483146 0.0157667
0.00680319 0.112459 0.00839142 0.154606 0.00943943 0.182625 0 0
0.0117134 0 0.0109755 0 0.0105623 0 0.0322555 0
2 3 4
) m G G X vg X 4 m 2Agz m m Rm tRmfg0 Umg0 cg5 þ m þ Rm rg Ug ¼ gg 0 Ug 0 þ H k h eff g 0 ¼1 0 g ¼1 g 0 –g o X 2 n m m m m þ 1 cm jgws g1 c g2 2cg3 2cg4 3H w ¼ x; u; v s ¼ l;r o 1 n m m m m jgzl þ jgzr 2Em þ m 1 Bm gz C gz gz dgz2 hz
ð33Þ
3. Implemented iterative approach in ACNECH Table 2 Comparison of ACNECH results for LWR-2D problem.
a b c d
Code
Keff
ek (%)
em (%)
ea (%)
Run time (s)
VENTURE a ACNECH b ACNECH c ACNECH d
1.05603 1.05397 1.05503 1.05522
– 0.195 0.095 0.077
– 12.8 8.40 6.85
– 4.63 2.20 1.98
– 0.15 0.20 0.20
Reference. Zeroth order solution. First order solution with Galerkin weighting. First order solution with Moments weighting.
In the formal scheme of nodal expansion methods, some divergence may be occurred due to missing some required properties of produced matrix. In these situations, convergence difficulties with the nodal fission source iteration scheme may be encountered,
Table 3 Fuel assembly cross-sections of VVER-440 problem. FA type
Group, g
Dg (cm)
RRg (cm1)
tRfg (cm1)
R1?2 (cm1)
1
1 2 1 2 1 2 1 2 1 2 1 2
1.3466 0.37169 1.3377 0.36918 1.3322 0.36502 1.1953 0.19313 1.4485 0.25176 1.3413 0.24871
0.025255 0.064277 0.024709 0.079361 0.02435 0.1001 0.035636 0.13498 0.033184 0.032839 0.029301 0.064655
0.004488 0.073753 0.0055337 0.10581 0.0070391 0.14964 0 0 0 0 0 0
0.016893 – 0.015912 – 0.014888 – 0.022264
2 m m m where the nodal coefficients Am gz Egz and c g1 cg7 are the rational
functions of
Dm g hm z
and
Dm g , H
respectively.
3 4
The final form of nodal equations is obtained by using the interface current equations (32) to eliminate the outgoing currents from the nodal balance equation (3). Thus, we have Eq. (33) in order to obtain the hexagonal-z node average flux:
5 6
Fig. 5. The radial shape of VVER-440 full core.
0.032262 0.027148
179
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
such as a breakdown of off-diagonal elements of the matrix. In this case, the produced linear system cannot be solved and the solution process will be diverged, Poursalehi et al. (2013a). Thus, in the ACNECH, we used an adopted iterative scheme for the solution of nodal equations in the hexagonal-z geometries. The flowchart of proposed iterative process is present in Fig. 3. In this approach, we do not use the formal pseudo finite difference, Putney (1984) and instead we use Eq. (33) for computing the node average flux for any node and energy group. Hence, in order to keep away from ineligible situations, we abscond from the solving formal linear system in which has been nominated to inner iteration solution in Putney (1984). According to Fig. 3, in any general iteration, Eqs. (18), (27), (32), and (33) are used for calculating the required parameters including higher order expansion coefficients in axial and radial directions, outgoing partial currents and node average fluxes, respectively, for all of nodes and energy groups in first step. In this step, in spite of using Gauss Seidel algorithm, the calculations do not repeat until average flux convergence of all nodes and energy groups are taken place, and the solution algorithm proceeds by obtaining the new eigenvalue based on last obtained average flux of nodes and previous eigenvalue. While the eigenvalue and flux convergence criteria are not met, the next general iteration is carried out for all of nodes and energy groups using the updated eigenvalue and other parameters until convergence criteria is met, then the calculation will be stopped. Moreover, it should be noted that in ACNECH, we implemented the fission source weighting in the power method scheme for the accelerating of eigenvalue convergence. For obtaining more details about fission source weighting acceleration method, one can revert to Poursalehi et al. (2012, 2013a).
4. Validation of ACNECH In order to calculate neutronic parameters of a nuclear reactor core with fuel assemblies by hexagonal shape, we developed a new variant of ACNEM implementing the first order accuracy in radial direction and the second order solution in axial direction using an adopted iterative approach for converging. For this purpose,
Fig. 7. The axial shape of three-dimensional VVER-440 problem.
Table 5 Comparison of ACNECH results for three-dimensional VVER-440 problem. Table 4 Comparison of ACNECH results for two-dimensional VVER-440 problem.
a b c d
Code
Keff
ek (%)
em (%)
ea (%)
Run time (s)
DIF3D-FD a ACNECH b ACNECH c ACNECH d
1.00970 1.00765 1.00900 1.00901
– 0.203 0.069 0.068
– 10.9 6.0 5.3
– 4.3 2.2 1.9
– 0.20 0.35 0.35
Reference. Zeroth order solution. First order solution with Galerkin weighting. First order solution with Moments weighting.
Code
Keff
ek (%)
em (%)
ea (%)
Run time (s)
DIF3D-FD a ACNECH b ACNECH c ACNECH d
1.01132 1.00977 1.01059 1.01054
– 0.153 0.072 0.077
– 15.4 5.8 4.9
– 4.2 1.6 1.4
– 9 16 16
a
Reference. Zeroth order solution in both radial and axial directions. c First order solution in radial direction with Galerkin weighting and second order solution in axial direction. d First order solution in radial direction with Moments weighting and second order solution in axial direction. b
Fig. 6. Fuel assembly relative power distribution of VVER-440 2D problem including reference solution and ACNECH results.
180
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Table 6 Three-dimensional node-wise power for VVER-440 problem obtained by ACNECH. Assembly identity
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
Axial zones 2 (Bottom)
3
4
5
6
7
8
9
10
11 (Top)
0.5120 0.4225 0.5448 0.5591 0.4509 0.4620 0.6019 0.6329 0.7588 0.4687 0.5402 0.4424 0.5535 0.5593 0.5798 0.4897 0.5088 0.6751 0.3894 0.5479 0.4340 0.4389 0.5716 0.6006 0.7026 0.4970 0.5286 0.4350 0.4574 0.4623 0.5605 0.5565 0.5664 0.5945 0.3578 0.5065 0.4077
1.0633 0.8784 1.1329 1.1632 0.9387 0.9626 1.2548 1.3213 1.5876 0.9823 1.1230 0.9206 1.1515 1.1640 1.2078 1.0220 1.0644 1.4143 0.8147 1.1398 0.9034 0.9140 1.1912 1.2542 1.4706 1.0410 1.0997 0.9066 0.9550 0.9679 1.1742 1.1613 1.1848 1.2456 0.7492 1.0620 0.8547
1.4428 1.1932 1.5400 1.5820 1.2776 1.3103 1.7087 1.8034 2.1692 1.3430 1.5263 1.2520 1.5668 1.5848 1.6454 1.3945 1.4546 1.9333 1.1145 1.5509 1.2305 1.2463 1.6264 1.7148 2.0115 1.4244 1.4998 1.2386 1.3069 1.3256 1.6082 1.5894 1.6235 1.7074 1.0274 1.4563 1.1725
1.5795 1.3085 1.6934 1.7422 1.4054 1.4396 1.8793 1.9932 2.4097 1.4968 1.6771 1.3787 1.7259 1.7438 1.8110 1.5397 1.6140 2.1539 1.2442 1.7097 1.3569 1.3753 1.7994 1.9053 2.2438 1.5924 1.6579 1.3742 1.4561 1.4824 1.8025 1.7737 1.8197 1.9185 1.1557 1.6379 1.3209
1.3551 1.1726 1.5732 1.6337 1.3014 1.2826 1.6146 1.8041 2.2723 1.4342 1.5480 1.2931 1.6171 1.6037 1.6189 1.3889 1.5113 2.0585 1.1989 1.6092 1.2727 1.2789 1.6723 1.7984 2.1515 1.5385 1.5622 1.3031 1.3938 1.4326 1.7519 1.7083 1.7698 1.8757 1.1326 1.6049 1.2984
0.0000 0.8224 1.2692 1.3462 1.0417 0.8775 0.0000 1.2837 1.8706 1.2195 1.2273 1.0656 1.3293 1.2578 1.1145 0.9797 1.2225 1.7415 1.0288 1.3351 1.0476 1.0291 1.3343 1.4831 1.8295 1.3254 1.2951 1.0899 1.1822 1.2338 1.5225 1.4650 1.5396 1.6441 0.9960 1.4110 1.1465
0.0000 0.5906 0.9510 1.0202 0.7776 0.6250 0.0000 0.9421 1.4381 0.9549 0.9124 0.8077 1.0067 0.9317 0.7997 0.7149 0.9323 1.3599 0.8109 1.0169 0.7951 0.7744 1.0062 1.1415 1.4345 1.0484 0.9899 0.8415 0.9245 0.9764 1.2131 1.1537 1.2273 1.3191 0.8013 1.1349 0.9257
0.0000 0.4175 0.6731 0.7232 0.5515 0.4444 0.0000 0.6803 1.0422 0.6941 0.6455 0.5725 0.7144 0.6622 0.5715 0.5151 0.6751 0.9881 0.5905 0.7221 0.5658 0.5534 0.7243 0.8273 1.0442 0.7649 0.7082 0.6066 0.6714 0.7129 0.8881 0.8390 0.8984 0.9688 0.5894 0.8346 0.6821
0.0000 0.2608 0.4215 0.4534 0.3462 0.2791 0.0000 0.4315 0.6635 0.4423 0.4040 0.3588 0.4484 0.4163 0.3601 0.3262 0.4294 0.6299 0.3765 0.4531 0.3558 0.3492 0.4589 0.5261 0.6658 0.4881 0.4470 0.3843 0.4271 0.4547 0.5673 0.5338 0.5732 0.6195 0.3771 0.5337 0.4367
0.0000 0.1151 0.1857 0.1999 0.1528 0.1237 0.0000 0.1918 0.2940 0.1958 0.1781 0.1580 0.1978 0.1841 0.1599 0.1449 0.1902 0.2789 0.1669 0.1999 0.1572 0.1546 0.2036 0.2334 0.2951 0.2163 0.1980 0.1704 0.1894 0.2015 0.2515 0.2370 0.2543 0.2747 0.1674 0.2366 0.1938
Fig. 8. Assemblies identification for three-dimensional VVER-440 problem.
Average Current Nodal Expansion Code for Hexagonal 3D geometry (ACNECH) has been prepared. In this section, we present the numerical results of ACNECH for three well-known benchmarks including LWR-2D, VVER-440-2D and VVER-440-3D and for validating of ACNECH; its results are compared with references. It should be noted that all of calculations are performed by a laptop computer with CPU 2.41 GHz.
4.1. Benchmark problem (A): Two dimensional LWR First problem is a tight-lattice LWR core configuration using uranium oxide (UO2) fuel assemblies, shown in Fig. 4. This configuration consist of 151 heterogeneous hexagonal fuel assemblies with six different types of homogenized pin cells and the flux zero boundary conditions are imposed along the outer surfaces of the
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
reflector. The two-group homogenized cross sections for assemblies are shown in Table 1. The calculation with ACNECH was carried out in two-dimensional geometry with a relatively large node size (one node per assembly). The results of the VENTURE code, which uses the triangular mesh-centered finite difference method with 1352 triangular meshes per assembly, serve as the reference solution, Cho and Noh (1995). This problem (LWR-2D) has been investigated by developed high order nodal method using zeroth and first order solutions with Galerkin and Moments weightings. Table 2 gives the ACNECH numerical results of LWR-2D problem including keff and its relative error (ek (%)), maximum and average relative power errors (em (%), ea (%)) in comparison to reference solution and the run time (s). According to Table 2, the highest accuracy of results is belonged to the first order solution using Moments weighting with less Keff and power errors relative to the first order expansion using Galerkin weighting and especially, the zeroth order solution by increasing the flux expansion degree. Moreover, very short execution time of ACNECH and its efficiency is found from Table 2 due to using coarse meshes and implementing the adopted iterative algorithm with fission source weighting acceleration method. The relative power distribution of fuel assemblies as calculated by ACNECH along with the reference solution is given in Fig. 4. Fig. 4 shows a relatively high difference in solution obtained from ACNECH and the reference solution in the assemblies adjacent to the water reflector, where strong material discontinuity is present. 4.2. Benchmark problem (B): Two-dimensional VVER-440 Second test case is a two-dimensional VVER-440 core type, with 25 hexagonal fuel elements across the diameter, as the full core is shown in Fig. 5. The core has seven control rods inserted and a layer of reflector at the boundary of the core. One point to be noted is that control rods in VVER-440 do not insert into fuel assemblies; rather, they form control rod assemblies that push out fuel assemblies underneath them. Once inserted in the core, control rod assemblies replace fuel assemblies, giving rise to steep flux gradients on the interfaces between a control rod assembly and the neighboring fuel assemblies. The core has one-twelfth reflective symmetry and assembly pitch is 14.7 cm. Vacuum boundary conditions are considered at the external boundary of the reflector. The neutronic characteristics including material cross sections for this problem are given in Table 3. The reference solution is extrapolated from DIF3D-FD runs with 600 and 864 triangle/hexagon subdivisions, Chao and Shatilla (1995). The ACNECH has been used to work out the diffusion equation with two energy groups in onetwelfth symmetry of the core with hexagonal nodes having a fuel assembly dimensions. Table 4 presents the relative error percent of multiplication factor; maximum and average fuel assembly power obtained by zeroth and first order solutions with Galerkin and Moments weightings in comparison to reference solution. According to results, the higher and sufficient accuracy of first order approach is found relative to zeroth order solution, particularly, the first order solution with Moments weighting. In addition, Table 4 gives the suitable runtime of ACNECH for zeroth and first order solutions using coarse meshes, one node per a fuel assembly. As we know, the relative short time of neutronic solver is very important for the nuclear problems in which the solution process is very consuming such as core pattern optimization. Therefore, using ACNECH in the loading pattern optimization (LPO) is purposive, Poursalehi et al. (2013b). However, the fuel assembly power distribution of VVER-440-2D obtained by ACNECH and reference are given in Fig. 6. The major power errors in periphery fuel assemblies and also the control rod fuel assemblies neighboring are present in Fig. 6 that these errors are created because of high flux gradient existence in these regions.
181
4.3. Benchmark problem (C): Three-dimensional VVER-440 Last test case is based on the same core for problem (B) in which the original two-dimensional core is extended to a three-dimensional core with 250.0 cm high. A reflector 25.0 cm thick is added to the top and bottom of the core. All control rods are halfway inserted in the core with the material number 4 and for another half, with the material number 2, as shown in Fig. 7. Vacuum boundary conditions are applied to the entire reflector outside boundaries. Also, the reference solution was extrapolated from DIF3D-FD runs with 216 and 294 triangle/hexagon subdivisions and 2.5 cm axial mesh spacing, Chao and Shatilla (1995). ACNECH has been ran for the three-dimensional VVER-440 problem using coarse nodes i.e. one node per a fuel assembly in radial shape and axial sections with 25.0 cm height. Table 5 gives the numerical results of ACNECH for 3D VVER-440 test case. The zeroth order solution in both radial and axial directions has been exploited and Table 5 gives the results. According to Table 5, implementing the first order expansion in the radial direction and second order accuracy in the axial direction, the accuracy of results has been improved three times relative to zeroth order solution by decreasing the maximum power error from 15.4% to 4.9%. Of course, similar to previous problems, the first order accuracy using Moments weighting has the better results in comparison to Galerkin weighting. Furthermore, Table 5 indicates the relatively appropriate execution time of ACNECH for a 3D problem due to coarse mesh calculation. At last, Table 6 presents the three-dimensional node-wise powers obtained by ACNECH based on the assembly identification numbers as shown in Fig. 8. 5. Conclusions In this work, we developed a high order Average Current Nodal Expansion Code, ACNECH, for hexagonal-z geometries using the first order accuracy in radial directions and the second order solution in axial direction exploiting an adopted iterative scheme. ACNECH has been tested against three popular benchmark problems. Obtained results present that using high order solution, the accuracy of results enhances greatly relative to zeroth order solution. Actually, results confirm that the proposed high order methodology can obtain the adequate accuracy for the neutronic treatment analysis of nuclear reactor core with fuel assemblies by hexagonal shape. Of course, the numerical errors can be reduced by increasing the order of solution, particularly in radial directions. But this accuracy is achieved in the suitable computational cost and time. Thus for some problems such as LPO, the obtained accuracy can be sufficient because the run time of neutronic solver is very significant in the LPO. In addition, we found the higher accuracy of results using the Moments weighting for obtaining the first order expansion coefficients in regard of Galerkin weighting. Totally, the numerical results demonstrate the ability of ACNECH to calculate the stationary neutronic parameters of reactor core with hexagonal assembly shape. Furthermore, ACNECH has the capability to use in the LPO problem and can be extended for solving the time dependent nuclear problems in order to simulate the realistic reactor core neutronic treatments. References Chao, Y.A., Shatilla, Y.A., 1995. Conformal mapping and hexagonal nodal methods— II: Implementation in the ANC-H code. Nucl. Sci. Eng. 121, 210–225. Chao, Y.A., Tsolfanidis, N., 1995. Conformal mapping and hexagonal nodal methods—I: Mathematical foundation. Nucl. Sci. Eng. 121, 202–209. 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 multi-group problems in hexagonal-z geometry. Nucl. Sci. Eng. 126, 35.
182
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Cho, N.Z., Lee, J., Kim, D.S., Yang, C.Y., Jo, J.C., 2005. Improved Features and Verification of an AFEN Method Code in Hexagonal-z 3-D Geometry for Neutron Diffusion Calculation submitted to American Nuclear Society Annual Meeting, San Diego, CA. Finnemann, H., 1975. A consistent nodal method for the analysis of space-time effects in large LWRs. In: Proc. of the Joint NEACRP/CSNI Specialists Meeting on New Developments in Three-Dimensional Neutron Kinetics and Review of Kinetics Benchmark Calculations. MRR 145, January 22–24, Munich, p. 131. Finnemann, H., Bennewitz, F., Wagner, M.R., 1977. Interface current techniques for multidimensional reactor calculations. Atomkernenergie 30, 123–128. Kuo, W.S., 1993. The General Evaluation of the Nodal Synthesis Method in Nuclear Reactor Transient. Submitted to the Department of Nuclear Engineering in partial fulfillment of the requirements for the degree of Doctor of Philosophy at The Massachusetts Institute of Technology. Lawrence, R.D., 1979. A nodal Green’s Function Method for Multidimensional Neutron Diffusion Calculations. Thesis, Nuclear Engineering Program, Univ. of Illinois, Urbana. Noh, J.M., Cho, N.Z., 1993. A new diffusion nodal method based on analytic basis function expansion. Trans. Am. Nucl. Soc. 69, 462–463. Noh, J.M., Cho, N.Z., 1994. A new approach of analytic basis function expansion to neutron diffusion calculation. Nucl. Sci. Eng. 116, 165. Poursalehi, N., Zolfaghari, A., Minuchehr, A., 2012. Performance comparison of zeroth order nodal expansion methods in 3D rectangular geometry. Nucl. Eng. Des. 252, 248–266. Poursalehi, N., Zolfaghari, A., Minuchehr, A., 2013a. Development of a high order and multi-dimensional nodal code, ACNEC3D, for reactor core analysis. Ann. Nucl. Energy 55, 211–224.
Poursalehi, N., Zolfaghari, A., Minuchehr, A., Valavi, K., 2013b. Self-adaptive global best harmony search algorithm applied to reactor core fuel management optimization. Ann. Nucl. Energy 62, 86–102. Putney, J.M., 1984. Nodal Methods for Solving the Diffusion Equation for Fast Reactor Analysis. Ph.D. Thesis. University of London. Smith, K.S., 1979. An Analytic Nodal Method for Solving The Two-Group, Multidimensional, Static and Transient Neutron Diffusion Equations. Submitted in Partial Fulfillment of The Requirements for The Degrees Of Nuclear Engineer and Master of Science at The Massachusetts Institute of Technology. Smith, K.S., 1985. QPANDA: an advanced nodal method for LWR analyses. Trans. Am. Nucl. Soc. 50, 532. Wagner, M.R., Finnemann, K.H., Koebke, K., Winter, J., 1977. Validation of the nodal expansion method and the depletion program MEDIUM-2 by bechmark calculations and direct comparisons with experiments. Atomkernenergie 30, 129–135. 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. Zerkle, M.L., 1992. Development of a Polynomial Nodal Methods with Flux and Current Discontinuity Factors, Ph.D. Thesis, Department of Nuclear Engineering, Massachusetts Institute of Technology, Cambridge, MA. Zimin, V.G., Baturin, D.M., 2002. Polynomial nodal method for solving neutron diffusion equations in hexagonal-z geometry. Ann. Nucl. Energy 29, 1105– 1117.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Three-dimensional high order nodal code, ACNECH, for the neutronic modeling of hexagonal-z geometry N. Poursalehi ⇑, A. Zolfaghari, A. Minuchehr Nuclear Engineering Department, Shahid Beheshti University, GC, P.O. Box 1983963113, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 10 July 2013 Received in revised form 21 November 2013 Accepted 8 January 2014 Available online 5 February 2014 Keywords: Diffusion equation Average current nodal expansion method Hexagonal geometry High order calculation
a b s t r a c t In this work, we developed a new high order nodal code for the neutronic analysis of hexagonal-z geometry using the first order accuracy of Average Current Nodal Expansion Method (ACNEM) for radial direction and the second order solution of ACNEM for axial direction. For this purpose, we prepared Average Current Nodal Expansion Code for three-dimensional Hexagonal geometries (ACNECH) in which it calculates with coarse meshes i.e. one node per hexagon assembly. In this code, we performed an adopted iterative approach for the solution of neutron diffusion equation in order to overbear some divergence situations. The results of ACNECH are validated against some well-known benchmark problems. Results show the accuracy of high order ACNEM calculation consuming relatively suitable computational time. In addition, for obtaining high order expansion coefficients, the computations are done using Galerkin and Moments weightings for the radial direction in which results confirm the higher accuracy of Moments weighting. At last, we can conclude that ACNECH can be used in the neutronic analysis of problems where the run time is important such as loading pattern optimization. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The fast and accurate determination of the instantaneous, threedimensional, spatial core power distribution and the total core reactivity can contribute significantly to the safety and control of a large commercial nuclear power plant. The more accurately the power distribution is obtained, the more optimally a plant can be operated. Also, the faster the power distribution and the reactivity are determined, the easier the power plant can be controlled. Conventionally, the power distribution is calculated by the diffusion theory approximation to the rigorous Boltzmann transport theory, which is impractical for use on a routine basis. Before the maturation of the nodal diffusion models, finite-difference methods for solving the diffusion equations were used. Although they can generate very accurate results, the requirement of very fine meshes (on the order of 1 cm) renders three-dimensional computations very expensive. Today, modern nodal methods are applied in the nuclear industry where use of relatively large meshes or nodes (about the size of a fuel assembly) is common and the computing speed is about two orders of magnitude faster than for the finite difference methods. This improvement in computing efficiency is due to a substantial reduction in the size of the system of equations to be solved. The basic strategy of nodal methods is ⇑ Corresponding author. Tel.: +98 21 22431595; fax: +98 21 29902546. E-mail address: [email protected] (N. Poursalehi). http://dx.doi.org/10.1016/j.anucene.2014.01.007 0306-4549/Ó 2014 Elsevier Ltd. All rights reserved.
to derive nodal balance equations and nodal coupling equations, and then solve these two sets of equations iteratively. The nodal balance equations involve node-averaged fluxes and surface-averaged currents. The nodal coupling equations provide the relationships between the node-averaged fluxes and the surface-averaged currents. In solving the coupling equations, the spatial shapes of the transverse leakages are not known, and various assumptions can be made. Procedures used to solve these coupling equations lead to the classification of the modern nodal method into two separate branches: analytic nodal methods, Smith (1979), and the polynomial nodal methods, Smith (1985) and Zerkle (1992). In the analytic nodal method, a quadratic shape of the transverse leakage is generally used and the coupling equations are solved analytically. On the other hand, in the modern polynomial nodal method, a special polynomial expansion technique is used to solve the coupling equations, Kuo (1993). Some advanced nodal methods have been developed for the reactor core analysis such as Nodal Expansion Method (NEM), Finnemann (1975), Nodal Green’s Function Method (NGFM), Lawrence (1979) and Analytic Function Expansion Method (AFEN), Noh and Cho (1993, 1994) in which AFEN was refined with transverse gradient basis functions and interface flux moments, Woo et al. (2001). Extension of AFEN for hexagonal-z geometry was also employed by Cho et al. (1997) and improved features and verification of this in hexagonal-z 3-D geometry presented by Cho et al. (2005). Furthermore, Chao and Shatilla (1995) and Chao and Tsolfanidis (1995) have been
173
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
provided a mathematical basis for an extension of nodal methods to hexagonal geometry using conformal mapping of the hexagon into a rectangle. In this work, a transverse integrated procedure is applied in the new coordinate system in which the methods developed for Cartesian geometry can be used, Zimin and Baturin (2002). Nodal Expansion Method (NEM) plays an important role in nodal diffusion calculations, Finnemann et al. (1977) and Wagner et al. (1977). NEM, which is the one of polynomial nodal methods, originated from the Nodal Synthesis Method by elimination of the finemesh calculations and improvement of the calculation of inter-node coupling coefficients. Furthermore, the NEM presents very suitable features, such as simplicity, computational efficiency and acceptable accuracy for most realistic nuclear reactor simulations. Recently, Poursalehi et al. (2012) exploited the zeroth order nodal expansion methods for three-dimensional Cartesian geometries in order to solve the neutron balance equation and then, compared their results together as respect to accuracy and convergence treatment. Moreover, they provided a high order Average Current Nodal Expansion Code (ACNEC) for 3D rectangular geometries in which this code can use the zeroth to second orders of ACNEM for performing the neutronic calculation of reactor core with fuel assemblies by rectangular shape, Poursalehi et al. (2013a). But in this research, we developed a high order ACNEM for hexagonal-z geometries. This method exploits the first order ACNEM in the radial direction using Galerkin and Moments weightings and the second order ACNEM in the axial direction using Moments weighting. For this purpose, we provided a 3D nodal code i.e. Average Current Nodal Expansion Code for Hexagonal-z geometries (ACNECH). In this code, an iterative scheme is also used for the diffusion equation solution. Numerical results of ACNECH demonstrate the efficiency of method because of gaining the relatively good accuracy consuming the suitable run time. The organization of the remainder of this paper is as follows. Section 2 presents the summary of high order average current nodal expansion method in 3D hexagonal geometry, in Section 3, the proposed iterative solution algorithm is given, Section 4 presents the results of ACNECH against some popular benchmarks and at last, conclusions are given in Section 5.
must be sufficient to provide the necessary additional equations to determine the values. Thus, the fundamental problem in deriving a hexagonal geometry nodal expansion method is to select a nodal expansion which satisfies the above requirements, and which leads to a system of nodal equations which are consistent with the inherent properties of a hexagon. The nodal equations for a 3D hexagonal geometry average current method, which has first order accuracy in the radial direction and second order accuracy in the axial direction, are developed below. 2.1. Nodal coordinate systems A typical node in 3D hexagonal geometry is shown in Fig. 1. The fundamental point to note is that, despite its name, the node only possesses hexagonal geometry characteristics in the radial directions; in the axial direction its characteristics are of an essentially rectangular geometry nature. Thus, the derivation of a nodal expansion method for this geometry requires the combination of a 2D hexagonal geometry analysis in the radial directions with a 1D rectangular geometry analysis in the axial direction. Evidently, this necessitates the use of a nodal expansion which, in the radial directions, reflects the inherent symmetry of a hexagon, but in the axial direction models the properties of a 1D slab. Irrespective of the geometry, in order to expand the neutron flux in each node it is necessary to introduce a set, or sets, of axes in its dependent variables. For our treatment of hexagonal geometry in radial direction, we shall introduce into each node both the local ‘‘rectangular Cartesian axes’’ (x, y) and local ‘‘hexagonal axes’’ (x, u, v) shown in Fig. 2. Thus, for 3D hexagonal geometry, we have the local rectangular-z Cartesian axes (x, y, z) and the local hexagonal-z axes (x, u, v, z). On each axis we define the local dimensionless variables:
nw ¼
w ; H
nz ¼
z m hz
w ¼ x; u; v ; y ð1Þ m
2. Average current nodal expansion method for hexagonal-z geometry In this section, we present the high order ACNEM for threedimensional hexagonal geometries, using the fundamentals of ACNEM in rectangular geometry, Putney (1984), Poursalehi et al. (2012, 2013a). However, we can point out that ACNEM can be interpreted as the solution of the (continuous) neutron diffusion equation via the incomplete determination of the coefficients of a polynomial flux expansion in each node, using a combination of weighted residual techniques and continuity conditions. Indeed, sufficient conditions to define the coefficients were obtained through the imposition of the neutron balance and continuity conditions then led to the necessary additional nodal equations to determine the nodal values. It is clearly unreasonable to expect all characteristics of the rectangular geometry nodal expansion methods to carry over to other geometries. Nevertheless, it is apparent that the basic requirements for a method to be classed as a nodal expansion method are as follows: (a) In its zeroth order variant, when the nodal expansions are fitted to the nodal values, the neutron balance and continuity conditions must be sufficient to determine them. (b) In the higher order variants, the higher order polynomials must be chosen to retain the fit of the zeroth order expansion and, after the local weighted diffusion equations have been applied, the neutron balance and continuity conditions
where H is the across flats width of a hexagon node and hz is the thickness of node Pm in the z-direction. Also, rectangular and hexagonal variables in radial directions are related by:
u¼
pffiffiffi 3y x ; 2
v¼
pffiffiffi 3y x 2
ð2Þ
The local hexagonal-z axes form the basis of our notation: m Cm ws , left (s = l)/right (s = r) w-surface of node P , w = x, u, v, z, m Um g , average flux for group g in P , m Wm gws , average flux for group g at Cws ,
Fig. 1. A three-dimensional hexagonal geometry.
174
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
1 1 n2w ; h5 ðnw Þ ¼ nw n2w 4 12
w ¼ x; u; v
ð6Þ
The various order polynomials gi(nz) for the axial direction also are obtained as:
g 1 ðnz Þ ¼ nz g 2 ðnz Þ ¼ n2z
1 12
1 g 3 ðnz Þ ¼ nz n2z 4 1 1 n2z n2z 4 20
g 4 ðnz Þ ¼ Fig. 2. Cartesian and hexagonal axes systems for the radial direction of a hexagon node.
wm gz , one dimensional spatially averaged flux in the z-direction of node Pm,
þm
m
jgws ; jgws , average outgoing (+) and incoming () partial currents for group g at Cm ws , km gws ,
value of boundary condition at C
Am ws ,
area of Cm ws ,
m ws ,
em ws , unit vector in the direction of the outward normal to C
m ws ,
m Pmws, node adjacent to surface Cm ws (of node P ).
2.2. A high order ACNEM formulation for three-dimensional hexagonal geometry The starting point for the derivation is the 3D multi group neutron continuity equation integrated over each hexagon node Pm:
n o 1 Xn o X 2 þm m þm m m jgws jgws þ m jgzs jgzs þ Rm rg Ug 3H w ¼ x; u; v hz s¼l;r ¼
s ¼ l; r G X
m Rm gg 0 Ug 0 þ
0
g ¼1 g 0 –g m ¼ 1; 2; . . . ; M;
Pm
/g dPm ¼ Um g
pffiffiffi Z 3 m /g dCm m ws ¼ Wgws ; Hhz Cmws 2 pffiffiffi 2 3H
Z Cm zs
w ¼ x; u; v
m /g dCm zs ¼ Wgzs ;
; s ¼ 1; r
s ¼ l; r
ð8Þ
For determining the higher order expansion coefficients of axial direction in Eq. (4) i.e. dgz1 and dgz2, the nodal expansions are forced to satisfy the weighted residual equations: 8 > > > > > <
Z P
m
9 > > > > > =
G G X v X W k ðzÞ div Dg r/g þ Rrg /g Rgg0 /g0 g tRfg0 /g0 dPm ¼ 0; > > k eff g 0 ¼1 > > 0 > > g ¼1 > > > > : ; g 0 –g
g ¼ 1; 2; . . . ; G;
k ¼ 1; 2
ð9Þ
where using the Moments weighting scheme, we have:
W 1 ¼ g1
ð10Þ
W 2 ¼ g2
g 0 ¼1
g ¼ 1; 2; . . . ; G
ð3Þ
In order to solve Eq. (3), we define the neutron flux distribution in each node and energy group g as the polynomial expansion i.e. Eq. (4) using the first order accuracy in the radial directions and the second order accuracy in the axial direction:
/g ðnx ; nu ; nv ; nz Þ ¼ Ag h0 þ agx h1 ðnx Þ þ bgx h2 ðnx Þ þ agu h1 ðnu Þ
2
Wgz ¼ pffiffiffi ð4Þ
where hi(nw) and gi(nz) are polynomials of degree i in nw(w = x, u, v) and nz, respectively. Using orthogonal characteristics of hi(nw) polynomials, they are obtained as:
h1 ðnw Þ ¼ nw ;
h2 ðnw Þ ¼
1 12
ð11Þ m
where if (P )2 denotes a horizontal plane through P
þ dgu h5 ðnu Þ þ dgv h5 ðnv Þ þ agz g 1 ðnz Þ
n2w
9 8 > > > > > > > > > > Z 2 = < G G X d Wm vg X gz m m m m m m m 0W 0 W ðzÞ D þ R W þ L R W t R dz ¼ 0; 0 0 rg k gz gz gz g gz gg fg 2 hm > > z k eff g 0 ¼1 2 d z > > > > g0 ¼ 1 > > > > ; : g 0 –g hm z 2
m
þ cg h1 ðnx Þh1 ðnu Þh1 ðnv Þ þ dgx h5 ðnx Þ þ bgz g 2 ðnz Þ þ dgz1 g 3 ðnz Þ þ dgz2 g 4 ðnz Þ
Due to the orthogonality properties of the expansion functions, Eq. (9) takes on the form of the one-dimensional weighted diffusion equation
m ¼ 1;2; .. .;M; g ¼ 1; 2;.. .; G; k ¼ 1;2
þ bgu h2 ðnu Þ þ agv h1 ðnv Þ þ bgv h2 ðnv Þ
h0 ¼ 1;
Z
1
m ¼ 1; 2; . . . ; M;
G vg X tRmfg0 Umg0 ;
keff
In order to determine the zeroth expansion coefficients, the polynomial of Eq. (4) is fitted to the average flux in the node and the average fluxes on its surfaces. That is, we impose the conditions:
vm
m Lm gz , z-dependent transverse leakage in node P ,
ð7Þ
ð5Þ
But in Eq. (4), dgw (w = x, u, v) are the first order expansion coefficients and the suitable polynomials h5(nw) are considered as:
Z
3H 2
ðPm Þ2
/g dðPm Þ2
ð12Þ
is the one-dimensional spatially averaged flux in the z-direction, and
2 Lgz ¼ pffiffiffi 2 3H
Z ðPm Þ 2
" Dg
2
d /g 2
d x
2
þ Dg
d /g 2
d y
# dðPm 2Þ
ð13Þ
is the z-dependent transverse leakage. The one dimensional spatially averaged flux for second order solution in the z-direction is also defined as:
175
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Wgz ¼ Um g þ agz g 1 ðnz Þ þ bgz g 2 ðnz Þ þ dgz1 g 3 ðnz Þ þ dgz2 g 4 ðnz Þ
ð14Þ
For a nodal expansion of up to 4th degree in the z-direction, the z-dependent transverse leakage variation will be of second degree, so that a consistent approximation is:
Lgz ðnz Þ ¼ Lgz0 þ Lgz1 g 1 ðnz Þ þ Lgz2 g 2 ðnz Þ
ð15Þ
In Eq. (15), Lm gz0 can be identified with the average z-dependent transverse leakage in the node, Lm gz , with the following relation:
" # Z hmz Z 2 2 2 d /g d /g m 1 2 ffiffiffi p D þ D d P2 dz g g m m 2 2 hz 3H2 h2z ðPm2 Þ d x d y o 2 X n þm m þm m j jgwr þ jgwl jgwl ¼ 3H w¼x;u;v gwr
Lgz ¼
ð16Þ
The most convenient means of determining the first and second m transverse leakage moments Lm gz1 and Lgz2 , is to fit the expansion (15) to the values of the average z-dependent; thereby enabling them to be calculated from the interface currents in these nodes. The resulting expressions for the transverse leakage moments take the form:
"
Lm gz1
#
Lm gz2
br 1 ¼ al br ar bl ar
bl al
"
m Lmzl gz Lgz
1 1 gx ¼ agx agu agv a 2 2
Dm g
m2 hz
)
m
þ Ak R m rg dgzk ¼
G X
n
8 > > > > > <
Pm
9 > > > > > =
G G X v X h5 ðnw Þ Dg r2 /g þ Rrg /g Rgg0 /g0 g tRfg0 /g0 dPm ¼ 0; > > k eff g 0 ¼1 > > 0 > > g ¼ 1 > > > > : ; 0 g –g
m ¼ 1; 2; . . . ; M;
g ¼ 1; 2; . . . ; G;
w ¼ x; u; v
ð19Þ
After evaluating the weighted integrals of Eq. (19) for the x-direction, we have:
61 13 13 d d d gx gu g v 9 36 36 H2
233 1 179 49 49 m gx cg þ dgx þ dgu þ dgv þ Rrg a 36 8 4752 19; 008 19; 008
G X 233 1 179 49 49 g 0 x cg 0 þ ¼ Rm dg 0 x þ dg0 u þ dg0 v a gg 0 36 8 4752 19; 008 19; 008 g0 ¼ 1 0 g –g
G vg X 233 1 179 49 49 g 0 x c g 0 þ þ tRmfg0 dg0 x þ dg 0 u þ dg0 v ; a 36 8 4752 19;008 19; 008 keff g0 ¼1
Dm g
m ¼ 1; 2;. . . ; M; g ¼ 1; 2; . . . ; G
ð24Þ
ð18Þ
m where Ak, Bk are constant coefficients. Note that em gz1 and egz2 are the m zeroth order nodal expansion coefficients Am and B , respectively, gz gz Poursalehi et al. (2013a). Here we need to calculate the first order expansion coefficients of radial directions in Eq. (4) i.e. dgx, dgu and dgv. Using Galerkin weighting, these coefficients are the solutions of the weighted diffusion equations
Z
ð23Þ
( ) ( ) m m 61 Dg 179 m 13 Dg 49 m þ R dgx þ þ R dgu 9 H2 4752 rg 36 H2 19; 008 rg ( ) m 13 Dg 49 þ þ Rm dgv ¼ Q m gx 36 H2 19; 008 rg
m ¼ 1; 2; .. .; M; g ¼ 1;2;. .. ;
G; k ¼ 1; 2
will be known if the groups are treated in the order of decreasing energy (‘‘-’’ denotes the value from the previous iteration). Rearranging gives:
g0 ¼ 1 g 0 –g G n o vg X þ tRm0 Ak dmg0 zk Bk emg0 zk þ Bk Rmrg emgzk keff g0 ¼1 fg þ Bk Lm gzk ;
ð22Þ
233 ~ 1 g0 x þ ~cg a 36 8
G X 233 ~ 1 179 ~ 49 ~ 49 ~ g0 x ~cg0 þ þ Rm dg0 x þ dg 0 u þ dg0 v a gg 0 36 8 4752 19; 008 19; 008 g0 ¼ 1 g 0 –g
G vg X 233 ~ 1 179 ~ 49 ~ 49 ~ g0 x ~cg0 þ þ dg0 x þ dg 0 u þ dg0 v tRm0 a 36 8 4752 19; 008 19; 008 keff g0 ¼1 fg
m Qm gx ¼ Rrg
ð17Þ
o
g ¼ 1; 2; . . . ; G
in which
mzs
m m Rm gg 0 Ak dg 0 zk Bk eg 0 zk
61 13 13 d d d gx gu gv 9 36 36 H2
179 49 49 ¼ Qm d þ d þ d þ Rm gx gu g v gx ; rg 4752 19; 008 19; 008
Dm g
m ¼ 1; 2; . . . ; M;
which as and bs are the function of hz and hz . Using Moments weighting in Eq. (11), the resulting system for the 1st and 2st order expansion coefficients in the z-direction takes the form:
(
ð21Þ
The first order nodal equations can be solved using a nodal fission source iteration scheme analogous to that adopted for the rectangular geometry methods. In the fission source iteration scheme effectively the fission source is frozen at each outer iteration and similarly for Eq. (20), the weighted integral of fission source can be frozen and computed using values of expansion coefficients from the previous iteration. Thus, Eq. (20) can be written:
#
m Lmzr gz Lgz m
where
Making the same assumptions for the weighted residual equations using the weighting functions h5(nu) and h5(nv), leads to a system of the form:
2
32
2
3 ð25Þ
where m 61 Dg 179 m þ R 9 H2 4752 rg m 13 Dg 49 bm þ Rm g ¼ 36 H2 19; 008 rg
amg ¼
ð26Þ
m Qm gu and Q g v follow by symmetry from Eq. (23). Solving the system of Eq. (25), we obtain:
dgx ¼
ð20Þ
3
Qm amg bmg bmg dgx gx 6 6 bm am bm 76 7 6 m 7 7 4 g g g 54 dgu 5 ¼ 4 Q gu 5 m m m m dgv bg bg ag Q gv
amg þ bmg Q mgx bmg Q mgu bmg Q mgv
amg bmg
amg þ 2bmg
m m m m m m bm g Q gx þ ag þ bg Q gu bg Q g v dgu ¼ amg bmg ðamg þ 2bmg Þ m m m m m m bm g Q gx bg Q gu þ ag þ bg Q g v dgv ¼ amg bmg amg þ 2bmg
ð27Þ
176
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Moments weighting also can be used in Eq. (19). In this work, we select the weighting function n3w as a Moments weighting. Therefore, the first order x-directed weighted diffusion equations arising from these Moments weighting scheme is given below, this for u- and v-directions follow by symmetry.
283 251 251 dgx dgu dgv 84 336 336 H G G X v X m 0m m 0m Rgg0 Ug0 þ g tRm0 U0m0 ; þ Rrg Ug ¼ keff g0 ¼1 fg g 0 g ¼1 g 0 –g m ¼ 1; 2; . . . ; M; g ¼ 1; 2; . . . ; G Dm g 2
Fig. 3. The flowchart of implemented adopted iterative scheme in ACNECH.
ð28Þ
177
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Similarly at surface Cm xl :
where
7 3
gx U0m a g ¼
1 1 1 1 cg þ dgx dgu dgv 32 112 896 896
ð29Þ
Having determined the high order expansion coefficients for the radial and axial directions, we can now proceed to calculate the interface currents in the integral neutron balance equation (3) for a hexagon node analogous to process in rectangular geometry treatment, Putney (1984) and Poursalehi et al. (2013a). Considering surface Cm xr we have, from Fick’s law: pffiffiffi Z pffiffiffi Z hmz
Z pHffi 2 @/g
2 3 3 3 þm m m jgxr jgxr ¼ m Dg r /g e m d C ¼ Dg dy m dz m xr xr h Hffi @x x¼H Hhz Cmxr Hhz 2z p 2 2 3 7 29 Dm g 47 þm m þm m þm m ¼ j þ jgxr þ jgxl þ jgxl þ j þj 5 10 gur gur H 5 gxr 11 þm m 29 þm m 11 þm m 36 m 1 jgul þ jgul þ jg v r þ jgv r jg v l þ jgv l Ug þ dgx 10 10 10 5 12 ð30Þ
þm
m
jgxl jgxl ¼
7 29 Dm g 47 þm m þm m þm m jgxl þ jgxl þ jgxr þ jgxr þ jgul þ jgul 5 10 H 5
29 11 þm m þm m jgur þ jgur þ jg v l þ jg v l 10 10 36 11 þm 1 m jg v r þ jg v r Um dgx g 10 5 12
Alike expressions for the other six surfaces (in the radial and axial directions) of node Pm follow by symmetry. The complete system, after separating the incoming and outgoing currents, takes the form of Eq. (32) as the interface current equations:
2
2
þm
3
2 m cg1 6 þm 7 6 jgxl 7 6 cm 7 6 g2 6 6 þm 7 6 m 6 jgur 7 6 cg4 7 6 6 6 þm 7 6 cm 6 jgul 7 6 g3 7 6 6 6 jþm 7 ¼ 6 cm 6 gv r 7 6 g4 6 þm 7 6 cm 6 j 7 6 g3 6 gv l 7 6 6 þm 7 6 0 6j 7 4 4 gzr 5 0 þm jgzl jgxr
cm g2
cm g4
cm g3
cm g4
cm g3
cm g1 cm g3 cm g4 cm g3 cm g4
cm g3 cm g1 cm g2 cm g4 cm g3
cm g4 cm g2 cm g1 cm g3 cm g4
cm g3 cm g4 cm g3 cm g1 cm g2
cm g4 cm g3 cm g4 cm g2 cm g1
0
cm g5
cm g6
cm g7
cm g7
0
0
0
cm g7
0
0
0
cm g7 cm g6 cm g6 cm g7 cm g7
0
0
cm g6 cm g7 cm g7 cm g7 cm g7
cm g7
0
cm g5 cm g5 cm g5 cm g5 cm g5 Am gz Am gz
cm g7
0
cm g6 cm g6
0 0
0 0
0 0
0
0
0
0
0
Cm gz
Bm gz
0
0
0
Dm gz
0
0
0
0
0
Bm gz
Cm gz
0
0
0
Dm gz
0
ð31Þ
m
jgxr
3
6 m 7 6 jgxl 7 7 6 6 m 7 6 jgur 7 36 m 7 7 0 6 6 jgul 7 6 m 7 7 0 76 j vr 7 7 76 7 6 gm 0 7 7 76 j 6 gv l 7 7 0 76 m 7 6 76 jgzl 7 7 0 7 7 76 7 6 jm 7 gzr 7 6 0 76 m 7 7 U 7 76 6 g 7 Em gz 56 7 6 dgx 7 Em 7 gz 6 6 dgu 7 7 6 6 dgv 7 7 6 7 6 4 dgz1 5 dgz2
Fig. 4. Fuel assembly relative power distribution of LWR-2D problem including reference solution and ACNECH results.
ð32Þ
178
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
(
Table 1 Cross-sections of LWR-2D problem. Region
Group, g
Dg (cm)
Rag (cm1)
tRag (cm1)
R1?2 (cm1)
1
1 2 1 2 1 2 1 2
1.230991 0.43165 1.23996 0.426453 1.24486 0.421773 1.24567 0.254215
0.0108167 0.0793297 0.011486 0.0100011 0.011943 0.11373 0.000483146 0.0157667
0.00680319 0.112459 0.00839142 0.154606 0.00943943 0.182625 0 0
0.0117134 0 0.0109755 0 0.0105623 0 0.0322555 0
2 3 4
) m G G X vg X 4 m 2Agz m m Rm tRmfg0 Umg0 cg5 þ m þ Rm rg Ug ¼ gg 0 Ug 0 þ H k h eff g 0 ¼1 0 g ¼1 g 0 –g o X 2 n m m m m þ 1 cm jgws g1 c g2 2cg3 2cg4 3H w ¼ x; u; v s ¼ l;r o 1 n m m m m jgzl þ jgzr 2Em þ m 1 Bm gz C gz gz dgz2 hz
ð33Þ
3. Implemented iterative approach in ACNECH Table 2 Comparison of ACNECH results for LWR-2D problem.
a b c d
Code
Keff
ek (%)
em (%)
ea (%)
Run time (s)
VENTURE a ACNECH b ACNECH c ACNECH d
1.05603 1.05397 1.05503 1.05522
– 0.195 0.095 0.077
– 12.8 8.40 6.85
– 4.63 2.20 1.98
– 0.15 0.20 0.20
Reference. Zeroth order solution. First order solution with Galerkin weighting. First order solution with Moments weighting.
In the formal scheme of nodal expansion methods, some divergence may be occurred due to missing some required properties of produced matrix. In these situations, convergence difficulties with the nodal fission source iteration scheme may be encountered,
Table 3 Fuel assembly cross-sections of VVER-440 problem. FA type
Group, g
Dg (cm)
RRg (cm1)
tRfg (cm1)
R1?2 (cm1)
1
1 2 1 2 1 2 1 2 1 2 1 2
1.3466 0.37169 1.3377 0.36918 1.3322 0.36502 1.1953 0.19313 1.4485 0.25176 1.3413 0.24871
0.025255 0.064277 0.024709 0.079361 0.02435 0.1001 0.035636 0.13498 0.033184 0.032839 0.029301 0.064655
0.004488 0.073753 0.0055337 0.10581 0.0070391 0.14964 0 0 0 0 0 0
0.016893 – 0.015912 – 0.014888 – 0.022264
2 m m m where the nodal coefficients Am gz Egz and c g1 cg7 are the rational
functions of
Dm g hm z
and
Dm g , H
respectively.
3 4
The final form of nodal equations is obtained by using the interface current equations (32) to eliminate the outgoing currents from the nodal balance equation (3). Thus, we have Eq. (33) in order to obtain the hexagonal-z node average flux:
5 6
Fig. 5. The radial shape of VVER-440 full core.
0.032262 0.027148
179
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
such as a breakdown of off-diagonal elements of the matrix. In this case, the produced linear system cannot be solved and the solution process will be diverged, Poursalehi et al. (2013a). Thus, in the ACNECH, we used an adopted iterative scheme for the solution of nodal equations in the hexagonal-z geometries. The flowchart of proposed iterative process is present in Fig. 3. In this approach, we do not use the formal pseudo finite difference, Putney (1984) and instead we use Eq. (33) for computing the node average flux for any node and energy group. Hence, in order to keep away from ineligible situations, we abscond from the solving formal linear system in which has been nominated to inner iteration solution in Putney (1984). According to Fig. 3, in any general iteration, Eqs. (18), (27), (32), and (33) are used for calculating the required parameters including higher order expansion coefficients in axial and radial directions, outgoing partial currents and node average fluxes, respectively, for all of nodes and energy groups in first step. In this step, in spite of using Gauss Seidel algorithm, the calculations do not repeat until average flux convergence of all nodes and energy groups are taken place, and the solution algorithm proceeds by obtaining the new eigenvalue based on last obtained average flux of nodes and previous eigenvalue. While the eigenvalue and flux convergence criteria are not met, the next general iteration is carried out for all of nodes and energy groups using the updated eigenvalue and other parameters until convergence criteria is met, then the calculation will be stopped. Moreover, it should be noted that in ACNECH, we implemented the fission source weighting in the power method scheme for the accelerating of eigenvalue convergence. For obtaining more details about fission source weighting acceleration method, one can revert to Poursalehi et al. (2012, 2013a).
4. Validation of ACNECH In order to calculate neutronic parameters of a nuclear reactor core with fuel assemblies by hexagonal shape, we developed a new variant of ACNEM implementing the first order accuracy in radial direction and the second order solution in axial direction using an adopted iterative approach for converging. For this purpose,
Fig. 7. The axial shape of three-dimensional VVER-440 problem.
Table 5 Comparison of ACNECH results for three-dimensional VVER-440 problem. Table 4 Comparison of ACNECH results for two-dimensional VVER-440 problem.
a b c d
Code
Keff
ek (%)
em (%)
ea (%)
Run time (s)
DIF3D-FD a ACNECH b ACNECH c ACNECH d
1.00970 1.00765 1.00900 1.00901
– 0.203 0.069 0.068
– 10.9 6.0 5.3
– 4.3 2.2 1.9
– 0.20 0.35 0.35
Reference. Zeroth order solution. First order solution with Galerkin weighting. First order solution with Moments weighting.
Code
Keff
ek (%)
em (%)
ea (%)
Run time (s)
DIF3D-FD a ACNECH b ACNECH c ACNECH d
1.01132 1.00977 1.01059 1.01054
– 0.153 0.072 0.077
– 15.4 5.8 4.9
– 4.2 1.6 1.4
– 9 16 16
a
Reference. Zeroth order solution in both radial and axial directions. c First order solution in radial direction with Galerkin weighting and second order solution in axial direction. d First order solution in radial direction with Moments weighting and second order solution in axial direction. b
Fig. 6. Fuel assembly relative power distribution of VVER-440 2D problem including reference solution and ACNECH results.
180
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Table 6 Three-dimensional node-wise power for VVER-440 problem obtained by ACNECH. Assembly identity
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
Axial zones 2 (Bottom)
3
4
5
6
7
8
9
10
11 (Top)
0.5120 0.4225 0.5448 0.5591 0.4509 0.4620 0.6019 0.6329 0.7588 0.4687 0.5402 0.4424 0.5535 0.5593 0.5798 0.4897 0.5088 0.6751 0.3894 0.5479 0.4340 0.4389 0.5716 0.6006 0.7026 0.4970 0.5286 0.4350 0.4574 0.4623 0.5605 0.5565 0.5664 0.5945 0.3578 0.5065 0.4077
1.0633 0.8784 1.1329 1.1632 0.9387 0.9626 1.2548 1.3213 1.5876 0.9823 1.1230 0.9206 1.1515 1.1640 1.2078 1.0220 1.0644 1.4143 0.8147 1.1398 0.9034 0.9140 1.1912 1.2542 1.4706 1.0410 1.0997 0.9066 0.9550 0.9679 1.1742 1.1613 1.1848 1.2456 0.7492 1.0620 0.8547
1.4428 1.1932 1.5400 1.5820 1.2776 1.3103 1.7087 1.8034 2.1692 1.3430 1.5263 1.2520 1.5668 1.5848 1.6454 1.3945 1.4546 1.9333 1.1145 1.5509 1.2305 1.2463 1.6264 1.7148 2.0115 1.4244 1.4998 1.2386 1.3069 1.3256 1.6082 1.5894 1.6235 1.7074 1.0274 1.4563 1.1725
1.5795 1.3085 1.6934 1.7422 1.4054 1.4396 1.8793 1.9932 2.4097 1.4968 1.6771 1.3787 1.7259 1.7438 1.8110 1.5397 1.6140 2.1539 1.2442 1.7097 1.3569 1.3753 1.7994 1.9053 2.2438 1.5924 1.6579 1.3742 1.4561 1.4824 1.8025 1.7737 1.8197 1.9185 1.1557 1.6379 1.3209
1.3551 1.1726 1.5732 1.6337 1.3014 1.2826 1.6146 1.8041 2.2723 1.4342 1.5480 1.2931 1.6171 1.6037 1.6189 1.3889 1.5113 2.0585 1.1989 1.6092 1.2727 1.2789 1.6723 1.7984 2.1515 1.5385 1.5622 1.3031 1.3938 1.4326 1.7519 1.7083 1.7698 1.8757 1.1326 1.6049 1.2984
0.0000 0.8224 1.2692 1.3462 1.0417 0.8775 0.0000 1.2837 1.8706 1.2195 1.2273 1.0656 1.3293 1.2578 1.1145 0.9797 1.2225 1.7415 1.0288 1.3351 1.0476 1.0291 1.3343 1.4831 1.8295 1.3254 1.2951 1.0899 1.1822 1.2338 1.5225 1.4650 1.5396 1.6441 0.9960 1.4110 1.1465
0.0000 0.5906 0.9510 1.0202 0.7776 0.6250 0.0000 0.9421 1.4381 0.9549 0.9124 0.8077 1.0067 0.9317 0.7997 0.7149 0.9323 1.3599 0.8109 1.0169 0.7951 0.7744 1.0062 1.1415 1.4345 1.0484 0.9899 0.8415 0.9245 0.9764 1.2131 1.1537 1.2273 1.3191 0.8013 1.1349 0.9257
0.0000 0.4175 0.6731 0.7232 0.5515 0.4444 0.0000 0.6803 1.0422 0.6941 0.6455 0.5725 0.7144 0.6622 0.5715 0.5151 0.6751 0.9881 0.5905 0.7221 0.5658 0.5534 0.7243 0.8273 1.0442 0.7649 0.7082 0.6066 0.6714 0.7129 0.8881 0.8390 0.8984 0.9688 0.5894 0.8346 0.6821
0.0000 0.2608 0.4215 0.4534 0.3462 0.2791 0.0000 0.4315 0.6635 0.4423 0.4040 0.3588 0.4484 0.4163 0.3601 0.3262 0.4294 0.6299 0.3765 0.4531 0.3558 0.3492 0.4589 0.5261 0.6658 0.4881 0.4470 0.3843 0.4271 0.4547 0.5673 0.5338 0.5732 0.6195 0.3771 0.5337 0.4367
0.0000 0.1151 0.1857 0.1999 0.1528 0.1237 0.0000 0.1918 0.2940 0.1958 0.1781 0.1580 0.1978 0.1841 0.1599 0.1449 0.1902 0.2789 0.1669 0.1999 0.1572 0.1546 0.2036 0.2334 0.2951 0.2163 0.1980 0.1704 0.1894 0.2015 0.2515 0.2370 0.2543 0.2747 0.1674 0.2366 0.1938
Fig. 8. Assemblies identification for three-dimensional VVER-440 problem.
Average Current Nodal Expansion Code for Hexagonal 3D geometry (ACNECH) has been prepared. In this section, we present the numerical results of ACNECH for three well-known benchmarks including LWR-2D, VVER-440-2D and VVER-440-3D and for validating of ACNECH; its results are compared with references. It should be noted that all of calculations are performed by a laptop computer with CPU 2.41 GHz.
4.1. Benchmark problem (A): Two dimensional LWR First problem is a tight-lattice LWR core configuration using uranium oxide (UO2) fuel assemblies, shown in Fig. 4. This configuration consist of 151 heterogeneous hexagonal fuel assemblies with six different types of homogenized pin cells and the flux zero boundary conditions are imposed along the outer surfaces of the
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
reflector. The two-group homogenized cross sections for assemblies are shown in Table 1. The calculation with ACNECH was carried out in two-dimensional geometry with a relatively large node size (one node per assembly). The results of the VENTURE code, which uses the triangular mesh-centered finite difference method with 1352 triangular meshes per assembly, serve as the reference solution, Cho and Noh (1995). This problem (LWR-2D) has been investigated by developed high order nodal method using zeroth and first order solutions with Galerkin and Moments weightings. Table 2 gives the ACNECH numerical results of LWR-2D problem including keff and its relative error (ek (%)), maximum and average relative power errors (em (%), ea (%)) in comparison to reference solution and the run time (s). According to Table 2, the highest accuracy of results is belonged to the first order solution using Moments weighting with less Keff and power errors relative to the first order expansion using Galerkin weighting and especially, the zeroth order solution by increasing the flux expansion degree. Moreover, very short execution time of ACNECH and its efficiency is found from Table 2 due to using coarse meshes and implementing the adopted iterative algorithm with fission source weighting acceleration method. The relative power distribution of fuel assemblies as calculated by ACNECH along with the reference solution is given in Fig. 4. Fig. 4 shows a relatively high difference in solution obtained from ACNECH and the reference solution in the assemblies adjacent to the water reflector, where strong material discontinuity is present. 4.2. Benchmark problem (B): Two-dimensional VVER-440 Second test case is a two-dimensional VVER-440 core type, with 25 hexagonal fuel elements across the diameter, as the full core is shown in Fig. 5. The core has seven control rods inserted and a layer of reflector at the boundary of the core. One point to be noted is that control rods in VVER-440 do not insert into fuel assemblies; rather, they form control rod assemblies that push out fuel assemblies underneath them. Once inserted in the core, control rod assemblies replace fuel assemblies, giving rise to steep flux gradients on the interfaces between a control rod assembly and the neighboring fuel assemblies. The core has one-twelfth reflective symmetry and assembly pitch is 14.7 cm. Vacuum boundary conditions are considered at the external boundary of the reflector. The neutronic characteristics including material cross sections for this problem are given in Table 3. The reference solution is extrapolated from DIF3D-FD runs with 600 and 864 triangle/hexagon subdivisions, Chao and Shatilla (1995). The ACNECH has been used to work out the diffusion equation with two energy groups in onetwelfth symmetry of the core with hexagonal nodes having a fuel assembly dimensions. Table 4 presents the relative error percent of multiplication factor; maximum and average fuel assembly power obtained by zeroth and first order solutions with Galerkin and Moments weightings in comparison to reference solution. According to results, the higher and sufficient accuracy of first order approach is found relative to zeroth order solution, particularly, the first order solution with Moments weighting. In addition, Table 4 gives the suitable runtime of ACNECH for zeroth and first order solutions using coarse meshes, one node per a fuel assembly. As we know, the relative short time of neutronic solver is very important for the nuclear problems in which the solution process is very consuming such as core pattern optimization. Therefore, using ACNECH in the loading pattern optimization (LPO) is purposive, Poursalehi et al. (2013b). However, the fuel assembly power distribution of VVER-440-2D obtained by ACNECH and reference are given in Fig. 6. The major power errors in periphery fuel assemblies and also the control rod fuel assemblies neighboring are present in Fig. 6 that these errors are created because of high flux gradient existence in these regions.
181
4.3. Benchmark problem (C): Three-dimensional VVER-440 Last test case is based on the same core for problem (B) in which the original two-dimensional core is extended to a three-dimensional core with 250.0 cm high. A reflector 25.0 cm thick is added to the top and bottom of the core. All control rods are halfway inserted in the core with the material number 4 and for another half, with the material number 2, as shown in Fig. 7. Vacuum boundary conditions are applied to the entire reflector outside boundaries. Also, the reference solution was extrapolated from DIF3D-FD runs with 216 and 294 triangle/hexagon subdivisions and 2.5 cm axial mesh spacing, Chao and Shatilla (1995). ACNECH has been ran for the three-dimensional VVER-440 problem using coarse nodes i.e. one node per a fuel assembly in radial shape and axial sections with 25.0 cm height. Table 5 gives the numerical results of ACNECH for 3D VVER-440 test case. The zeroth order solution in both radial and axial directions has been exploited and Table 5 gives the results. According to Table 5, implementing the first order expansion in the radial direction and second order accuracy in the axial direction, the accuracy of results has been improved three times relative to zeroth order solution by decreasing the maximum power error from 15.4% to 4.9%. Of course, similar to previous problems, the first order accuracy using Moments weighting has the better results in comparison to Galerkin weighting. Furthermore, Table 5 indicates the relatively appropriate execution time of ACNECH for a 3D problem due to coarse mesh calculation. At last, Table 6 presents the three-dimensional node-wise powers obtained by ACNECH based on the assembly identification numbers as shown in Fig. 8. 5. Conclusions In this work, we developed a high order Average Current Nodal Expansion Code, ACNECH, for hexagonal-z geometries using the first order accuracy in radial directions and the second order solution in axial direction exploiting an adopted iterative scheme. ACNECH has been tested against three popular benchmark problems. Obtained results present that using high order solution, the accuracy of results enhances greatly relative to zeroth order solution. Actually, results confirm that the proposed high order methodology can obtain the adequate accuracy for the neutronic treatment analysis of nuclear reactor core with fuel assemblies by hexagonal shape. Of course, the numerical errors can be reduced by increasing the order of solution, particularly in radial directions. But this accuracy is achieved in the suitable computational cost and time. Thus for some problems such as LPO, the obtained accuracy can be sufficient because the run time of neutronic solver is very significant in the LPO. In addition, we found the higher accuracy of results using the Moments weighting for obtaining the first order expansion coefficients in regard of Galerkin weighting. Totally, the numerical results demonstrate the ability of ACNECH to calculate the stationary neutronic parameters of reactor core with hexagonal assembly shape. Furthermore, ACNECH has the capability to use in the LPO problem and can be extended for solving the time dependent nuclear problems in order to simulate the realistic reactor core neutronic treatments. References Chao, Y.A., Shatilla, Y.A., 1995. Conformal mapping and hexagonal nodal methods— II: Implementation in the ANC-H code. Nucl. Sci. Eng. 121, 210–225. Chao, Y.A., Tsolfanidis, N., 1995. Conformal mapping and hexagonal nodal methods—I: Mathematical foundation. Nucl. Sci. Eng. 121, 202–209. 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 multi-group problems in hexagonal-z geometry. Nucl. Sci. Eng. 126, 35.
182
N. Poursalehi et al. / Annals of Nuclear Energy 68 (2014) 172–182
Cho, N.Z., Lee, J., Kim, D.S., Yang, C.Y., Jo, J.C., 2005. Improved Features and Verification of an AFEN Method Code in Hexagonal-z 3-D Geometry for Neutron Diffusion Calculation submitted to American Nuclear Society Annual Meeting, San Diego, CA. Finnemann, H., 1975. A consistent nodal method for the analysis of space-time effects in large LWRs. In: Proc. of the Joint NEACRP/CSNI Specialists Meeting on New Developments in Three-Dimensional Neutron Kinetics and Review of Kinetics Benchmark Calculations. MRR 145, January 22–24, Munich, p. 131. Finnemann, H., Bennewitz, F., Wagner, M.R., 1977. Interface current techniques for multidimensional reactor calculations. Atomkernenergie 30, 123–128. Kuo, W.S., 1993. The General Evaluation of the Nodal Synthesis Method in Nuclear Reactor Transient. Submitted to the Department of Nuclear Engineering in partial fulfillment of the requirements for the degree of Doctor of Philosophy at The Massachusetts Institute of Technology. Lawrence, R.D., 1979. A nodal Green’s Function Method for Multidimensional Neutron Diffusion Calculations. Thesis, Nuclear Engineering Program, Univ. of Illinois, Urbana. Noh, J.M., Cho, N.Z., 1993. A new diffusion nodal method based on analytic basis function expansion. Trans. Am. Nucl. Soc. 69, 462–463. Noh, J.M., Cho, N.Z., 1994. A new approach of analytic basis function expansion to neutron diffusion calculation. Nucl. Sci. Eng. 116, 165. Poursalehi, N., Zolfaghari, A., Minuchehr, A., 2012. Performance comparison of zeroth order nodal expansion methods in 3D rectangular geometry. Nucl. Eng. Des. 252, 248–266. Poursalehi, N., Zolfaghari, A., Minuchehr, A., 2013a. Development of a high order and multi-dimensional nodal code, ACNEC3D, for reactor core analysis. Ann. Nucl. Energy 55, 211–224.
Poursalehi, N., Zolfaghari, A., Minuchehr, A., Valavi, K., 2013b. Self-adaptive global best harmony search algorithm applied to reactor core fuel management optimization. Ann. Nucl. Energy 62, 86–102. Putney, J.M., 1984. Nodal Methods for Solving the Diffusion Equation for Fast Reactor Analysis. Ph.D. Thesis. University of London. Smith, K.S., 1979. An Analytic Nodal Method for Solving The Two-Group, Multidimensional, Static and Transient Neutron Diffusion Equations. Submitted in Partial Fulfillment of The Requirements for The Degrees Of Nuclear Engineer and Master of Science at The Massachusetts Institute of Technology. Smith, K.S., 1985. QPANDA: an advanced nodal method for LWR analyses. Trans. Am. Nucl. Soc. 50, 532. Wagner, M.R., Finnemann, K.H., Koebke, K., Winter, J., 1977. Validation of the nodal expansion method and the depletion program MEDIUM-2 by bechmark calculations and direct comparisons with experiments. Atomkernenergie 30, 129–135. 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. Zerkle, M.L., 1992. Development of a Polynomial Nodal Methods with Flux and Current Discontinuity Factors, Ph.D. Thesis, Department of Nuclear Engineering, Massachusetts Institute of Technology, Cambridge, MA. Zimin, V.G., Baturin, D.M., 2002. Polynomial nodal method for solving neutron diffusion equations in hexagonal-z geometry. Ann. Nucl. Energy 29, 1105– 1117.