材料和求积组文件：material_and_quadrature

• 一维多群MOC输运方程为：

  $$ \mu \frac{\mathrm{d} \Phi_g(x, \mu)}{\mathrm{d} x}+ \Sigma_{g,t}(x)\Phi_g(x, \mu)=S_g(x, \mu),\quad x \in(0, X), \mu \in(-1,1) $$

  式子中：

  • $x$ ：沿x轴的位置， $(0, X)$为堆芯边界坐标
  • $\mu = \cos \theta$：中子飞行方向与 $x$ 轴正方向夹角的余弦。
  • $\Phi_g(x, \mu)$：表示能群 $g$ 的角通量。
  • $\Sigma_{g,t}(x)$：能群 $g$的总宏观截面。
  • $S_g(x, \mu)$：总源项。

• 将其应用于一维板型堆（即仅x方向有变化，其余方向可看作无限均匀介质）

  

  • 以如下所示的建模及网格划分为例：

  

  """ 上图对应的建模参数如下 """
  # 设置材料区域的边界和对应的材料名称以及边界条件
  region_widths = np.array([5.0, 4.0, 5.0, 4.0])
  material_names = ["c5g7_UO2_7G", "c5g7_MOX43_7G", "c5g7_MOX7_7G", "c5g7_MOX87_7G"]
  boundary_conditions = [1, 1]
  
  mesh_size = 1 # 期望的网格尺寸
  num_polars = 4  # 离散的角度个数
  num_groups = 2 # 离散的能群数量
  

  • 材料分配及网格划分initialize_grids函数的代码如下

• 对于一维多群 MOC 输运方程，其总源项 $S_g(x, \mu)$ 由散射项，裂变项和固定源组成:

  1. 散射源:

     来自其他能群中子的散射。这个部分通常是通过对所有能群和所有方向的积分来计算的，表示为:

     $$ S_{\text {scatter }, g}(x, \mu)=\frac{1}{2} \sum_{g^{\prime}} \int_{-1}^1 \Sigma_{s, g^{\prime} \rightarrow g}\left(x, \mu^{\prime} \rightarrow \mu\right) \Phi_{g^{\prime}}\left(x, \mu^{\prime}\right) d \mu^{\prime} $$

     其中：

     • $\Sigma_{s, g^{\prime} \rightarrow g}\left(x, \mu^{\prime} \rightarrow \mu\right)$是从能群 $g^{\prime}$ 到能群 $g$ 的方向依赖散射截面
     • $\Phi_{g^{\prime}}\left(x, \mu^{\prime}\right)$ 是在位置 $x$ 和方向 $\mu^{\prime}$的中子角通量。

  2. 裂变源

     如果系统内部包含可裂变材料，那么裂变事件也会产生中子，这些中子形成了源项的一部分。裂变源通常计算为：

     $$

     S_{\text{fission}, g}(x, \mu) = \frac{1}{2} \sum_{g'} \frac{\chi_{g^{\prime} \rightarrow g}(\mu)}{k_{\text{eff}}}\nu_{g^{\prime}} \Sigma_{f, g'}(x) \int_{-1}^1 \Phi_{g'}(x, \mu') \, d\mu'

     $$

     其中：

     • $\chi_{g^{\prime} \rightarrow g}(\mu)$是裂变中子能谱（表示裂变中子在方向 $\mu$，从能群 $g^{\prime}$ 到能群 $g$的转移概率）
     • $\nu_{g^{\prime}} \Sigma_{f, g'}(x)$ 是能群 $g'$的裂变产额截面
     • $\Phi_{g'}(x, \mu')$是相应的中子通量

  3. 外部源:

     这是由外界施加的中子源，可能是由实验或其他中子产生机制提供的，与物质内部的核反应无关。通常表示为 $Q_g(x, \mu)$。

• 综合以上三部分，源项 $S_g(x, \mu)$可以表示为:

  $$ S_g(x, \mu)=S_{\text {scatter }, g}(x, \mu)+S_{\text {fission }, g}(x, \mu)+Q_g(x, \mu) $$

  • 为了简化，我们在这里忽略掉不同方向之间的转化 $\mu^{\prime} \rightarrow \mu$, 讨论平几何带各向同性广义源 $S_g(x)$ 的输运方程

    $$ \mu \frac{\mathrm{d} \Phi_g(x, \mu)}{\mathrm{d} x}+ \Sigma_{g,t}(x)\Phi_g(x, \mu)=S_g(x),\quad x \in(0, X), \mu \in(-1,1) $$

  • 此时，广义源 $S_g(x)$可表示为:

    $$ \begin{split} S_g(x)&= S_{\text {fission }, g}(x)+S_{\text {scatter }, g}(x)+Q_g(x) \\ &= \frac{1}{2}\sum_{g'} \frac{\chi_{g^{\prime} \rightarrow g}}{k_{\text{eff}}}\nu_{g^{\prime}} \Sigma_{f, g'}(x) \int_{-1}^1 \Phi_{g'}(x, \mu') \, d\mu'+\frac{1}{2}\sum_{g^{\prime}} \Sigma_{s, g^{\prime} \rightarrow g}\left(x\right) \int_{-1}^1\Phi_{g^{\prime}}\left(x, \mu^{\prime}\right) d \mu^{\prime} +Q_g(x) \\ &= \frac{1}{2} \sum_{g'} \left( \frac{\chi_{g^{\prime} \rightarrow g}}{k_{\text{eff}}}\nu_{g^{\prime}} \Sigma_{f, g'}(x) + \Sigma_{s, g^{\prime} \rightarrow g}\left(x\right) \right)\int_{-1}^1 \Phi_{g'}(x, \mu') \, d\mu'+ Q_g(x)\end{split} $$

  • 源中的标通量积分可以通过角度离散方法求得，即利用特别选定的一组方向余弦值 $\left\{\mu_p\right\}$采用高斯求积公式近似求解积分。假定在区间 $(-1,1)$内已选定了一组数量为 $N$ 的方向余弦值 $\left\{\mu_p\right\}$ 及其相应的求积权重 $\left\{w_p\right\}, p=1,2, \cdots, N$。则标通量 $\phi_g$ 可近似表示为：

    $$ \phi_g(x) = \int_{-1}^1 \Phi_g\left(x, \mu^{\prime}\right) d \mu^{\prime} \simeq \sum_{p=1}^N w_p \Phi_g \left(x, \mu_p\right)

    $$

  • 而后，我们采用平源近似，即假设平源区网格 $i$ 内各源项为常数，同时将裂变中子能谱 $\chi_{g^{\prime} \rightarrow g}$简化为产生的裂变中子在每个能量群 $g$ 中的平均概率的 $\chi_g$ , 此时，网格 $i$ 中广义源 $S_g(i)$可被简化为:

    $$ S_g(i)=\frac{1}{2}\left(\frac{\chi_{g}}{k_{\text{eff}}} \nu_g \Sigma_{f, g}(i) +\sum_{g'} \Sigma_{s, g^{\prime} \rightarrow g}\left(i \right)\right)\phi_{g}(i)+ Q_g(i)

    $$

• 以角度$\mu = \cos \theta$，网格 $i$ 的特征线段为例，输运扫描过程中，入射与出射角通量如下所示：

  

  • $i$ 为网格编号
  • $x_i$ 为网格宽度。
  • $L_(μ,i)$为角度μ的特征线段在该网格上的长度，$L_(μ,i)=x_i/(cos⁡θ)$。
  • $phi\_left$ 为从左射入该网格的入射角通量，$phi\_left^{\prime}$ 为从左射出该网格的出射角通量，也即下一个网格 $i+1$ 在 $μ$ 方向的入射角通量。

• 沿该特征线段，角通量满足式：

  $$ \Phi_{\mu,g,i}(s) =\Phi_{\mu,g,i}^{in} \times  e^{-s\Sigma_{g,t}}+\frac{Q_{g,i} }{\Sigma_{g,t}}\left(1-e^{-s\Sigma_{g,t}}\right),\quad  s \in(0, L_{\mu,i}) $$

  其中：

  • $s$ 为沿着 $μ$ 方向前进，距离网格左边界的距离
  • $Q_(g,i)$ 网格 $i$ 能群 $g$ 的总源项，也即前文的 $S_g (i)$ 。

• 故而，沿着该特征线段角通量的变化量del_phi为:

  $$ \begin{split}\Delta\Phi_{\mu,g,i} &= \Phi_{\mu,g,i}^{out} - \Phi_{\mu,g,i}^{in} \\&=\Phi_{\mu,g,i}(L_i) - \Phi_{\mu,g,i}^{in}\\&=\left( \frac{Q_{g,i} }{\Sigma_{g,t}} - \Phi_{\mu,g,i}^{in} \right)\left(1-e^{-s\Sigma_{g,t}}\right)\end{split} $$

• 故而根据平源近似，我们可以求得该特征线段平均角通量phi_avg：

  $$ \bar{\Phi}{\mu,g,i}=\frac{Q{g,i}}{\Sigma_{g,t}}+\frac{\Delta\Phi}{\Sigma_{g,t} L_{\mu,i}} $$

• 沿着该方向不断重复这一过程直到边界，并在扫描结束后，根据边界条件存储下一次扫描的起始边界角通量，就完成了一条特征线的扫描。

• 具体的，根据对称关系，在反射边界上，下一次扫描时编号 $-μ,g,i$ 的入射角通量为此次扫描编号 $μ,g,i$ 的出射角通量；而真空边界条件下，出射角通量不会返回，故而下一次扫描时真空边界上的起始角通量均为0：

  $$ \begin{aligned} \left.\Phi_{-\mu,g,i}(s)\right|{s=refl} &= \left.\Phi{\mu,g,i}(s)\right|{s=refl} \\ \left.\Phi{-\mu,g,i}(s)\right|_{s=vacu} &= 0\end{aligned} $$

• 求得所有角通量后便可根据高斯求积公式得到标通量信息，进而根据源项计算公式求得此次扫描堆芯内的总中子数目，除以与上一次扫描对应的中子数目，便得到了此次迭代的有效增殖因子$k_{eff}$。不断重复此过程直到$k_{eff}$不再变化，我们便完成了稳态下的堆芯输运计算。

• 根据以上推导，我们可以完成一维MOC输运计算程序的开发。我们将程序分为三大模块，建模及网格划分模块、单条特征线输运扫描模块、有效增殖因子求解模块，分别对应前文中给出的函数initialize_grids函数，还有sweep_1d_single_group函数及eigen_solver函数。

  • 特征线输运扫描sweep_1d_single_group函数的代码如下
  • 源迭代eigen_solver函数代码如下

  函数计算所需的不会变的参数，如边界条件、每个网格的长度及对应的材料等信息，我们将其存储为全局变量。而计算过程中会变化的参数，如扫描中源项、边界角通量、能群编号，则作为函数的输入参数。

• 我们使用无限大均匀介质基准题及一维C5G7-MOX基准题验证程序开发的准确性。其中无限大均匀介质中子输运问题为两群问题，可以通过给有限空间赋予全反射边界条件进行计算。其增殖因子可以通过双群截面解析求解，等于：

  $$ k_{e f f}=\frac{v \Sigma_{f ,1}\left(\Sigma_{t ,2}-\Sigma_{s, 2\to2}\right)+v \Sigma_{f ,2} \Sigma_{s, 1\to2}}{\left(\Sigma_{t ,2}-\Sigma_{s ,2\to2}\right)\left(\Sigma_{t ,1}-\Sigma_{s, 1\to1}\right)-\Sigma_{s ,2\to1} \Sigma_{s, 1\to 2}}=1.17952 $$

  • 无限大均匀介质题代码

• 一维C5G7-MOX基准题几何布置、尺寸及边界条件如下图所示。该问题的参考解为1.02157(来自程序PEACH-1D) 和1.02156±0∙00005(来自程序MCMG)。

  

  • 一维C5G7-MOX基准题代码
  • 功率密度绘图代码