CN103150424B - 一种获取反应堆堆芯三维中子通量密度精细分布的方法 - Google Patents

Abstract

一种获取反应堆堆芯三维中子通量密度精细分布的方法，1、对堆芯几何建模，划分计算区域，离散角度空间，生成特征线信息，指定各计算区域的材料，获取其宏观截面参数，为计算区域的中子通量密度通量、反应堆边界条件、特征值设置初始值；2、计算每一个子区的矩阵方法需要的系数矩阵和右端项不随着迭代计算变化的部分；3、求粗网中子通量密度，修正一维、二维的细网中子通量密度；4、求每个栅元的一维中子通量密度；5、迭代求解每一层的二维中子通量密度；6、更新三维粗网参数，判断特征值和三维中子通量密度是否收敛，如果不收敛，则转到步骤3继续迭代，直至收敛，即得到三维中子通量密度精细分布；能够快速得到反应堆堆芯三维中子通量密度的精细分布。

Description

一种获取反应堆堆芯三维中子通量密度精细分布的方法

技术领域

本发明涉及核反应堆堆芯设计和反应堆物理计算领域，具体涉及一种获取反应堆堆芯三维中子通量密度精细分布的方法。

背景技术

为了保证反应堆堆芯设计安全和运行安全，需要准确快速地计算出反应堆及相关的设备内三维中子通量密度分布的情况。

目前广泛采用的传统反应堆物理分析计算方法“两步法”是建立在等效均匀化理论基础之上的。随着反应堆堆芯设计越来越复杂，安全要求越来越高，“两步法”逐渐不能满足工程计算要求。

所谓“两步法”，第一步是在全反射边界条件下对各种非均匀组件进行多群中子输运计算，得到组件内的中子通量密度分布从而归并出组件的少群等效均匀化群常数及不连续因子等物理量，以备堆芯计算使用；第二步，是根据前一步生成的均匀化参数，采用粗网节块方法对堆芯进行少群中子扩散或输运计算，通过功率重构计算，便可获得三维堆芯中子通量密度分布。

然而，“两步法”由于只针对有限种孤立的组件来计算等效均匀化参数，导致它无法考虑在真实堆芯布置中组件位置不同受到的影响，因此在反应堆精细中子通量密度分布计算要求不断提高的今天，该方法越来越受到局限。尽管诸多改进方法被相继提出，但依然不能跳出等效均匀化的思想框架。

能够从根本上解决上述问题的方法只有进行“一步”全堆芯非均匀中子输运计算。然而就目前的计算条件而言，直接对大型核反应堆堆芯进行三维中子输运方程求解还不现实。所以需要发明一种高效、准确的获取反应堆堆芯三维中子通量密度的方法。

发明内容

为解决上述现有技术存在的问题，本发明提供一种获取反应堆堆芯三维中子通量密度精细分布的方法，能够快速得到反应堆堆芯三维中子通量密度的精细分布，给反应堆堆芯设计和安全提供可靠的信息。

为了达到上述目的，本发明采用如下技术方案：

一种获取反应堆堆芯三维中子通量密度精细分布的方法，包括如下步骤：

步骤1：对反应堆堆芯进行几何建模，将反应堆堆芯按照轴向划分成多层，用Auto-CAD分别描述反应堆堆芯每一层几何尺寸，将每一层沿径向划分成多个子区，再用Auto-CAD对每个子区划分计算区域，根据计算需求对每一层计算区域的三维角度空间进行离散，利用VBA二次开发的Auto-CAD的工程，在每一个子区内生成特征线信息，子区内边界通过插值进行耦合；根据反应堆堆芯材料信息，读取各层中每一个计算区域的宏观截面参数；对二维、一维计算区域的中子通量密度、反应堆边界条件和特征值赋初值；

步骤2：在步骤1所得几何尺寸、材料和特征线信息的基础上，计算每一层二维计算区域中每一个子区的矩阵MOC方法需要的系数矩阵和右端项不随着迭代计算变化的部分，系数矩阵不随迭代计算变化；

步骤3：用步骤1中子通量密度初值计算三维粗网有限差分CMFD需要的净流修正系数，进行一次全堆芯三维粗网有限差分CMFD计算，求得粗网中子通量密度后，修正一维、二维的细网中子通量密度，所述二维的细网中子通量密度修正如公式（25）所示：

${\mathrm{\φ}}_g^{i,(l+1/2)}=\frac{{\stackrel{\‾}{\mathrm{\φ}}}_g^{n,(l+1/2)}}{{\stackrel{\‾}{\mathrm{\φ}}}_g^{n,\left(l\right)}}{\mathrm{\φ}}_g^{i,\left(l\right)},$ i∈n公式(25)

式中：

——修正后的平源细区i第g群的平均标通量，用以更新第l+1次矩阵MOC计算的源项，中子通量密度的单位是cm^-2·s^-1，

——第l次矩阵MOC计算后的三维CMFD计算得到的粗网n的第g群平均标通量，

——由第l次矩阵MOC计算得到的粗网n的第g群均匀化标通量，

——第l次矩阵MOC计算所得到的平源细区i第g群的平均标通量；

对于一维的细网中子通量密度修正如公式（26）所示：

$a_{\mathrm{gi}}^{k,(l+1/2)}=\frac{{\stackrel{\‾}{\mathrm{\φ}}}_g^{n,(l+1/2)}}{a_{g0}^{k,\left(l\right)}}{\mathrm{a\φ}}_{\mathrm{gi}}^{k,\left(l\right)},$ k∈n,i＝0,1,2,3,4    公式(26)

式中：

——修正后的一维节块k第g群的第i阶通量展开系数，用以更新第l+1次NEM计算的源项；

——第l次一维NEM计算后的三维CMFD计算得到的粗网n的第g群平均标通量；

——第l次一维NEM计算得到节块k内第g群的第i阶中子通量密度展开系数；

步骤4：更新一维节块通量展开需要的展开系数和径向泄漏项，对每一个栅元逐层、每一能群地进行一维节块展开NEM固定源计算，将每个节块的中子通量密度和中子源项用多项式展开，通过求解多项式的展开系数，求得每个栅元的一群中子通量密度，具体的迭代过程如下所示：

1）更新径向泄漏项展开系数

由以下公式求得：

$l_{g0}^k={\mathrm{TL}}_{g,k}^{\mathrm{Radial}}$ 公式(31)

$l_{g1}^k=\frac{{\mathrm{\Δz}}_k}{d}[z_{k-1}({\mathrm{TL}}_{g,k+1}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k}^{\mathrm{Radial}})+z_{k+1}({\mathrm{TL}}_{g,k}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k-1}^{\mathrm{Radial}})]$ 公式(32)

$l_{g2}^k=\frac{{\left({\mathrm{\Δz}}_k\right)}^2}{d}[({\mathrm{\Δz}}_{k-1}+{\mathrm{\Δz}}_k)({\mathrm{TL}}_{g,k+1}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k}^{\mathrm{Radial}})-({\mathrm{\Δz}}_{k+1}+{\mathrm{\Δz}}_k)({\mathrm{TL}}_{g,k}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k-1}^{\mathrm{Radial}})]$

公式(33)

其中：

z_i=(Δz_k+Δz_i)(Δz_k+2Δz_i)d＝(Δz_k-1+Δz_k)(Δz_k+Δz_k+1)(Δz_k-1+Δz_k+Δz_k+1)

${\mathrm{TL}}_{g,k}^{\mathrm{Radial}}=\frac{{\mathrm{\ξ}}_m}{\mathrm{\Δx}}({\mathrm{\Ψ}}_{g,x+}-{\mathrm{\Ψ}}_{g,x-})+\frac{{\mathrm{\η}}_m}{\mathrm{\Δy}}\left({\mathrm{\Ψ}}_{g,y+}{-\mathrm{\Ψ}}_{g,y-}\right)$

Δzk、Δx、Δy分别是x、y、z方向网格间距，

${\mathrm{\Ψ}}_{g,x\±}^p\left(z\right)=\frac{1}{\mathrm{\Δy}}{\∫}_{y_{j-1/2}}^{y_{j+1/2}}{\mathrm{\Ψ}}_g(x_{i\±1/2},y,z)\mathrm{dy}$

${\mathrm{\Ψ}}_{g,y\±}^p\left(z\right)=\frac{1}{\mathrm{\Δy}}{\∫}_{x_{j-1/2}}^{x_{j+1/2}}{\mathrm{\Ψ}}_g({x,y}_{j\±1/2},z)\mathrm{dx}$

ξ_m、η_m分别是x、y方向的方向余弦

2）假定初值：（g＝1,K,G;k=1,K,K;n=0,K,4）；

其中：

表示节块k第g能群轴向的上表面和下表面净流值，净流的单位

是cm^-2·s^-1，

表示节块k第g能群中子通量密度展开多项式第一项系数，

中子通量密度多项式展开如下所示：

${\mathrm{\φ}}_g^k\left(z\right)=\underset{n=0}{\overset{4}{\mathrm{\Σ}}}{a}_{\mathrm{gn}}^k{H}_n\left(z\right)$ 公式(34)

其中：

$H_n\left(z\right)=\{1,\mathrm{\γ},{3\mathrm{\γ}}^2-\frac{1}{4},{\mathrm{\γ}}^4-\frac{{3\mathrm{\γ}}^2}{10}+\frac{1}{80}\}$

z∈[z_k-1/2,z_k+1/2

k表示节块编号，g表示能群编号，n表示展开阶数；

3）构造源项展开系数

一维节块展开计算NEM中源项多项式展开如下所示：

$Q_g^k\left(z\right)=\underset{n=0}{\overset{4}{\mathrm{\Σ}}}{q}_{\mathrm{gn}}^k{H}_n\left(z\right)$ 公式(35)

其中：

$q_{\mathrm{gn}}^k=\frac{x_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left({\mathrm{v\Σ}}_f\right)}_{g^{\′},k}{a}_{g^{\′}n}^k+\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s,g^{\′}-g,k}{a}_{g^{\′}n}^k,$ n＝0,1,2,3,4

是节块k第g能群中子通量密度多项式展开的系数，求解公式如下所示：

$a_{g0}^k={\stackrel{\‾}{\mathrm{\φ}}}_g^k$ 公式(36)

$a_{g1}^k={\mathrm{\φ}}_{g+}^k-{\mathrm{\φ}}_{g-}^k$ 公式(37)

$a_{g2}^k={\mathrm{\φ}}_{g+}^k+{\mathrm{\φ}}_{g-}^k-2{\stackrel{\‾}{\mathrm{\φ}}}_g^k$ 公式(38)

$a_{g3}^k=b_{g3}^k{a}_{g1}^k+c_{g3}^k$ 公式(39)

公式(39)中：

$b_{g3}^k=\frac{{\mathrm{\Σ}}_{r,g,k}}{\frac{{6D}_g^k}{{\left({\mathrm{\Δz}}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{10}}$

$c_{g3}^k=\frac{\frac{q_{g3}^k}{10}-q_{g1}^k}{\frac{{6D}_g^k}{{\left(\mathrm{\Δ}{z}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{10}}$

$a_{g4}^k=b_{g4}^k{a}_{g2}^k+c_{g4}^k$ 公式(40)公式(40)中：

$b_{g4}^k=\frac{{\mathrm{\Σ}}_{r,g,k}}{\frac{{4D}_g^k}{{\left({\mathrm{\Δz}}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{35}}$

$c_{g4}^k=\frac{\frac{q_{g4}^k}{35}-q_{g2}^k}{\frac{{4D}_g^k}{{\left({\mathrm{\Δz}}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{35}}$

其中：

是节块k第g群内的中子通量密度，

是节块k第g群上下界面的中子通量密度，

∑_r，g,k是节块k第g群移除截面，单位是cm^-1

节块k第g群扩散系数，单位是cm

4）求解由两个边界条件方程以及净中子流耦合方程构成的线性代数方程组，得到

三个方程如下所示：

$A_{g2}^k{J}_{\mathrm{zg}-}^k-(A_{g1}^k+A_{g1}^{k+1})J_{\mathrm{zg}+}^k+A_{g2}^{k+1}{J}_{\mathrm{zg}+}^{k+1}=A_{g0}^k-A_{g3}^{k+1},$ k＝1,K-1公式(41) $A_{g2}^K{J}_{\mathrm{zg}-}^K+[4(\frac{{\mathrm{\β}}^{+}}{{1-\mathrm{\β}}^{+}}+\frac{1}{2})-A_{g1}^K]J_{\mathrm{zg}+}^K=A_{g0}^K$ 公式(42) $[4(\frac{{\mathrm{\β}}^{-}}{{\mathrm{\β}}^{-}-1}-\frac{1}{2})+A_{g1}^1]J_{\mathrm{zg}-}^1-A_{g2}^1{J}_{\mathrm{zg}+}^1=A_{g3}^1$ 公式(43)其中：

表示节块k第g能群轴向的上表面和下表面净中子流值

$\left\{\begin{array}{c}{A}_{g0}^k=\frac{{\hat{a}}_{g2}^k{\hat{r}}_{g2}^k-{\hat{a}}_{g1}^k{\hat{r}}_{g1}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\\ A_{g1}^k\frac{\frac{{\mathrm{\Δz}}_k}{D_g^k}{\hat{a}}_{g1}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\\ A_{g2}^k=\frac{\frac{{\mathrm{\Δz}}_k}{D_g^k}{\hat{a}}_{g2}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\\ A_{g3}^k=\frac{{\hat{a}}_{g1}^k{\hat{r}}_{g2}^k-{\stackrel{k}{a}}_{g2}^k{\hat{r}}_{g1}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\end{array}\right.$

$\left\{\begin{array}{c}{\hat{a}}_{g1}^k=\frac{1}{2}{b}_{g3}^k+\frac{1}{5}{b}_{g4}^k+4\\ {\hat{a}}_{g2}^k=\frac{1}{2}{b}_{g3}^k-\frac{1}{5}{b}_{g4}^k-2\\ {\hat{r}}_{g1}^k=(6+\frac{2}{5}{b}_{g4}^k){\stackrel{\‾}{\mathrm{\φ}}}_g^k-\frac{1}{2}{c}_{g3}^k-\frac{1}{5}{c}_{g4}^k\\ {\hat{r}}_{g2}^k=-{\hat{r}}_{g1}^k-c_{g3}^k\end{array}\right.$

其中

Δz_k由上一步计算求得；

5）求解下述中子平衡方程，得到节块平均通量

$\frac{1}{{\mathrm{\Δz}}_k}(J_{\mathrm{zg}+}^k-J_{\mathrm{zg}-}^k)+{\mathrm{\Σ}}_{r,g,k}{\stackrel{\‾}{\mathrm{\φ}}}_g^k=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left({\mathrm{v\Σ}}_f\right)}_{g^{\′},k}{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^k+\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s,g^{\′}-g,k}{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^k+{\mathrm{TL}}_{g,k}^{\mathrm{Radial}}\left(z\right)$ 公式(44)

6）由节块表面中子通量与净中子流关系式，求得节块界面中子通量

${\mathrm{\φ}}_{g\±}^{k,\left(l\right);}$

${\mathrm{\φ}}_{g+}^k=A_{g1}^k{J}_{\mathrm{zg}+}^k-A_{g2}^k{J}_{\mathrm{zg}-}^k+A_{g0}^k$ 公式(45)

${\mathrm{\φ}}_{g-}^k=A_{g2}^k{J}_{\mathrm{zg}+}^k-A_{g1}^k{J}_{\mathrm{zg}-}^k+A_{g3}^k$ 公式(46)

7）更新通量高阶展开系数

如步骤3）中所示；

8）对每一个栅元的每一个能群进行3）到6）迭代计算，若通量收敛则结束计算，收敛准则如下：

$\left|\right|\frac{{\stackrel{\‾}{\mathrm{\φ}}}_g^k-{\stackrel{\‾}{\mathrm{\φ}}}_g^{k-1}}{{\stackrel{\‾}{\mathrm{\φ}}}_g^k}\left|\right|\≤{\mathrm{\ϵ}}_1$ 公式(47)

ε₁为很小的正数，也就是收敛条件；

步骤5：更新轴向泄漏项，根据步骤1划分的几何区域，对轴向的每一层进行二维输运计算，根据步骤2所得到的每个子区的系数矩阵，采用GMRES算法，并行求解每个子区每一群的中子通量密度，在边界处通过插值互相耦合，迭代求解每一层的二维中子通量密度，具体步骤如下：

1）更新轴向泄漏项

该泄漏项由一维计算提供；

2）读取步骤2产生的子区系数矩阵和右端项不随迭代改变部分，设定子区内、外边界入射角通量和子区内中子通量密度初始值；

3）生成子区群外源，即右端项随着迭代变化的部分

$Q_{\mathrm{\ρi}}=\left\{\begin{array}{cc}{Q}_i^f+Q_i^{\mathrm{iso},g^{\′}-g},& \mathrm{\ρ}=0\\ \underset{\stackrel{g^{\′}=1}{g^{\′}\≠g}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{\mathrm{sl},g^{\′}-g,i}{\mathrm{\φ}}_{\mathrm{\ρ},i,g^{\′}},& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.;$

4）采用GMRES方法求解子区群内方程；

5）与相邻子区进行通信，交换内边界信息，即交换内边界出射角通量信息；

6）判断子区通量是否收敛，若不收敛，继续迭代求解；

7）判断所有能群是否都已经计算，对各个能群进行迭代，直到所有区域通量收敛；

步骤6：用步骤4得到的一维通量密度和步骤5得到的二维通量密度更新三维粗网参数，判断特征值和三维中子通量密度是否收敛，如果不收敛，则转到步骤3继续迭代，并且计算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收敛，就可以得到三维中子通量密度精细分布；

所述收敛条件为：

$\left|\right|\frac{{\mathrm{\φ}}_g^{i,\left(l\right)}{-\mathrm{\φ}}_g^{i,(l-1)}}{{\mathrm{\φ}}_g^{i,\left(l\right)}}\left|\right|{\≤\mathrm{\ϵ}}_1,$ $\left|\right|\frac{k_{\mathrm{eff}}^{\left(l\right)}{-k}_{\mathrm{eff}}^{(l-1)}}{k_{\mathrm{eff}}^{\left(l\right)}}\left|\right|{\≤\mathrm{\ϵ}}_2$ 公式(48)

其中：

为第l次迭代第i区第g群通量

为第l次迭代特征值

ε₁，ε₂为两个很小的正数，也就是收敛条件。

步骤2所述计算每一层二维计算区域中每一个子区的矩阵MOC方法需要的系数矩阵和右端项不随着迭代计算变化的部分的计算方法为：

任意子区域sub上的矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程如下式所示：

公式(16)

ρ＝0,1,2,3

公式(17)

公式(18)其中：

φ_ρ，i是中子通量密度球谐函数矩，

是子区外边界入射中子角通量密度，

是子区内边界入射中子角通量密度，

$a_{0,i}=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{s0,g-g,i}}{{\mathrm{\Σ}}_{t,i}},& \mathrm{\ρ}=0\\ 0,& \mathrm{\ρ}=1,\mathrm{2,3}\end{array}\right.$

$a_{{\mathrm{\ρ}}^{\′},i}=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{s0,g-g,i}}{{\mathrm{\Σ}}_{t,i}},& {\mathrm{\ρ}}^{\′}=\mathrm{\ρ}\\ 0,& {\mathrm{\ρ}}^{\′}\≠1,\mathrm{2,3}\end{array}\right.$

${\mathrm{\τ}}_{\mathrm{\ρ},{\mathrm{\ρ}}^{\′},j,i}^l={\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l{\mathrm{\σ}}_{{\mathrm{\ρ}}^{\′},j}^l$

${\mathrm{\ϵ}}_{\mathrm{\ρ},i}^l={\mathrm{\λ}}_{\mathrm{\ρ},i}^l{\mathrm{\δe}}_{1,i}^l$

$q_{\mathrm{\ρ},i}^{\mathrm{loc}}=\left\{\begin{array}{cc}\frac{Q_i^f+Q_i^{\mathrm{iso},g^{\′}-g}}{{\mathrm{\Σ}}_{t,i}},& \mathrm{\ρ}=0\\ \frac{\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s1,g^{\′}-g,i}{\mathrm{\φ}}_{\mathrm{\ρ},i,g^{\′}}}{{\mathrm{\Σ}}_{t,i}},& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.$

$q_{\mathrm{\ρ},i}^{\mathrm{trans}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l{Q}_j^{\mathrm{outgrp}}\left({\mathrm{\Ω}}_m\right)$

δ_k0克罗内特尔向量

$e_{j,p}^l=\left\{\begin{array}{cc}\exp (-\underset{n=j}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{t,n}{\hat{S}}_{\mathrm{nl}}),& j\≤p\\ 1,& j>p\end{array}\right.$

${\mathrm{\σ}}_{{\mathrm{\ρ}}^{\′},j}^l=\left\{\begin{array}{cc}{\mathrm{\Σ}}_{s0,g-g,j},& {\mathrm{\ρ}}^{\′}=0\\ {3A}_{{\mathrm{\ρ}}^{\′}}\left({\mathrm{\Ω}}_m\right){\mathrm{\Σ}}_{s1,g-g,j},& {\mathrm{\ρ}}^{\′}=\mathrm{1,2,3}\end{array}\right.$

$c_{j,i}^l=\frac{{\mathrm{\δe}}_{j+1,i}^l{-\mathrm{\δe}}_{j,i}^l}{{\mathrm{\Σ}}_{t,j}}$

${\mathrm{\δe}}_{j,p}^l=e_{j,p-1}^l-e_{j,p}^l$

任意子区域sub上的线性代数方程组便可写为如下形式：

$A_{\mathrm{sub}}\left[\begin{array}{c}{\stackrel{v}{\mathrm{\φ}}}_{\mathrm{sub}}\\ {\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{ob},-}\\ {\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},+}\end{array}\right]={\stackrel{v}{b}}_{\mathrm{sub}}\left({\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},-}\right)$ 公式(19)

式中：

A_sub——子区域sub的方程组系数矩阵；

——子区域sub的角通量密度矩向量；

——子区域sub的外边界入射角通量向量（若子区域sub不存在外边界，则不含此向量）；

——子区域sub的内边界出射角通量向量；

——子区域sub的方程组右端项，它与内边界入射角通量向量ψ

相关；

将右端项拆解成以下形式：

${\stackrel{v}{b}}_{(N{\mathrm{\φ}}_{\mathrm{sub}}+N{\mathrm{\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×1}=R_{(N_{{\mathrm{\φ}}_{\mathrm{sub}}+}{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{\φ}}_{\mathrm{sub}}}{Q}_{N{\mathrm{\φ}}_{\mathrm{sub}}\×1}+E_{({\mathrm{N\φ}}_{\mathrm{sub}}+N_{}{\mathrm{\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+N{\mathrm{\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},n})\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}{\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},-}$

公式(20)

或者，

$\left[\begin{array}{c}{B}_{{\mathrm{N\φ}}_{\mathrm{dub}}\×1}^1\\ B_{{\mathrm{N\Ψ}}_{\mathrm{sub}}\×1}^2\\ B_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}{\×1}_{-}}^3\end{array}\right]=\left[\begin{array}{c}{R}_{{\mathrm{N\φ}}_{\mathrm{sub}}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^1\\ R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^2\\ R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^3\end{array}\right]Q_{N{\mathrm{\φ}}_{\mathrm{sub}}\×1}+\left[\begin{array}{c}{E}_{{\mathrm{N\φ}}_{\mathrm{sub}}\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}^1\\ E_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}^2\\ E_{N_{\mathrm{sub}}^{\mathrm{inb},+}\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}^3\end{array}\right]{\stackrel{v}{\mathrm{\ψ}}}_{\mathrm{sub}}^{\mathrm{inb},-}$ 公式(21)

式中：

Nφ_sub——子区域角通量密度矩向量

的长度；

——子区域外边界入射角通量向量ψ

的长度；

——子区域内边界出射角通量向量ψ

的长度；

$Q_{N{\mathrm{\φ}}_{\mathrm{sub}}\×1}={\left[Q_{\mathrm{\ρi}}\right]}_{\mathrm{\ρ}\∈\left[\mathrm{0,3}\right]i\∈[1,\mathrm{Nr}]}$

$Q_{\mathrm{\ρi}}=\left\{\begin{array}{cc}{Q}_i^f+Q_i^{\mathrm{iso},g^{\′}-g},& \mathrm{\ρ}=0\\ \underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s1,g^{\′}-g,i}{\mathrm{\φ}}_{\mathrm{\ρ},i,g^{\′}},& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.$ 右端项随迭代变化部分，

步骤2中求解的就是每一层中每个子区的系数矩阵A_sub和不随迭代变化的模板矩阵 $R_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}){\×\mathrm{\φ}}_{\mathrm{sub}}},$ $E_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}};$

A_sub具体求解由前述的子区矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程得到；

$R_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}){\×\mathrm{\φ}}_{\mathrm{sub}}}$ 求解如下：

$R_{{\mathrm{N\φ}}_{\mathrm{sub}}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^1=\left[\begin{array}{ccccccccccc}L& M& M& M& M& L& M& M& M& M& L\\ L& {\hat{r}}_i^1+r_{i,i,\mathrm{0,0}}^1& r_{i,i,\mathrm{0,1}}^1& r_{i,i,\mathrm{0,2}}^1& r_{i,i,\mathrm{0,3}}^1& L& r_{i,j,\mathrm{0,0}}^1& r_{i,j,0,1}^1& r_{i,j,\mathrm{0,2}}^1& r_{i,j,\mathrm{0,3}}^1& L\\ L& r_{i,i,\mathrm{1,0}}^1& {\hat{r}}_i^1+r_{i,i,\mathrm{1,1}}^1& r_{i,i,\mathrm{1,2}}^1& r_{i,i,\mathrm{1,3}}^1& L& r_{i,j,\mathrm{1,0}}^1& r_{i,j,\mathrm{1,1}}^1& r_{i,j,\mathrm{1,2}}^1& r_{i,j,\mathrm{1,3}}^1& L\\ L& r_{i,i,\mathrm{2,0}}^1& r_{i,i,\mathrm{2,1}}^1& {\hat{r}}_i^1+r_{i,i,\mathrm{2,2}}^1& r_{i,i,\mathrm{2,3}}^1& L& r_{i,j,\mathrm{2,0}}^1& r_{i,j,\mathrm{2,1}}^1& r_{i,j,\mathrm{2,2}}^1& r_{i,j,\mathrm{2,3}}^1& L\\ L& r_{i,i,\mathrm{3,0}}^1& r_{i,i,\mathrm{3,1}}^1& r_{i,i,\mathrm{3,2}}^1& {\hat{r}}_i^1+r_{i,i,\mathrm{3,3}}^1& L& r_{i,j,\mathrm{3,0}}^1& r_{i,j,\mathrm{3,1}}^1& r_{i,j,\mathrm{3,2}}^1& r_{i,j,\mathrm{3,3}}^1& L\\ L& M& M& M& M& L& M& M& M& M& L\\ L& r_{j,i,\mathrm{0,0}}^1& r_{j,i,\mathrm{0,1}}^1& r_{j,i,\mathrm{0,2}}^1& r_{j,i,\mathrm{0,3}}^1& L& {\hat{r}}_j^1+r_{j,j,\mathrm{0,0}}^1& r_{j,j,\mathrm{0,1}}^1& r_{j,j,\mathrm{0,2}}^1& r_{j,j,\mathrm{0,3}}^1& L\\ L& r_{j,i,\mathrm{1,0}}^1& r_{j,i,\mathrm{1,1}}^1& r_{j,i,\mathrm{1,2}}^1& r_{j,i,\mathrm{1,3}}^1& L& r_{j,j,\mathrm{1,0}}^1& {\hat{r}}_j^1+r_{j,j,\mathrm{1,1}}^1& r_{j,j,\mathrm{1,2}}^1& r_{j,j,\mathrm{1,3}}^1& L\\ L& r_{j,i,\mathrm{2,0}}^1& r_{j,i,\mathrm{2,1}}^1& r_{j,i,\mathrm{2,2}}^1& r_{j,i,\mathrm{2,3}}^1& L& r_{j,j,\mathrm{2,0}}^1& r_{j,j,2,1}^1& {\hat{r}}_j^1+r_{j,j,\mathrm{2,2}}^1& r_{j,j,\mathrm{2,3}}^1& L\\ L& r_{j,i,\mathrm{3,0}}^1& r_{j,i,\mathrm{3,1}}^1& r_{j,i,\mathrm{3,2}}^1& r_{j,i,\mathrm{3,3}}^1& L& r_{j,j,\mathrm{3,0}}^1& r_{j,j,\mathrm{3,1}}^1& r_{j,j,\mathrm{3,2}}^1& {\hat{r}}_j^1+r_{j,j,\mathrm{3,3}}^1& L\\ L& M& M& M& M& L& M& M& M& M& L\end{array}\right]$

公式(22)

其中的元素表达形式为：

${\hat{r}}_i^1=\frac{1}{{\mathrm{\Σ}}_{t,i}}$

$r_{i,j,\mathrm{\ρ},{\mathrm{\ρ}}^{\′}}^1=\left\{\begin{array}{cc}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}\left({\mathrm{\λ}}_{\mathrm{\ρ},i}^l\right),& {\mathrm{\ρ}}^{\′}=0\\ \underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{3A}_{{\mathrm{\ρ}}^{\′}}\left({\mathrm{\Ω}}_m\right){\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l,& {\mathrm{\ρ}}^{\′}=\mathrm{1,2,3}\end{array}\right.$

$R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^{2L}=\left[\begin{array}{cccccc}L& M& M& M& M& L\\ L& r_{l_r,j,0}^2& r_{l_r,j,1}^2& r_{l_r,j,2}^2& r_{l_r,j,3}^2& L\\ L& M& M& M& M& L\end{array}\right]$ 公式(23)

其中的元素表达形式为：

$R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^{2L}=\left[\begin{array}{cccccc}L& M& M& M& M& L\\ L& r_{l,j,0}^2& r_{l,j,1}^2& r_{l,j,2}^2& r_{l,j,3}^2& L\\ L& M& M& M& M& L\end{array}\right]$ 公式(24)

其中的元素表达形式为：

则由矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程右端第二项的系数求的。

步骤3所述计算三维粗网有限差分CMFD需要的净流修正系数的公式为：

公式(27)

公式(28)

${\stackrel{\‾}{D}}_g^n=\frac{1}{3{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{tr},g}^n}$ 公式(29)

式中：

——高阶计算得到的粗网n的界面净流，

——第g群的宏观移除截面，

——网格间距。

步骤3所述的全堆芯三维粗网有限差分CMFD计算公式为：

$d_{\mathrm{zg}-}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n-1_z}+d_{\mathrm{yg}-}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n{-1}_y}+d_{\mathrm{xg}-}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n-1_x}+({\stackrel{\‾}{\mathrm{\Σ}}}_{r,g}^n+d_{\mathrm{xg}}^n+d_{\mathrm{zg}}^n){\stackrel{\‾}{\mathrm{\φ}}}_g^n$

公式(30)

$-d_{\mathrm{xg}+}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n{+1}_x}-d_{\mathrm{yg}+}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n+1_y}{-d}_{\mathrm{zg}+}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n{+1}_z}=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left(v{\stackrel{\‾}{\mathrm{\Σ}}}_f\right)}_{g^{\′}}^n{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^n+\underset{\underset{g^{\′}\≠g}{g^{\′}}=1}{\overset{G}{\mathrm{\Σ}}}{\stackrel{\‾}{\mathrm{\Σ}}}_{s,g^{\′}g}^n{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^n$

式中：

是中子通量密度，

是扩散系数，

是中子裂变产额截面，

是中子散射截面，

x_g是第g群裂变谱，

k_eff是特征值，也称为有效增殖因数，

其中：

${\stackrel{\‾}{D}}_g^n=\frac{1}{3{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{tr},g}^n}$

本发明和现有技术相比，具有如下优点：

1.一步法输运计算考虑了组件在堆芯不同位置所受到的影响，这在传统两步法中不能精细考虑；

2.避免传统MOC方法中特征线反复迭代计算，使用预先求解的系数矩阵，提高计算效率；

3.引入区域分解并行算法后可以直接获得大型反应堆的三维中子通量密度分布，并且具有较高的计算效率；

4.改进一维NEM方法，提高计算效率；

5.几何处理采用AutoCAD的二次开发技术，高效快速地处理二维非均匀几何，生成特征线信息。

附图说明

图1为二维/一维耦合计算流程图。

图2为反应堆堆芯几何描述。

图3为计算区域角度空间离散说明。

图4为计算区域特定方向特征线示意图。

图5为一维NEM方法计算流程示意图。

图6为区域分解矩阵MOC方法计算流程图。

图7为三维非均匀压水堆堆芯相对棒功率分布。

具体实施方式

下面结合附图和具体实施方式对本发明结构进行详细说明。

如图1所示，本发明一种获取反应堆堆芯三维中子通量密度精细分布的方法，包括如下步骤：

步骤1：对反应堆堆芯进行几何建模，将反应堆堆芯按照轴向划分成多层，用Auto-CAD分别描述反应堆堆芯每一层几何尺寸，将每一层沿径向划分成多个子区，再用Auto-CAD对每个子区划分计算区域，根据计算需求对每一层计算区域的三维角度空间进行离散，利用VBA二次开发的Auto-CAD的工程，在每一个子区内生成特征线信息，子区内边界通过插值进行耦合；根据反应堆堆芯材料信息，读取各层中每一个计算区域的宏观截面参数；对二维、一维计算区域的中子通量密度、反应堆边界条件和特征值赋初值；

步骤1中的几何描述、计算区域划分包括轴向分层、径向划分子区、子区内划分计算区域、所有计算区域的角度空间离散、子区内特征线生成均需要根据不同的反应堆堆芯和不同的计算条件选择合适的方案。如图2所示，是一个压水堆堆芯布置示意图，一般压水堆堆芯轴向分为10到15层，径向每一层分为9到25个子区，每个子区中几百个计算区域不等，通常使用48到80个角度方向去离散角度空间，角度空间的离散方案确定了特征线的方向，特征线线宽一般采用0.001cm到0.05cm。

如图3所示，的角度空间离散，图中用了9个方向离散八分之一角度空间，θ代表极角，

代表辐角。

如图4所示的某区域中某一个方向的所有特征线，图中实线为特征线，虚线为两条相邻特征线的间距中心线，

是出射中子角通量，是入射中子角通量密度，s_i,k是特征线的长度，δA_k是特征线线宽。

步骤2：在步骤1所得几何尺寸、材料和特征线信息的基础上，计算每一层二维计算区域中每一个子区的矩阵MOC方法需要的系数矩阵和右端项不随着迭代计算变化的部分，系数矩阵不随迭代计算变化；

步骤2中系数矩阵的求解，需要通过对传统MOC方法求解中子输运方程的过程做一些转化。

三维直角坐标系中，稳态多群中子输运方程如公式(1)所示，

$\mathrm{\Ω}\·\▿{\mathrm{\Ψ}}_g(r,\mathrm{\Ω})+{\mathrm{\Σ}}_{t,g}\left(r\right){\mathrm{\Ψ}}_g(r,\mathrm{\Ω})=Q_g(r,\mathrm{\Ω}),$ g＝1,K,G    公式(1)

式中：

y_g——第g能群中子角通量密度

g——能群标识；

G——能群总数；

Q——中子输运方程源项。

以特征值问题为例，第g能群源项的具体形式为：

$Q_g(r,\mathrm{\Ω})=Q_g^f\left(r\right)+Q_g^s(r,\mathrm{\Ω})$

$=\frac{{\mathrm{\χ}}_g}{4{\mathrm{\πk}}_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left({\mathrm{v\Σ}}_f\left(r\right)\right)}_{g^{\′}}{\mathrm{\φ}}_{g^{\′}}\left(r\right)+\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\∫}_{{\mathrm{\Ω}}^{\′}}{\mathrm{\Σ}}_{s,g^{\′}-g}(r,{\mathrm{\Ω}}^{\′}\→\mathrm{\Ω}){\mathrm{\Ψ}}_{g^{\′}}(r,{\mathrm{\Ω}}^{\′})d{\mathrm{\Ω}}^{\′}$

公式(2)

式中：

φ_g——第g能群中子通量密度

——第g能群的裂变源项；

——第g能群的散射源项。

经过角度离散后，在第g群、第m离散方向上的形式如公式(2)所示。

${\mathrm{\ξ}}_m\frac{{\∂\mathrm{\Ψ}}_g(r,{\mathrm{\Ω}}_m)}{\∂x}{+\mathrm{\η}}_m\frac{{\∂\mathrm{\Ψ}}_g(r,{\mathrm{\Ω}}_m)}{\∂y}+{\mathrm{\μ}}_m\frac{{\∂\mathrm{\Ψ}}_g\left({r,\mathrm{\Ω}}_m\right)}{\∂z}+{\mathrm{\Σ}}_{t,g}\left(r\right){\mathrm{\Ψ}}_g\left({r,\mathrm{\Ω}}_m\right)$

$=Q_g(r,{\mathrm{\Ω}}_m)$

公式(3)

根据步骤1对于反应堆堆芯进行的几何划分，对轴向的每一层进行积分，并处以网格宽度，最终整理得：

${\mathrm{\ξ}}_m\frac{{\∂\mathrm{\Ψ}}_{g,m,k}(x,y)}{\∂x}{+\mathrm{\η}}_m\frac{{\∂\mathrm{\Ψ}}_{g,m,k}(x,y)}{\∂y}+{\mathrm{\Σ}}_{t,g,k}(x,y){\mathrm{\Ψ}}_{g,m,k}(x,y)$

$=Q_{g,k}(x,y)-{\mathrm{TL}}_{g,m,k}^{\mathrm{Axial}}(x,y)$ 公式(4)

在该方程中，定义轴向泄漏项为：

${\mathrm{TL}}_{g,m,k}^{\mathrm{Axial}}(x,y)=\frac{{\mathrm{\μ}}_m}{{\mathrm{\Δz}}_k}[{\mathrm{\Ψ}}_{g,m,k+1/2}(x,y)-{\mathrm{\Ψ}}_{g,m,k-1/2}(x,y)]$ 公式(5)

其物理意义是第k层沿轴向在离散方向m上的泄漏。在二维/一维耦合算法中，轴向泄漏项是由一维轴向求解计算提供的，因此在求解径向二维方程时，认为该项是已知源项。

MOC方法是将公式（2）转换成沿着某条特征线的微分方程，如下所示：

$\frac{{\mathrm{d\φ}}_{\mathrm{ik},g}\left(s\right)}{\mathrm{ds}}+{\mathrm{\Σ}}_{t,g}{\mathrm{\φ}}_{\mathrm{ik},g}\left(s\right)=Q_g$ 公式(6)

认为∑_t，g,Q_g在某个区域内是常数，则公式（6）是一个常微分方程，可以有解析解。步骤1中二维的计算区域在合适的范围内，MOC方法可以得到非常精确的输运解。沿着某一段特征线的平均中子角通量密度表达式如下所示：

${\stackrel{\‾}{\mathrm{\φ}}}_{\mathrm{ik},g}{s}_{i,k}=\frac{Q_g}{{\mathrm{\Σ}}_{t,g}}{s}_{i,k}+\frac{{\mathrm{\φ}}_{\mathrm{ik},g}^{\mathrm{in}}-{\mathrm{\φ}}_{\mathrm{ik},g}^{\mathrm{out}}}{{\mathrm{\Σ}}_{t,g}}$ 公式(7)

对同一方向所有特征线的平均中子角通量密度进行体积加权平均得到该计算区域这个方向的平均中子角通量密度，对该计算区域所有方向的中子角通量密度进行加权平均就可以得到该计算区域的平均中子通量密度。

由于公式（7）的右端含有源项，该源项需要通过中子通量密度进行求解，所以一般采用迭代计算求解MOC方程。

本发明中采用的区域分解矩阵MOC方法，是根据MOC方程的解，推导出以区域中子通量密度和边界入射中子角通量为未知量的线性方程组的系数矩阵和右端项；利用求解线性方程组代替传统方法中的内迭代计算。步骤2是生成矩阵MOC需要使用的系数矩阵和右端项不随迭代变化的部分。

根据球谐函数的定义，定义第(k,l)阶中子角通量密度球谐函数矩为：

${\mathrm{\φ}}_{k,l}(r,E)=\frac{1}{4\mathrm{\π}}\underset{4\mathrm{\π}}{\∫}{A}_{k,l}\left(\mathrm{\Ω}\right)\mathrm{\ψ}(r,E,\mathrm{\Ω})\mathrm{d\Ω},$ k∈N;|l|≤k    公式(8)

第(0,0)阶中子角通量密度矩具有中子标通量密度的物理含义，其值是中子标通量密度的

而第(1,1)、(1,-1)和(1,0)阶矩则分别具有沿x、y和z轴方向的净中子流密度的物理含义，同样也与对应的净中子流密度相差倍数

二维问题的z向净流恒为零，所以(1,0)阶矩为0。球谐函数A_0，0、A_1，1、A_1,-1及A_1,0也分别简记为A₀、A₁、A₂及A₃。MOC的离散格式可以转化成：

${\mathrm{\φ}}_{k,\mathrm{li},g}=\frac{1}{4\mathrm{\π}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ω}}_m{A}_{k,l}\left({\mathrm{\Ω}}_m\right){\stackrel{\‾}{\mathrm{\Ψ}}}_{i,g}\left({\mathrm{\Ω}}_m\right),$ k∈N;|l|≤k    公式(9)

由于输运方程中的源是通过中子标通量密度来求得的，所以也可以用球谐函数中子通量密度矩来表示。

裂变源：

$Q_i^f=\frac{{\mathrm{\χ}}_g}{{4\mathrm{\πk}}_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left(v{\mathrm{\Σ}}_f\right)}_{g^{\′},i}{\stackrel{\‾}{\mathrm{\φ}}}_{i,g^{\′}}=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left(v{\mathrm{\Σ}}_f\right)}_{g^{\′}i}\left(\frac{1}{4\mathrm{\π}}{\stackrel{\‾}{\mathrm{\φ}}}_{i,g^{\′}}\right)$ 公式(10)

$=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left({\mathrm{v\Σ}}_f\right)}_{g^{\′},i}{\mathrm{\φ}}_{0,i,g^{\′}}$

各项同性散射源：

$Q_i^{\mathrm{iso}}=\frac{1}{4\mathrm{\π}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s0,g^{\′}-g,i}{\stackrel{\‾}{\mathrm{\φ}}}_{i,g^{\′}}=\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s0,g^{\′}-g,i}\left(\frac{1}{4\mathrm{\π}}{\stackrel{\‾}{\mathrm{\φ}}}_{i,g^{\′}}\right)=\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s0,g^{\′}-g,i}{\mathrm{\φ}}_{0,i,g^{\′}}$ 公式(11)

各向异性散射源：

$Q_i^{\mathrm{aniso}}\left({\mathrm{\Ω}}_m\right)$

公式（12）

$=\frac{3}{\mathrm{4\π}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s1,g^{\′}-g,i}4\mathrm{\π}[A_1\left({\mathrm{\Ω}}_m\right){\mathrm{\φ}}_{1,i,g^{\′}}+A_2\left({\mathrm{\Ω}}_m\right){\mathrm{\φ}}_{2,i,g^{\′}}+A_3\left({\mathrm{\Ω}}_m\right){\mathrm{\φ}}_{3,i,g^{\′}}]$

$=3\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s1,g^{\′}-g,i}\underset{\mathrm{\ρ}=1}{\overset{3}{\mathrm{\Σ}}}{A}_{\mathrm{\ρ}}\left({\mathrm{\Ω}}_m\right){\mathrm{\φ}}_{\mathrm{\ρ},i,g^{\′}}$

根据输运方程数值求解过程的源迭代策略，按照源项是否需要在内迭代中随本群通量（或矩）更新，将以上源项分为群内源项和群外源项两类，即：

$Q_i\left({\mathrm{\Ω}}_m\right)=Q_i^{\mathrm{ingrp}}\left({\mathrm{\Ω}}_m\right)+Q_i^{\mathrm{outgrp}}\left({\mathrm{\Ω}}_m\right)$ 公式(13)

其中，群内源项定义为：

$Q_i^{\mathrm{ingrp}}\left({\mathrm{\Ω}}_m\right)=Q_i^{\mathrm{iso},g-g}+Q_i^{\mathrm{aniso},g-g}\left({\mathrm{\Ω}}_m\right)$ 公式(14)

群外源项则定义为：

$Q_i^{\mathrm{outgrp}}\left({\mathrm{\Ω}}_m\right)=Q_i^f+Q_i^{\mathrm{iso},g^{\′}-g}+Q_i^{\mathrm{aniso},g^{\′}-g}\left({\mathrm{\Ω}}_m\right)$ 公式(15)

从角通量密度矩的定义出发，根据传统MOC的基本理论，得到以平源网格区各阶角通量密度矩和外边界入射角通量为未知量的矩相关方程。

公式(16)

ρ＝0,1,2,3

外边界入射角通量方程的建立过程及最终形式是与外边界条件密切相关的，求出出射角通量，根据出射角通量和入射角通量的关系得到外边界入射角通量方程。区域分解算法中将子区内边界和子区外边界做了区分，得到子区内边界入射方程和子区外边界入射方程，如下所示：

公式(17)

公式(18)

其中：

φ_ρ，i——中子通量密度球谐函数矩，

——子区外边界入射中子角通量密度，

——子区内边界入射中子角通量密度，

$a_{0,i}=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{s0,g-g,i}}{{\mathrm{\Σ}}_{t,i}},& \mathrm{\ρ}=0\\ 0,& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.$

$a_{{\mathrm{\ρ}}^{\′},i}=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{s1,g-g,i}}{{\mathrm{\Σ}}_{t,i}}& {\mathrm{\ρ}}^{\′}\\ 0,& {\mathrm{\ρ}}^{\′}\≠\mathrm{\ρ}\end{array}\right.$

${\mathrm{\τ}}_{\mathrm{\ρ},{\mathrm{\ρ}}^{\′},j,i}^l={\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l{\mathrm{\σ}}_{{\mathrm{\ρ}}^{\′},j}^l$

${\mathrm{\ϵ}}_{\mathrm{\ρ},i}^l={\mathrm{\λ}}_{\mathrm{\ρ},i}^l{\mathrm{\δe}}_{1,i}^l$

$q_{\mathrm{\ρ},i}^{\mathrm{loc}}=\left\{\begin{array}{cc}\frac{Q_i^f+Q_i^{\mathrm{iso},g^{\′}-g}}{{\mathrm{\Σ}}_{t,i}},& \mathrm{\ρ}=0\\ \frac{\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}}{{\mathrm{\Σ}}_{t,i}}& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.$

$q_{\mathrm{\ρ},i}^{\mathrm{trans}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l{Q}_j^{\mathrm{outgrp}}\left({\mathrm{\Ω}}_m\right)$

δ_k0——克罗内特尔向量

$e_{j,p}^l=\left\{\begin{array}{cc}\exp (-\underset{n=j}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{t,n}{\hat{s}}_{\mathrm{nl}})& j\≤p\\ 1,& j>p\end{array}\right.$

${\mathrm{\σ}}_{{\mathrm{\ρ}}^{\′},j}^l=\left\{\begin{array}{cc}{\mathrm{\Σ}}_{s0,g-g,j},& {\mathrm{\ρ}}^{\′}=0\\ {3A}_{{\mathrm{\ρ}}^{\′}}\left({\mathrm{\Ω}}_m\right){\mathrm{\Σ}}_{s1,g-g,j},& {\mathrm{\ρ}}^{\′}=\mathrm{1,2,3}\end{array}\right.$

$c_{j,i}^l=\frac{{\mathrm{\δe}}_{j+1,i}^l-{\mathrm{\δe}}_{j,i}^l}{{\mathrm{\Σ}}_{t,j}}$

${\mathrm{\δe}}_{j,p}^l=e_{j,p-1}^l-e_{j,p}^l$

任意子区域sub上的线性代数方程组便可写为如下形式：

$A_{\mathrm{sub}}\left[\begin{array}{c}{\stackrel{v}{\mathrm{\φ}}}_{\mathrm{sub}}\\ {\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{ob},-}\\ {\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},+}\end{array}\right]={\stackrel{v}{b}}_{\mathrm{sub}}\left({\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},-}\right)$ 公式(19)

式中：

A_sub——子区域sub的方程组系数矩阵；

——子区域sub的角通量密度矩向量；

——子区域sub的外边界入射角通量向量（若子区域sub不存在外边界，则不含此向量）；

——子区域sub的内边界出射角通量向量；

——子区域sub的方程组右端项，它与内边界入射角通量向量ψ

相关。

将右端项拆解成以下形式：

${\stackrel{v}{b}}_{(N{\mathrm{\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×1}=R_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{N\φ}}_{\mathrm{sub}}}{Q}_{{\mathrm{N\φ}}_{\mathrm{sub}}\×1}+E_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}{+\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}{\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},-}$

公式(20)

或者，

$\left[\begin{array}{c}{B}_{N{\mathrm{\φ}}_{\mathrm{sub}}\×1}^1\\ B_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×1}^2\\ B_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}\×1}^3\end{array}\right]=\left[\begin{array}{c}{R}_{{\mathrm{N\φ}}_{\mathrm{sub}}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^1\\ R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^2\\ R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^3\end{array}\right]Q_{{\mathrm{N\φ}}_{\mathrm{sub}}\×1}+\left[\begin{array}{c}{E}_{{\mathrm{N\φ}}_{\mathrm{sub}}\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}^1\\ E_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+}\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}^2\\ E_{N_{\mathrm{sub}}^{\mathrm{inb},+}\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}^2\end{array}\right]{\stackrel{v}{\mathrm{\Ψ}}}_{\mathrm{sub}}^{\mathrm{inb},-}$ 公式(21)

式中：

Nδ_sub——子区域角通量密度矩向量

的长度；

——子区域外边界入射角通量向量ψ

的长度；

——子区域内边界出射角通量向量ψ

的长度；

$Q_{{\mathrm{N\φ}}_{\mathrm{sub}}\×1}={\left[Q_{\mathrm{\ρi}}\right]}_{\mathrm{\ρ}\∈\left[\mathrm{0,3}\right],i\∈[1,\mathrm{Nr}]}$

$Q_{\mathrm{\ρ}}=\left\{\begin{array}{cc}{Q}_i^f+Q_i^{\mathrm{iso},g^{\′}-g}& \mathrm{\ρ}=0\\ \underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s1,g^{\′}-g,i}{\mathrm{\φ}}_{\mathrm{\ρ},i,g^{\′}},& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.$ 右端项随着迭代变化的部分

步骤2中求解的就是每一层中每个子区的系数矩阵A_sub和不随迭代变化的模板矩阵 $R_{({\mathrm{N\φ}}_{\mathrm{dub}}+{\mathrm{N\ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{N\φ}}_{\mathrm{sub}}},$ $E_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},-}}.$

A_sub具体求解由前面的子区矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程得到。

$R_{({\mathrm{N\φ}}_{\mathrm{sub}}+{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{inb},+})\×{\mathrm{N\φ}}_{\mathrm{sub}}}$ 求解如下：

$R_{{\mathrm{N\φ}}_{\mathrm{sub}}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^1=\left[\begin{array}{ccccccccccc}L& M& M& M& M& L& M& M& M& M& L\\ L& {\hat{r}}_i^1+r_{i,i,\mathrm{0,0}}^1& r_{i,i,\mathrm{0,1}}^1& r_{i,i,\mathrm{0,2}}^1& r_{i,i,\mathrm{0,3}}^1& L& r_{i,j,\mathrm{0,0}}^1& r_{i,j,0,1}^1& r_{i,j,\mathrm{0,2}}^1& r_{i,j,\mathrm{0,3}}^1& L\\ L& r_{i,i,\mathrm{1,0}}^1& {\hat{r}}_i^1+r_{i,i,\mathrm{1,1}}^1& r_{i,i,\mathrm{1,2}}^1& r_{i,i,\mathrm{1,3}}^1& L& r_{i,j,\mathrm{1,0}}^1& r_{i,j,\mathrm{1,1}}^1& r_{i,j,\mathrm{1,2}}^1& r_{i,j,\mathrm{1,3}}^1& L\\ L& r_{i,i,\mathrm{2,0}}^1& r_{i,i,\mathrm{2,1}}^1& {\hat{r}}_i^1+r_{i,i,\mathrm{2,2}}^1& r_{i,i,\mathrm{2,3}}^1& L& r_{i,j,\mathrm{2,0}}^1& r_{i,j,\mathrm{2,1}}^1& r_{i,j,\mathrm{2,2}}^1& r_{i,j,\mathrm{2,3}}^1& L\\ L& r_{i,i,\mathrm{3,0}}^1& r_{i,i,\mathrm{3,1}}^1& r_{i,i,\mathrm{3,2}}^1& {\hat{r}}_i^1+r_{i,i,\mathrm{3,3}}^1& L& r_{i,j,\mathrm{3,0}}^1& r_{i,j,\mathrm{3,1}}^1& r_{i,j,\mathrm{3,2}}^1& r_{i,j,\mathrm{3,3}}^1& L\\ L& M& M& M& M& L& M& M& M& M& L\\ L& r_{j,i,\mathrm{0,0}}^1& r_{j,i,\mathrm{0,1}}^1& r_{j,i,\mathrm{0,2}}^1& r_{j,i,\mathrm{0,3}}^1& L& {\hat{r}}_j^1+r_{j,j,\mathrm{0,0}}^1& r_{j,j,\mathrm{0,1}}^1& r_{j,j,\mathrm{0,2}}^1& r_{j,j,\mathrm{0,3}}^1& L\\ L& r_{j,i,\mathrm{1,0}}^1& r_{j,i,\mathrm{1,1}}^1& r_{j,i,\mathrm{1,2}}^1& r_{j,i,\mathrm{1,3}}^1& L& r_{j,j,\mathrm{1,0}}^1& {\hat{r}}_j^1+r_{j,j,\mathrm{1,1}}^1& r_{j,j,\mathrm{1,2}}^1& r_{j,j,\mathrm{1,3}}^1& L\\ L& r_{j,i,\mathrm{2,0}}^1& r_{j,i,\mathrm{2,1}}^1& r_{j,i,\mathrm{2,2}}^1& r_{j,i,\mathrm{2,3}}^1& L& r_{j,j,\mathrm{2,0}}^1& r_{j,j,2,1}^1& {\hat{r}}_j^1+r_{j,j,\mathrm{2,2}}^1& r_{j,j,\mathrm{2,3}}^1& L\\ L& r_{j,i,\mathrm{3,0}}^1& r_{j,i,\mathrm{3,1}}^1& r_{j,i,\mathrm{3,2}}^1& r_{j,i,\mathrm{3,3}}^1& L& r_{j,j,\mathrm{3,0}}^1& r_{j,j,\mathrm{3,1}}^1& r_{j,j,\mathrm{3,2}}^1& {\hat{r}}_j^1+r_{j,j,\mathrm{3,3}}^1& L\\ L& M& M& M& M& L& M& M& M& M& L\end{array}\right]$

公式(22)其中的元素表达形式为：

${\hat{r}}_i^1=\frac{1}{{\mathrm{\Σ}}_{t,i}}$

$r_{i,j,\mathrm{\ρ},{\mathrm{\ρ}}^{\′}}^1=\left\{\begin{array}{cc}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}\left({\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l\right),& {\mathrm{\ρ}}^{\′}=0\\ \underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{3A}_{{\mathrm{\ρ}}^{\′}}\left({\mathrm{\Ω}}_m\right){\mathrm{\λ}}_{\mathrm{\ρ},i}^l{c}_{j,i}^l,& {\mathrm{\ρ}}^{\′}=\mathrm{1,2,3}\end{array}\right.$

$R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^2=\left[\begin{array}{cccccc}L& M& M& M& M& L\\ L& r_{l_r,j,0}^2& r_{l_r,j,1}^2& r_{l_r,j,2}^2& r_{l_r,j,3}^2& L\\ L& M& M& M& M& L\end{array}\right]$ 公式(23)其中的元素表达形式为：

$R_{{\mathrm{N\Ψ}}_{\mathrm{sub}}^{\mathrm{ob},-}\×{\mathrm{N\φ}}_{\mathrm{sub}}}^3=\left[\begin{array}{cccccc}L& M& M& M& M& L\\ L& r_{l_r,j,0}^3& r_{l_r,j,1}^3& r_{l_r,j,2}^3& r_{l_r,j,3}^3& L\\ L& M& M& M& M& L\end{array}\right]$ 公式(24)其中的元素表达形式为：

则由矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程右端第二项的系数求的。

由前面的推导可知，系数矩阵A只与特征线几何、材料等信息有关，而这些信息在中子输运方程的数值求解过程中是不需要随着迭代而更新的。因此，可预先对没个子区每一能群仅做一次全局特征线扫描，即可生成各个子区各群的系数矩阵和右端项不随迭代变化的部分。

步骤3：用步骤1中子通量密度初值计算三维粗网有限差分CMFD需要的净流修正系数，进行一次全堆芯三维粗网有限差分CMFD计算，求得粗网中子通量密度后，修正一维、二维的细网中子通量密度，所述二维的细网中子通量密度修正如公式（25）所示：

${\mathrm{\φ}}_g^{i,(l+1/2)}=\frac{{\stackrel{\‾}{\mathrm{\φ}}}_g^{n,(l+1/2)}}{{\stackrel{\‾}{\mathrm{\φ}}}_g^{n,\left(l\right)}}{\mathrm{\φ}}_g^{i,\left(l\right)},$ i∈n公式(25)

式中：

——修正后的平源细区i第g群的平均标通量，用以更新第l+1次矩阵MOC计算的源项，中子通量密度的单位是cm^-2·s^-1，

——第l次矩阵MOC计算后的三维CMFD计算得到的粗网n的第g群平均标通量，

——由第l次矩阵MOC计算得到的粗网n的第g群均匀化标通量，

——第l次矩阵MOC计算所得到的平源细区i第g群的平均标通量；

对于一维的细网中子通量密度修正如公式（26）所示：

$a_{\mathrm{gi}}^{k,(l+1/2)}=\frac{{\stackrel{\‾}{\mathrm{\φ}}}_g^{n,(l+1/2)}}{a_{g0}^{k,\left(l\right)}}{a}_{\mathrm{gi}}^{k,\left(l\right)},$ k∈n,i=0,1,2,3,4    公式(26)

式中：

——修正后的一维节块k第g群的第i阶通量展开系数，用以更新第l+1次NEM计算的源项；

——第l次一维NEM计算后的三维CMFD计算得到的粗网n的第g群平均标通量；

——第l次一维NEM计算得到节块k内第g群的第i阶中子通量密度展开系数；

计算三维粗网有限差分CMFD需要的净流修正系数的公式为：

公式(27)

公式(28)

${\stackrel{\‾}{D}}_g^n=\frac{1}{3{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{tr},g}^n}$ 公式(29)

式中：

——高阶计算得到的粗网n的界面净流。

——第g群的宏观移除截面

——网格间距

三维多群CMFD某一群方程如下所示：

$d_{\mathrm{zg}-}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n-1_z}+d_{\mathrm{yg}-}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n-1_y}+d_{\mathrm{xg}-}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n-1_x}+({\stackrel{\‾}{\mathrm{\Σ}}}_{r,g}^n+d_{\mathrm{xg}}^n+d_{\mathrm{yg}}^n+d_{\mathrm{zg}}^n){\stackrel{\‾}{\mathrm{\φ}}}_g^n$

公式(30)

$-d_{\mathrm{xg}+}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n+1_x}-d_{\mathrm{yg}+}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n+1_y}{-d}_{\mathrm{zg}+}^n{\stackrel{\‾}{\mathrm{\φ}}}_g^{n+1_z}=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left(v{\stackrel{\‾}{\mathrm{\Σ}}}_f\right)}_{g^{\′}}^n{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^n+\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\stackrel{\‾}{\mathrm{\Σ}}}_{s,g^{\′}g}^n{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^n$

其中：

${\stackrel{\‾}{D}}_g^n=\frac{1}{3{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{tr},g}^n}$

是中子通量密度

是扩散系数

是中子裂变产额截面

是中子散射截面

xg是第g群裂变谱

k_eff是特征值，也称为有效增殖因数

CMFD群内方程组的系数矩阵是一个七对角矩阵，并且规模往往较大，因此采用了具有广泛适用性的GMRES解法。依次求解每一群，可得到三维粗网通量和特征值，迭代收敛后，用三维粗网差分求得的粗网通量修正一维、二维通量，耦合一维、二维计算同时加速整个计算。

步骤4：更新一维节块通量展开需要的展开系数和径向泄漏项，对每一个栅元逐层、每一能群地进行一维节块展开NEM固定源计算，将每个节块的中子通量密度和中子源项用多项式展开，通过求解多项式的展开系数，求得每个栅元的一维中子通量密度，具体的迭代过程如图5所示：

1）更新径向泄漏项展开系数

由以下公式求得：

$l_{g0}^k={\mathrm{TL}}_{g,k}^{\mathrm{Radial}}$ 公式(31)

$l_{g1}^k=\frac{{\mathrm{\Δz}}_k}{d}\left[z_{k-1}({\mathrm{TL}}_{g,k+1}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k}^{\mathrm{Radial}}){+z}_{k+1}\left({\mathrm{TL}}_{g,k}^{\mathrm{Radial}}{-\mathrm{TL}}_{g,k-1}^{\mathrm{Radial}}\right)\right]$ 公式(32)

$l_{g2}^k=\frac{{\left({\mathrm{\Δz}}_k\right)}^2}{d}[({\mathrm{\Δz}}_{k-1}+{\mathrm{\Δz}}_k)({\mathrm{TL}}_{g,k+1}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k}^{\mathrm{Radial}})-({\mathrm{\Δz}}_{k+1}+{\mathrm{\Δz}}_k)({\mathrm{TL}}_{g,k}^{\mathrm{Radial}}-{\mathrm{TL}}_{g,k-1}^{\mathrm{Radial}})]$

公式(33)

其中：

z_i＝(Δz_k+Δz_i)(Δz_k+2Δz_i)

d＝(Δz_k-1+Δz_k)(Δz_k+Δz_k+1)(Δz_k-1+Δz_k+Δz_k+1)

${\mathrm{TL}}_{g,k}^{\mathrm{Radial}}=\frac{{\mathrm{\ξ}}_m}{\mathrm{\Δx}}({\mathrm{\Ψ}}_{g,x+}-{\mathrm{\Ψ}}_{g,x-})+\frac{{\mathrm{\η}}_m}{\mathrm{\Δy}}\left({\mathrm{\Ψ}}_{g,y+}{-\mathrm{\Ψ}}_{g,y-}\right)$

Δz_k、Δx、Δy分别是x、y、z方向网格间距，

${\mathrm{\Ψ}}_{g,x\±}^p\left(z\right)=\frac{1}{\mathrm{\Δy}}{\∫}_{y_{j-1/2}}^{y_{j+1/2}}{\mathrm{\Ψ}}_g(x_{i\±1/2},y,z)\mathrm{dy}$

${\mathrm{\Ψ}}_{g,y\±}^p\left(z\right)=\frac{1}{\mathrm{\Δx}}{\∫}_{x_{j-1/2}}^{x_{j+1/2}}{\mathrm{\Ψ}}_g(x_{i\±1/2},y,z)\mathrm{dx}$

ξm、ηm分别是x、y方向的方向余弦

2）假定初值：

（g=1,K,G；k=1,K,K;n=0,K,4）；

其中：

表示节块k第g能群轴向的上表面和下表面净流值，净流的单位是cm^-2·s^-1，

表示节块k第g能群中子通量密度展开多项式第一项系数，

中子通量密度多项式展开如下所示：

${\mathrm{\φ}}_g^k\left(z\right)=\underset{n=0}{\overset{4}{\mathrm{\Σ}}}{a}_{\mathrm{gn}}^k{H}_n\left(z\right)$ 公式(34)

其中：

$H_n\left(z\right)=\left\{\begin{array}{ccccc}1,& \mathrm{\γ},& {3\mathrm{\γ}}^2-\frac{1}{4},& {\mathrm{\γ}}^3-\frac{\mathrm{\γ}}{4}& {\mathrm{\γ}}^4-\frac{{3\mathrm{\γ}}^2}{10}+\frac{1}{80}\end{array}\right\}$

z∈[z_k-1/2,z_k+1/2]

k表示节块编号，g表示能群编号，n表示展开阶数；

3）构造源项展开系数

一维节块展开计算NEM中源项多项式展开如下所示：

$Q_g^k\left(z\right)=\underset{n=0}{\overset{4}{\mathrm{\Σ}}}{q}_{\mathrm{gn}}^k{H}_n\left(z\right)$ 公式(35)

其中：

$q_{\mathrm{gn}}^k=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left({\mathrm{v\Σ}}_g\right)}_{g^{\′},k}{a}_{g^{\′}n}^k+\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s,g^{\′}-g,k}{a}_{g^{\′}n}^k,$ n＝0,1,2,3,4

是节块k第g能群中子通量密度多项式展开的系数，求解公式如下所示：

公式(36)

$a_{g1}^k={\mathrm{\φ}}_{g+}^k-{\mathrm{\φ}}_{g-}^k$ 公式(37)

$a_{g2}^k={\mathrm{\φ}}_{g+}^k+{\mathrm{\φ}}_{g-}^k-2{\stackrel{\‾}{\mathrm{\φ}}}_g^k$ 公式(38)

$a_{g3}^k=b_{g3}^k{a}_{g1}^k+c_{g3}^k$ 公式(39)

公式(39)中：

$b_{g3}^k=\frac{{\mathrm{\Σ}}_{r,g,k}}{\frac{{6D}_g^k}{{\left({\mathrm{\Δz}}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{10}}$

$c_{g3}^k\frac{\frac{q_{g3}^k}{10}-q_{g1}^k}{\frac{{6D}_g^k}{{\left(\mathrm{\Δ}{z}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{10}}$

$a_{g4}^k=b_{g4}^k{a}_{g2}^k+c_{g4}^k$ 公式(40)公式(40)中：

$b_{g4}^k=\frac{{\mathrm{\Σ}}_{r,g,k}}{\frac{{4D}_g^k}{{\left({\mathrm{\Δz}}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{35}}$

$c_{g4}^k\frac{\frac{q_{g4}^k}{35}-q_{g2}^k}{\frac{{4D}_g^k}{{\left(\mathrm{\Δ}{z}_k\right)}^2}+\frac{{\mathrm{\Σ}}_{r,g,k}}{35}}$

其中：

是节块k第g群内的中子通量密度，

是节块k第g群上下界面的中子通量密度，

∑_r，g,k是节块k第g群移除截面，单位是cm^-1

节块k第g群扩散系数，单位是cm

4）求解由两个边界条件方程以及净中子流耦合方程构成的线性代数方程组，得到

三个方程如下所示：

$A_{g2}^k{J}_{\mathrm{zg}-}^k-(A_{g1}^k+A_{g1}^{k+1})J_{\mathrm{zg}+}^k+A_{g2}^{k+1}{J}_{\mathrm{zg}+}^{k+1}=A_{g0}^k-A_{g3}^{k+1},$ k＝1,K-1公式(41) $A_{g2}^K{J}_{\mathrm{zg}-}^K+[4(\frac{{\mathrm{\β}}^{+}}{1-{\mathrm{\β}}^{+}}+\frac{1}{2})-A_{g1}^K]J_{\mathrm{zg}+}^K=A_{g0}^K$ 公式(42) $[4(\frac{{\mathrm{\β}}^{-}}{{\mathrm{\β}}^{-}-1}-\frac{1}{2})+A_{g1}^1]J_{\mathrm{zg}-}^1-A_{g2}^1{J}_{\mathrm{zg}+}^1=A_{g3}^1$ 公式(43)其中：

表示节块k第g能群轴向的上表面和下表面净中子流值

$\left\{\begin{array}{c}{A}_{g0}^k=\frac{{\hat{a}}_{g2}^k{\hat{r}}_{g2}^{\mathrm{kg}}-{\hat{a}}_{g1}^k{\hat{r}}_{g1}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\\ A_{g1}^k=\frac{\frac{{\mathrm{\Δz}}_k}{D_g^k}{\hat{a}}_{g2}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g2}^k\right)}^2}\\ A_{g2}^k=\frac{\frac{{\mathrm{\Δz}}_k}{D_g^k}{\hat{a}}_{g2}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\\ A_{g3}^k=\frac{{\hat{a}}_{g1}^k{\hat{r}}_{g2}^k-{\hat{a}}_{g2}^k{\hat{r}}_{g1}^k}{{\left({\hat{a}}_{g2}^k\right)}^2-{\left({\hat{a}}_{g1}^k\right)}^2}\end{array}\right.$

$\left\{\begin{array}{c}{\hat{a}}_{g1}^k=\frac{1}{2}{b}_{g3}^k+\frac{1}{5}{b}_{g4}^k+4\\ {\hat{a}}_{g2}^k=\frac{1}{2}{b}_{g3}^k-\frac{1}{5}{b}_{g4}^k-2\\ {\hat{r}}_{g1}^k=(6+\frac{2}{5}{b}_{g4}^k){\stackrel{\‾}{\mathrm{\φ}}}_g^k-\frac{1}{2}{c}_{g3}^k-\frac{1}{5}{c}_{g4}^k\\ {\hat{r}}_{g2}^k=-{\hat{r}}_{g1}^k{-c}_{g3}^k\end{array}\right.$

其中

Δz_k由上一步计算求得；

5）求解下述中子平衡方程，得到节块平均通量

$\frac{1}{{\mathrm{\Δz}}_k}(J_{\mathrm{zg}+}^k-J_{\mathrm{zg}-}^k)+{\mathrm{\Σ}}_{r,g,k}{\stackrel{\‾}{\mathrm{\φ}}}_g^k=\frac{{\mathrm{\χ}}_g}{k_{\mathrm{eff}}}\underset{g^{\′}=1}{\overset{G}{\mathrm{\Σ}}}{\left({\mathrm{v\Σ}}_f\right)}_{g^{\′},k}{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^k+\underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s,g^{\′}-g,k}{\stackrel{\‾}{\mathrm{\φ}}}_{g^{\′}}^k+{\mathrm{TL}}_{g,k}^{\mathrm{Radial}}\left(z\right)$ 公式(44)

6）由节块表面中子通量与净中子流关系式，求得节块界面中子通量

${\mathrm{\φ}}_{g+}^k=A_{g1}^k{J}_{\mathrm{zg}+}^k{J}_{\mathrm{zg}+}^k-A_{g2}^k{J}_{\mathrm{zg}-}^k+A_{g0}^k$ 公式(45)

${\mathrm{\φ}}_{g-}^k=A_{g2}^k{J}_{\mathrm{zg}+}^k-A_{g1}^k{J}_{\mathrm{zg}-}^k{+A}_{g3}^k$ 公式(46)

7）更新通量高阶展开系数

如步骤3）中所示；

8）对每一个栅元的每一个能群进行3）到6）迭代计算，若通量收敛则结束计算，收敛准则如下：

$\left|\right|\frac{{\stackrel{\‾}{\mathrm{\φ}}}_g^k-{\stackrel{\‾}{\mathrm{\φ}}}_g^{k-1}g}{{\stackrel{\‾}{\mathrm{\φ}}}_g^k}\left|\right|\≤{\mathrm{\ϵ}}_1$ 公式(47)

ε₁为很小的正数，也就是收敛条件；

步骤5：更新轴向泄漏项，根据步骤1划分的几何区域，对轴向的每一层进行二维输运计算，根据步骤2所得到的每个子区的系数矩阵，采用GMRES算法，并行求解每个子区每一群的中子通量密度，在边界处通过插值互相耦合，迭代求解每一层的二维中子通量密度，具体步骤如图6所示：

1）更新轴向泄漏项

该泄漏项由一维计算提供；

2）读取步骤2产生的子区系数矩阵和右端项不随迭代改变部分，设定子区内、外边界入射角通量和子区内中子通量密度初始值；

3）生成子区群外源，即右端项随着迭代变化的部分

$Q_{\mathrm{\ρi}}=\left\{\begin{array}{cc}{Q}_i^f+Q_i^{\mathrm{iso},g^{\′}-g},& \mathrm{\ρ}=0\\ \underset{\underset{g^{\′}\≠g}{g^{\′}=1}}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{s1,g^{\′}-g,i}{\mathrm{\φ}}_{\mathrm{\ρ},i,g^{\′}},& \mathrm{\ρ}=\mathrm{1,2,3}\end{array}\right.$

4）采用GMRES方法求解子区群内方程；

5）与相邻子区进行通信，交换内边界信息，即交换内边界出射角通量信息；

6）判断子区通量是否收敛，若不收敛，继续迭代求解；

7）判断所有能群是否都已经计算，对各个能群进行迭代，直到所有区域通量收敛；

步骤6：用步骤4得到的一维通量密度和步骤5得到的二维通量密度更新三维粗网参数，判断特征值和三维中子通量密度是否收敛，如果不收敛，则转到步骤3继续迭代，并且计算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收敛，就可以得到三维中子通量密度精细分布；

所述收敛条件为：

$\left|\right|\frac{{\mathrm{\φ}}_g^{i,\left(l\right)}{-\mathrm{\φ}}_g^{i,(l-1)}}{{\mathrm{\φ}}_g^{i,\left(l\right)}}\left|\right|{\≤\mathrm{\ϵ}}_1,$ $\left|\right|\frac{k_{\mathrm{eff}}^{\left(l\right)}{-k}_{\mathrm{eff}}^{(l-1)}}{k_{\mathrm{eff}}^{\left(l\right)}}\left|\right|{\≤\mathrm{\ϵ}}_2$ 公式(48)

其中：

为第l次迭代第i区第g群通量

为第l次迭代特征值

ε₁，ε₂为两个很小的正数，也就是收敛条件。

大量数值验证结果显示，本发明具有可靠的精度、很好的效率和很好的几何适应性，适应工程实际中的计算要求。图7是三维非均匀压水堆堆芯相对棒功率分布，由于外围的反射层和组件内的控制棒没有功率，所以图中反射层对应的是外围的黑色区域，控制棒对应的是组件内黑色的棒状。从图中可以看出径向和轴向的网格划分。每个区域内的颜色深浅代表了功率大小，图7就是利用本发明得到的反应堆堆芯三维中子通量密度精细分布简单计算得来的。

Claims (4)

1.一种获取反应堆堆芯三维中子通量密度精细分布的方法，其特征在于：包括如下步骤： 

步骤1：对反应堆堆芯进行几何建模，将反应堆堆芯按照轴向划分成多层，用Auto-CAD分别描述反应堆堆芯每一层几何尺寸，将每一层沿径向划分成多个子区，再用Auto-CAD对每个子区划分计算区域，根据计算需求对每一层计算区域的三维角度空间进行离散，利用VBA二次开发的Auto-CAD的工程，在每一个子区内生成特征线信息，子区内边界通过插值进行耦合；根据反应堆堆芯材料信息，读取各层中每一个计算区域的宏观截面参数；对二维、一维计算区域的中子通量密度、反应堆边界条件和特征值赋初值； 

步骤2：在步骤1所得几何尺寸、材料和特征线信息的基础上，计算每一层二维计算区域中每一个子区的矩阵MOC方法需要的系数矩阵和右端项不随着迭代计算变化的部分，系数矩阵不随迭代计算变化； 

步骤3：用步骤1中子通量密度初值计算三维粗网有限差分CMFD需要的净流修正系数，进行一次全堆芯三维粗网有限差分CMFD计算，求得粗网中子通量密度后，修正一维、二维的细网中子通量密度，所述二维的细网中子通量密度修正如公式（25）所示： 

i∈n    公式(25) 

式中： 

——修正后的平源细区i第g群的平均标通量，用以更新第l+1次矩阵MOC计算的源项，中子通量密度的单位是cm^-2·s^-1， 

——第l次矩阵MOC计算后的三维CMFD计算得到的粗网n的第g群平均标通量， 

——由第l次矩阵MOC计算得到的粗网n的第g群均匀化标通量， 

——第l次矩阵MOC计算所得到的平源细区i第g群的平均标通量； 

对于一维的细网中子通量密度修正如公式（26）所示： 

k∈n,i＝0,1,2,3,4    公式(26) 

式中： 

——修正后的一维节块k第g群的第i阶通量展开系数，用以更新第l+1次NEM计算的源项； 

——第l次一维NEM计算后的三维CMFD计算得到的粗网n的第g群平均标通量； 

——第l次一维NEM计算得到节块k内第g群的第i阶中子通量密度展开系数； 

步骤4：更新一维节块通量展开需要的展开系数和径向泄漏项，对每一个栅元逐层、每一能群地进行一维节块展开NEM固定源计算，将每个节块的中子通量密度和中子源项用多项式展开，通过求解多项式的展开系数，求得每个栅元的一维中子通量密度，具体的迭代过程如下所示： 

1）更新径向泄漏项展开系数

由以下公式求得： 

公式(31) 

公式(32) 

公式(33) 

其中： 

z_i＝(Δz_k+Δz_i)(Δz_k+2Δz_i) 

d=(Δz_k-1+Δz_k)(Δz_k+Δz_k+1)(Δz_k-1+Δz_k+Δz_k+1) 

Δz_k、Δx、Δy分别是x、y、z方向网格间距， 

ξ_m、η_m分别是x、y方向的方向余弦 

2）假定初值：

（g＝1,K,G；k＝1,K,K;n=0,K,4）； 

其中： 

表示节块k第g能群轴向的上表面和下表面净流值，净流的单 

位是cm^-2·s^-1， 

表示节块k第g能群中子通量密度展开多项式第一项系数， 

中子通量密度多项式展开如下所示： 

公式(34) 

其中： 

z∈[zk-1/2,zk+1/2] 

k表示节块编号，g表示能群编号，n表示展开阶数； 

3）构造源项展开系数

一维节块展开计算NEM中源项多项式展开如下所示： 

公式(35) 

其中： 

n=0,1,2,3,4 

是节块k第g能群中子通量密度多项式展开的系数，求解公式如下所示： 

公式(36) 

公式(37) 

公式(38) 

公式(39) 

公式(1)中： 

公式(40)公式(2)中： 

其中： 

是节块k第g群内的中子通量密度， 

是节块k第g群上下界面的中子通量密度， 

∑_r，g,k是节块k第g群移除截面，单位是cm^-1

节块k第g群扩散系数，单位是cm 

4）求解由两个边界条件方程以及净中子流耦合方程构成的线性代数方程组，得到三个方程如下所示： 

k＝1,K-1公式(41) 

公式(42) 

公式(43)其中： 

表示节块k第g能群轴向的上表面和下表面净中子流值 

其中

Δz_k由上一步计算求得； 

5）求解下述中子平衡方程，得到节块平均通量

公式(44) 

6）由节块表面中子通量与净中子流关系式，求得节块界面中子通量 

公式(45) 

公式(46) 

7）更新通量高阶展开系数如步骤3）中所示； 

8）对每一个栅元的每一个能群进行3）到6）迭代计算，若通量收敛则结束计算，收敛准则如下： 

公式(47) 

ε₁为很小的正数，也就是收敛条件； 

步骤5：更新轴向泄漏项，根据步骤1划分的几何区域，对轴向的每一层进行二维输运计算，根据步骤2所得到的每个子区的系数矩阵，采用GMRES算法，并行求解每个子区每一群的中子通量密度，在边界处通过插值互相耦合，迭代求解每一层的二维中子通量密度，具体步骤如下： 

1）更新轴向泄漏项该泄漏项由一维计算提供； 

2）读取步骤2产生的子区系数矩阵和右端项不随迭代改变部分，设定子区内、外边界入射角通量和子区内中子通量密度初始值； 

3）生成子区群外源，即右端项随着迭代变化的部分 

4）采用GMRES方法求解子区群内方程； 

5）与相邻子区进行通信，交换内边界信息，即交换内边界出射角通量信息； 

6）判断子区通量是否收敛，若不收敛，继续迭代求解； 

7）判断所有能群是否都已经计算，对各个能群进行迭代，直到所有区域通量收敛； 

步骤6：用步骤4得到的一维通量密度和步骤5得到的二维通量密度更新三维粗网参数，判断特征值和三维中子通量密度是否收敛，如果不收敛，则转到步骤3继续迭代，并且计算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收敛，就可以得到三维中子通量密度精细分布； 

所述收敛条件为： 

公式(48) 

其中： 

为第l次迭代第i区第g群通量 

为第l次迭代特征值 

ε₁，ε₂为两个很小的正数，也就是收敛条件。 

2.根据权利要求1所述的一种获取反应堆堆芯三维中子通量密度精细分布的方法，其特征在于：步骤2所述计算每一层二维计算区域中每一个子区的矩阵MOC方法需要的系数矩阵和右端项不随着迭代计算变化的部分的计算方法为： 

任意子区域sub上的矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程如下式所示： 

公式(16) 

ρ＝0,1,2,3 

公式(17) 

公式(18) 

其中： 

φ_ρ，i是中子通量密度球谐函数矩， 

是子区外边界入射中子角通量密度， 

是子区内边界入射中子角通量密度， 

δ_k0克罗内特尔向量 

任意子区域sub上的线性代数方程组便可写为如下形式： 

公式(19) 

式中： 

A_sub——子区域sub的方程组系数矩阵； 

——子区域sub的角通量密度矩向量； 

——子区域sub的外边界入射角通量向量（若子区域sub不存在外边界，则不含此向量）； 

——子区域sub的内边界出射角通量向量； 

——子区域sub的方程组右端项，它与内边界入射角通量向量 

相关； 

将右端项拆解成以下形式： 

公式(20) 

或者， 

公式(21) 

式中： 

Nφ_sub——子区域角通量密度矩向量

的长度； 

——子区域外边界入射角通量向量ψ

的长度； 

——子区域内边界出射角通量向量ψ的长度； 

右端项随迭代变化部分， 

步骤2中求解的就是每一层中每个子区的系数矩阵A_sub和不随迭代变化的模板矩阵

A_sub具体求解由前述的子区矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程得到； 

求解如下： 

公式(22) 

其中的元素表达形式为： 

公式(23) 

其中的元素表达形式为： 

公式(24) 

其中的元素表达形式为： 

则由矩相关方程、子区内边界出射角通量方程和子区外边界出射角通量方程右端第二项的系数求的。 

3.根据权利要求1所述的一种获取反应堆堆芯三维中子通量密度精细分布的方法，其特征在于：步骤3所述计算三维粗网有限差分CMFD需要的净流修正系数的公式为： 

公式(27) 

公式(28) 

公式(29) 

式中： 

——高阶计算得到的粗网n的界面净流， 

——第g群的宏观移除截面， 

——网格间距。 

4.根据权利要求1所述的一种获取反应堆堆芯三维中子通量密度精细分布的方法，其特征在于：步骤3所述的全堆芯三维粗网有限差分CMFD计算公式为： 

公式(30) 

式中： 

是中子通量密度， 

是扩散系数， 

是中子裂变产额截面， 

是中子散射截面， 

x_g是第g群裂变谱， 

k_eff是特征值，也称为有效增殖因数， 

其中： 

。

Images