# Development and Validation of a Three-Dimensional Diffusion Code Based on a High Order Nodal Expansion Method for Hexagonal-z Geometry

Authors: Daogang Lu; Chao Guo
Journal: Science and Technology of Nuclear Installations
 (2016)

---

## Abstract

A three-dimensional, multigroup, diffusion code based on a high order nodal expansion method for hexagonal-z geometry (HNHEX) was developed to perform the neutronic analysis of hexagonal-z geometry. In this method, one-dimensional radial and axial spatially flux of each node and energy group are defined as quadratic polynomial expansion and four-order polynomial expansion, respectively. The approximations for one-dimensional radial and axial spatially flux both have second-order accuracy. Moment weighting is used to obtain high order expansion coefficients of the polynomials of one-dimensional radial and axial spatially flux. The partially integrated radial and axial leakages are both approximated by the quadratic polynomial. The coarse-mesh rebalance method with the asymptotic source extrapolation is applied to accelerate the calculation. This code is used for calculation of effective multiplication factor, neutron flux distribution, and power distribution. The numerical calculation in this paper for three-dimensional SNR and VVER 440 benchmark problems demonstrates the accuracy of the code. In addition, the results show that the accuracy of the code is improved by applying quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux in comparison to flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux.


---

## Body

## 1. Introduction

It has been nearly 30 years since the initial implementation of finite-difference method (FDM) [1] and finite element method (FEM) [2] in computer codes designed to solve the few-group neutron diffusion equations in more than one spatial dimension. A large number of mesh points were required for solving few-group neutron diffusion equations by traditional finite-difference techniques in order to represent accurately the spatial variation of the neutron flux. Due to the large number of unknowns involved, these calculations were very expensive, particularly for fuel management studies which require repeated solution of the diffusion equation. In order to solve this problem, the development of nodal methods had enabled computing times for LWR analysis to be significantly reduced. The nodal methods consisted of two main approaches which are polynomial nodal method [3] and Analytic Nodal Method (ANM) [4]. Two nodal methods have one feature in common which is that they are based on approximations to one-dimensional equations derived by integrating the three-dimensional equation over the two directions transverse to each coordinate axis. The essential difference of the two methods suggests classification of the methods developed for the solution on the basis of whether information obtained from an analytic solution of the diffusion equation within the node is incorporated into the numerical scheme. In the polynomial nodal method, the basic scheme is that the one-dimensional fluxes are approximated by a polynomial without the use of analytic information. A well-known example of the polynomial methods is the nodal expansion method (NEM) [5, 6]. Its basis is the solution of equations for the average flux on a coarse spatial mesh, combined with a high order local polynomial flux expansion, which is used to determine a relationship between the mean interface partial currents and the mesh average fluxes. Accuracies of both ANM and NEM are comparable. However, the ANM method, due to its complexity, is restricted to not more than two energy groups. NEM is an example of first-generation polynomial nodal method. One of its main advantages is that there is no restriction on number of energy groups.Recently, the capacity to compute accurate numerical solutions to the neutron diffusion equation in hexagonal-z geometry is required for the physics and safety analysis of hexagonal fuel assembly (FA) designs. However, the NEM, which is initially developed for LWR, is formulated in rectangular geometry. The analysis of hexagonal-z geometry requires an extension to hexagonal geometry of the transverse integration procedure widely used in the development of Cartesian-geometry nodal schemes. The success of these Cartesian-geometry nodal methods has prompted the development of analogous techniques [7, 8] for fast reactor calculations in hexagonal geometry. Extension of AFEN for hexagonal-z geometry was also employed by Cho et al. [9]. Besides, a mathematical basis for an extension of nodal methods to hexagonal geometry using conformal mapping of the hexagon into a rectangle has been provided by Chao and Shatilla [10]. In the traditional nodal expansion method for hexagonal-z geometry [11], partially integrated axial flux is approximated by quadratic polynomial and flat approximation is used for partially integrated axial leakage. The accuracy of this traditional method can be improved by using high order approximation.In this paper, we developed a 3D nodal code based on a high order nodal expansion method for hexagonal-z geometry (HNHEX). This method exploits quadratic polynomial expansion for partially integrated radial flux using which is based on the NEM [11] and four-order polynomial expansion for partially integrated axial flux. Moment weighting is used for radial direction and axial direction to obtain high order expansion coefficients of the polynomials. The partially integrated radial and axial leakages are both approximated by the quadratic polynomial. Numerical results of HNHEX demonstrate the efficiency and good accuracy of method.The organization of this paper is as follows. Section2 presents the summary of high order nodal expansion method for hexagonal-z geometry. Section 3 presents calculation of the leakage moments. In the next section, the proposed iterative solution algorithm is given and the results of HNHEX against three-dimensional SNR and VVER 440 benchmarks are presented. Section 6 presents the conclusion.


## 2. A High Order Nodal Expansion Method for Hexagonal-z Geometry

In this section, we present a high order nodal expansion method for hexagonal-z geometry. The multigroup neutron diffusion equation is solved by using a nodal scheme with one node per hexagonal assembly. One-dimensional radial spatially flux of each node and energy group is defined as the quadratic polynomial expansion. One-dimensional axial spatially flux of each node and energy group is defined as the four-order polynomial expansion. The polynomial for the radial spatially flux is equivalent to a quadratic approximation and has second-order accuracy. The polynomial for the axial spatially flux is a four-order approximation and has second-order accuracy. The partially integrated radial and axial leakages are both approximated by the quadratic polynomial. The final equations, which are cast in response matrix form, involve spatial moments of the node flux distribution plus surface-averaged partial currents across the faces of the node. The coarse-mesh rebalance method with the asymptotic source extrapolation is applied to accelerate the calculation.


### 2.1. Nodal Coordinate Systems

The hexagonal-z geometry only has difference with rectangular geometry in the radial direction. There is no difference in the axial direction. Thus, a combination of a hexagonal geometry analysis in the radial direction with a 1D rectangular geometry analysis in the axial direction is required to the derivate of a nodal expansion method for hexagonal-z geometry.The hexagonal axes which consist ofx-direction,u-direction, andv-direction are employed in the hexagonal geometry as shown in Figure 1. Theu-direction is rotated 60° counterclockwise with respect to the x-axis and thev-direction is rotated 120° from the x-axis. Using the orientation shown in Figure 2, with the origin (in local coordinates) taken as the center of the hexagon, a homogeneous node volume ofkth node can be defined as shown in (1) in the local coordinate:(1)Vk:x,y,zx∈-h2,h2,y∈-yx,yx,z∈-Δzk2,Δzk2,where(2)yx=13h-x.h is the lattice pitch and the x-axis is taken as perpendicular to one pair of opposite faces of the hexagon. Analogous partially integrated fluxes and currents are defined for the two opposite faces.Figure 1
Hexagonal coordinate system.Figure 2
A three-dimensional hexagonal-z geometry.


### 2.2. High Order Nodal Expansion Method for Hexagonal-z Geometry

In this method, reactor is divided intok number of homogeneous hexagonal nodes. Cross sections within the nodes are assumed to be constant. Let us consider a homogeneous node k shown in Figure 2, the multigroup neutron diffusion equation for kth node, and gth group can then be written in the form shown in(3)-∇Dgk∇ϕgkx,y,z+Σgr,kϕgkx,y,z=Qgkx,y,z,k=1,…,K,g=1,…,G,where(4)Qgkx,y,z=χgkeff∑g′=1GυΣg′f,kϕg′kx,y,z+∑g′=1g′<gGΣg′→gs,kϕ_g′kx,y,z.The solution of (3) requires additional relationships between the surface-averaged leakages and the nodal fluxes in thekth node and its immediate neighbors. The additional relationships are given by approximations to one-dimensional equations derived by integrating the three-dimensional equation over the two directions transverse to each coordinate axis. The one-dimensional radial equation can be obtained by integrating the three-dimensional equation overy andz directions and take the form shown in(5)ddwJgwkw+Σgr,kϕgwkw=Qgwkw-23Jgwykw,yw-Jgwykw,-yw-Lgzkw,w=x,u,v,where(6)ϕgwkw=∫-Δzk/2Δzk/2dz∫-ywywdyϕgkx,y,z,Jgwkw=∫-Δzk/2Δzk/2dz∫-ywywdy-Dgk∂∂wϕgkx,y,z,Qgwkw=∫-Δzk/2Δzk/2dz∫-ywywdyQgkx,y,z,Lgzkx,y=∫-Δzk/2Δzk/2dz-Dgk∂2∂z2ϕgkx,y,z,w=x,u,v.The one-dimensional axial equation can be obtained by integrating the three-dimensional equation overx andy directions and take the form shown in(7)ddzJgzkz+Σgr,kϕgzkz=Qgzkz-Lgxykz,where(8)ϕgzkz=∫-h/2h/2dx∫-yxyxdyϕgkx,y,z,Jgzkz=∫-h/2h/2dx∫-yxyxdy-Dgk∂∂zϕgkx,y,z,Qgzkz=∫-h/2h/2dx∫-yxyxdyQgkx,y,z,Lgxykz=∫-h/2h/2dx∫-yxyxdy-Dgk∂2∂x2+∂2∂y2ϕgkx,y,z.In order to solve (5) and (7), the partially integrated radial and axial fluxes are defined as quadratic polynomial expansion and four-order polynomial expansion shown in (9) and (11), respectively. The approximation defined in (9) for the partially integrated radial fluxes has second-order accuracy due to the fact that this latter procedure is equivalent to enforcing a neutron balance over each of the half-nodes. The approximation defined in (11) for the axial direction has second-order accuracy. The approximation of the partially integrated radial fluxes takes the form shown in(9)ϕ¯gwkw=ϕgwkw2yw=ϕ¯gk+agw1kf1w+agw2kf2w+agw3kf3w+agw4kf4w,w=x,u,v,where(10)ϕ¯gk=1Vk∫-Δzk/2Δzk/2dz∫-h/2h/2dx∫-yxyxdyϕgkx,y,z,f1w=wh,f2w=3613wh2-526,f3w=1013wh2-12wh+352,f4w=whwh-12,w=x,u,v.f1w and f2w are obtained by using orthogonal characteristics of fnw polynomials. Higher order approximations are obtained by first adding a basis function f3w which has a first derivative discontinuity at w=0 and then adding an additional basis function f4w which provides a quadratic approximation over the half-intervals -h/2<w<0 and 0<w<h/2.The approximation of the partially integrated axial fluxes takes the form shown in(11)ϕ¯gzkz=ΔzkϕgzkzVk=ϕ¯gk+agz1kg1z+agz2kg2z+agz3kg3z+agz4kg4z,where(12)g1z=zΔzk,g2z=zΔzk2-112,g3z=zΔzkzΔzk2-14,g4z=zΔzk2-14zΔzk2-120.In order to determine the zeroth expansion coefficients, the polynomials of (9) and (11) are fitted to the average flux in the node and the average fluxes on its surfaces. Thus, we impose the conditions shown in(13)ϕ¯gwkw=+h2=ϕ¯gk+agw1kf1h2+agw2kf2h2+agw3kf3h2+agw4kf4h2=ϕ¯gw+k⟹agw1k+agw2k=2ϕ¯gw+k-2ϕ¯gk,w=x,u,v,ϕ¯gwkw=-h2=ϕ¯gk+agw1kf1-h2+agw2kf2-h2+agw3kf3-h2+agw4kf4-h2=ϕ¯gw-k⟹agw1k-agw2k=2ϕ¯gk-2ϕ¯gw-k,w=x,u,v,ϕ¯gzkz=+Δzk2=ϕ¯gk+agz1kg1Δzk2+agz2kg2Δzk2+agz3kg3Δzk2+agz4kg4Δzk2=ϕ¯gz+k⟹agz1k+agz2k=2ϕ¯gz+k-2ϕ¯gk,ϕ¯gzkz=-Δzk2=ϕ¯gk+agz1kg1-Δzk2+agz2kg2-Δzk2+agz3kg3-Δzk2+agz4kg4-Δzk2=ϕ¯gz-k⟹agz1k-agz2k=2ϕ¯gk-2ϕ¯gz-k.We solve (13) to get the zeroth expansion coefficients shown in(14)agw1k=ϕ¯gw+k-ϕ¯gw-k,agw2k=ϕ¯gw+k+ϕ¯gw-k-2ϕ¯gk,agz1k=ϕ¯gz+k-ϕ¯gz-k,agz2k=ϕ¯gz+k+ϕ¯gz-k-2ϕ¯gk,w=x,u,v.The coefficientsagw3k(w=x,u,v) are determined by requiring they-integrated net current Jgwkw(w=x,u,v) to be continuous at w=0 and the imposed condition is shown in(15)limε→0⁡-Dgkddxϕgxkxx=-εx=ε=-Dgk2yxagx3kf3′x+2y′x2yxϕgxkxx=-εx=ε=2Dgk3agx3k+2ϕ¯gxkxx=0.Apply the two-step leakage approximation and the final form ofagw3k(w=x,u,v) is given by(16)agx3k=-26189hDgkJ¯gu+k-J¯gu-k+J¯gv+k-J¯gv-k-821ϕ¯gx+k+ϕ¯gx-k-2ϕ¯gk,agu3k=-26189hDgkJ¯gx-k-J¯gx+k+J¯gv+k-J¯gv-k-821ϕ¯gu+k+ϕ¯gu-k-2ϕ¯gk,agv3k=-26189hDgkJ¯gx-k-J¯gx+k+J¯gu-k-J¯gu+k-821ϕ¯gv+k+ϕ¯gv-k-2ϕ¯gk,where(17)J¯gw±k=J¯gwk±h2,w=x,u,v.Higher order coefficientsagw4k(w=x,u,v) are determined by applying a moment weighting scheme to the one-dimensional balance equation (5). The weighted residual equation is obtained by weighting equation (5) with a weight function hiw and require the result to be zero when integrated over the interval w∈-h/2,h/2(w=x,u,v) as shown in(18)1Vk∫-h/2+h/2hiwddwJgwkw+Σgr,kϕgwkw-Qgwkw+Lgzwkw+23Jgkw,yw-Jgkw,-ywdw=0,w=x,u,v.For moment weighting,h1w=f1w=w/h(w=x,u,v); integrate and rewrite (18) to get flux moment equation in three radial directions shown in(19)Σgr,k+32Dgkh2ϕgx1k+23hΓ¯gxk+Γ¯guk-Γ¯gvk-40Dgk9h2agx1k=Qgx1k-1ΔzkLgzx1k,Σgr,k+32Dgkh2ϕgu1k+23hΓ¯gxk+Γ¯guk+Γ¯gvk-40Dgk9h2agu1k=Qgu1k-1ΔzkLgzu1k,Σgr,k+32Dgkh2ϕgv1k+23h-Γ¯gxk+Γ¯guk+Γ¯gvk-40Dgk9h2agv1k=Qgv1k-1ΔzkLgzv1k,where(20)ϕgw1k=1Vk∫-h/2h/2dwf1wϕgwkw,Qgw1k=1Vk∫-h/2h/2dwf1wQgwkw,Lgzw1k=ΔzkVk∫-h/2h/2dwf1w∫-ywywdyLgzkw,y,w=x,u,v,(21)Γ¯gxk=J¯gx+k+J¯gx-k,Γ¯guk=J¯gu+k+J¯gu-k,Γ¯gvk=J¯gv+k+J¯gv-k.Now substituteϕgwkw in ϕgw1k with (9) and then integrate to get agw4k as shown in(22)ϕgw1k=1Vk∫-h/2h/2dxf1wϕgwkw=29agw1k-124agw4k⟹agw4k=-24ϕgw1k+163agw1k=-24ϕgw1k+163ϕ¯gw+k-ϕ¯gw-k,w=x,u,v.Higher order coefficientsagz3k and agz4k are determined by applying a moment weighting scheme to the one-dimensional balance equation (7). The weighted residual equation is obtained by weighting equation (7) with a weight function hziz(i=1,2) and requiring the result to be zero when integrated over the interval z∈-Δzk/2,Δzk/2(23)1Vk∫-Δzk/2Δzk/2hzizddzJgzkz+Σgr,kϕgzkz-Qgzkz+Lgxykzdz=0,i=1,2.For moment weighting,i=1, hz1z=g1z=z/Δzk, integrate and rewrite (23) to get one-order flux moment equation in axial direction shown in(24)Σgr,kϕgz1k+12J¯gz+k+J¯gz-kΔzk+DgkΔzk2agz1k=Qgz1k-23hLgxyz1k,where(25)ϕgz1k=1Vk∫-Δzk/2Δzk/2dzg1zϕgzkz,Qgz1k=1Vk∫-Δzk/2Δzk/2dzg1zQgzkz,Lgxyz1k=3h21Vk∫-Δzk/2Δzk/2dzg1zLgxykz,J¯gz±k=J¯gzk±Δzk2.For moment weighting,i=2, hz2z=g2z=z/Δzk2-1/12, integrate and rewrite (23) to get second-order flux moment equation in axial direction shown in(26)Σgr,kϕgz2k+16ΔzkJ¯gz+k-J¯gz-k+Dgk3Δzk2agz2k=Qgz2k-23hLgxyz2k,where(27)ϕgz2k=1Vk∫-Δzk/2Δzk/2dzg2zϕgzkz,Qgz2k=1Vk∫-Δzk/2Δzk/2dzg2zQgzkz,Lgxyz2k=3h21Vk∫-Δzk/2Δzk/2dzg2zLgxykz.Substituteϕgzkz into ϕgz1k and ϕgz2k with (11) and then integrate to get (28a) and (28b):(28a)ϕgz1k=1Vk∫-Δzk/2Δzk/2g1zϕgzkzdz=1Vk∫-Δzk/2Δzk/2zΔzkVkΔzkϕ¯gk+agz1kfz1z+agz2kfz2z+agz3kfz3z+agz4kfz4zdz=112agz1k-1120agz3k,(28b)ϕgz2k=1Vk∫-Δzk/2Δzk/2g2zϕgzkzdz=1Vk∫-Δzk/2Δzk/2zΔzk2-112ϕ¯gk+agz1kfz1z+agz2kfz2z+agz3kfz3z+agz4kfz4zdz=1180agz2k-12100agz4k.Thus,agz3k and agz4k are shown in(29)agz3k=-120ϕgz1k+10agz1k=-120ϕgz1k+10ϕ¯gz+k-ϕ¯gz-k,agz4k=-2100ϕgz2k+353agz2k=-2100ϕgz2k+353ϕ¯gz+k+ϕ¯gz-k-2ϕ¯gk.The 3D multigroup nodal balance equation is obtained by integrating over a homogeneous node volumeVk(30)23hJ¯gx+k-J¯gx-k+J¯gu+k-J¯gu-k+J¯gv+k-J¯gv-k+1ΔzkJ¯gz+k-J¯gz-k+Σgr,kϕ¯gk=Q¯gk,where(31)Q¯gk=1Vk∫-Δzk/2Δzk/2dz∫-h/2h/2dx∫-ysxysxdyQgkx,y,z.ϕ¯gk is derived from (30) and shown in (32). ϕgx1k, ϕgu1k, and ϕgv1k are derived from (19) and shown in (33), (34), and (35). ϕgz1k and ϕgz2k are derived from (24) and (26) and shown in (36) and (37):(32)ϕ¯gk=1Σgr,kQ¯gk-23hΣgr,kJ¯gx+k-J¯gx-k+J¯gu+k-J¯gu-k+J¯gv+k-J¯gv-k-1ΔzkΣgr,kJ¯gz+k-J¯gz-k,(33)ϕgx1k=1Σgr,k+32Dgk/h2Qgx1k-1ΔzkLgzx1k-23hΣgr,k+32Dgk/h2Γ¯gxk+Γ¯guk-Γ¯gvk+40Dgk9h2Σgr,k+32Dgk/h2agx1k,(34)ϕgu1k=1Σgr,k+32Dgk/h2Qgu1k-1ΔzkLgzu1k-23hΣgr,k+32Dgk/h2Γ¯gxk+Γ¯guk+Γ¯gvk+40Dgk9h2Σgr,k+32Dgk/h2agu1k,(35)ϕgv1k=1Σgr,k+32Dgk/h2Qgv1k-1ΔzkLgzv1k-23hΣgr,k+32Dgk/h2-Γ¯gxk+Γ¯guk+Γ¯gvk+40Dgk9h2Σgr,k+32Dgk/h2agv1k,(36)ϕgz1k=1Σgr,kQgz1k-23hLgxyz1k-12ΔzkΣgr,kJ¯gz+k+J¯gz-k-DgkΔzk2Σgr,kagz1k,(37)ϕgz2k=1Σgr,kQgz2k-23hLgxyz2k-16ΔzkΣgr,kJ¯gz+k-J¯gz-k-Dgk3Δzk2Σgr,kagz2k.The partial net current and one-dimensional average spatially flux can be written in the form of outgoing and incoming partial currents as shown in(38)±J¯gw±k=J¯gw±out,k-J¯gw±in,k,ϕ¯gw±k=2J¯gw±out,k+J¯gw±in,k,w=x,u,v,z.By using (9) and (38), outgoing partial currents through radial-directed surfaces ofkth node in terms of flux expansion coefficients are shown in(39)J¯gw+out,k=J¯gw+k+J¯gw+in,k=-Dgkddwϕ¯gwkw+Dgky′w2ywEgwkwx=h/2+J¯gw+in,k=-Dgkhagw1k+3613agw2k+726agw3k+12agw4k+Egw+k+J¯gw+in,k,w=x,u,v,J¯gw-out,k=-J¯gw-k+J¯gw-in,k=--Dgkddwϕ¯gwkw+Dgky′w2ywEgwkwx=-h/2+J¯gw-in,k=Dgkhagw1k-3613agw2k-726agw3k+12agw4k-Egw-k+J¯gw-in,k,w=x,u,v,where(40)Egx+k=-1285hDgk17J¯gu+k-J¯gv-k-2J¯gv+k-J¯gu-k-11235179ϕ¯gx+k+49ϕ¯gx-k-228ϕ¯gk+1130agx3k,Egu+k=-1285hDgk17J¯gv+k-J¯gx+k-2J¯gx-k-J¯gv-k-11235179ϕ¯gu+k+49ϕ¯gu-k-228ϕ¯gk+1130agu3k,Egv+k=-1285hDgk17J¯gx-k-J¯gu+k-2J¯gu-k-J¯gx+k-11235179ϕ¯gv+k+49ϕ¯gv-k-228ϕ¯gk+1130agv3k,Egx-k=-1285hDgk17J¯gu+k-J¯gv-k-2J¯gv+k-J¯gu-k-1123549ϕ¯gx+k+179ϕ¯gx-k-228ϕ¯gk+1130agx3k,Egu-k=-1285hDgk17J¯gv+k-J¯gx+k-2J¯gx-k-J¯gv-k-1123549ϕ¯gu+k+179ϕ¯gu-k-228ϕ¯gk+1130agu3k,Egv-k=-1285hDgk17J¯gx-k-J¯gu+k-2J¯gu-k-J¯gx+k-1123549ϕ¯gv+k+179ϕ¯gv-k-228ϕ¯gk+1130agv3k.By using (11) and (38), outgoing partial currents through axial-directed surfaces ofkth node in terms of flux expansion coefficients are shown in(41)J¯gz+out,k=J¯gz+k+J¯gz+in,k=-Dgkddzϕ¯gzkzz=Δzk/2+J¯gz+in,k=-DgkΔzkVkddzϕ^gzkzz=Δzk/2+J¯gz+in,k=-DgkΔzkagz1k+agz2k+12agz3k+1740agz4k+J¯gz+in,k,J¯gz-out,k=-J¯gz-k+J¯gz-in,k=--Dgkddzϕ¯gzkzz=-Δzk/2+J¯gz-in,k=--DgkΔzkVkddzϕ^gzkzz=-Δzk/2+J¯gz-in,k=DgkΔzkagz1k-agz2k+12agz3k-1740agz4k+J¯gz-in,k.We use (14), (16), (22), and (29) for expansion coefficients along with (32), (33), (34), (35), (36), and (37) for node averaged flux and its moments and (49), (50), (51), (56), and (58) for node leakage moments. One-dimensional average spatially flux is expressed by partial currents with (38). Eight outgoing partial currents can be expressed in terms of eight incoming partial currents, source, and its moments and leakage moments.(42)
a
11
a
12
a
13
a
14
a
15
a
16
a
17
a
18
a
21
a
22
a
23
a
24
a
25
a
26
a
27
a
28
a
31
a
32
a
33
a
34
a
35
a
36
a
37
a
38
a
41
a
42
a
43
a
44
a
45
a
46
a
47
a
48
a
51
a
52
a
53
a
54
a
55
a
56
a
57
a
58
a
61
a
62
a
63
a
64
a
65
a
66
a
67
a
68
a
71
a
72
a
73
a
74
a
75
a
76
a
77
a
78
a
81
a
82
a
83
a
84
a
85
a
86
a
87
a
88
J
¯
g
x
+
out
,
k
J
¯
g
u
+
out
,
k
J
¯
g
v
+
out
,
k
J
¯
g
x
-
out
,
k
J
¯
g
u
-
out
,
k
J
¯
g
v
-
out
,
k
J
¯
g
z
+
out
,
k
J
¯
g
z
-
out
,
k
=
b
11
b
12
0
0
0
0
b
21
0
b
23
0
0
0
b
31
0
0
b
34
0
0
b
41
b
42
0
0
0
0
b
51
0
b
53
0
0
0
b
61
0
0
b
64
0
0
b
71
0
0
0
b
75
b
76
b
81
0
0
0
b
85
b
86
Q
¯
g
k
Q
g
x
1
k
-
1
Δ
z
k
L
g
z
x
1
k
Q
g
u
1
k
-
1
Δ
z
k
L
g
z
u
1
k
Q
g
v
1
k
-
1
Δ
z
k
L
g
z
v
1
k
Q
g
z
1
k
-
2
3
h
L
g
x
y
z
1
k
Q
g
z
2
k
-
2
3
h
L
g
x
y
z
2
k
+
c
11
c
12
c
13
c
14
c
15
c
16
c
17
c
18
c
21
c
22
c
23
c
24
c
25
c
26
c
27
c
28
c
31
c
32
c
33
c
34
c
35
c
36
c
37
c
38
c
41
c
42
c
43
c
44
c
45
c
46
c
47
c
48
c
51
c
52
c
53
c
54
c
55
c
56
c
57
c
58
c
61
c
62
c
63
c
64
c
65
c
66
c
67
c
68
c
71
c
72
c
73
c
74
c
75
c
76
c
77
c
78
c
81
c
82
c
83
c
84
c
85
c
86
c
87
c
88
J
¯
g
x
+
in
,
k
J
¯
g
u
+
in
,
k
J
¯
g
v
+
in
,
k
J
¯
g
x
-
in
,
k
J
¯
g
u
-
in
,
k
J
¯
g
v
-
in
,
k
J
¯
g
z
+
in
,
k
J
¯
g
z
-
in
,
k
⟹
A
g
k
J
_
g
out
,
k
=
B
g
k
Q
_
g
k
-
L
_
g
k
+
C
g
k
J
_
g
in
,
k
⟹
J
_
g
out
,
k
=
A
g
k
-
1
B
g
k
Q
_
g
k
-
L
_
g
k
+
A
g
k
-
1
C
g
k
J
_
g
in
,
k
⟹
J
_
g
out
,
k
=
P
g
k
Q
_
g
k
-
L
_
g
k
+
R
g
k
J
_
g
in
,
k
,where(43)Q_gk=col⁡Q¯gk,Qgx1k,Qgu1k,Qgv1k,Qgz1k,Qgz2k,L_gk=col⁡0,1ΔzkLgzx1k,1ΔzkLgzu1k,1ΔzkLgzv1k,23hLgxyz1k,23hLgxyz2k.Matrix elements (aij,bij,cij) are given below. These elements are characterized by node dimensions and homogenized material of the node.(44)a11=a22=a33=a44=a55=a66=1+2Dgkh468133+83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk,a12=a23=a34=a45=a56=a21=a32=a43=a54=a65=a16=a61=24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk-13133,a13=a24=a35=a46=a31=a42=a53=a64=a15=a26=a51=a62=24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk-6133,a14=a25=a36=a41=a52=a63=2Dgkh216133-83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk,a17=a27=a37=a47=a57=a67=a18=a28=a38=a48=a58=a68=36Dgk7ΔzkhΣgr,k,a71=a72=a73=a74=a75=a76=a81=a82=a83=a84=a85=a86=2Dgk3ΔzkhΣgr,k14312+595DgkΔzk2Σgr,k,a77=a88=DgkΔzk28712+715DgkΔzk2Σgr,k+5723ΔzkΣgr,k+595DgkΔzk3Σgr,k2+1,a78=a87=DgkΔzk-112+475DgkΔzk2Σgr,k+3923ΔzkΣgr,k+595DgkΔzk3Σgr,k2,b11=b21=b31=b41=b51=b61=36Dgk7hΣgr,k,b12=b23=b34=12hDgkh2Σgr,k+32Dgk,b42=b53=b64=-12hDgkh2Σgr,k+32Dgk,b71=b81=DgkΔzk14312Σgr,k+595DgkΔzk2Σgr,k2,b75=60DgkΔzkΣgr,k,b85=-60DgkΔzkΣgr,k,b76=b86=1785Dgk2ΔzkΣgr,k,c11=c22=c33=c44=c55=c66=1-2Dgkh468133+83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk,c12=c23=c34=c45=c56=c21=c32=c43=c54=c65=c16=c61=24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk-13133,c13=c24=c35=c46=c31=c42=c53=c64=c15=c26=c51=c62=24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk-6133,c14=c25=c36=c41=c52=c63=-2Dgkh216133-83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk,c17=c27=c37=c47=c57=c67=c18=c28=c38=c48=c58=c68=36Dgk7ΔzkhΣgr,k,c71=c72=c73=c74=c75=c76=c81=c82=c83=c84=c85=c86=2Dgk3ΔzkhΣgr,k14312+595DgkΔzk2Σgr,k,c77=c88=DgkΔzk5723ΔzkΣgr,k+595Dgk3Δzk3Σgr,k2-715DgkΔzk2Σgr,k-28712+1,c78=c87=DgkΔzk3923ΔzkΣgr,k+595DgkΔzk3Σgr,k2+112-475DgkΔzk2Σgr,k.


## 2.1. Nodal Coordinate Systems

The hexagonal-z geometry only has difference with rectangular geometry in the radial direction. There is no difference in the axial direction. Thus, a combination of a hexagonal geometry analysis in the radial direction with a 1D rectangular geometry analysis in the axial direction is required to the derivate of a nodal expansion method for hexagonal-z geometry.The hexagonal axes which consist ofx-direction,u-direction, andv-direction are employed in the hexagonal geometry as shown in Figure 1. Theu-direction is rotated 60° counterclockwise with respect to the x-axis and thev-direction is rotated 120° from the x-axis. Using the orientation shown in Figure 2, with the origin (in local coordinates) taken as the center of the hexagon, a homogeneous node volume ofkth node can be defined as shown in (1) in the local coordinate:(1)Vk:x,y,zx∈-h2,h2,y∈-yx,yx,z∈-Δzk2,Δzk2,where(2)yx=13h-x.h is the lattice pitch and the x-axis is taken as perpendicular to one pair of opposite faces of the hexagon. Analogous partially integrated fluxes and currents are defined for the two opposite faces.Figure 1
Hexagonal coordinate system.Figure 2
A three-dimensional hexagonal-z geometry.


## 2.2. High Order Nodal Expansion Method for Hexagonal-z Geometry

In this method, reactor is divided intok number of homogeneous hexagonal nodes. Cross sections within the nodes are assumed to be constant. Let us consider a homogeneous node k shown in Figure 2, the multigroup neutron diffusion equation for kth node, and gth group can then be written in the form shown in(3)-∇Dgk∇ϕgkx,y,z+Σgr,kϕgkx,y,z=Qgkx,y,z,k=1,…,K,g=1,…,G,where(4)Qgkx,y,z=χgkeff∑g′=1GυΣg′f,kϕg′kx,y,z+∑g′=1g′<gGΣg′→gs,kϕ_g′kx,y,z.The solution of (3) requires additional relationships between the surface-averaged leakages and the nodal fluxes in thekth node and its immediate neighbors. The additional relationships are given by approximations to one-dimensional equations derived by integrating the three-dimensional equation over the two directions transverse to each coordinate axis. The one-dimensional radial equation can be obtained by integrating the three-dimensional equation overy andz directions and take the form shown in(5)ddwJgwkw+Σgr,kϕgwkw=Qgwkw-23Jgwykw,yw-Jgwykw,-yw-Lgzkw,w=x,u,v,where(6)ϕgwkw=∫-Δzk/2Δzk/2dz∫-ywywdyϕgkx,y,z,Jgwkw=∫-Δzk/2Δzk/2dz∫-ywywdy-Dgk∂∂wϕgkx,y,z,Qgwkw=∫-Δzk/2Δzk/2dz∫-ywywdyQgkx,y,z,Lgzkx,y=∫-Δzk/2Δzk/2dz-Dgk∂2∂z2ϕgkx,y,z,w=x,u,v.The one-dimensional axial equation can be obtained by integrating the three-dimensional equation overx andy directions and take the form shown in(7)ddzJgzkz+Σgr,kϕgzkz=Qgzkz-Lgxykz,where(8)ϕgzkz=∫-h/2h/2dx∫-yxyxdyϕgkx,y,z,Jgzkz=∫-h/2h/2dx∫-yxyxdy-Dgk∂∂zϕgkx,y,z,Qgzkz=∫-h/2h/2dx∫-yxyxdyQgkx,y,z,Lgxykz=∫-h/2h/2dx∫-yxyxdy-Dgk∂2∂x2+∂2∂y2ϕgkx,y,z.In order to solve (5) and (7), the partially integrated radial and axial fluxes are defined as quadratic polynomial expansion and four-order polynomial expansion shown in (9) and (11), respectively. The approximation defined in (9) for the partially integrated radial fluxes has second-order accuracy due to the fact that this latter procedure is equivalent to enforcing a neutron balance over each of the half-nodes. The approximation defined in (11) for the axial direction has second-order accuracy. The approximation of the partially integrated radial fluxes takes the form shown in(9)ϕ¯gwkw=ϕgwkw2yw=ϕ¯gk+agw1kf1w+agw2kf2w+agw3kf3w+agw4kf4w,w=x,u,v,where(10)ϕ¯gk=1Vk∫-Δzk/2Δzk/2dz∫-h/2h/2dx∫-yxyxdyϕgkx,y,z,f1w=wh,f2w=3613wh2-526,f3w=1013wh2-12wh+352,f4w=whwh-12,w=x,u,v.f1w and f2w are obtained by using orthogonal characteristics of fnw polynomials. Higher order approximations are obtained by first adding a basis function f3w which has a first derivative discontinuity at w=0 and then adding an additional basis function f4w which provides a quadratic approximation over the half-intervals -h/2<w<0 and 0<w<h/2.The approximation of the partially integrated axial fluxes takes the form shown in(11)ϕ¯gzkz=ΔzkϕgzkzVk=ϕ¯gk+agz1kg1z+agz2kg2z+agz3kg3z+agz4kg4z,where(12)g1z=zΔzk,g2z=zΔzk2-112,g3z=zΔzkzΔzk2-14,g4z=zΔzk2-14zΔzk2-120.In order to determine the zeroth expansion coefficients, the polynomials of (9) and (11) are fitted to the average flux in the node and the average fluxes on its surfaces. Thus, we impose the conditions shown in(13)ϕ¯gwkw=+h2=ϕ¯gk+agw1kf1h2+agw2kf2h2+agw3kf3h2+agw4kf4h2=ϕ¯gw+k⟹agw1k+agw2k=2ϕ¯gw+k-2ϕ¯gk,w=x,u,v,ϕ¯gwkw=-h2=ϕ¯gk+agw1kf1-h2+agw2kf2-h2+agw3kf3-h2+agw4kf4-h2=ϕ¯gw-k⟹agw1k-agw2k=2ϕ¯gk-2ϕ¯gw-k,w=x,u,v,ϕ¯gzkz=+Δzk2=ϕ¯gk+agz1kg1Δzk2+agz2kg2Δzk2+agz3kg3Δzk2+agz4kg4Δzk2=ϕ¯gz+k⟹agz1k+agz2k=2ϕ¯gz+k-2ϕ¯gk,ϕ¯gzkz=-Δzk2=ϕ¯gk+agz1kg1-Δzk2+agz2kg2-Δzk2+agz3kg3-Δzk2+agz4kg4-Δzk2=ϕ¯gz-k⟹agz1k-agz2k=2ϕ¯gk-2ϕ¯gz-k.We solve (13) to get the zeroth expansion coefficients shown in(14)agw1k=ϕ¯gw+k-ϕ¯gw-k,agw2k=ϕ¯gw+k+ϕ¯gw-k-2ϕ¯gk,agz1k=ϕ¯gz+k-ϕ¯gz-k,agz2k=ϕ¯gz+k+ϕ¯gz-k-2ϕ¯gk,w=x,u,v.The coefficientsagw3k(w=x,u,v) are determined by requiring they-integrated net current Jgwkw(w=x,u,v) to be continuous at w=0 and the imposed condition is shown in(15)limε→0⁡-Dgkddxϕgxkxx=-εx=ε=-Dgk2yxagx3kf3′x+2y′x2yxϕgxkxx=-εx=ε=2Dgk3agx3k+2ϕ¯gxkxx=0.Apply the two-step leakage approximation and the final form ofagw3k(w=x,u,v) is given by(16)agx3k=-26189hDgkJ¯gu+k-J¯gu-k+J¯gv+k-J¯gv-k-821ϕ¯gx+k+ϕ¯gx-k-2ϕ¯gk,agu3k=-26189hDgkJ¯gx-k-J¯gx+k+J¯gv+k-J¯gv-k-821ϕ¯gu+k+ϕ¯gu-k-2ϕ¯gk,agv3k=-26189hDgkJ¯gx-k-J¯gx+k+J¯gu-k-J¯gu+k-821ϕ¯gv+k+ϕ¯gv-k-2ϕ¯gk,where(17)J¯gw±k=J¯gwk±h2,w=x,u,v.Higher order coefficientsagw4k(w=x,u,v) are determined by applying a moment weighting scheme to the one-dimensional balance equation (5). The weighted residual equation is obtained by weighting equation (5) with a weight function hiw and require the result to be zero when integrated over the interval w∈-h/2,h/2(w=x,u,v) as shown in(18)1Vk∫-h/2+h/2hiwddwJgwkw+Σgr,kϕgwkw-Qgwkw+Lgzwkw+23Jgkw,yw-Jgkw,-ywdw=0,w=x,u,v.For moment weighting,h1w=f1w=w/h(w=x,u,v); integrate and rewrite (18) to get flux moment equation in three radial directions shown in(19)Σgr,k+32Dgkh2ϕgx1k+23hΓ¯gxk+Γ¯guk-Γ¯gvk-40Dgk9h2agx1k=Qgx1k-1ΔzkLgzx1k,Σgr,k+32Dgkh2ϕgu1k+23hΓ¯gxk+Γ¯guk+Γ¯gvk-40Dgk9h2agu1k=Qgu1k-1ΔzkLgzu1k,Σgr,k+32Dgkh2ϕgv1k+23h-Γ¯gxk+Γ¯guk+Γ¯gvk-40Dgk9h2agv1k=Qgv1k-1ΔzkLgzv1k,where(20)ϕgw1k=1Vk∫-h/2h/2dwf1wϕgwkw,Qgw1k=1Vk∫-h/2h/2dwf1wQgwkw,Lgzw1k=ΔzkVk∫-h/2h/2dwf1w∫-ywywdyLgzkw,y,w=x,u,v,(21)Γ¯gxk=J¯gx+k+J¯gx-k,Γ¯guk=J¯gu+k+J¯gu-k,Γ¯gvk=J¯gv+k+J¯gv-k.Now substituteϕgwkw in ϕgw1k with (9) and then integrate to get agw4k as shown in(22)ϕgw1k=1Vk∫-h/2h/2dxf1wϕgwkw=29agw1k-124agw4k⟹agw4k=-24ϕgw1k+163agw1k=-24ϕgw1k+163ϕ¯gw+k-ϕ¯gw-k,w=x,u,v.Higher order coefficientsagz3k and agz4k are determined by applying a moment weighting scheme to the one-dimensional balance equation (7). The weighted residual equation is obtained by weighting equation (7) with a weight function hziz(i=1,2) and requiring the result to be zero when integrated over the interval z∈-Δzk/2,Δzk/2(23)1Vk∫-Δzk/2Δzk/2hzizddzJgzkz+Σgr,kϕgzkz-Qgzkz+Lgxykzdz=0,i=1,2.For moment weighting,i=1, hz1z=g1z=z/Δzk, integrate and rewrite (23) to get one-order flux moment equation in axial direction shown in(24)Σgr,kϕgz1k+12J¯gz+k+J¯gz-kΔzk+DgkΔzk2agz1k=Qgz1k-23hLgxyz1k,where(25)ϕgz1k=1Vk∫-Δzk/2Δzk/2dzg1zϕgzkz,Qgz1k=1Vk∫-Δzk/2Δzk/2dzg1zQgzkz,Lgxyz1k=3h21Vk∫-Δzk/2Δzk/2dzg1zLgxykz,J¯gz±k=J¯gzk±Δzk2.For moment weighting,i=2, hz2z=g2z=z/Δzk2-1/12, integrate and rewrite (23) to get second-order flux moment equation in axial direction shown in(26)Σgr,kϕgz2k+16ΔzkJ¯gz+k-J¯gz-k+Dgk3Δzk2agz2k=Qgz2k-23hLgxyz2k,where(27)ϕgz2k=1Vk∫-Δzk/2Δzk/2dzg2zϕgzkz,Qgz2k=1Vk∫-Δzk/2Δzk/2dzg2zQgzkz,Lgxyz2k=3h21Vk∫-Δzk/2Δzk/2dzg2zLgxykz.Substituteϕgzkz into ϕgz1k and ϕgz2k with (11) and then integrate to get (28a) and (28b):(28a)ϕgz1k=1Vk∫-Δzk/2Δzk/2g1zϕgzkzdz=1Vk∫-Δzk/2Δzk/2zΔzkVkΔzkϕ¯gk+agz1kfz1z+agz2kfz2z+agz3kfz3z+agz4kfz4zdz=112agz1k-1120agz3k,(28b)ϕgz2k=1Vk∫-Δzk/2Δzk/2g2zϕgzkzdz=1Vk∫-Δzk/2Δzk/2zΔzk2-112ϕ¯gk+agz1kfz1z+agz2kfz2z+agz3kfz3z+agz4kfz4zdz=1180agz2k-12100agz4k.Thus,agz3k and agz4k are shown in(29)agz3k=-120ϕgz1k+10agz1k=-120ϕgz1k+10ϕ¯gz+k-ϕ¯gz-k,agz4k=-2100ϕgz2k+353agz2k=-2100ϕgz2k+353ϕ¯gz+k+ϕ¯gz-k-2ϕ¯gk.The 3D multigroup nodal balance equation is obtained by integrating over a homogeneous node volumeVk(30)23hJ¯gx+k-J¯gx-k+J¯gu+k-J¯gu-k+J¯gv+k-J¯gv-k+1ΔzkJ¯gz+k-J¯gz-k+Σgr,kϕ¯gk=Q¯gk,where(31)Q¯gk=1Vk∫-Δzk/2Δzk/2dz∫-h/2h/2dx∫-ysxysxdyQgkx,y,z.ϕ¯gk is derived from (30) and shown in (32). ϕgx1k, ϕgu1k, and ϕgv1k are derived from (19) and shown in (33), (34), and (35). ϕgz1k and ϕgz2k are derived from (24) and (26) and shown in (36) and (37):(32)ϕ¯gk=1Σgr,kQ¯gk-23hΣgr,kJ¯gx+k-J¯gx-k+J¯gu+k-J¯gu-k+J¯gv+k-J¯gv-k-1ΔzkΣgr,kJ¯gz+k-J¯gz-k,(33)ϕgx1k=1Σgr,k+32Dgk/h2Qgx1k-1ΔzkLgzx1k-23hΣgr,k+32Dgk/h2Γ¯gxk+Γ¯guk-Γ¯gvk+40Dgk9h2Σgr,k+32Dgk/h2agx1k,(34)ϕgu1k=1Σgr,k+32Dgk/h2Qgu1k-1ΔzkLgzu1k-23hΣgr,k+32Dgk/h2Γ¯gxk+Γ¯guk+Γ¯gvk+40Dgk9h2Σgr,k+32Dgk/h2agu1k,(35)ϕgv1k=1Σgr,k+32Dgk/h2Qgv1k-1ΔzkLgzv1k-23hΣgr,k+32Dgk/h2-Γ¯gxk+Γ¯guk+Γ¯gvk+40Dgk9h2Σgr,k+32Dgk/h2agv1k,(36)ϕgz1k=1Σgr,kQgz1k-23hLgxyz1k-12ΔzkΣgr,kJ¯gz+k+J¯gz-k-DgkΔzk2Σgr,kagz1k,(37)ϕgz2k=1Σgr,kQgz2k-23hLgxyz2k-16ΔzkΣgr,kJ¯gz+k-J¯gz-k-Dgk3Δzk2Σgr,kagz2k.The partial net current and one-dimensional average spatially flux can be written in the form of outgoing and incoming partial currents as shown in(38)±J¯gw±k=J¯gw±out,k-J¯gw±in,k,ϕ¯gw±k=2J¯gw±out,k+J¯gw±in,k,w=x,u,v,z.By using (9) and (38), outgoing partial currents through radial-directed surfaces ofkth node in terms of flux expansion coefficients are shown in(39)J¯gw+out,k=J¯gw+k+J¯gw+in,k=-Dgkddwϕ¯gwkw+Dgky′w2ywEgwkwx=h/2+J¯gw+in,k=-Dgkhagw1k+3613agw2k+726agw3k+12agw4k+Egw+k+J¯gw+in,k,w=x,u,v,J¯gw-out,k=-J¯gw-k+J¯gw-in,k=--Dgkddwϕ¯gwkw+Dgky′w2ywEgwkwx=-h/2+J¯gw-in,k=Dgkhagw1k-3613agw2k-726agw3k+12agw4k-Egw-k+J¯gw-in,k,w=x,u,v,where(40)Egx+k=-1285hDgk17J¯gu+k-J¯gv-k-2J¯gv+k-J¯gu-k-11235179ϕ¯gx+k+49ϕ¯gx-k-228ϕ¯gk+1130agx3k,Egu+k=-1285hDgk17J¯gv+k-J¯gx+k-2J¯gx-k-J¯gv-k-11235179ϕ¯gu+k+49ϕ¯gu-k-228ϕ¯gk+1130agu3k,Egv+k=-1285hDgk17J¯gx-k-J¯gu+k-2J¯gu-k-J¯gx+k-11235179ϕ¯gv+k+49ϕ¯gv-k-228ϕ¯gk+1130agv3k,Egx-k=-1285hDgk17J¯gu+k-J¯gv-k-2J¯gv+k-J¯gu-k-1123549ϕ¯gx+k+179ϕ¯gx-k-228ϕ¯gk+1130agx3k,Egu-k=-1285hDgk17J¯gv+k-J¯gx+k-2J¯gx-k-J¯gv-k-1123549ϕ¯gu+k+179ϕ¯gu-k-228ϕ¯gk+1130agu3k,Egv-k=-1285hDgk17J¯gx-k-J¯gu+k-2J¯gu-k-J¯gx+k-1123549ϕ¯gv+k+179ϕ¯gv-k-228ϕ¯gk+1130agv3k.By using (11) and (38), outgoing partial currents through axial-directed surfaces ofkth node in terms of flux expansion coefficients are shown in(41)J¯gz+out,k=J¯gz+k+J¯gz+in,k=-Dgkddzϕ¯gzkzz=Δzk/2+J¯gz+in,k=-DgkΔzkVkddzϕ^gzkzz=Δzk/2+J¯gz+in,k=-DgkΔzkagz1k+agz2k+12agz3k+1740agz4k+J¯gz+in,k,J¯gz-out,k=-J¯gz-k+J¯gz-in,k=--Dgkddzϕ¯gzkzz=-Δzk/2+J¯gz-in,k=--DgkΔzkVkddzϕ^gzkzz=-Δzk/2+J¯gz-in,k=DgkΔzkagz1k-agz2k+12agz3k-1740agz4k+J¯gz-in,k.We use (14), (16), (22), and (29) for expansion coefficients along with (32), (33), (34), (35), (36), and (37) for node averaged flux and its moments and (49), (50), (51), (56), and (58) for node leakage moments. One-dimensional average spatially flux is expressed by partial currents with (38). Eight outgoing partial currents can be expressed in terms of eight incoming partial currents, source, and its moments and leakage moments.(42)
a
11
a
12
a
13
a
14
a
15
a
16
a
17
a
18
a
21
a
22
a
23
a
24
a
25
a
26
a
27
a
28
a
31
a
32
a
33
a
34
a
35
a
36
a
37
a
38
a
41
a
42
a
43
a
44
a
45
a
46
a
47
a
48
a
51
a
52
a
53
a
54
a
55
a
56
a
57
a
58
a
61
a
62
a
63
a
64
a
65
a
66
a
67
a
68
a
71
a
72
a
73
a
74
a
75
a
76
a
77
a
78
a
81
a
82
a
83
a
84
a
85
a
86
a
87
a
88
J
¯
g
x
+
out
,
k
J
¯
g
u
+
out
,
k
J
¯
g
v
+
out
,
k
J
¯
g
x
-
out
,
k
J
¯
g
u
-
out
,
k
J
¯
g
v
-
out
,
k
J
¯
g
z
+
out
,
k
J
¯
g
z
-
out
,
k
=
b
11
b
12
0
0
0
0
b
21
0
b
23
0
0
0
b
31
0
0
b
34
0
0
b
41
b
42
0
0
0
0
b
51
0
b
53
0
0
0
b
61
0
0
b
64
0
0
b
71
0
0
0
b
75
b
76
b
81
0
0
0
b
85
b
86
Q
¯
g
k
Q
g
x
1
k
-
1
Δ
z
k
L
g
z
x
1
k
Q
g
u
1
k
-
1
Δ
z
k
L
g
z
u
1
k
Q
g
v
1
k
-
1
Δ
z
k
L
g
z
v
1
k
Q
g
z
1
k
-
2
3
h
L
g
x
y
z
1
k
Q
g
z
2
k
-
2
3
h
L
g
x
y
z
2
k
+
c
11
c
12
c
13
c
14
c
15
c
16
c
17
c
18
c
21
c
22
c
23
c
24
c
25
c
26
c
27
c
28
c
31
c
32
c
33
c
34
c
35
c
36
c
37
c
38
c
41
c
42
c
43
c
44
c
45
c
46
c
47
c
48
c
51
c
52
c
53
c
54
c
55
c
56
c
57
c
58
c
61
c
62
c
63
c
64
c
65
c
66
c
67
c
68
c
71
c
72
c
73
c
74
c
75
c
76
c
77
c
78
c
81
c
82
c
83
c
84
c
85
c
86
c
87
c
88
J
¯
g
x
+
in
,
k
J
¯
g
u
+
in
,
k
J
¯
g
v
+
in
,
k
J
¯
g
x
-
in
,
k
J
¯
g
u
-
in
,
k
J
¯
g
v
-
in
,
k
J
¯
g
z
+
in
,
k
J
¯
g
z
-
in
,
k
⟹
A
g
k
J
_
g
out
,
k
=
B
g
k
Q
_
g
k
-
L
_
g
k
+
C
g
k
J
_
g
in
,
k
⟹
J
_
g
out
,
k
=
A
g
k
-
1
B
g
k
Q
_
g
k
-
L
_
g
k
+
A
g
k
-
1
C
g
k
J
_
g
in
,
k
⟹
J
_
g
out
,
k
=
P
g
k
Q
_
g
k
-
L
_
g
k
+
R
g
k
J
_
g
in
,
k
,where(43)Q_gk=col⁡Q¯gk,Qgx1k,Qgu1k,Qgv1k,Qgz1k,Qgz2k,L_gk=col⁡0,1ΔzkLgzx1k,1ΔzkLgzu1k,1ΔzkLgzv1k,23hLgxyz1k,23hLgxyz2k.Matrix elements (aij,bij,cij) are given below. These elements are characterized by node dimensions and homogenized material of the node.(44)a11=a22=a33=a44=a55=a66=1+2Dgkh468133+83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk,a12=a23=a34=a45=a56=a21=a32=a43=a54=a65=a16=a61=24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk-13133,a13=a24=a35=a46=a31=a42=a53=a64=a15=a26=a51=a62=24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk-6133,a14=a25=a36=a41=a52=a63=2Dgkh216133-83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk,a17=a27=a37=a47=a57=a67=a18=a28=a38=a48=a58=a68=36Dgk7ΔzkhΣgr,k,a71=a72=a73=a74=a75=a76=a81=a82=a83=a84=a85=a86=2Dgk3ΔzkhΣgr,k14312+595DgkΔzk2Σgr,k,a77=a88=DgkΔzk28712+715DgkΔzk2Σgr,k+5723ΔzkΣgr,k+595DgkΔzk3Σgr,k2+1,a78=a87=DgkΔzk-112+475DgkΔzk2Σgr,k+3923ΔzkΣgr,k+595DgkΔzk3Σgr,k2,b11=b21=b31=b41=b51=b61=36Dgk7hΣgr,k,b12=b23=b34=12hDgkh2Σgr,k+32Dgk,b42=b53=b64=-12hDgkh2Σgr,k+32Dgk,b71=b81=DgkΔzk14312Σgr,k+595DgkΔzk2Σgr,k2,b75=60DgkΔzkΣgr,k,b85=-60DgkΔzkΣgr,k,b76=b86=1785Dgk2ΔzkΣgr,k,c11=c22=c33=c44=c55=c66=1-2Dgkh468133+83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk,c12=c23=c34=c45=c56=c21=c32=c43=c54=c65=c16=c61=24Dgk7h2Σgr,k+8Dgkh2Σgr,k+32Dgk-13133,c13=c24=c35=c46=c31=c42=c53=c64=c15=c26=c51=c62=24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk-6133,c14=c25=c36=c41=c52=c63=-2Dgkh216133-83-160Dgk3h2Σgr,k+32Dgk+24Dgk7h2Σgr,k-8Dgkh2Σgr,k+32Dgk,c17=c27=c37=c47=c57=c67=c18=c28=c38=c48=c58=c68=36Dgk7ΔzkhΣgr,k,c71=c72=c73=c74=c75=c76=c81=c82=c83=c84=c85=c86=2Dgk3ΔzkhΣgr,k14312+595DgkΔzk2Σgr,k,c77=c88=DgkΔzk5723ΔzkΣgr,k+595Dgk3Δzk3Σgr,k2-715DgkΔzk2Σgr,k-28712+1,c78=c87=DgkΔzk3923ΔzkΣgr,k+595DgkΔzk3Σgr,k2+112-475DgkΔzk2Σgr,k.


## 3. Calculation of the Leakage Moments

Leakage moments in the response matrix are derived in this section. The partially integrated axial leakageLgzkx and partially integrated radial leakage Lgxykz are expanded in terms of polynomials, but here it is restricted up to second order as the accuracy achieved in calculation is insignificant with respect to computational effort involved in higher order expansion.


### 3.1. The Radial Moment of the Partially Integrated Axial Leakage

Thex-direction moment of the partially integrated axial leakage Lgzx1k is calculated by using the approximation shown in (45) in which f1x and f2x are as given in (10). The partially integrated axial leakage δgzkx is given by (45) and the total axial leakage is shown in (46):(45)∫-ysxysxdyLgzkx,y=Lgzkx=L¯gzk+δgz1kf1x+δgz2kf2x,where(46)L¯gzk=ΔzkVk∫-h/2h/2dx∫-ysxysxLgzkx,y=ΔzkVk∫-h/2h/2dxδgzkx,y=J¯gz+k-J¯gz-k.The coefficientsδgz1k and δgz2k are calculated in the following manner. Let kx- and kx+ denote the neighboring nodes in the minus and plusx-directions and the quadratic polynomial extends over the three nodes kx-, k, and kx+. The expansion coefficients in (45) are calculated such that the total axial leakages shown in (46) in the nodes kx- and kx+ are preserved as shown in (47). The expansion coefficients δgz1k and δgz2k are calculated by using the constraints given in (47) and the results are shown in (48):(47)L¯gzkx-=ΔzkVk∫-h/2-h-h/2dxδgzkx,L¯gzkx+=ΔzkVk∫h/2h/2+hdxδgzkx,(48)δgz1k=VkL¯gzkx+-L¯gzkx-2Δzkh,δgz2k=1373VkL¯gzkx++L¯gzkx-Δzkh-2L¯gzk.The required leakage momentLgzx1k is calculated by substituting (45) into (20) and performing the necessary integrations and the total axial leakages L¯gzkx+ and L¯gzkx- shown in (46) can be computed by axial direction net currents. The final form of thex-direction moment of the partially integrated axial leakage Lgzx1k is shown in (49). Similar equation for Lgzu1k and Lgzv1k can be obtained as shown in (50) and (51):(49)Lgzx1k=ΔzkVk∫-h/2h/2dxf1x∫-ysxysxdyLgzkx,y=ΔzkhVkL¯gzk+δgz2k26=ΔzkhVkL¯gzk+1146VkL¯gzkx++L¯gzkx-Δzkh-2L¯gzk,(50)Lgzu1k=ΔzkhVkL¯gzk+1146VkL¯gzku++L¯gzku-Δzkh-2L¯gzk,(51)Lgzv1k=ΔzkhVkL¯gzk+1146VkL¯gzkv++L¯gzkv-Δzkh-2L¯gzk.


### 3.2. The Axial Moment of the Partially Integrated Radial Leakage

The axial moment of the partially integrated radial leakageLgxyz1k is calculated by using the approximation equation (52). g1z and g2z are given in (12). The total radial leakage L¯gxyk can be written in terms of the average leakages in the three radial directions as shown in (53):(52)Lgxykz=L¯gxyk+δgxy1kg1z+δgxy2kg2z,where(53)L¯gxyk=3h21Vk∫-Δzk/2Δzk/2Lgxykz=L¯gxk+L¯guk+L¯gvk.The coefficientsδgxy1k and δgxy2k are calculated in the following manner. Let kz- and kz+ denote the neighboring nodes in the minus and plusz-directions and the quadratic polynomial extends over the three nodes kz-, k, and kz+. Integrating (54) over kz- and kz+, respectively, and the expansion coefficients δgxy1k and δgxy2k is calculated by using the constraints given in (54) and the results are shown in (55).(54)L¯gxykz-=3h21Vkz-∫-Δzk/2-Δzkz--Δzk/2dzLgxykz,L¯gxykz+=3h21Vkz+∫Δzk/2Δzk/2+Δzkz+dzLgxykz,(55)δgxy1k=23hVkd2Δzkz-+ΔzkΔzkz-+ΔzkL¯gxykz+-L¯gxyk+Δzk+Δzkz+Δzk+2Δzkz+L¯gxyk-L¯gxykz-,δgxy2k=2Δzk3hVkdΔzkz-+ΔzkL¯gxykz+-L¯gxyk+Δzk+Δzkz+L¯gxyk-L¯gxykz-,d=Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk.The required leakage momentLgxyz1k is calculated by substituting (52) into (25) and performing the necessary integrations to get the final form for Lgxyz1k shown in(56)Lgxyz1k=3h21Vk∫-Δzk/2Δzk/2dzg1zLgxykz=h8ΔzkVkδgxy1k=ξ1kz-L¯gxyk-L¯gxykz--ξ1kz+L¯gxyk-L¯gxykz+,where(57)ξ1kz-=Δzk2Δzkz++ΔzkΔzkz++Δzk12Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk,ξ1kz+=Δzk2Δzkz-+ΔzkΔzkz-+Δzk12Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk.The required leakage momentLgxyz2k is calculated by substituting (52) into (27) and performing the necessary integrations to get the final form for Lgxyz2k shown in(58)Lgxyz2k=3h21Vk∫-Δzk/2Δzk/2dzg2zLgxykz=h8ΔzkVkδgxy1k+115δgxy2k=ξ1kz-+ξ2kz-L¯gxyk-L¯gxykz--ξ1kz++ξ2kz+L¯gxyk-L¯gxykz+,where(59)ξ2kz-=Δzk2Δzkz++Δzk180Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk,ξ2kz+=Δzk2Δzkz-+Δzk180Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk.


## 3.1. The Radial Moment of the Partially Integrated Axial Leakage

Thex-direction moment of the partially integrated axial leakage Lgzx1k is calculated by using the approximation shown in (45) in which f1x and f2x are as given in (10). The partially integrated axial leakage δgzkx is given by (45) and the total axial leakage is shown in (46):(45)∫-ysxysxdyLgzkx,y=Lgzkx=L¯gzk+δgz1kf1x+δgz2kf2x,where(46)L¯gzk=ΔzkVk∫-h/2h/2dx∫-ysxysxLgzkx,y=ΔzkVk∫-h/2h/2dxδgzkx,y=J¯gz+k-J¯gz-k.The coefficientsδgz1k and δgz2k are calculated in the following manner. Let kx- and kx+ denote the neighboring nodes in the minus and plusx-directions and the quadratic polynomial extends over the three nodes kx-, k, and kx+. The expansion coefficients in (45) are calculated such that the total axial leakages shown in (46) in the nodes kx- and kx+ are preserved as shown in (47). The expansion coefficients δgz1k and δgz2k are calculated by using the constraints given in (47) and the results are shown in (48):(47)L¯gzkx-=ΔzkVk∫-h/2-h-h/2dxδgzkx,L¯gzkx+=ΔzkVk∫h/2h/2+hdxδgzkx,(48)δgz1k=VkL¯gzkx+-L¯gzkx-2Δzkh,δgz2k=1373VkL¯gzkx++L¯gzkx-Δzkh-2L¯gzk.The required leakage momentLgzx1k is calculated by substituting (45) into (20) and performing the necessary integrations and the total axial leakages L¯gzkx+ and L¯gzkx- shown in (46) can be computed by axial direction net currents. The final form of thex-direction moment of the partially integrated axial leakage Lgzx1k is shown in (49). Similar equation for Lgzu1k and Lgzv1k can be obtained as shown in (50) and (51):(49)Lgzx1k=ΔzkVk∫-h/2h/2dxf1x∫-ysxysxdyLgzkx,y=ΔzkhVkL¯gzk+δgz2k26=ΔzkhVkL¯gzk+1146VkL¯gzkx++L¯gzkx-Δzkh-2L¯gzk,(50)Lgzu1k=ΔzkhVkL¯gzk+1146VkL¯gzku++L¯gzku-Δzkh-2L¯gzk,(51)Lgzv1k=ΔzkhVkL¯gzk+1146VkL¯gzkv++L¯gzkv-Δzkh-2L¯gzk.


## 3.2. The Axial Moment of the Partially Integrated Radial Leakage

The axial moment of the partially integrated radial leakageLgxyz1k is calculated by using the approximation equation (52). g1z and g2z are given in (12). The total radial leakage L¯gxyk can be written in terms of the average leakages in the three radial directions as shown in (53):(52)Lgxykz=L¯gxyk+δgxy1kg1z+δgxy2kg2z,where(53)L¯gxyk=3h21Vk∫-Δzk/2Δzk/2Lgxykz=L¯gxk+L¯guk+L¯gvk.The coefficientsδgxy1k and δgxy2k are calculated in the following manner. Let kz- and kz+ denote the neighboring nodes in the minus and plusz-directions and the quadratic polynomial extends over the three nodes kz-, k, and kz+. Integrating (54) over kz- and kz+, respectively, and the expansion coefficients δgxy1k and δgxy2k is calculated by using the constraints given in (54) and the results are shown in (55).(54)L¯gxykz-=3h21Vkz-∫-Δzk/2-Δzkz--Δzk/2dzLgxykz,L¯gxykz+=3h21Vkz+∫Δzk/2Δzk/2+Δzkz+dzLgxykz,(55)δgxy1k=23hVkd2Δzkz-+ΔzkΔzkz-+ΔzkL¯gxykz+-L¯gxyk+Δzk+Δzkz+Δzk+2Δzkz+L¯gxyk-L¯gxykz-,δgxy2k=2Δzk3hVkdΔzkz-+ΔzkL¯gxykz+-L¯gxyk+Δzk+Δzkz+L¯gxyk-L¯gxykz-,d=Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk.The required leakage momentLgxyz1k is calculated by substituting (52) into (25) and performing the necessary integrations to get the final form for Lgxyz1k shown in(56)Lgxyz1k=3h21Vk∫-Δzk/2Δzk/2dzg1zLgxykz=h8ΔzkVkδgxy1k=ξ1kz-L¯gxyk-L¯gxykz--ξ1kz+L¯gxyk-L¯gxykz+,where(57)ξ1kz-=Δzk2Δzkz++ΔzkΔzkz++Δzk12Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk,ξ1kz+=Δzk2Δzkz-+ΔzkΔzkz-+Δzk12Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk.The required leakage momentLgxyz2k is calculated by substituting (52) into (27) and performing the necessary integrations to get the final form for Lgxyz2k shown in(58)Lgxyz2k=3h21Vk∫-Δzk/2Δzk/2dzg2zLgxykz=h8ΔzkVkδgxy1k+115δgxy2k=ξ1kz-+ξ2kz-L¯gxyk-L¯gxykz--ξ1kz++ξ2kz+L¯gxyk-L¯gxykz+,where(59)ξ2kz-=Δzk2Δzkz++Δzk180Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk,ξ2kz+=Δzk2Δzkz-+Δzk180Δzkz-+ΔzkΔzkz-+Δzk+Δzkz+Δzkz++Δzk.


## 4. Numerical Solution of the Nodal Equations

### 4.1. Overview of the Solution Procedure

The effective multiplication factor, neutron flux, and fission power distribution are obtained by using iteration technique to solve the neutron diffusion equation. The nodal equations are solved by using a conventional fission source iteration procedure accelerated by coarse-mesh rebalance in combination with the asymptotic source extrapolation technique. At each fission source iteration, the solution of the interface partial currents for each group is accomplished via a series of sweeps through the spatial mesh.The algorithm used to solve the nodal equations is shown in Figure3 and steps are briefly described here:(1)
We calculate nodal coupling coefficients shown in (42). This calculation is performed only once before the iteration starts.(2)
initialize the node flux moment vectorϕ_gk, eigenvalue keff0, rebalance factors f_m and partial currents vector J_gout,k and J_gin,k to initial value (assuming ϕ_gk=1.0,keff0=1.0, J_gout,k=J_gin,k=0.0) for all energy groups g and all nodes k. The initial total fission source vector forkth node is calculated by(60)F_k0=∑g′=1GνΣg′f,kϕ_g′k0.(3)
At the beginning of outer iterationn, the group partial currents J_gout,k are scaled by rebalance factors f_m calculated by coarse-mesh rebalance procedure as shown in(61)J_^gout,kn=f_mJ_gout,kn.(4)
Fission source vector forkth nodes and gth group, which is defined in (62), is calculated from flux moment vector and eigenvalue obtained in previous outer iteration (n-1) (or, it will be guess value if it is the first outer iteration)(62)QF_gkn-1=1keffn-1χg∑g′=1GνΣg′f,kϕ_g′kn-1=1keffn-1χgFgkn-1.(5)
From this step, inner iteration (group iteration) starts. Scattering source vector, as given in (63), is calculated with the most recently computed flux moment vector(63)QS_gkn=∑g′=1g′<gGΣg′→gs,kϕ_g′kn.Then, construct group source vector forkth node and gth group due to fission source vector and scattering source vector by(64)Q_gkn=QF_gkn-1+QS_gkn.(6)
Leakage moment vector forkth nodes and gth group defined in (49), (50), (51), (56), and (58) is calculated by using the most recently computed incoming and outgoing partial currents.(7)
The response matrix defined in (42) forkth nodes and gth group is calculated by using latest incoming currents (these outgoing partial currents are assigned to the incoming partial currents of neighboring nodes through respective surfaces).(8)
The flux moment vector defined in (32), (33), (34), (35), (36), and (37) is calculated by using the most recently computed partial currents, total source vector calculated in step (5), and leakage moment vector calculated in step (6).(9)
Steps (5)–(8) are repeated for all nodes and all energy groups with keeping fission source vector calculated in step (4) constant.(10)
Once all energy groups have been processed, the coarse-mesh rebalance equations are solved to get coarse-mesh rebalance factors and eigenvalue. The total fission source vector forkth nodes is scaled by rebalance factors prior to checking the convergence of the fission source as shown in(65)F_^kn=f_mF_kn,k∈Vm.(11)
A new eigenvaluekeffn is calculated from keffn-1, the total fission source Fkn-1 obtained in previous outer iteration (n-1), and the total fission source Fkn obtained in present outer iteration n as shown in(66)keffn=keffn-1FnFn-1,where(67)Fkn=∑g′=1GνΣg′f,kϕ¯g′kn.The outer iterations are terminated when three convergence criteria shown in (68) are satisfied:(68)keffn-keffn-1<ε1=1.0×10-7,maxk⁡Fkn-Fkn-1Fkn<ε2=1.0×10-5,1K∑k=1KFkn-Fkn-1Fkn21/2<ε3=1.0×10-5.Figure 3
Overview of the nodal solution algorithm.


### 4.2. An Improved Coarse-Mesh Rebalance Method

The outer (fission source) iterations are accelerated by using an improved coarse-mesh rebalance method in combination with the asymptotic source extrapolation technique. The basic idea of the method is to scale the fluxes calculated at each outer iteration on the “fine-mesh” by rebalance factors such that a neutron balance is enforced over each cell (region) of a “coarse-mesh” superimposed on the fine mesh. This approach is nonlinear since the fine mesh fluxes are used to compute the coefficients of the coarse-mesh equations. In this section, the traditional coarse-mesh rebalance method improved by Wielandt method and vector normalization procedure is described.The outer iterations are accelerated by using coarse-mesh rebalancing in which each ring of hexagons corresponds to a hex-plane coarse-mesh rebalance zone. The coarse-mesh equations can be combined in the form shown in(69)Mf_m=1λPf_m,where f_m is coarse-mesh rebalance factor vector. Equation (69) is constructed and solved following each outer iteration. The solution to this eigenvalue problem can be obtained using the Wielandt method. For problems in which the [M] matrix can be inverted directly, the Wielandt method is often more efficient. This approach is based on the application of the power method to the “shifted” eigenvalue problem obtained by rewriting (69) to(70)M~f_m=1λ~P~f_m,where(71)M~=M-1λeP,1λ~=1λ-1λe.The convergence rate of the power method is determined by the dominance ratio of the matrixM-1P; the closer this ratio is to 1, the slower the convergence rate will be. It can be shown that, for λe>λ0, the dominance ratio of M~-1P is smaller than that of M-1P. Hence the Wielandt method, which is obtained by applying the power method to (70), will converge faster than the power method applied directly to (69). However, the application of power method to calculate eigenvalues of a matrix and corresponding eigenvectors may lead to divergence of calculation. Thus, it is necessary to improve the power method by normalizing vector. The Wielandt method improved by vector normalization procedure, which is obtained by applying the power method to (70), will converge faster and be more stable than the power method applied directly to (69). We solve (70) by using the following iterative procedure as shown in(72)S_0=f_0m≠0,S_k=M~-1Pf_k-1m,λ~k=maxS_k,f_km=S_kλ~k,k=1,2,….Numerical calculations to date have demonstrated that only two power iterations are required at each power iteration (k in (72) is equal to 2) due to the efficiency of the Wielandt method. The computed rebalance factors are used to scale the partial currents and fission source moments in accordance with (62) and (65). The fission source is scaled prior to checking the convergence of the fission source, while the group partial currents are scaled at the beginning of the loop over groups in the next outer iteration in order to avoid an additional group loop following the rebalance procedure. It is not necessary to scale the flux moments ϕ_gk(n) because, as shown in (65), they do not enter into the calculation of the group source term Q_gk(n+1) due to the assumption that there is no upscatter. The final estimate for the eigenvalue at the nth outer iteration λn is obtained from (73). The coarse-mesh rebalance is performed once every five outer iterations:(73)1λn=1λ~+1λe.


## 4.1. Overview of the Solution Procedure

The effective multiplication factor, neutron flux, and fission power distribution are obtained by using iteration technique to solve the neutron diffusion equation. The nodal equations are solved by using a conventional fission source iteration procedure accelerated by coarse-mesh rebalance in combination with the asymptotic source extrapolation technique. At each fission source iteration, the solution of the interface partial currents for each group is accomplished via a series of sweeps through the spatial mesh.The algorithm used to solve the nodal equations is shown in Figure3 and steps are briefly described here:(1)
We calculate nodal coupling coefficients shown in (42). This calculation is performed only once before the iteration starts.(2)
initialize the node flux moment vectorϕ_gk, eigenvalue keff0, rebalance factors f_m and partial currents vector J_gout,k and J_gin,k to initial value (assuming ϕ_gk=1.0,keff0=1.0, J_gout,k=J_gin,k=0.0) for all energy groups g and all nodes k. The initial total fission source vector forkth node is calculated by(60)F_k0=∑g′=1GνΣg′f,kϕ_g′k0.(3)
At the beginning of outer iterationn, the group partial currents J_gout,k are scaled by rebalance factors f_m calculated by coarse-mesh rebalance procedure as shown in(61)J_^gout,kn=f_mJ_gout,kn.(4)
Fission source vector forkth nodes and gth group, which is defined in (62), is calculated from flux moment vector and eigenvalue obtained in previous outer iteration (n-1) (or, it will be guess value if it is the first outer iteration)(62)QF_gkn-1=1keffn-1χg∑g′=1GνΣg′f,kϕ_g′kn-1=1keffn-1χgFgkn-1.(5)
From this step, inner iteration (group iteration) starts. Scattering source vector, as given in (63), is calculated with the most recently computed flux moment vector(63)QS_gkn=∑g′=1g′<gGΣg′→gs,kϕ_g′kn.Then, construct group source vector forkth node and gth group due to fission source vector and scattering source vector by(64)Q_gkn=QF_gkn-1+QS_gkn.(6)
Leakage moment vector forkth nodes and gth group defined in (49), (50), (51), (56), and (58) is calculated by using the most recently computed incoming and outgoing partial currents.(7)
The response matrix defined in (42) forkth nodes and gth group is calculated by using latest incoming currents (these outgoing partial currents are assigned to the incoming partial currents of neighboring nodes through respective surfaces).(8)
The flux moment vector defined in (32), (33), (34), (35), (36), and (37) is calculated by using the most recently computed partial currents, total source vector calculated in step (5), and leakage moment vector calculated in step (6).(9)
Steps (5)–(8) are repeated for all nodes and all energy groups with keeping fission source vector calculated in step (4) constant.(10)
Once all energy groups have been processed, the coarse-mesh rebalance equations are solved to get coarse-mesh rebalance factors and eigenvalue. The total fission source vector forkth nodes is scaled by rebalance factors prior to checking the convergence of the fission source as shown in(65)F_^kn=f_mF_kn,k∈Vm.(11)
A new eigenvaluekeffn is calculated from keffn-1, the total fission source Fkn-1 obtained in previous outer iteration (n-1), and the total fission source Fkn obtained in present outer iteration n as shown in(66)keffn=keffn-1FnFn-1,where(67)Fkn=∑g′=1GνΣg′f,kϕ¯g′kn.The outer iterations are terminated when three convergence criteria shown in (68) are satisfied:(68)keffn-keffn-1<ε1=1.0×10-7,maxk⁡Fkn-Fkn-1Fkn<ε2=1.0×10-5,1K∑k=1KFkn-Fkn-1Fkn21/2<ε3=1.0×10-5.Figure 3
Overview of the nodal solution algorithm.


## 4.2. An Improved Coarse-Mesh Rebalance Method

The outer (fission source) iterations are accelerated by using an improved coarse-mesh rebalance method in combination with the asymptotic source extrapolation technique. The basic idea of the method is to scale the fluxes calculated at each outer iteration on the “fine-mesh” by rebalance factors such that a neutron balance is enforced over each cell (region) of a “coarse-mesh” superimposed on the fine mesh. This approach is nonlinear since the fine mesh fluxes are used to compute the coefficients of the coarse-mesh equations. In this section, the traditional coarse-mesh rebalance method improved by Wielandt method and vector normalization procedure is described.The outer iterations are accelerated by using coarse-mesh rebalancing in which each ring of hexagons corresponds to a hex-plane coarse-mesh rebalance zone. The coarse-mesh equations can be combined in the form shown in(69)Mf_m=1λPf_m,where f_m is coarse-mesh rebalance factor vector. Equation (69) is constructed and solved following each outer iteration. The solution to this eigenvalue problem can be obtained using the Wielandt method. For problems in which the [M] matrix can be inverted directly, the Wielandt method is often more efficient. This approach is based on the application of the power method to the “shifted” eigenvalue problem obtained by rewriting (69) to(70)M~f_m=1λ~P~f_m,where(71)M~=M-1λeP,1λ~=1λ-1λe.The convergence rate of the power method is determined by the dominance ratio of the matrixM-1P; the closer this ratio is to 1, the slower the convergence rate will be. It can be shown that, for λe>λ0, the dominance ratio of M~-1P is smaller than that of M-1P. Hence the Wielandt method, which is obtained by applying the power method to (70), will converge faster than the power method applied directly to (69). However, the application of power method to calculate eigenvalues of a matrix and corresponding eigenvectors may lead to divergence of calculation. Thus, it is necessary to improve the power method by normalizing vector. The Wielandt method improved by vector normalization procedure, which is obtained by applying the power method to (70), will converge faster and be more stable than the power method applied directly to (69). We solve (70) by using the following iterative procedure as shown in(72)S_0=f_0m≠0,S_k=M~-1Pf_k-1m,λ~k=maxS_k,f_km=S_kλ~k,k=1,2,….Numerical calculations to date have demonstrated that only two power iterations are required at each power iteration (k in (72) is equal to 2) due to the efficiency of the Wielandt method. The computed rebalance factors are used to scale the partial currents and fission source moments in accordance with (62) and (65). The fission source is scaled prior to checking the convergence of the fission source, while the group partial currents are scaled at the beginning of the loop over groups in the next outer iteration in order to avoid an additional group loop following the rebalance procedure. It is not necessary to scale the flux moments ϕ_gk(n) because, as shown in (65), they do not enter into the calculation of the group source term Q_gk(n+1) due to the assumption that there is no upscatter. The final estimate for the eigenvalue at the nth outer iteration λn is obtained from (73). The coarse-mesh rebalance is performed once every five outer iterations:(73)1λn=1λ~+1λe.


## 5. Validation of HNHEX

In order to validate the effective multiplication factor, neutron flux, and fission power distribution calculation of code HNHEX, two 3D benchmark problems are selected. First problem is a three-dimensional LMFBR benchmark problem which is a simplified model of MARK-I core design of SNR 300 prototype LMFBR. The second problem is a three-dimensional VVER-440 benchmark problem. It should be noted that all calculations are performed by a desktop computer with CPU 3.1 GHz.


### 5.1. The Three-Dimensional SNR Benchmark Problem

The 3D SNR benchmark problem was introduced by Lawrence in ANL-7416 [12]. This benchmark problem is derived from the triangular-z model by altering the outer radial boundary (while preserving the total volume of the core) to allow imposition of boundary conditions on surfaces of hexagons. This modification thus permits solution by both hexagonal-geometry and triangular-geometry codes. Figures 4 and 5 show the three-dimensional layout of the model reactor. The model consists of a two-zone core surrounded by radial and axial blankets without a reflector. The height of the active core is 95 cm and each axial blanket is 40 cm thick. A total of 11 rings of hexagons (including the central hexagon) are included in the model with a lattice pitch of 11.2 cm. Vacuum boundary conditions are imposed on the outer surfaces of the blankets. Materials 1 and 2 are the inner and outer zone consisting of a PuO2–UO2 mixture as fuel whereas materials 3 and 4 represent axial and radial blanket consisting of U235-depleted uranium. Finally, materials 5 and 6 stand for absorber-material and follower. The whole core model includes a total of 18 control rods. In order to simulate a realistic three-dimensional problem, each of the absorber rods on the inner ring is withdrawn to the upper core/blanket boundary whereas the outer ring of control rods is half-inserted. More details about four-group homogenized cross sections for assemblies are from [13].Figure 4
Cross section of SNR benchmark calculation model for hexagonal-z geometry.Figure 5
Vertical section of SNR benchmark calculation model for hexagonal-z geometry.The HNHEX of traditional method and high order method both have been used to solve the diffusion equation with four energy groups for whole core calculation. The hexagonal node size of the radial direction is a fuel assembly dimension. The reference solution is obtained by Richardson extrapolation of DIF3D (6 triangular meshes per assembly) and DIF3D (24 triangular meshes per assembly) solutions.Several codes were used to calculate whole core of 3D SNR problem and Table1 shows the numerical results of 3D SNR problem including effective multiplication factor and its relative percentage error εk; relative percentage error of region-averaged fluxes in comparison to reference solution; and the CPU times. εIC, εOC, εRB, εAB, and εAO are the errors of the group- and region-averaged fluxes for the inner core, outer core, radial blanket, axial blanket, and absorber regions, respectively. DIF3D nodal conducts diffusion calculation with using nodal expansion method. DIF3D6 and DIF3D(24) perform finite-difference calculations with 6 and 24 triangular mesh cells per hexagonal assembly, respectively. It should be noted that calculations of HNHEX are performed by a desktop computer with CPU 3.1 GHz. The calculations of DIF3D Nodal and DIF(6,24) were performed using the IBM 370/195 computer. The calculations of AFEN were performed using the SUN Sparc/2 computer [9]. According to Table 1, the higher accuracy of results is belonging to the solution with using quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux with less effective multiplication factor and region-averaged fluxes errors relative to the flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux. The results calculated by HNHEX with 18 axial planes are more accurate than that with 8 axial planes. Table 2 shows the region-averaged fluxes calculated by HNHEX with a high order nodal expansion method. These fluxes are normalized to a total power of 3 watts over the whole-core model using a power conversion factor of 3.1×1010 fissions/watt-sec. It can be seen from Table 2 that the major flux errors appear in the absorber regions and blanket assemblies and it is due to high flux gradient existence in these regions.Table 1
HNHEX results for whole core 3D SNR benchmark in hexagonal-z geometry.
Code Number of axial planes k e f f k e f f error (pcm) ε I C (%) ε O C (%) ε R B (%) ε A B (%) ε C R (%) CPU time (sec) DIFa — 1.00989 — — — — — — — HNHEXb 8 1.01143 154 0.087 0.071 0.74 0.25 1.52 24.3 HNHEXc 8 1.01056 67 0.075 0.061 0.51 0.23 0.83 27.1 HNHEXc 18 1.01036 47 0.072 0.06 0.49 0.21 0.8 60.1 DIF3D Nodal 8 1.01150 161 −0.17 0.23 0.95 −0.30 −0.6 86.4 DIF3D Nodal 18 1.01125 136 −0.18 0.22 0.96 −0.11 −0.44 230.4 DIF3D(6) 36 1.01280 291 −0.27 0.42 0.47 −0.60 −1.72 576.0 DIF3D(24) 36 1.01118 129 −0.04 0.13 0.05 −0.64 −0.64 2304.0 AFEN 8 1.01334 345 0.24 0.05 −1.46 −0.81 −2.11 3207.6 AFEN 16 1.01134 145 0.08 −0.05 −1.29 −0.41 −0.92 7160.4 aReference. bApplying flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux. cApplying quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux.Table 2
Region-averaged fluxes calculated by HNHEX for 3D SNR benchmark in hexagonal-z geometry.
Group Region Average fluxes Calculated Reference Relative error (%) 1 Inner core 4.0504 E + 06 4.0569 E + 06 −0.16 Outer core 2.4888 E + 06 2.4870 E + 06 0.07 Radial blanket 3.0776 E + 05 3.0649 E + 05 0.41 Axial blanket 3.6263 E + 05 3.6273 E + 05 −0.03 Rods 1.2732 E + 06 1.2640 E + 06 0.73 Rod followers 2.3874 E + 06 2.3688 E + 06 0.79  2 Inner core 1.9106 E + 07 1.9027 E + 07 0.42 Outer core 1.0584 E + 07 1.0594 E + 07 −0.09 Radial blanket 2.2665 E + 06 2.2209 E + 06 2.05 Axial blanket 3.4752 E + 06 3.4705 E + 06 0.14 Rods 6.2578 E + 06 6.1898 E + 06 1.10 Rod followers 1.2544 E + 07 1.2375 E + 07 1.37  3 Inner core 1.7613 E + 06 1.7579 E + 06 0.19 Outer core 9.2373 E + 05 9.2860 E + 05 −0.52 Radial blanket 3.1861 E + 05 3.0960 E + 05 2.91 Axial blanket 5.2894 E + 05 5.2703 E + 05 0.36 Rods 4.6395 E + 05 4.6361 E + 05 0.07 Rod followers 1.4109 E + 06 1.3897 E + 06 1.53  4 Inner core 3.2760 E + 05 3.2706 E + 05 0.17 Outer core 1.5468 E + 05 1.5343 E + 05 0.81 Radial blanket 9.1536 E + 04 9.0518 E + 04 1.12 Axial blanket 1.9719 E + 05 1.9564 E + 05 0.79 Rods 5.3574 E + 04 5.3702 E + 04 −0.24 Rod followers 3.4212 E + 05 3.4029 E + 05 0.54


### 5.2. The Three-Dimensional VVER 440 Benchmark Problem

Three-dimensional VVER 440 benchmark [14] models a VVER-440 core in 30-degree symmetry using two-group diffusion approximation and given cross sections. The core is 250 cm high and covered with axial and radial reflectors and contains fuel elements of 3 different enrichments. The radial fuel assembly pitch is 14.7 cm. The control rods are large control elements with fuel followers. Rod Bank 6 is partially inserted, as shown in Figure 6. Vacuum boundary conditions are applied to the entire reflector outside boundaries. The two-group diffusion parameters are given in Table 3. The task is to calculate the effective multiplication factor and the 3D power distribution. The power distribution should be normalized to core average power density of unity and presented for axial node size of 25 cm. The reference solution was extrapolated from DIF3D-FD runs with 216 and 294 triangle/hexagon subdivisions and 2.5 cm axial mesh spacing. The node size used in HNHEX for 3D VVER 440 benchmark is one node per a fuel assembly in radial shape and axial sections with 25.0 cm height. Several codes were used to calculate 3D VVER 440 problem and Table 4 gives the numerical results of 3D VVER-440 test case. It should be noted that calculations of HNHEX are performed by a desktop computer with CPU 3.1 GHz. The calculations of ACNECH [15] were performed using a laptop computer with CPU 2.41 GHz. The CPU times are referred to HP9000/735-125 for ANC-H and AFEN.Table 3
Fuel assembly cross sections of three-dimensional VVER-440 benchmark.
Parameter FA type 1 2 3 4 5 6 D 1 (cm) 1.3466 1.3377 1.3322 1.1953 1.4485 1.3413 D 2 (cm) 0.37169 0.36918 0.36502 0.19313 0.25176 0.24871 Σ R , 1 (cm−1) 2.5255 E - 02 2.4709 E - 02 2.4350 E - 02 2.5636 E - 02 3.3184 E - 02 2.9301 E - 02 Σ R , 2 (cm−1) 6.4277 E - 02 7.9361 E - 02 1.0010 E - 01 1.3498 E - 01 3.2839 E - 02 6.4655 E - 02 Σ 1 → 2 (cm−1) 1.6893 E - 02 1.5912 E - 02 1.4888 E - 02 2.2264 E - 02 3.2262 E - 02 2.7148 E - 02 Σ f , 1 (cm−1) 2.21676 E - 03 2.79212 E - 03 3.59068 E - 03 0.0 0.0 0.0 Σ f , 2 (cm−1) 3.94368 E - 02 5.65720 E - 02 8.00000 E - 02 0.0 0.0 0.0 υ Σ f , 1 (cm−1) 4.4488 E - 03 5.5337 E - 03 7.0391 E - 03 0.0 0.0 0.0 υ Σ f , 2 (cm−1) 7.3753 E - 02 1.0581 E - 01 1.4964 E - 01 0.0 0.0 0.0Table 4
Numerical results of three-dimensional VVER 440 benchmark.
Code Number of axial planes Eigenvalue k e f f error (pcm) ε m a x (%) CPU time (sec) DIFa — 1.01132 — — — HNHEXb 12 1.01243 111 −8.6 5.16 HNHEXc 12 1.01186 54 −3.72 6.12 HNHEXc 24 1.01161 29 −3.32 12.3 ACNECH 12 1.01054 −78 −4.9 16 AFEN 12 1.01162 30 3.2 592.8 ANC-H 12 1.01157 25 1.28 124.8 aReference. bApplying flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux. cApplying quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux.Figure 6
Axial core cut for three-dimensional VVER 440 benchmark problem.The HNHEX of traditional method and high order method both have been used to solve the diffusion equation with two energy groups for whole-core calculation. According to Table4, implementing quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux, the accuracy of results has been improved relative to flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux by decreasing the maximum power error from −8.6% to −3.72%. There is no obvious increase in execution time with quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux. The keff error of 24 axial plane decreases by 25 pcm comparing to that of 12 axial planes for VVER 440. At last, Table 5 presents the three-dimensional power distribution obtained by HNHEX of high order method based on the assembly identification numbers as shown in Figure 7. Table 6 shows relative power error of the HNHEX solution of high order method to DIF3D reference solution for 3D VVER-440 problem.Table 5
Power distribution obtained by HNHEX for three-dimensional VVER-440 problem.
Assembly identity Axial layer 1 (bottom) 2 3 4 5 6 7 8 9 10 (top) 1 0.5153 1.0686 1.4489 1.5830 1.3531 0.0000 0.0000 0.0000 0.0000 0.0000 2 0.4212 0.8744 1.1866 1.2991 1.1610 0.8110 0.5808 0.4098 0.2560 0.1130 3 0.5467 1.1349 1.5412 1.6929 1.5700 1.2626 0.9446 0.6681 0.4180 0.1843 4 0.5605 1.1640 1.5816 1.7402 1.6293 1.3390 1.0129 0.7179 0.4498 0.1985 5 0.4500 0.9354 1.2718 1.3977 1.2922 1.0313 0.7690 0.5449 0.3419 0.1510 6 0.4616 0.9603 1.3061 1.4328 1.2736 0.8692 0.6182 0.4384 0.2754 0.1221 7 0.6063 1.2617 1.7174 1.8857 1.6150 0.0000 0.0000 0.0000 0.0000 0.0000 8 0.6383 1.3302 1.8146 2.0031 1.8099 1.2860 0.9433 0.6796 0.4311 0.1918 9 0.7672 1.6013 2.1873 2.4284 2.2882 1.8810 1.4465 1.0472 0.6661 0.2957 10 0.4705 0.9837 1.3447 1.4982 1.4348 1.2189 0.9545 0.6935 0.4415 0.1958 11 0.5420 1.1249 1.5273 1.6762 1.5443 1.2200 0.9054 0.6400 0.4003 0.1766 12 0.4410 0.9162 1.2448 1.3696 1.2824 1.0540 0.7974 0.5651 0.3539 0.1560 13 0.5561 1.1549 1.5699 1.7279 1.6163 1.3251 1.0020 0.7108 0.4458 0.1969 14 0.5624 1.1686 1.5896 1.7472 1.6042 1.2543 0.9282 0.6590 0.4141 0.1833 15 0.5835 1.2135 1.6520 1.8155 1.6193 1.1125 0.7970 0.5683 0.3581 0.1591 16 0.4902 1.0215 1.3929 1.5360 1.3829 0.9737 0.7097 0.5104 0.3233 0.1437 17 0.5088 1.0622 1.4509 1.6088 1.5051 1.2152 0.9270 0.6705 0.4262 0.1890 18 0.6799 1.4206 1.9417 2.1623 2.0653 1.7453 1.3631 0.9898 0.6304 0.2797 19 0.3970 0.8290 1.1339 1.2654 1.2187 1.0450 0.8241 0.5996 0.3820 0.1696 20 0.5503 1.1430 1.5625 1.7114 1.6083 1.3309 1.0121 0.7178 0.4506 0.1990 21 0.4336 0.9013 1.2264 1.3513 1.2656 1.0392 0.7879 0.5604 0.3522 0.1557 22 0.4392 0.9134 1.2443 1.3719 1.2741 1.0229 0.7694 0.5493 0.3465 0.1535 23 0.5766 1.1997 1.6367 1.8092 1.6795 1.3370 1.0081 0.7249 0.4591 0.2039 24 0.6063 1.2640 1.7271 1.9178 1.8087 1.4890 1.1461 0.8301 0.5276 0.2344 25 0.7107 1.4840 2.0294 2.2628 2.1683 1.8416 1.4443 1.0506 0.6694 0.2972 26 0.5042 1.0536 1.4416 1.6112 1.5558 1.3393 1.0598 0.7727 0.4927 0.2187 27 0.5334 1.1083 1.5100 1.6680 1.5696 1.2982 0.9911 0.7090 0.4473 0.1983 28 0.4360 0.9074 1.2387 1.3735 1.3010 1.0861 0.8381 0.6039 0.3825 0.1698 29 0.4586 0.9560 1.3074 1.4560 1.3925 1.1794 0.9219 0.6694 0.4257 0.1890 30 0.4619 0.9651 1.3212 1.4770 1.4264 1.2272 0.9710 0.7088 0.4519 0.2005 31 0.5678 1.1867 1.6252 1.8210 1.7690 1.5362 1.2243 0.8958 0.5718 0.2540 32 0.5621 1.1712 1.6019 1.7869 1.7195 1.4725 1.1590 0.8427 0.5360 0.2382 33 0.5708 1.1918 1.6326 1.8292 1.7778 1.5450 1.2316 0.9010 0.5747 0.2553 34 0.6013 1.2568 1.7227 1.9351 1.8909 1.6561 1.3287 0.9754 0.6234 0.2770 35 0.3622 0.7631 1.0465 1.1768 1.1527 1.0130 0.8156 0.5992 0.3832 0.1704 36 0.5077 1.0620 1.4563 1.6374 1.6036 1.4089 1.1333 0.8330 0.5325 0.2365 37 0.4136 0.8650 1.1867 1.3365 1.3132 1.1590 0.9362 0.6893 0.4411 0.1961Table 6
Relative power error of the HNHEX power distribution solution to DIF3D reference solution for three-dimensional VVER-440 problem (unit: %).
Assembly identity Axial layer 1 (bottom) 2 3 4 5 6 7 8 9 10 (top) 1 −0.96 −0.74 −0.63 −0.33 0.23 0.00 0.00 0.00 0.00 0.00 2 0.46 0.69 0.85 1.10 1.52 2.15 2.61 2.89 2.92 2.84 3 −0.52 −0.27 −0.12 0.04 0.31 0.79 1.03 1.14 1.26 1.12 4 −0.36 −0.10 0.03 0.17 0.41 0.82 1.09 1.12 1.23 1.10 5 0.31 0.54 0.69 0.84 1.07 1.54 1.71 1.85 1.93 1.82 6 0.12 0.36 0.48 0.71 1.08 1.46 1.68 2.10 2.07 2.01 7 −1.08 −0.82 −0.76 −0.50 −0.04 0.00 0.00 0.00 0.00 0.00 8 −1.25 −1.00 −0.91 −0.73 −0.48 −0.27 −0.18 0.15 0.15 0.03 9 −1.62 −1.27 −1.23 −1.14 −1.03 −0.82 −0.86 −0.72 −0.59 −0.85 10 −0.57 −0.21 −0.19 −0.14 −0.06 0.08 0.06 0.14 0.26 0.00 11 −0.50 −0.25 −0.10 0.08 0.36 0.90 1.17 1.30 1.40 1.27 12 0.48 0.72 0.87 1.01 1.27 1.68 1.97 1.99 2.11 1.96 13 −0.69 −0.44 −0.30 −0.17 0.07 0.47 0.71 0.76 0.87 0.71 14 −0.83 −0.59 −0.45 −0.29 −0.05 0.42 0.57 0.73 0.80 0.69 15 −0.94 −0.70 −0.59 −0.37 −0.04 0.27 0.51 0.84 0.83 0.72 16 −0.16 0.07 0.17 0.36 0.65 0.93 1.10 1.39 1.37 1.27 17 0.01 0.31 0.38 0.48 0.63 0.91 0.87 1.05 1.14 0.93 18 −1.06 −0.66 −0.64 −0.58 −0.49 −0.33 −0.36 −0.26 −0.11 −0.43 19 −2.81 −2.52 −2.50 −2.46 −2.38 −2.27 −2.35 −2.22 −2.11 −2.32 20 −0.66 −0.41 −1.11 −0.15 0.09 0.48 0.71 0.66 0.85 0.70 21 0.12 0.35 0.50 0.62 0.85 1.23 1.40 1.48 1.56 1.45 22 −0.11 0.10 0.24 0.37 0.56 0.92 0.98 1.13 1.19 1.07 23 −1.27 −1.05 −0.93 −0.81 −0.64 −0.30 −0.28 −0.13 −0.08 −0.23 24 −1.40 −1.15 −1.06 −0.97 −0.85 −0.59 −0.60 −0.50 −0.42 −0.61 25 −1.68 −1.33 −1.30 −1.25 −1.15 −0.98 −1.01 −0.91 −0.79 −1.06 26 −2.11 −1.77 −1.76 −1.72 −1.64 −1.53 −1.59 −1.48 −1.37 −1.62 27 −1.34 −1.15 −1.00 −0.90 −0.70 −0.35 −0.19 −0.18 −0.11 −0.21 28 −0.34 −0.14 −0.01 0.08 0.24 0.53 0.61 0.67 0.71 0.57 29 −0.38 −0.16 −0.06 0.01 0.14 0.36 0.42 0.46 0.51 0.35 30 0.14 0.44 0.50 0.55 0.65 0.81 0.83 0.88 0.95 0.72 31 −1.90 −1.55 −1.55 −1.50 −1.43 −1.32 −1.35 −1.27 −1.17 −1.44 32 −1.48 −1.25 −1.16 −1.10 −0.96 −0.76 −0.68 −0.66 −0.61 −0.77 33 −1.15 −0.88 −0.83 −0.77 −0.67 −0.52 −0.52 −0.44 −0.40 −0.61 34 −1.68 −1.32 −1.31 −1.27 −1.19 −1.08 −1.07 −1.01 −0.92 −1.22 35 −3.72 −2.66 −2.66 −2.62 −2.55 −2.46 −2.56 −2.40 −2.32 −2.54 36 −0.36 −0.01 0.00 0.04 0.12 0.22 0.21 0.29 0.35 0.08 37 −2.10 −1.76 −1.77 −1.73 −1.66 −1.59 −1.66 −1.54 −1.47 −1.71Figure 7
Radial core map and assemblies identification for three-dimensional VVER 440 benchmark.


## 5.1. The Three-Dimensional SNR Benchmark Problem

The 3D SNR benchmark problem was introduced by Lawrence in ANL-7416 [12]. This benchmark problem is derived from the triangular-z model by altering the outer radial boundary (while preserving the total volume of the core) to allow imposition of boundary conditions on surfaces of hexagons. This modification thus permits solution by both hexagonal-geometry and triangular-geometry codes. Figures 4 and 5 show the three-dimensional layout of the model reactor. The model consists of a two-zone core surrounded by radial and axial blankets without a reflector. The height of the active core is 95 cm and each axial blanket is 40 cm thick. A total of 11 rings of hexagons (including the central hexagon) are included in the model with a lattice pitch of 11.2 cm. Vacuum boundary conditions are imposed on the outer surfaces of the blankets. Materials 1 and 2 are the inner and outer zone consisting of a PuO2–UO2 mixture as fuel whereas materials 3 and 4 represent axial and radial blanket consisting of U235-depleted uranium. Finally, materials 5 and 6 stand for absorber-material and follower. The whole core model includes a total of 18 control rods. In order to simulate a realistic three-dimensional problem, each of the absorber rods on the inner ring is withdrawn to the upper core/blanket boundary whereas the outer ring of control rods is half-inserted. More details about four-group homogenized cross sections for assemblies are from [13].Figure 4
Cross section of SNR benchmark calculation model for hexagonal-z geometry.Figure 5
Vertical section of SNR benchmark calculation model for hexagonal-z geometry.The HNHEX of traditional method and high order method both have been used to solve the diffusion equation with four energy groups for whole core calculation. The hexagonal node size of the radial direction is a fuel assembly dimension. The reference solution is obtained by Richardson extrapolation of DIF3D (6 triangular meshes per assembly) and DIF3D (24 triangular meshes per assembly) solutions.Several codes were used to calculate whole core of 3D SNR problem and Table1 shows the numerical results of 3D SNR problem including effective multiplication factor and its relative percentage error εk; relative percentage error of region-averaged fluxes in comparison to reference solution; and the CPU times. εIC, εOC, εRB, εAB, and εAO are the errors of the group- and region-averaged fluxes for the inner core, outer core, radial blanket, axial blanket, and absorber regions, respectively. DIF3D nodal conducts diffusion calculation with using nodal expansion method. DIF3D6 and DIF3D(24) perform finite-difference calculations with 6 and 24 triangular mesh cells per hexagonal assembly, respectively. It should be noted that calculations of HNHEX are performed by a desktop computer with CPU 3.1 GHz. The calculations of DIF3D Nodal and DIF(6,24) were performed using the IBM 370/195 computer. The calculations of AFEN were performed using the SUN Sparc/2 computer [9]. According to Table 1, the higher accuracy of results is belonging to the solution with using quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux with less effective multiplication factor and region-averaged fluxes errors relative to the flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux. The results calculated by HNHEX with 18 axial planes are more accurate than that with 8 axial planes. Table 2 shows the region-averaged fluxes calculated by HNHEX with a high order nodal expansion method. These fluxes are normalized to a total power of 3 watts over the whole-core model using a power conversion factor of 3.1×1010 fissions/watt-sec. It can be seen from Table 2 that the major flux errors appear in the absorber regions and blanket assemblies and it is due to high flux gradient existence in these regions.Table 1
HNHEX results for whole core 3D SNR benchmark in hexagonal-z geometry.
Code Number of axial planes k e f f k e f f error (pcm) ε I C (%) ε O C (%) ε R B (%) ε A B (%) ε C R (%) CPU time (sec) DIFa — 1.00989 — — — — — — — HNHEXb 8 1.01143 154 0.087 0.071 0.74 0.25 1.52 24.3 HNHEXc 8 1.01056 67 0.075 0.061 0.51 0.23 0.83 27.1 HNHEXc 18 1.01036 47 0.072 0.06 0.49 0.21 0.8 60.1 DIF3D Nodal 8 1.01150 161 −0.17 0.23 0.95 −0.30 −0.6 86.4 DIF3D Nodal 18 1.01125 136 −0.18 0.22 0.96 −0.11 −0.44 230.4 DIF3D(6) 36 1.01280 291 −0.27 0.42 0.47 −0.60 −1.72 576.0 DIF3D(24) 36 1.01118 129 −0.04 0.13 0.05 −0.64 −0.64 2304.0 AFEN 8 1.01334 345 0.24 0.05 −1.46 −0.81 −2.11 3207.6 AFEN 16 1.01134 145 0.08 −0.05 −1.29 −0.41 −0.92 7160.4 aReference. bApplying flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux. cApplying quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux.Table 2
Region-averaged fluxes calculated by HNHEX for 3D SNR benchmark in hexagonal-z geometry.
Group Region Average fluxes Calculated Reference Relative error (%) 1 Inner core 4.0504 E + 06 4.0569 E + 06 −0.16 Outer core 2.4888 E + 06 2.4870 E + 06 0.07 Radial blanket 3.0776 E + 05 3.0649 E + 05 0.41 Axial blanket 3.6263 E + 05 3.6273 E + 05 −0.03 Rods 1.2732 E + 06 1.2640 E + 06 0.73 Rod followers 2.3874 E + 06 2.3688 E + 06 0.79  2 Inner core 1.9106 E + 07 1.9027 E + 07 0.42 Outer core 1.0584 E + 07 1.0594 E + 07 −0.09 Radial blanket 2.2665 E + 06 2.2209 E + 06 2.05 Axial blanket 3.4752 E + 06 3.4705 E + 06 0.14 Rods 6.2578 E + 06 6.1898 E + 06 1.10 Rod followers 1.2544 E + 07 1.2375 E + 07 1.37  3 Inner core 1.7613 E + 06 1.7579 E + 06 0.19 Outer core 9.2373 E + 05 9.2860 E + 05 −0.52 Radial blanket 3.1861 E + 05 3.0960 E + 05 2.91 Axial blanket 5.2894 E + 05 5.2703 E + 05 0.36 Rods 4.6395 E + 05 4.6361 E + 05 0.07 Rod followers 1.4109 E + 06 1.3897 E + 06 1.53  4 Inner core 3.2760 E + 05 3.2706 E + 05 0.17 Outer core 1.5468 E + 05 1.5343 E + 05 0.81 Radial blanket 9.1536 E + 04 9.0518 E + 04 1.12 Axial blanket 1.9719 E + 05 1.9564 E + 05 0.79 Rods 5.3574 E + 04 5.3702 E + 04 −0.24 Rod followers 3.4212 E + 05 3.4029 E + 05 0.54


## 5.2. The Three-Dimensional VVER 440 Benchmark Problem

Three-dimensional VVER 440 benchmark [14] models a VVER-440 core in 30-degree symmetry using two-group diffusion approximation and given cross sections. The core is 250 cm high and covered with axial and radial reflectors and contains fuel elements of 3 different enrichments. The radial fuel assembly pitch is 14.7 cm. The control rods are large control elements with fuel followers. Rod Bank 6 is partially inserted, as shown in Figure 6. Vacuum boundary conditions are applied to the entire reflector outside boundaries. The two-group diffusion parameters are given in Table 3. The task is to calculate the effective multiplication factor and the 3D power distribution. The power distribution should be normalized to core average power density of unity and presented for axial node size of 25 cm. The reference solution was extrapolated from DIF3D-FD runs with 216 and 294 triangle/hexagon subdivisions and 2.5 cm axial mesh spacing. The node size used in HNHEX for 3D VVER 440 benchmark is one node per a fuel assembly in radial shape and axial sections with 25.0 cm height. Several codes were used to calculate 3D VVER 440 problem and Table 4 gives the numerical results of 3D VVER-440 test case. It should be noted that calculations of HNHEX are performed by a desktop computer with CPU 3.1 GHz. The calculations of ACNECH [15] were performed using a laptop computer with CPU 2.41 GHz. The CPU times are referred to HP9000/735-125 for ANC-H and AFEN.Table 3
Fuel assembly cross sections of three-dimensional VVER-440 benchmark.
Parameter FA type 1 2 3 4 5 6 D 1 (cm) 1.3466 1.3377 1.3322 1.1953 1.4485 1.3413 D 2 (cm) 0.37169 0.36918 0.36502 0.19313 0.25176 0.24871 Σ R , 1 (cm−1) 2.5255 E - 02 2.4709 E - 02 2.4350 E - 02 2.5636 E - 02 3.3184 E - 02 2.9301 E - 02 Σ R , 2 (cm−1) 6.4277 E - 02 7.9361 E - 02 1.0010 E - 01 1.3498 E - 01 3.2839 E - 02 6.4655 E - 02 Σ 1 → 2 (cm−1) 1.6893 E - 02 1.5912 E - 02 1.4888 E - 02 2.2264 E - 02 3.2262 E - 02 2.7148 E - 02 Σ f , 1 (cm−1) 2.21676 E - 03 2.79212 E - 03 3.59068 E - 03 0.0 0.0 0.0 Σ f , 2 (cm−1) 3.94368 E - 02 5.65720 E - 02 8.00000 E - 02 0.0 0.0 0.0 υ Σ f , 1 (cm−1) 4.4488 E - 03 5.5337 E - 03 7.0391 E - 03 0.0 0.0 0.0 υ Σ f , 2 (cm−1) 7.3753 E - 02 1.0581 E - 01 1.4964 E - 01 0.0 0.0 0.0Table 4
Numerical results of three-dimensional VVER 440 benchmark.
Code Number of axial planes Eigenvalue k e f f error (pcm) ε m a x (%) CPU time (sec) DIFa — 1.01132 — — — HNHEXb 12 1.01243 111 −8.6 5.16 HNHEXc 12 1.01186 54 −3.72 6.12 HNHEXc 24 1.01161 29 −3.32 12.3 ACNECH 12 1.01054 −78 −4.9 16 AFEN 12 1.01162 30 3.2 592.8 ANC-H 12 1.01157 25 1.28 124.8 aReference. bApplying flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux. cApplying quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux.Figure 6
Axial core cut for three-dimensional VVER 440 benchmark problem.The HNHEX of traditional method and high order method both have been used to solve the diffusion equation with two energy groups for whole-core calculation. According to Table4, implementing quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux, the accuracy of results has been improved relative to flat approximation for partially integrated axial leakage and quadratic approximation for one-dimensional axial spatially flux by decreasing the maximum power error from −8.6% to −3.72%. There is no obvious increase in execution time with quadratic approximation for partially integrated axial leakage and four-order approximation for one-dimensional axial spatially flux. The keff error of 24 axial plane decreases by 25 pcm comparing to that of 12 axial planes for VVER 440. At last, Table 5 presents the three-dimensional power distribution obtained by HNHEX of high order method based on the assembly identification numbers as shown in Figure 7. Table 6 shows relative power error of the HNHEX solution of high order method to DIF3D reference solution for 3D VVER-440 problem.Table 5
Power distribution obtained by HNHEX for three-dimensional VVER-440 problem.
Assembly identity Axial layer 1 (bottom) 2 3 4 5 6 7 8 9 10 (top) 1 0.5153 1.0686 1.4489 1.5830 1.3531 0.0000 0.0000 0.0000 0.0000 0.0000 2 0.4212 0.8744 1.1866 1.2991 1.1610 0.8110 0.5808 0.4098 0.2560 0.1130 3 0.5467 1.1349 1.5412 1.6929 1.5700 1.2626 0.9446 0.6681 0.4180 0.1843 4 0.5605 1.1640 1.5816 1.7402 1.6293 1.3390 1.0129 0.7179 0.4498 0.1985 5 0.4500 0.9354 1.2718 1.3977 1.2922 1.0313 0.7690 0.5449 0.3419 0.1510 6 0.4616 0.9603 1.3061 1.4328 1.2736 0.8692 0.6182 0.4384 0.2754 0.1221 7 0.6063 1.2617 1.7174 1.8857 1.6150 0.0000 0.0000 0.0000 0.0000 0.0000 8 0.6383 1.3302 1.8146 2.0031 1.8099 1.2860 0.9433 0.6796 0.4311 0.1918 9 0.7672 1.6013 2.1873 2.4284 2.2882 1.8810 1.4465 1.0472 0.6661 0.2957 10 0.4705 0.9837 1.3447 1.4982 1.4348 1.2189 0.9545 0.6935 0.4415 0.1958 11 0.5420 1.1249 1.5273 1.6762 1.5443 1.2200 0.9054 0.6400 0.4003 0.1766 12 0.4410 0.9162 1.2448 1.3696 1.2824 1.0540 0.7974 0.5651 0.3539 0.1560 13 0.5561 1.1549 1.5699 1.7279 1.6163 1.3251 1.0020 0.7108 0.4458 0.1969 14 0.5624 1.1686 1.5896 1.7472 1.6042 1.2543 0.9282 0.6590 0.4141 0.1833 15 0.5835 1.2135 1.6520 1.8155 1.6193 1.1125 0.7970 0.5683 0.3581 0.1591 16 0.4902 1.0215 1.3929 1.5360 1.3829 0.9737 0.7097 0.5104 0.3233 0.1437 17 0.5088 1.0622 1.4509 1.6088 1.5051 1.2152 0.9270 0.6705 0.4262 0.1890 18 0.6799 1.4206 1.9417 2.1623 2.0653 1.7453 1.3631 0.9898 0.6304 0.2797 19 0.3970 0.8290 1.1339 1.2654 1.2187 1.0450 0.8241 0.5996 0.3820 0.1696 20 0.5503 1.1430 1.5625 1.7114 1.6083 1.3309 1.0121 0.7178 0.4506 0.1990 21 0.4336 0.9013 1.2264 1.3513 1.2656 1.0392 0.7879 0.5604 0.3522 0.1557 22 0.4392 0.9134 1.2443 1.3719 1.2741 1.0229 0.7694 0.5493 0.3465 0.1535 23 0.5766 1.1997 1.6367 1.8092 1.6795 1.3370 1.0081 0.7249 0.4591 0.2039 24 0.6063 1.2640 1.7271 1.9178 1.8087 1.4890 1.1461 0.8301 0.5276 0.2344 25 0.7107 1.4840 2.0294 2.2628 2.1683 1.8416 1.4443 1.0506 0.6694 0.2972 26 0.5042 1.0536 1.4416 1.6112 1.5558 1.3393 1.0598 0.7727 0.4927 0.2187 27 0.5334 1.1083 1.5100 1.6680 1.5696 1.2982 0.9911 0.7090 0.4473 0.1983 28 0.4360 0.9074 1.2387 1.3735 1.3010 1.0861 0.8381 0.6039 0.3825 0.1698 29 0.4586 0.9560 1.3074 1.4560 1.3925 1.1794 0.9219 0.6694 0.4257 0.1890 30 0.4619 0.9651 1.3212 1.4770 1.4264 1.2272 0.9710 0.7088 0.4519 0.2005 31 0.5678 1.1867 1.6252 1.8210 1.7690 1.5362 1.2243 0.8958 0.5718 0.2540 32 0.5621 1.1712 1.6019 1.7869 1.7195 1.4725 1.1590 0.8427 0.5360 0.2382 33 0.5708 1.1918 1.6326 1.8292 1.7778 1.5450 1.2316 0.9010 0.5747 0.2553 34 0.6013 1.2568 1.7227 1.9351 1.8909 1.6561 1.3287 0.9754 0.6234 0.2770 35 0.3622 0.7631 1.0465 1.1768 1.1527 1.0130 0.8156 0.5992 0.3832 0.1704 36 0.5077 1.0620 1.4563 1.6374 1.6036 1.4089 1.1333 0.8330 0.5325 0.2365 37 0.4136 0.8650 1.1867 1.3365 1.3132 1.1590 0.9362 0.6893 0.4411 0.1961Table 6
Relative power error of the HNHEX power distribution solution to DIF3D reference solution for three-dimensional VVER-440 problem (unit: %).
Assembly identity Axial layer 1 (bottom) 2 3 4 5 6 7 8 9 10 (top) 1 −0.96 −0.74 −0.63 −0.33 0.23 0.00 0.00 0.00 0.00 0.00 2 0.46 0.69 0.85 1.10 1.52 2.15 2.61 2.89 2.92 2.84 3 −0.52 −0.27 −0.12 0.04 0.31 0.79 1.03 1.14 1.26 1.12 4 −0.36 −0.10 0.03 0.17 0.41 0.82 1.09 1.12 1.23 1.10 5 0.31 0.54 0.69 0.84 1.07 1.54 1.71 1.85 1.93 1.82 6 0.12 0.36 0.48 0.71 1.08 1.46 1.68 2.10 2.07 2.01 7 −1.08 −0.82 −0.76 −0.50 −0.04 0.00 0.00 0.00 0.00 0.00 8 −1.25 −1.00 −0.91 −0.73 −0.48 −0.27 −0.18 0.15 0.15 0.03 9 −1.62 −1.27 −1.23 −1.14 −1.03 −0.82 −0.86 −0.72 −0.59 −0.85 10 −0.57 −0.21 −0.19 −0.14 −0.06 0.08 0.06 0.14 0.26 0.00 11 −0.50 −0.25 −0.10 0.08 0.36 0.90 1.17 1.30 1.40 1.27 12 0.48 0.72 0.87 1.01 1.27 1.68 1.97 1.99 2.11 1.96 13 −0.69 −0.44 −0.30 −0.17 0.07 0.47 0.71 0.76 0.87 0.71 14 −0.83 −0.59 −0.45 −0.29 −0.05 0.42 0.57 0.73 0.80 0.69 15 −0.94 −0.70 −0.59 −0.37 −0.04 0.27 0.51 0.84 0.83 0.72 16 −0.16 0.07 0.17 0.36 0.65 0.93 1.10 1.39 1.37 1.27 17 0.01 0.31 0.38 0.48 0.63 0.91 0.87 1.05 1.14 0.93 18 −1.06 −0.66 −0.64 −0.58 −0.49 −0.33 −0.36 −0.26 −0.11 −0.43 19 −2.81 −2.52 −2.50 −2.46 −2.38 −2.27 −2.35 −2.22 −2.11 −2.32 20 −0.66 −0.41 −1.11 −0.15 0.09 0.48 0.71 0.66 0.85 0.70 21 0.12 0.35 0.50 0.62 0.85 1.23 1.40 1.48 1.56 1.45 22 −0.11 0.10 0.24 0.37 0.56 0.92 0.98 1.13 1.19 1.07 23 −1.27 −1.05 −0.93 −0.81 −0.64 −0.30 −0.28 −0.13 −0.08 −0.23 24 −1.40 −1.15 −1.06 −0.97 −0.85 −0.59 −0.60 −0.50 −0.42 −0.61 25 −1.68 −1.33 −1.30 −1.25 −1.15 −0.98 −1.01 −0.91 −0.79 −1.06 26 −2.11 −1.77 −1.76 −1.72 −1.64 −1.53 −1.59 −1.48 −1.37 −1.62 27 −1.34 −1.15 −1.00 −0.90 −0.70 −0.35 −0.19 −0.18 −0.11 −0.21 28 −0.34 −0.14 −0.01 0.08 0.24 0.53 0.61 0.67 0.71 0.57 29 −0.38 −0.16 −0.06 0.01 0.14 0.36 0.42 0.46 0.51 0.35 30 0.14 0.44 0.50 0.55 0.65 0.81 0.83 0.88 0.95 0.72 31 −1.90 −1.55 −1.55 −1.50 −1.43 −1.32 −1.35 −1.27 −1.17 −1.44 32 −1.48 −1.25 −1.16 −1.10 −0.96 −0.76 −0.68 −0.66 −0.61 −0.77 33 −1.15 −0.88 −0.83 −0.77 −0.67 −0.52 −0.52 −0.44 −0.40 −0.61 34 −1.68 −1.32 −1.31 −1.27 −1.19 −1.08 −1.07 −1.01 −0.92 −1.22 35 −3.72 −2.66 −2.66 −2.62 −2.55 −2.46 −2.56 −2.40 −2.32 −2.54 36 −0.36 −0.01 0.00 0.04 0.12 0.22 0.21 0.29 0.35 0.08 37 −2.10 −1.76 −1.77 −1.73 −1.66 −1.59 −1.66 −1.54 −1.47 −1.71Figure 7
Radial core map and assemblies identification for three-dimensional VVER 440 benchmark.


## 6. Conclusion

In this work, a three-dimensional, multigroup, diffusion code based on a high order nodal expansion method for hexagonal-z geometry (HNHEX) is developed. In this method, one-dimensional radial and axial spatially flux of each node and energy group are defined as quadratic polynomial expansion and four-order polynomial expansion, respectively. The partially integrated axial leakage is approximated by the quadratic polynomial instead of flat approximation to improve the accuracy for three-dimensional problem. HNHEX has been tested against two popular benchmark problems. The numerical results confirm that the high order nodal expansion method for hexagonal-z geometry can obtain the adequate accuracy for the neutronic calculation of nuclear reactor core with fuel assemblies of hexagonal shape. Besides, obtained results present that the accuracy of four-order polynomial expansion for one-dimensional axial spatially flux and quadratic approximation for partially integrated axial leakage is enhanced relative to quadratic approximation for one-dimensional axial spatially flux and flat approximation for partially integrated axial leakage. The improved coarse-mesh rebalance method effectively accelerates the computation. In future work, HNHEX will continue to be improved to increase the accuracy and calculation speed with our efforts. Furthermore, HNHEX can be extended for calculation of the time dependent nuclear problems in order to conduct simulation of the realistic reactor core neutron kinetic treatments.



---

*Source: 6340652-2016-10-17.xml*