Abstract
A three-dimensional, multigroup, diffusion code based on a high order nodal expansion method for hexagonal- geometry (HNHEX) was developed to perform the neutronic analysis of hexagonal- 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.
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- 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- 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- 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- 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- 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. Section 2 presents the summary of high order nodal expansion method for hexagonal- 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- Geometry
In this section, we present a high order nodal expansion method for hexagonal- 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- 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- geometry.
The hexagonal axes which consist of x-direction, u-direction, and v-direction are employed in the hexagonal geometry as shown in Figure 1. The u-direction is rotated 60° counterclockwise with respect to the -axis and the v-direction is rotated 120° from the -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 of kth node can be defined as shown in (1) in the local coordinate:where is the lattice pitch and the -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.
2.2. High Order Nodal Expansion Method for Hexagonal- Geometry
In this method, reactor is divided into number of homogeneous hexagonal nodes. Cross sections within the nodes are assumed to be constant. Let us consider a homogeneous node shown in Figure 2, the multigroup neutron diffusion equation for th node, and th group can then be written in the form shown inwhere
The solution of (3) requires additional relationships between the surface-averaged leakages and the nodal fluxes in the kth 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 over y and z directions and take the form shown inwhere
The one-dimensional axial equation can be obtained by integrating the three-dimensional equation over x and y directions and take the form shown inwhere
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 inwhere and are obtained by using orthogonal characteristics of polynomials. Higher order approximations are obtained by first adding a basis function which has a first derivative discontinuity at and then adding an additional basis function which provides a quadratic approximation over the half-intervals and .
The approximation of the partially integrated axial fluxes takes the form shown inwhere
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
We solve (13) to get the zeroth expansion coefficients shown in
The coefficients are determined by requiring the y-integrated net current to be continuous at and the imposed condition is shown in
Apply the two-step leakage approximation and the final form of is given bywhere
Higher order coefficients 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 and require the result to be zero when integrated over the interval as shown in
For moment weighting, ; integrate and rewrite (18) to get flux moment equation in three radial directions shown inwhere
Now substitute in with (9) and then integrate to get as shown in
Higher order coefficients and 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 and requiring the result to be zero when integrated over the interval
For moment weighting, , , integrate and rewrite (23) to get one-order flux moment equation in axial direction shown inwhere
For moment weighting, , , integrate and rewrite (23) to get second-order flux moment equation in axial direction shown inwhere
Substitute into and with (11) and then integrate to get (28a) and (28b):
Thus, and are shown in
The 3D multigroup nodal balance equation is obtained by integrating over a homogeneous node volume where is derived from (30) and shown in (32). , , and are derived from (19) and shown in (33), (34), and (35). and are derived from (24) and (26) and shown in (36) and (37):
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
By using (9) and (38), outgoing partial currents through radial-directed surfaces of kth node in terms of flux expansion coefficients are shown inwhere
By using (11) and (38), outgoing partial currents through axial-directed surfaces of kth node in terms of flux expansion coefficients are shown in
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.whereMatrix elements are given below. These elements are characterized by node dimensions and homogenized material of the node.
3. Calculation of the Leakage Moments
Leakage moments in the response matrix are derived in this section. The partially integrated axial leakage and partially integrated radial leakage 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
The x-direction moment of the partially integrated axial leakage is calculated by using the approximation shown in (45) in which and are as given in (10). The partially integrated axial leakage is given by (45) and the total axial leakage is shown in (46):where
The coefficients and are calculated in the following manner. Let and denote the neighboring nodes in the minus and plus x-directions and the quadratic polynomial extends over the three nodes , , and . The expansion coefficients in (45) are calculated such that the total axial leakages shown in (46) in the nodes and are preserved as shown in (47). The expansion coefficients and are calculated by using the constraints given in (47) and the results are shown in (48):
The required leakage moment is calculated by substituting (45) into (20) and performing the necessary integrations and the total axial leakages and shown in (46) can be computed by axial direction net currents. The final form of the x-direction moment of the partially integrated axial leakage is shown in (49). Similar equation for and can be obtained as shown in (50) and (51):
3.2. The Axial Moment of the Partially Integrated Radial Leakage
The axial moment of the partially integrated radial leakage is calculated by using the approximation equation (52). and are given in (12). The total radial leakage can be written in terms of the average leakages in the three radial directions as shown in (53):where
The coefficients and are calculated in the following manner. Let and denote the neighboring nodes in the minus and plus z-directions and the quadratic polynomial extends over the three nodes , , and . Integrating (54) over and , respectively, and the expansion coefficients and is calculated by using the constraints given in (54) and the results are shown in (55).
The required leakage moment is calculated by substituting (52) into (25) and performing the necessary integrations to get the final form for shown inwhere
The required leakage moment is calculated by substituting (52) into (27) and performing the necessary integrations to get the final form for shown inwhere
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 Figure 3 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 , eigenvalue , rebalance factors and partial currents vector and to initial value (assuming , ) for all energy groups and all nodes . The initial total fission source vector for kth node is calculated by(3)At the beginning of outer iteration n, the group partial currents are scaled by rebalance factors calculated by coarse-mesh rebalance procedure as shown in(4)Fission source vector for kth nodes and th group, which is defined in (62), is calculated from flux moment vector and eigenvalue obtained in previous outer iteration (or, it will be guess value if it is the first outer iteration)(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 Then, construct group source vector for kth node and th group due to fission source vector and scattering source vector by(6)Leakage moment vector for kth nodes and th 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) for kth nodes and th 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 (), and leakage moment vector calculated in step ().(9) Steps ()–() are repeated for all nodes and all energy groups with keeping fission source vector calculated in step () 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 for kth nodes is scaled by rebalance factors prior to checking the convergence of the fission source as shown in(11)A new eigenvalue is calculated from , the total fission source obtained in previous outer iteration , and the total fission source obtained in present outer iteration n as shown in where The outer iterations are terminated when three convergence criteria shown in (68) are satisfied:
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 inwhere 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 [] 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) towhere
The convergence rate of the power method is determined by the dominance ratio of the matrix ; the closer this ratio is to 1, the slower the convergence rate will be. It can be shown that, for , the dominance ratio of is smaller than that of . 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
Numerical calculations to date have demonstrated that only two power iterations are required at each power iteration ( 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 because, as shown in (65), they do not enter into the calculation of the group source term due to the assumption that there is no upscatter. The final estimate for the eigenvalue at the th outer iteration is obtained from (73). The coarse-mesh rebalance is performed once every five outer iterations:
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- 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].
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 Table 1 shows the numerical results of 3D SNR problem including effective multiplication factor and its relative percentage error ; relative percentage error of region-averaged fluxes in comparison to reference solution; and the CPU times. , , , , and 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. DIF3D and DIF3D 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 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 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.
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.
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 Table 4, 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.
6. Conclusion
In this work, a three-dimensional, multigroup, diffusion code based on a high order nodal expansion method for hexagonal- 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- 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.
Nomenclature
| One-dimensional radial spatially flux for th nodes and th group | |
| One-dimensional axial spatially flux for th nodes and th group | |
| Lattice pitch of the hexagonal assembly | |
| Axial height of th node | |
| Macroscopic fission cross section of th node and th group | |
| Macroscopic scattering cross section from group to group | |
| Macroscopic removal cross section of th node and th group | |
| Diffusion coefficient for th node and th group | |
| Partially integrated axial leakage for th nodes and th group | |
| Partially integrated radial leakage for th nodes and th group | |
| One-dimensional face-averaged spatially flux on the -directed face and for th node and th energy group | |
| Left /right -surface of node , , (the left surface is , , , and ; the right surface is , , , and ) | |
| Face-averaged net current on the -directed face and for th node and th energy group | |
| Total fission source vector for th node | |
| Fission source vector for th node and th group | |
| Scattering source vector for th node and th group | |
| Group source vector for th node and th group | |
| Leakage moment vector for th nodes and th group | |
| Volume-averaged flux for th nodes and th group. |
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by the Fundamental Research Funds for the Central Universities (JB2016163).