A Revisit to CMFD Schemes: Fourier Analysis and Enhancement

Abstract

The coarse-mesh finite difference (CMFD) scheme is a very effective nonlinear diffusion acceleration method for neutron transport calculations. CMFD can become unstable and fail to converge when the computational cell optical thickness is relatively large in k-eigenvalue problems or diffusive fixed-source problems. Some variants and fixups have been developed to enhance the stability of CMFD, including the partial current-based CMFD (pCMFD), optimally diffusive CMFD (odCMFD), and linear prolongation-based CMFD (lpCMFD). Linearized Fourier analysis has proven to be a very reliable and accurate tool to investigate the convergence rate and stability of such coupled high-order transport/low-order diffusion iterative schemes. It is shown in this paper that the use of different transport solvers in Fourier analysis may have some potential implications on the development of stabilizing techniques, which is exemplified by the odCMFD scheme. A modification to the artificial diffusion coefficients of odCMFD is proposed to improve its stability. In addition, two explicit expressions are presented to calculate local optimal successive overrelaxation (SOR) factors for lpCMFD to further enhance its acceleration performance for fixed-source problems and k-eigenvalue problems, respectively.

1. Introduction

The coarse-mesh finite difference (CMFD) scheme, originally developed by Smith [1,2], is being widely used for accelerating neutron transport calculations. A well-known issue with CMFD is that the scheme can become unstable and even fail when the optical thickness becomes large, which has been reported in various numerical and theoretical studies. Some stabilization techniques have been developed to improve the stability of CMFD [3,4,5,6]. One simple method is the under-relaxation approach, which arbitrarily applies some damping on the drift flux coefficient or the scalar flux [4,5,6]. However, a judicious choice of under-relaxation has to be problem-dependent.

A variant of CMFD, called partial current-based CMFD (pCMFD), is found to be stable, but it is slower than the standard CMFD at an intermediate optical thickness [3]. Zhu et al. recently proposed an optimally diffusive CMFD (odCMFD) scheme, which generalizes pCMFD by adding some artificial diffusion to the diffusion coefficient of CMFD [7]. It is noted that caution should be exercised in determining adequate artificial diffusion. If the added diffusion is more than necessary, it would hurt the convergence rate; conversely, it would be insufficient for stability. In this paper, Fourier analysis and numerical results are presented to show that odCMFD can become unstable and even fail at a relatively large optical thickness because of insufficient diffusion. In this paper, it is shown that the use of different transport solvers in Fourier analysis can have some effect on the prediction of the convergence rate. A modification to the original artificial diffusion coefficients of odCMFD is proposed to improve its stability.

The linear prolongation-based CMFD (lpCMFD) scheme is the latest development in the CMFD family [8]. It replaces the conventional flat flux ratio-based scaling approach in the standard CMFD method with a linear interpolation of the scalar flux differences at the coarse-mesh cell edges between the high-order transport and low-order diffusion calculations. Various studies have demonstrated that lpCMFD has superior convergence performance in terms of stability and effectiveness, as compared to other CMFD schemes [9,10,11,12,13,14]. It is worth noting that lpCMFD does not use any “ad hoc” fixups, such as under-relaxation or artificial diffusion. Most recently, Wang has proposed to enhance the convergence performance of lpCMFD by updating the scalar flux (or the scattering source) with successive overrelaxation (SOR) [15,16]. It was found that lpCMFD with SOR can increase the convergence rate by almost 40% for neutron transport problems. However, the global optimal SOR value is problem-dependent. Instead of finding the global optimal SOR factor, the local SOR can be determined for each coarse-mesh cell based on its optical thickness and scattering ratio. In this paper, two explicit expressions are developed for the calculation of local optimal SOR factors—one is for fixed-source problems, and the other for k-eigenvalue problems, which can be easily implemented in a transport code.

The remainder of this paper is organized as follows: Section 2 presents the Fourier analysis of the CMFD schemes, such as CMFD, pCMFD, odCMFD, and lpCMFD, for both fixed-source problems and k-eigenvalue problems. The implications on the convergence rate of the step characteristic (SC) method and the diamond difference (DD) method for the transport solution are presented and discussed. In Section 3, the modification of odCMFD is discussed. Section 4 presents the two explicit expressions for the calculation of local SOR factors for lpCMFD. The numerical results of the one- and two-dimensional test problems are given in Section 5. Final conclusions are given in Section 6. A brief summary of the CMFD algorithm and formulation is presented in Appendix A.

2. Fourier Analysis of CMFD Schemes

The general algorithm and formulation of the coupled S_N and CMFD iterative procedure is summarized in Appendix A. For simplicity, only the final Fourier analysis results of the CMFD schemes are presented for both fixed-source problems and k-eigenvalue problems. Detailed derivation and formulation can be found in the literature [6,7,8,15,16].

2.1. Spectral Radius for Fixed-Source and K-Eigenvalue Problems

Fourier analysis is typically performed based on a 1D homogeneous problem with a uniform mesh and a specified number $\rho$ of fine-mesh cells per coarse cell. A single coarse-mesh cell with periodic boundary conditions is employed. The coarse-mesh spacing is $\Delta$. The fine mesh spacing is $h=\frac{\mathsf{\Delta }}{p}$. The monoenergetic S_N transport equation is first linearized near the exact solution. Then, the error terms in the linearized equations are expressed as Fourier modes. The spectral radius can be obtained from the error iteration matrix.

For fixed-source problems, the final expression for the fine-mesh scalar flux iterative error modes is given as

$$\mathit{T}\left(\omega \right)=\mathit{H}+\mathit{\theta }\frac{∆\left(e^{i{\mathsf{\Sigma }}_t\mathsf{\Delta }\omega }-1\right)\mathit{G}+{\mathsf{\Sigma }}_a{\mathsf{\Delta }}^2\mathit{C}}{D\left[2\cos \left({\mathsf{\Sigma }}_t\mathsf{\Delta }\omega \right)-2\right]-{\mathsf{\Sigma }}_a{\mathsf{\Delta }}^2},$$

(1)

where

$$\mathit{H}={\displaystyle \sum }_{n=1}^N{w}_n{\mathit{P}}_n{Y}_n^{-1}\frac{{\mathsf{\Sigma }}_s}{2}$$

(2)

$$\mathit{C}=\frac{1}{p}\times {\left[\begin{array}{ccc}1& \cdots & 1\end{array}\right]}_{1\times p}\mathit{H}$$

(3)

$$\mathit{G}={\left[\begin{array}{cc}\begin{array}{cc}1& 0\end{array}& \begin{array}{cc}\cdots & 0\end{array}\end{array}\right]}_{1\times p}\mathit{H}$$

(4)

$${\mathit{P}}_n={\left[\begin{array}{ccccc}\frac{1-{\alpha }_n}{2}& \frac{1+{\alpha }_n}{2}& 0& 0& 0\\ 0& \frac{1-{\alpha }_n}{2}& \frac{1+{\alpha }_n}{2}& 0& 0\\ 0& 0& 0& \ddots & 0\\ \frac{1+{\alpha }_n}{2}{e}^{i{\mathsf{\Sigma }}_t\mathsf{\Delta }\mathsf{\omega }}& 0& 0& 0& \frac{1-{\alpha }_n}{2}\end{array}\right]}_{p\times p}$$

(5)

$${\mathit{Y}}_n={\left[\begin{array}{ccccc}a& b& 0& 0& 0\\ 0& a& b& 0& 0\\ 0& 0& 0& \ddots & 0\\ b{e}^{i{\mathsf{\Sigma }}_t\mathsf{\Delta }\mathsf{\omega }}& 0& 0& 0& a\end{array}\right]}_{p\times p}$$

(6)

$$a=-\frac{{\mu }_n}{h}+{\mathsf{\Sigma }}_t\left(\frac{1-{\alpha }_n}{2}\right)$$

(7)

$$b=\frac{{\mu }_n}{h}+{\mathsf{\Sigma }}_t\left(\frac{1+{\alpha }_n}{2}\right)$$

(8)

$${\alpha }_n=\left\{\begin{array}{ll}\frac{1+e^{-\frac{{\mathsf{\Sigma }}_{t}h}{{\mu }_n}}}{1-e^{-\frac{{\mathsf{\Sigma }}_{t}h}{{\mu }_n}}}-\frac{2{\mu }_n}{{\mathsf{\Sigma }}_{t}h}& \mathrm{for} \mathrm{SC}\\ 0& \mathrm{for} \mathrm{DD}\end{array}\right.$$

(9)

$$h=\frac{\mathsf{\Delta }}{p}$$

(10)

$$D=\frac{1}{3{\mathsf{\Sigma }}_t}$$

(11)

Note that ${\mu }_n$ and $w_n$ are the nodes and weights of a quadrature set (e.g., Gauss–Legendre quadrature).

The vector $\mathit{\theta }$ is defined as

$$\mathit{\theta }={({\theta }_1 {\theta }_2\cdots {\theta }_p )}^T$$

(12)

The following different iterative schemes can be obtained by defining the vector $\mathit{\theta }$ as

$$\mathbf{SI}: {\theta }_k=0.$$

(13)

$$\mathbf{CMFD}:{\theta }_k=1.$$

(14)

$$\mathbf{pCMFD}:{\theta }_k=\frac{1}{1+\frac{3}{4}\tau }.$$

(15)

$$\mathbf{odCMFD}: {\theta }_k=\theta \left(\tau \right).$$

(16)

$$\mathbf{lpCMFD}: {\theta }_k= \frac{2p+\left(2p-2k+1\right)e^{-i\mathsf{\tau }\omega }+\left(2k-1\right)e^{i\mathsf{\tau }\omega }}{4p}.$$

(17)

where $1\le k\le p$, and the coarse-mesh optical thickness $\tau ={\mathsf{\Sigma }}_t\mathsf{\Delta }$.

It can be seen that when ${\theta }_k=0$, it reverts to the transport source iteration (SI) method with no CMFD acceleration. For the standard CMFD scheme, ${\theta }_k=1$, which means a flat prolongation approach is used. The function of ${\theta }_k$ for pCMFD was derived by Jarret et al. [6], and it implies that pCMFD tends to the CMFD scheme when the coarse-mesh cell optical thickness $\tau$ ($={\mathsf{\Sigma }}_t\mathsf{\Delta }$) becomes large. The ${\theta }_k$ of odCMFD is a function of $\tau$, described by Zhu et al. [7]. The ${\theta }_k$ values of lpCMFD are reported by Wang and Xiao [8].

The spectral radius of the CMFD schemes as a function of Fourier frequency $\omega$ is given as

$$\rho \left(\omega \right)=\max (\mathrm{abs}(eig(\mathit{T}\left(\omega \right)))).$$

(18)

where Fourier frequency is expressed as follows

$$\omega =\frac{2\pi s}{{\mathsf{\Sigma }}_{t}L},s=1, 2, \dots , J,$$

(19)

$$J=\frac{L}{\mathsf{\Delta }}.$$

(20)

$L$ is the domain size. Finally, the spectral radius is defined as

$$\rho =\max \left(\rho \left(\omega \right)\right).$$

(21)

For k-eigenvalue problems, the final expression for the fine-mesh scalar flux iterative error modes is given as

$$\mathit{T}\left(\omega \right)=\mathit{H}+\mathit{\theta }\frac{∆\left(e^{i{\mathsf{\Sigma }}_t\mathsf{\Delta }\omega }-1\right)\mathit{G}}{D\left[2\cos \left({\mathsf{\Sigma }}_t\mathsf{\Delta }\omega \right)-2\right]},$$

(22)

where

$$\mathit{H}={{\displaystyle \sum }}_{n=1}^N{w}_n{\mathit{P}}_n{Y}_n^{-1}\frac{{\mathsf{\Sigma }}_t}{2},$$

(23)

$$\mathit{G}={\left[\begin{array}{cc}\begin{array}{cc}1& 0\end{array}& \begin{array}{cc}\cdots & 0\end{array}\end{array}\right]}_{1\times p}\mathit{H}.$$

(24)

Compared to the expression for the fixed-source case, only matrices $\mathit{H}$ and $\mathit{G}$ are different, in which the scattering cross-section is replaced with the total cross-section. Therefore, the spectral radius is independent of scattering for k-eigenvalue problems. It is interesting to note that the spectral radius of fixed-source problems will be limited to that of k-eigenvalue problems, as the scattering ratio $c\to 0$ (pure scattering).

2.2. The Implication of Transport Solver

In Fourier analysis, it is customary to use the SC method to discretize the S_N transport equation, in which the transport angular flux on each computational cell is weighted as

$${\psi }_{n,k}=\left(\frac{1-{\alpha }_n}{2}\right){\psi }_{n,k-\frac{1}{2}}+\left(\frac{1+{\alpha }_n}{2}\right){\psi }_{n,k+\frac{1}{2}},$$

(25)

where $1\le k\le p$. ${\psi }_{n,k-\frac{1}{2}}$ and ${\psi }_{n,k+\frac{1}{2}}$ are the left and right edge fluxes of the $k-th$ cell. Then,

$${\alpha }_n=\frac{1+e^{-\frac{{\mathsf{\Sigma }}_{t}h}{{\mu }_n}}}{1-e^{-\frac{{\mathsf{\Sigma }}_{t}h}{{\mu }_n}}}-\frac{2{\mu }_n}{{\mathsf{\Sigma }}_{t}h}.$$

(26)

Remark 1.

If the cell optical thickness${\Sigma }_{t}h\to 0$, then${\alpha }_n\to 0$, thus the SC solution tends to the DD solution; if${\Sigma }_{t}h\to \infty$, then it tends to the first-order upwind method (i.e., the step difference method) because of${\alpha }_n\to 1$for${\mu }_n>0$, and${\alpha }_n\to -1$for${\mu }_n<0$. Therefore, there should be no difference in the spectral radius predicted by Fourier analysis when using the SC or the DD method at a thin optical thickness. However, the difference will become appreciable when the coarse-mesh optical thickness$\tau >1$. This is confirmed by the Fourier analysis as shown in Figure 1. Note that the results are obtained with the domain size of$L=200$cm, the coarse-mesh size$\mathsf{\Delta }=1$cm, and the number of fine cells in each coarse cell$p=10$.

The following observations can be drawn from the above Fourier analyses:

• For CMFD, the results of the SC and DD are very close (the black and blue curves are almost coincident), and both solvers predict failure when $\mathit{\tau }>1.$
• For pCMFD, the results of the SC and DD agree very well within the differences of less than 0.006, except that the DD predicts failure when $\mathit{\tau }>90.$
• For odCMFD, significant differences show up when $\mathit{\tau }>2.5.$ The DD predicts failure at $\mathit{\tau }>5.5$, while the SC results show that the scheme is more stable and becomes unstable until $\mathit{\tau }>20.$ It is worth mentioning that the SC solution is somewhat sensitive to the domain size. If a smaller domain is used, e.g., $\mathit{L}=100$ cm, then there is no occurrence of a failure in the SC solution.
• For lpCMFD, the SC results show that the scheme becomes unstable when $\mathit{\tau }>27$, but interestingly the DD does not predict any failure. Numerical testing indicates that the lpCMFD scheme is stable even up to a much larger optical thickness, which implies that the SC prediction is less reliable.

Now that the spectral radius of the DD is slightly larger than that of the SC, the question arises: which transport solver (DD or SC) should be used in Fourier analysis? It is less of a problem for CMFD, pCMFD, and lpCMFD, since both transport solvers agree very well, except for very large optical thicknesses. However, this is not the case for the odCMFD, which shows that the DD results are quite different from those of the SC. The reason is that the original optimal artificial diffusion coefficients were determined based on the SC solution, and they are slightly insufficient for the DD solution. In the next section, it will be shown that the scheme can be stabilized by using the DD-based results.

In addition, it should be noted that when the optical thickness becomes large, the spectral radius gradually reaches the spectral radius of the power (or source) iteration scheme, $\rho =\frac{ta{n}^{-1}\left(\frac{\pi }{{\mathsf{\Sigma }}_{t}L}\right)}{\frac{\pi }{{\mathsf{\Sigma }}_{t}L}}$, which becomes almost 1 for the relatively large domain. Therefore, the convergence or stability of CMFD at a large optical thickness has the least practical importance, in particular for neutron transport problems.

3. Modification of odCMFD

Zhu et al. optimized the artificial diffusion coefficients based on Fourier analysis with the SC for the transport solution [7]. The odCMFD generalizes the CMFD by defining the modified diffusion coefficient as

$$\tilde{D}=\frac{1}{3{\mathsf{\Sigma }}_{t,CM}}+\theta \mathsf{\Delta },$$

(27)

where the artificial diffusion coefficient is defined as

$$\theta =\left\{\begin{array}{ll}0,& \tau <1\\ {{\displaystyle \sum }}_{i=0}^6{a}_i{\tau }^i,& 1\le \tau <14\\ 0.127,& \tau \ge 14\end{array}\right.,$$

(28)

which is a function of the coarse-mesh cell optical thickness $\tau ={\mathsf{\Sigma }}_{t,CM}\mathsf{\Delta }$. The polynomial coefficients $a_i$ are given in

Table 1. Coefficients for the polynomial fitting of optimal $\theta$.

$\mathit{i}$	Original ${\mathit{a}}_{\mathit{i}}$ [7]	Modified ${\mathit{a}}_{\mathit{i}}$
0	−5.543 × 10−2	−1.764 × 10−2
1	8.741 × 10−2	7.156 × 10−2
2	−2.153 × 10−2	−1.409 × 10−2
3	3.146 × 10−3	1.591 × 10−3
4	−2.684 × 10−4	−1.012 × 10−4
5	1.223 × 10−5	3.372 × 10−6
6	−2.285 × 10−7	−4.557 × 10−8

[7]. The Fourier analysis in the previous section shows that odCMFD can lose its effectiveness and even become unstable at a relatively large optical thickness. The following modified artificial diffusion coefficients are obtained by using the DD instead of the SC in the Fourier analysis.

$${\theta }_{mod}=\left\{\begin{array}{ll}0,\hfill & \tau <0.8\\ {{\displaystyle \sum }}_{i=0}^6{a}_i{\tau }^i,\hfill & 0.8\le \tau <14\\ 0.1685,\hfill & 14\le \tau <28\\ 0.1685+0.002067\left(\tau -28\right),\hfill & 28\le \tau <140\\ >0.4,\hfill & \tau >140\end{array}\right.,$$

(29)

where the modified coefficients $a_i$ are given in Table 1. A polynomial function is used for $0.8\le \tau <14$, which is of practical importance for neutron transport problems. A simple linear function is used for $28\le \tau <140$, and a larger value is suggested for $\tau >140$.

Figure 2a shows that the modified optimal $\theta$ values are larger than the original ones, implying that more artificial diffusion should be added to the CMFD scheme. Herein, the odCMFD scheme with the new $\theta$ function is called “modCMFD”. Figure 2b illustrates the spectral radius curves of modCMFD for the SC and DD, respectively. It can be seen that both transport solvers predict almost the same convergence rate with the modified artificial diffusion coefficients, and the scheme is stable for any optical thickness.

4. Optimal SOR for lpCMFD

It has been illustrated that the lpCMFD can be enhanced with SOR. For the fixed-source case, the Fourier analysis indicates that the optimal SOR factor depends on not only the optical thickness, but also the scattering ratio [15], while for the k-eigenvalue case, it only depends on the optical thickness [16]. However, it is found in numerical testing that for better performance, some damping (relaxation) may be necessary for the SOR factors predicted by the Fourier analysis, in particular, for large scattering ratios. In this section, two explicit expressions are developed to analytically determine the local optimal SOR factors for fixed-source problems and k-eigenvalue problems, respectively.

4.1. SOR for Fixed-Source Problems

The spectral radius of the lpCMFD scheme for fixed-source problems is a function of the coarse-mesh cell optical thickness and scattering ratio. Therefore, the optimal SOR factor should also depend on optical thickness and scattering.

Based on a comprehensive Fourier analysis, optimal SOR factors are obtained for the optical thickness in the range of $0.1\le \tau \le 20$, and the scattering ratio in the range of $0.1\le c\le 0.998$. The SOR distribution in terms of the optical thickness and scattering ratio is plotted in Figure 3.

The fitted SOR function is given as

$${\mathrm{SOR}}_{\mathrm{FA}}={{\displaystyle \sum }}_{i=0}^5{{\displaystyle \sum }}_{j=0}^3{a}_{i,j}{\tau }^i{c}^j, i+j\le 5.$$

(30)

Note that $\mathrm{SOR}=1$ when $c<0.1$, $\mathrm{SOR}=1.1$ when $c<0.998$, and $\mathrm{SOR}=\mathrm{SOR} \left(\tau =20\right)$ for $\tau <20$. The polynomial fitting coefficients are given in

Table 2. Polynomial fitting coefficients.

Order $\mathbf{(}\mathit{i},\mathit{j}\mathbf{)}$	Coefficient $\mathbf{\left(}{\mathit{a}}_{\mathit{i},\mathit{j}}\mathbf{\right)}$
(0, 0)	9.898 × 10−1
(1, 0)	1.241 × 10−2
(0, 1)	8.727 × 10−2
(2, 0)	−4.468 × 10−3
(1, 1)	1.579 × 10−1
(0, 2)	−3.770 × 10−2
(3, 0)	6.021 × 10−4
(2, 1)	−2.131 × 10−2
(1, 2)	4.672 × 10−2
(0, 3)	−4.223 × 10−3
(4, 0)	−3.473 × 10−5
(3, 1)	1.316 × 10−3
(2, 2)	−7.051 × 10−3
(1, 3)	6.083 × 10−2
(5, 0)	7.160 × 10−7
(4, 1)	−2.931 × 10−5
(3, 2)	2.328 × 10−4
(2, 3)	−1.871 × 10−3

.

It is seen that SOR increases with the optical thickness and scattering ratio. In practical implementation, some relaxation is needed for relatively large scattering ratios because of numerical effects, such as the treatment of boundary conditions. To apply the relaxation, the domain of the scattering ratio $\left(c\in \left[0.1, 0.998\right]\right)$ is divided into three regions:

• Region 1: $0.1\le c<0.9$;
• Region 2: $0.9\le c<0.99$;
• Region 3: $0.99\le c\le 0.998$.

The numerical results indicate that no relaxation is needed for Region 1. However, for both Regions 2 and 3, the relaxation is required. In order to obtain better fitting, the following scaling is used for Regions 2 and 3.

$$\left\{\begin{array}{ll}{c}_{scaled}=\frac{9\ast \left(c-0.9\right)}{0.098}+1,& 0.9\le c<0.99 \\ c_{scaled}=\frac{9\ast \left(c-0.99\right)}{0.008}+1,& 0.99<c\le 0.998\end{array}\right..$$

(31)

The relaxation functions for Regions 2 and 3 are given as

$$\text{Relax1}={{\displaystyle \sum }}_{i=0}^4{{\displaystyle \sum }}_{j=0}^5{b}_{i,j}{\tau }^i{c}_{scaled}^j, i+j\le 5, (0.9\le c<0.99),$$

(32)

$$\text{Relax2}={{\displaystyle \sum }}_{i=0}^4{{\displaystyle \sum }}_{j=0}^5{c}_{i,j}{\tau }^i{c}_{scaled}^j, i+j\le 5, (0.99\le c<0.998).$$

(33)

The polynomial fitting coefficients $b_{i,j}$ and $c_{i,j}$ are given in

Table 3. Polynomial fitting coefficients for relaxation factors.

Order $\mathbf{(}\mathit{i},\mathit{j}\mathbf{)}$	Coefficient $\mathbf{\left(}{\mathit{b}}_{\mathit{i},\mathit{j}}\mathbf{\right)}$	Coefficient $\mathbf{\left(}{\mathit{c}}_{\mathit{i},\mathit{j}}\mathbf{\right)}$
(0, 0)	1.025 × 10+0	9.851 × 10−1
(1, 0)	−5.911 × 10−3	−1.978 × 10−2
(0, 1)	−1.027 × 10−3	−1.975 × 10−2
(2, 0)	8.845 × 10−4	1.988 × 10−3
(1, 1)	7.404 × 10−4	−5.955 × 10−3
(0, 2)	5.243 × 10−2	1.422 × 10−2
(3, 0)	−7.166 × 10−5	−6.609 × 10−5
(2, 1)	−2.291 × 10−4	3.229 × 10−4
(1, 2)	−2.535 × 10−3	7.285 × 10−4
(0, 3)	−1.216 × 10−2	−3.725 × 10−3
(4, 0)	1.562 × 10−6	8.453 × 10−7
(3, 1)	6.605 × 10−6	−6.618 × 10−7
(2, 2)	2.695 × 10−5	−3.883 × 10−5
(1, 3)	3.953 × 10−4	−4.949 × 10−5
(0, 4)	1.251 × 10−3	3.988 × 10−4
(4, 1)	7.090 × 10−8	−2.618 × 10−7
(3, 2)	−1.480 × 10−6	9.277 × 10−7
(2, 3)	1.859 × 10−6	8.170 × 10−7
(1, 4)	−2.373 × 10−5	5.419 × 10−7
(0, 5)	−4.642 × 10−5	−1.515 × 10−5

.

The final function of SOR is given as

$$\mathrm{SOR}=\left\{\begin{array}{ll}{\mathrm{SOR}}_{\mathrm{FA}},& c<0.9\\ {\mathrm{SOR}}_{\mathrm{FA}}\ast \text{Relax1},& 0.9\le c<0.99\\ {\mathrm{SOR}}_{\mathrm{FA}}\ast \text{Relax2},& c\ge 0.99\end{array}\right..$$

(34)

The relaxed SOR distribution in terms of the optical thickness and scattering ratio is plotted in Figure 4.

4.2. SOR for K-Eigenvalue Problems

As discussed in Section 2.1, differing from fixed-source problems, the spectral radius of the CMFD schemes for k-eigenvalue problems is dependent on the coarse-mesh optical thickness, but independent of scattering.

Based on the Fourier analysis, optimal SOR factors are found for the optical thickness in the range of $0.1\le {\mathsf{\Sigma }}_{t,CM}\mathsf{\Delta }\le 20$. The fitted SOR function is given as

$${\mathrm{SOR}}_{\mathrm{FA}}={{\displaystyle \sum }}_{i=0}^6{a}_i{\tau }^i, 0.1\le \tau \le 20.$$

(35)

When $\tau <0.1$, $\mathrm{SOR}=1$, and when $\tau >20$, $\mathrm{SOR}=\mathrm{SOR} \left(\tau =20\right)$. The polynomial fitting coefficients are given in

Table 4. Polynomial fitting coefficients.

Order (i)	Coefficient (${\mathit{a}}_{\mathit{i}}$)
0	1.019 × 10+0
1	2.474 × 10−1
2	−2.030 × 10−2
3	−8.307 × 10−4
4	2.174 × 10−4
5	−1.163 × 10−5
6	2.049 × 10−7

, and the function above is plotted in Figure 5.

This shows that for k-eigenvalue problems, the SOR factor increases with $\tau$ significantly when $\tau <8$, and after that, the SOR is almost a constant. However, for numerical implementation, some relaxation should be applied to the SOR factor. A relaxation factor of 0.9 is recommended for the entire domain, and additional relaxation is needed for high scattering ratios.

The final expression of SOR for k-eigenvalue problems is summarized as follows:

$$\mathrm{SOR}=0.9\ast \left\{\begin{array}{ll}{\mathrm{SOR}}_{\mathrm{FA}},& c<0.8\\ {\mathrm{SOR}}_{\mathrm{FA}}\ast \left[1+\left(\mathrm{Relax}-1\right)\frac{c-0.8}{0.19}\right],& 0.8<c<0.99\\ {\mathrm{SOR}}_{\mathrm{FA}}\ast \mathrm{Relax},& c\ge 0.99\end{array}\right.,$$

(36)

where

$$\mathrm{Relax}=\left\{\begin{array}{ll}1,& \tau <4.35\\ 1.1\tau -0.0237,& \tau \ge 4.35\end{array}\right..$$

(37)

It should be noted that if the relaxed SOR is smaller than 1, then $\mathrm{SOR}=1$.

5. Numerical Results

The numerical results of the one- and two-dimensional test problems are presented to demonstrate the convergence performance of the modCMFD and the lpCMFD with SOR (lpCMFD-SOR). For comparison, the results of CMFD and pCMFD are presented as well.

5.1. One-Dimensional Iron-Water Test Problem

The 1D modified iron-water test problem is defined in Figure 6 [9,12,15]. The problem domain is uniformly divided into 30 coarse-mesh cells with a mesh size of $∆$ = 1.0 cm, and each coarse-mesh cell consists of 10 uniform fine-mesh cells. The high-order S_N neutron transport equation is discretized by using the Gauss–Legendre S₁₂ quadrature set for the angle, and the SC and DD methods for space, respectively. The S_N equation is solved on the fine-mesh grid. The low-order CMFD equation is discretized and solved on the coarse-mesh grid to accelerate the S_N source iteration. The numerical results are shown in Figure 7.

This shows that the source iteration (SI) scheme converges very slowly because the problem is relatively thick and diffusive ($\tau ={\mathsf{\Sigma }}_{t,CM}\mathsf{\Delta }=3.3333$ and $c= 0.9941$). The CMFD fails to converge. Note that the original odCMFD scheme is slower than the pCMFD for this problem. The modCMFD can improve the convergence rate significantly, and is faster than the pCMFD, as supposed. It can be seen that the original odCMFD is sensitive to the S_N transport solver. It converges more slowly for the DD than SC, which is consistent with the theoretical finding. The modCMFD is not sensitive to the transport solver and it yields almost the same performance for both the DD and SC. As expected, the lpCMFD scheme performs best among all the CMFD schemes. More interestingly, the lpCMFD-SOR scheme can further accelerate the convergence, in which the local optimal SOR factor for each coarse-mesh cell is precomputed by Equation (34). The results are similar to the case reported by Wang [15], where a global SOR is used.

5.2. Two-Dimensional Test Problem

The second test case is a 2D k-eigenvalue problem as shown in Figure 8. The problem domain is a $6\times 6$ cm square with the reflective boundary condition on each side. The domain is divided into $5\times 5$ uniform coarse-mesh cells, with each consisting of $10\times 10$ fine-mesh cells. The numerical solutions for the S_N transport are obtained on the fine-mesh ($50\times 50$) using the level-symmetric S₁₂ quadrature set for angular discretization and the DD method for spatial discretization. The CMFD solutions are obtained on the coarse-mesh grids ($5\times 5$). Two cases are presented for the problem: Case A has the scattering ratio of 0.9, and it is 0.7 in Case B. The problem is to investigate whether the scattering has any impact on the convergence, although the Fourier analysis shows that the spectral radius of the CMFD schemes should be independent of scattering.

The convergence performance is illustrated in Figure 9. It is interesting to notice that the scattering does have some effect on the convergence of the pCMFD and odCMFD. In general, the convergence slows down as the scattering ratio decreases. However, both the modCMFD and lpCMFD are much less sensitive to the scattering. Note that the local SOR factors of the coarse-mesh cells are calculated according to Equation (36).

5.3. Two-Dimensional Mini-Assembly Problem

In this section, the convergence comparison of the CMFD schemes for a 2D mini-assembly k-eigenvalue problem is shown in Figure 10. Similar to the previous problem, the problem domain is a $6\times 6$ cm square with the reflective boundary condition on each side. The domain is divided into $5\times 5$ uniform coarse-mesh cells, and each coarse cell consists of $10\times 10$ fine-mesh cells. In each coarse cell, the central $4\times 4$ fine cells represent the fuel region, and all the other cells are the water region. The numerical solutions for the S_N transport are obtained on the fine mesh ($50\times 50$) using the level-symmetric S₁₂ quadrature set for angular discretization and the DD method for spatial discretization. The CMFD solutions are obtained on the coarse-mesh grids ($5\times 5$). The homogenized cross-sections on the coarse-mesh cells are updated during each iteration by weighting with the volume and flux.

The normalized scalar flux is plotted in Figure 11a. The convergence comparison is illustrated in Figure 11b. In this problem, CMFD also fails to converge. The convergence of odCMFD is similar to that of modCMFD. Again, the lpCMFD scheme outperforms other schemes, and SOR can improve the convergence rate significantly.

6. Conclusions

Fourier analysis results have been presented for the acceleration schemes: CMFD, pCMFD, odCMFD, and lpCMFD. It has been shown that the DD transport solver predicts slightly larger spectral radii than the SC, which can have some implications on the development of diffusion-based stabilizing techniques, as exemplified by the odCMFD scheme. The odCMFD becomes less effective and unstable for a relatively large optical thickness in the k-eigenvalue problems. This is because the original optimal artificial diffusion coefficients determined by the SC solution are less sufficient for stability. To stabilize the odCMFD, these coefficients have been redeveloped based on the DD-based Fourier analysis. The numerical results for different test problems have confirmed that the modified odCMFD (modCMFD) scheme is more stable and effective.

The lpCMFD convergence can be enhanced with SOR. However, finding an optimal global SOR is not straightforward and is problem-dependent. In this paper, two explicit expressions have been proposed for the calculation of local optimal SOR factors—one is for fixed source problems, and the other is for k-eigenvalue problems. These local SOR factors can be precomputed or updated on-the-fly during the iteration.