Figure 1 illustrates the balance for neutron flux within a node of interest i, and its associated balance equation is expressed as follows: J ¯ i , i + 1 − J ¯ i − 1 , i + Σ r i Δ x i ϕ ¯ i = 1 k ν ∑ f i Δ x i ϕ ¯ i (1) where J ¯ i , i + 1 denotes the net neutron current from node i toward its adjacent node i+1, k represents the multiplication factor, and all the other notations are those of the convention. By approximating the gradient of the flux to be linear in the diffusion theory, one garners the following equation which is known as the finite difference method (FDM): J ¯ i , i + 1 = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) (2) D ˜ i , i + 1 = 2 ( D i Δ x i ) ( D i + 1 Δ x i + 1 ) D i Δ x i + D i + 1 Δ x i + 1   (3)

The correction factor D ^ i , i + 1 can be envisioned for Eq. 2 which alters the net current to be the reference current as follows: J ¯ i , i + 1 r e f = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) − D ^ i , i + 1 ( ϕ ¯ i + 1 + ϕ ¯ i ) , (4) D ^ i , i + 1 = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) − J ¯ i , i + 1 r e f ( ϕ ¯ i + 1 + ϕ ¯ i ) . (5)

The neutron balance equation, i.e., Eq. 1, is known to be a well-defined problem, in which a single solution (eigenvalue) exists for a certain current value. Hence, the inclusion of Eqs 4, 5 while implementing the FDM would provide an equivalent solution to that of the reference, which is exploited for acquiring J ¯ i , i + 1 r e f (Smith, 1983). The simplicity of the FDM with its coarse node size will manifest as a reduction in the computing burden (acceleration), hence attaining the name of coarse-mesh finite difference (CMFD) method.

It is apparent that the preservation of partial currents will guarantee retaining net current. From such a perspective, two correction factors for incoming and outgoing partial currents can be considered as depicted in Figure 2, attaining the name of partial current–based CMFD (pCMFD) method (Cho et al., 2003): J ¯ i , i + 1 + r e f = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) + 2 D ^ i , i + 1 + ϕ ¯ i 2 (6) J ¯ i , i + 1 − r e f = D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) + 2 D ^ i , i + 1 − ϕ ¯ i + 1 2 (7)

The net current can be expressed as J ¯ i , i + 1 r e f = J ¯ i , i + 1 + r e f − J ¯ i , i + 1 − r e f = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) − ( D ^ i , i + 1 + ϕ ¯ i + D ^ i , i + 1 − ϕ ¯ i + 1 ) (8) where D ^ i , i + 1 + = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) + 2 J ¯ i , i + 1 + r e f 2 ϕ ¯ i (9) D ^ i , i + 1 − = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) − 2 J ¯ i , i + 1 − r e f 2 ϕ ¯ i + 1 (10)

The correction factors for the preservation of net current and surface flux could be introduced separately for each node as shown in Figure 3, where incoming partial current acts as a boundary condition for invoking kernel calculation (Shin and Kim, 1999). Such an approach differs from the conventional CMFD and pCMFD methods which introduce correction alongside two contiguous nodes: J ¯ i , i + 1 = − 2 D i Δ x i ( ϕ ¯ s − ϕ ¯ i ) − 2 D ^ i R Δ x i ( ϕ ¯ s + ϕ ¯ i ) (11) J ¯ i , i + 1 = − 2 D i + 1 Δ x i + 1 ( ϕ ¯ i + 1 − ϕ ¯ s ) − 2 D ^ i + 1 L Δ x i + 1 ( ϕ ¯ s + ϕ ¯ i + 1 ) (12)

Equating Eqs 11, 12 yields the following expression for the surface flux: ϕ ¯ s = Δ x i + 1 ( D i − D ^ i R ) ϕ ¯ i + Δ x i ( D i + 1 + D ^ i + 1 L ) ϕ ¯ i + 1 Δ x i + 1 ( D i + D ^ i R ) + Δ x i ( D i + 1 − D ^ i + 1 L ) (13) where D ^ i R = − Δ x i   J ¯ i , i + 1 r e f + 2 D i ( ϕ ¯ s r e f − ϕ ¯ i ) 2 ( ϕ ¯ s r e f + ϕ ¯ i ) (14) D ^ i + 1 L = − Δ x i + 1   J ¯ i , i + 1 r e f + 2 D i + 1 ( ϕ ¯ i + 1 − ϕ ¯ s r e f ) 2 ( ϕ ¯ s r e f + ϕ ¯ i + 1 ) (15)

Originally devised for parallel acceleration, the aforementioned acceleration scheme, which is referred to as one-node CMFD, can still be applied in a similar manner to that of the conventional CMFD or pCMFD method. Note that, for such an implementation, preservation of net current becomes irrelevant to surface flux values.

The multigroup diffusion equation can be written as follows: M g ϕ g = χ g k ∑ g ′ = 1 G ν Σ f , g ′ ϕ g ′ + ∑ g ′ = 1 ( g ′ ≠ g ) G Σ s , g ′ → g ϕ g ′ (19) where G and k represent the number of energy groups and multiplication factor, M g denotes the migration operator for group g ( M g ϕ g : = − ∇ ⋅ D g ∇ ϕ g + Σ r , g ϕ g ), and all the other notations are those of the convention. Note that, through proper discretization, the given equation can be represented in a matrix form. The numerical adjoint flux can then be calculated by transposing Eq. 19: M g † ϕ g † = 1 k ∑ g ′ = 1 G χ g ′ ν Σ f , g ϕ g ′ † + ∑ g ′ = 1 ( g ′ ≠ g ) G Σ s , g → g ′ ϕ g ′ † (20) where superscript dagger ( † ) signifies the adjoint operation.

Mathematically, it could be shown that the continuous migration operator for a certain group is self-adjoint (Ott and Neuhold, 1985). Such a feature is retained for the FDM approach, however, but not for the CMFD accelerated nodal calculation due to the presence of correction factors. The discretized balance equation can be generalized as follows: a i , i − 1 ϕ ¯ i − 1 + a i , i ϕ ¯ i + a i , i + 1 ϕ ¯ i + 1 = s i (21) where a i , j represents the contribution from ϕ ¯ j for the neutron balance concerning ϕ ¯ i and s i denotes the source term for node i containing both scattering and fission. Table 1 summarizes the matrix entries for the FDM and CMFD methods, where only the FDM approach retains the self-adjoint of the migration matrix, i.e., a i − 1 , i =   a i , i − 1 and a i , i + 1 =   a i + 1 , i .

It could be observed that all the enumerated CMFD methods do not retain the self-adjoint property for the group-wise migration matrix. In addition, for the conventional CMFD and one-node CMFD, it could be easily shown that the absolute magnitude of corrections factors will dwindle with a decrease in the mesh size, i.e., numerators for Eqs 5, 14, and 15 converge to zero: − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) →   J ¯ i , i + 1 r e f (22) − D i Δ x i / 2 ( ϕ ¯ s r e f − ϕ ¯ i ) → J ¯ i , i + 1 r e f (23)

However, correction factors for pCMFD do not converge to zero like the other two CMFD-based acceleration schemes: D ^ i , i + 1 + = − D ˜ i , i + 1 ( ϕ ¯ i + 1 − ϕ ¯ i ) + 2 J ¯ i , i + 1 + r e f 2 ϕ ¯ i ≅ J i , i + 1 r e f − 2 J ¯ i , i + 1 + r e f 2 ϕ ¯ i = ( J ¯ i , i + 1 + r e f − J ¯ i , i + 1 − r e f ) − 2 J ¯ i , i + 1 + r e f 2 ϕ ¯ i = − J ¯ i , i + 1 + r e f + J ¯ i , i + 1 − r e f 2 ϕ ¯ i . (24)

The inclusion of correction factors can be expressed in the following manner: ( M + δ M 1 + δ M 2 ) ϕ C M F D = 1 k F ϕ C M F D (25) where ϕ C M F D represents the CMFD-based flux and δ M 1 and δ M 2 denote diagonal and off-diagonal correction entries, respectively. The transpose of Eq. 25 can be written as ( M † + δ M 1 † + δ M 2 † ) ϕ C M F D † = 1 k F † ϕ C M F D † (26)

Since the diagonal matrix is self-adjoint, i.e., δ M 1 = δ M 1 † , subtraction of the two equations above garners ( M + δ M 1 ) ( ϕ C M F D − ϕ C M F D † ) + ( δ M 2 ϕ C M F D − δ M 2 † ϕ C M F D † ) = 1 k F ( ϕ C M F D − ϕ C M F D † ) . (27)

If δ M 2 ≅ δ M 2 † is satisfied, ϕ C M F D = ϕ C M F D † becomes the solution for Eq. 27, which corresponds to the preservation of self-adjoint feature.

Since CMFD-induced and one-node CMFD–induced correction factors converge to zero with an increase in the number of nodes, their corresponding numerical adjoints would also converge to the reference, which is not expected for the pCMFD-based numerical adjoint flux.

The correction factors in the discretized balance equation could result in the occurrence of negative numerical adjoint flux values as pointed out in previous studies which implemented the analytic nodal method (ANM) while deducing correction factors (Müller, 2014). Such an anomaly ensues when the off-diagonal and its corresponding diagonal entry of the migration matrix attain the same sign, which cannot be prevented for the conventional CMFD method. To circumvent such an issue, a different formula for net current preservation can be partially utilized under certain conditions; however, such an approach deteriorates the consistency in the CMFD formulation, i.e., ad hoc up to a certain extent.

Recalling that the usage of one-node CMFD in a two-node manner, i.e., not parallelized, could preserve the net current regardless of its surface flux values, one could exclude the occurrence of negative adjoint flux values through proper adjustment of the surface flux values. It is noteworthy to articulate that consistent usage of the same surface flux value while formulating the correction factors is responsible for the preservation of the current.

From Eqs 14, 15, it could be recognized that the one-node CMFD correction factors have the same unit as the diffusion coefficient. Since it is unphysical for the correction factor–included diffusion coefficient to be negative, the following conditions can be envisioned: ϕ ¯ s r e f ≤ Δ x i 4 D i | J ¯ i , i + 1 r e f | (28) ϕ ¯ s r e f ≤ Δ x i + 1 4 D i + 1 | J ¯ i , i + 1 r e f | (29)

Eqs 28, 29 represent the criterion for correction factors D ^ i R and D ^ i + 1 L being less than their associated diffusion coefficient in magnitude, respectively. Through adjustment of surface flux to suffice Eqs 28, 29, the occurrence of negative numerical adjoint flux can be stifled.

To test the deviation in the self-adjoint feature for various CMFD-based numerical adjoint fluxes, a simple one-group two-dimensional reactor problem was considered as shown in Figure 4. Two types of assemblies with a length of 20 cm are considered as shown in the cartoon, where their cross-section (XS) values are enumerated in Table 2.

For the acquisition of numerical adjoints, the nodal expansion method (NEM) kernel was utilized, whereas the reference adjoint flux was obtained via the FDM while dividing each assembly into equally spaced 2,500 nodes (50 × 50 per assembly). Since the FDM-based numerical adjoint always retains the self-adjoint feature, the acquired fine node-based result was regarded as a reference after condensing into an assembly-wise value according to the following equation (Downar et al., 2009): ∫ 0 ∞ d E ∫ V d V ϕ † ( r , E ′ ) ϕ ( r , E ′ )   = 1.0 (30)

The acquired multiplication factors for both forward and adjoint calculations are given in Table 3, where all the cases exhibit the same value.

Figure 5 illustrates the acquired adjoint fluxes from each acceleration scheme, where normalization according to Eq. 30 was performed for comparison. It could be recognized that only the pCMFD-based numerical adjoint flux exhibits a different distribution. Figure 6 summarizes the calculation result where the absolute value of percentage error for each case is given for both forward and adjoint fluxes. Note that the numerical adjoint flux exhibits the most conspicuous error for the TYPE 2 assembly region due to its enlarged absorption XS value: E R R ( % ) = | ϕ r e f − ϕ C M F D ϕ r e f | × 100 . (31)

As aforementioned, the deviation from the self-adjoint feature weakens as the size of the node dwindles. A similar analysis was performed by imposing a node size of 10 cm, where reduction in the adjoint flux error is observed for both CMFD and one-node CMFD (1NCMFD) cases as depicted in Figure 7.

The KWU PWR 2D benchmark problem has been considered while excluding the presence of soluble boron under fully rodded conditions (Benchmark Source Situation, 1985). The configuration of the reactor problem is given in Figure 8 alongside four different positions for imposing localized perturbation in the absorption XS.

The attainment of reference adjoint flux was done in a similar fashion to that of one-group reactor problem while dividing each assembly into equally spaced 10,000 nodes (100 × 100 per assembly). Table 4 juxtaposes the calculated multiplication factors, where the same values are obtained regardless of the CMFD acceleration schemes as expected. Note that each assembly was taken as a single node during the CMFD accelerated nodal calculation.

The acquired numerical adjoint fluxes for both fast and thermal groups are shown in Figures 9, 10, where pCMFD-based results do not conform with the other results. In addition, negative adjoint flux values (red color) are observed for the thermal group adjoint flux in the peripheral regions as depicted in Figure 10, where the surface flux attained from the NEM kernel calculation was directly utilized for one-node CMFD acceleration, i.e., no correction was made for acquisition of correction factors.

In order to stifle the occurrence of negative adjoint flux as shown in Figure 10, the surface flux was adjusted as follows to alter the correction factor to be zero when one of Eqs 28, 29 is met during a one-node CMFD calculation. The resulting numerical adjoint is illustrated in Figures 11, 12, where no negative values are observed: ϕ ¯ s r e f = − 1 2 ⋅ Δ x i   J ¯ i , i + 1 r e f D i + ϕ ¯ i     f o r   D ^ i R (32) ϕ ¯ s r e f = 1 2 ⋅ Δ x i + 1   J ¯ i , i + 1 r e f D i + 1 + ϕ ¯ i + 1   f o r   D ^ i + 1 L (33)

For a systematic comparison between the numerical adjoint fluxes, the first-order perturbation theory was utilized, where assessment in the change of reactivity was made and compared with the reference value. Note that the reference reactivity change was evaluated via a direct solution of the perturbed system. The reactivity change can be estimated as follows: Δ ρ = ∫ 0 ∞ ∫ V d E d V ϕ 0 † ( r , E ) ( λ 0 Δ F − Δ M ) ϕ 0 ( r , E ) ∫ 0 ∞ ∫ V d E d V ϕ 0 † ( r , E ) F 0 ϕ 0 ( r , E ) (34) where λ denotes the reciprocal of multiplication factor, F and M represent augmented fission and migration operators, respectively, and all the other notations are those of the convention.

Four different perturbation scenarios are envisaged as shown in Figure 8, where a change in the absorption XS was locally imposed (yellow colored assemblies). The extent of variation in the XS compared to its original value was set to be 30% for cases 1 to 3 and 80% for case 4. Two different adjoint fluxes are considered for the one-node CMFD method, namely, the original and the adjusted one. Note that the former result is subjected to a negative value issue. The calculated results are enumerated from Tables 5–8.

Where the asterisk denotes the adjusted one-node CMFD numerical adjoint and UNIFORM indicates the usage of unit vector while estimating reactivity change according to the first-order perturbation theory. It could be observed that the exploitation of one-node CMFD–based adjoint flux renders the estimation to be more accurate compared to other approaches. Especially, for the perturbation in the reflector region, i.e., case 4, only the negative adjoint flux issue–resolved one-node CMFD exhibits a reliable result.