中子输运方程一维求解程序

开发该程序主要用于测试各种算法，算法的稳定性，算法的效率，以及学习并行算法的实施。目前程序角度上采用球谐函数方法，空间上采用有限元方法。

单能一维中子输运方程

$$μ\frac{∂\phi(z,μ)}{∂z}+ σ_t(z)ϕ(z,μ) = \frac{1}{2\pi}\int_{0}^{2\pi}\int_{-1}^{1}σ_s(z,\mu_0)\phi(z,\mu^{‘}))d\mu^{‘}d\varphi^{‘} + \frac{s(z)}{2}$$

推导过程

单能三维中子输运方程为
$$Ω\cdot ∇ϕ(r,Ω) + σ_t(r)ϕ(r,Ω) = \int_{Ω’}\sigma_s(r,Ω’\cdot Ω)\phi(r,Ω’)dΩ’ + s(r,Ω)$$
对于一维情况，根据对称性有
$$\frac{∂\phi(r,Ω)}{∂x} = 0$$

$$\frac{∂\phi(r,Ω)}{∂y} = 0$$

$$\phi(r,Ω) = \frac{1}{2\pi}\phi(z,\mu)$$

$$s(r,Ω) = \frac{1}{4\pi}s(z)$$

对三维的输运方程对 φ 在 $[0,2\pi]$ 积分，得
$$μ\frac{∂\phi(z,μ)}{∂z}+ σ_t(z)ϕ(z,μ) =\int_{0}^{2\pi}\int_{Ω’}\sigma_s(r,Ω’\cdot Ω)\phi(r,Ω’)dΩ’d\varphi + \frac{s(z)}{2}$$

定义
$$σ_s(z,\mu_0) = \int_{0}^{2\pi}σ_s(z,Ω’\cdot Ω)d\varphi$$

所以，一维中子输运方程为
$$μ\frac{∂\phi(z,μ)}{∂z}+ σ_t(z)ϕ(z,μ) = \frac{1}{2\pi}\int_{0}^{2\pi}\int_{-1}^{1}σ_s(z,\mu_0)\phi(z,\mu^{‘}))d\mu^{‘}d\varphi^{‘} + \frac{s(z)}{2}$$

边界条件

该方程是一阶双曲方程，需要给定如流边界条件，即
$$\phi(z,\mu) = g,n\cdot \mu <0$$
所以，真空边界条件为
$$\phi(z,\mu) = 0,n\cdot \mu <0$$
反射边界条件为
$$\phi(z,\mu) = \phi(z,-\mu),n\cdot \mu <0$$

Pn方法

用勒让得多项式展开角度变量

将中子角通量密度和散射截面都用勒让得多项式展开，得

$$\phi(z.μ) = \sum_{n=0}^{∞}\frac{2n+1}{2}\phi_n(z)P_n(\mu)$$

$$σ_s(z,\mu_0) = \sum_{n=0}^{∞}\frac{2n+1}{2}\sigma_{sn}(z)P_n(\mu_0)$$

由勒让得多项式加法定理及正交关系得散射源项为
$$\frac{1}{2\pi}\int_{0}^{2\pi}\int_{-1}^{1}σ_s(z,\mu_0)\phi(z,\mu^{‘}))d\mu^{‘}d\varphi^{‘} = \sum_{n=0}^{∞}\frac{2n+1}{2}\sigma_{sn}(z)\phi_n(z)P_n(\mu)$$

由勒让得多项式递推公式得

$$\sum_{n=0}^{∞}\frac{d\phi_n(z)}{dz}[(n+1)P_{n+1}(\mu)+n P_{n-1}(μ)]+ \sum_{n=0}^{∞}(2n+1)σ_t(z)\phi_n(z)P_n(\mu) = \sum_{n=0}^{∞}(2n+1)\sigma_{sn}(z)\phi_n(z)P_n(\mu) + s(z)$$

在方程两端同时乘以$P_{n}(\mu),n=0,1,2,…$,并对 $-1\le \mu \le 1$积分，得

$$\frac{n+1}{2n+1}\frac{d\phi_{n+1}(z)}{dz}+\frac{n}{2n+1}\frac{d\phi_{n-1}(z)}{dz}+\sigma_{n}\phi_n(z) = s(z)δ_{0n}$$

式中 $σ_n = σ_t - σ_{sn}$

在实际计算中，只取方程组前 $N+1$个方程，即认为$\frac{d\phi_{n+1}}{dz} = 0$以及$\phi_{n} = 0,n>N$，此时
$$\phi(z.μ) = \sum_{n=0}^{N}\frac{2n+1}{2}\phi_n(z)P_n(\mu)$$

将其写为矩阵形式为

$$\mathbf{A}\frac{d\vec{\phi}(z)}{dz}+\mathbf{\Sigma}\vec{\phi}(z) = \vec{s}$$

其中
$$\vec{\phi} = [\phi_0,\phi_1,\phi_2,…,\phi_N]^{T}$$
$$\vec{s} = [s_0,s_1,s_2,…,s_N]^{T}$$
$$\mathbf{\Sigma} = diag[\sigma_n]$$
$$\mathbf{A} = \begin{bmatrix} 0 & 1 & 0 &0 &…
\ \frac{1}{3}&0 &\frac{2}{3} & 0 & …
\ 0&\frac{2}{5} &0& \frac{3}{5} & …
\ …&… &…&… & …\end{bmatrix} $$

此时矩阵 $\mathbf{A}$ 不对称

用正交归一勒让得多项式展开

$$\phi(z.μ) = \sum_{n=0}^{∞}\phi_n^{‘}(z)P_n^{‘}(\mu)$$
$$σ_s(z,\mu_0) = \sum_{n=0}^{∞}\sigma_{sn}^{‘}(z)P_n^{‘}(\mu_0)$$
其中
$$P_n^{‘}(\mu) = \sqrt{\frac{2n+1}{2}}P_n(\mu)$$
$$\int_{-1}^{1}P_n^{‘}(\mu)P_m^{‘}(\mu)d\mu = δ_{nm}$$

递推公式
$$\mu P_n^{‘}(\mu) = \frac{n+1}{\sqrt{(2n+1)(2n+3)}}{P_{n+1}^{‘}(\mu)}+\frac{n}{\sqrt{(2n-1)(2n+1)}}{P_{n-1}^{‘}(\mu)}$$

将其带入输运方程，并利用勒让得多项式递推公式和加法定理得
$$\sum_{n=0}^{∞}\frac{d\phi_n^{‘}(z)}{dz}[\frac{n+1}{\sqrt{(2n+1)(2n+3)}}{P_{n+1}^{‘}(\mu)}+\frac{n}{\sqrt{(2n-1)(2n+1)}}{P_{n-1}^{‘}(\mu)}]+ \sum_{n=0}^{∞}σ_t(z)\phi_n^{‘}(z)P_n^{‘}(\mu) = \sum_{n=0}^{∞}\sigma_{sn}^{‘}(z)\phi_n^{‘}(z)P_n^{‘}(\mu) +\frac{s(z)}{2}$$

在方程两端同时乘以$P_{n}^{‘}(\mu),n=0,1,2,…$,并对 $-1\le \mu \le 1$积分，得
$$\mathbf{A}\frac{d\vec{\phi}(z)}{dz}+\mathbf{\Sigma}\vec{\phi}(z) = \vec{s}$$

其中
$$\vec{\phi} = [\phi_0,\phi_1,\phi_2,…,\phi_N]^{T}= \sqrt{2}[\phi_0^{‘},\phi_1^{‘},\phi_2^{‘},…,\phi_N^{‘}]^{T}$$
$$\vec{s} = [s_0,s_1,s_2,…,s_N]^{T}$$
$$\mathbf{\Sigma} = diag[\sigma_n]$$
$$\mathbf{A} = \begin{bmatrix} 0 & \frac{1}{\sqrt{3}} & 0 &0 &…
\ \frac{1}{\sqrt{3}}&0 &\frac{2}{\sqrt{15}} & 0 & …
\ 0&\frac{2}{\sqrt{15}} &0& \frac{3}{\sqrt{35}} & …
\ 0&0&\frac{3}{\sqrt{35}}&0 & …
\ …&…&…&… & …\end{bmatrix} $$

此时矩阵 $\mathbf{A}$ 是对称的，其中
$$a_{i+1,i} = a_{i,i+1} = \frac{i}{\sqrt{4i^2-1}}$$
其中0阶矩$\phi_0 = \sqrt{2}\phi_{0}^{‘} = \int_{-1}^{1}\phi(z,\mu)d\mu$ 表示标通量。外中子源是各向同性的，所以 $s_1,s_2,…,s_N$都为0。

有限元方法

对以下一阶方程组使用有限元方法

将$\vec{\phi}(z)$中的每一个通量矩 $\phi_{j}(z)$用有限元基函数展开得
$$\phi_{j}(z) = \sum_{i=1}^I\phi_{ji}N_{i}(z) = \vec{\phi}_{j}^{T}\vec{N}(z)$$

将离散的未知量写成如下矩阵形式
$$\mathbf{P} = \begin{bmatrix}
\phi_{11}& \phi_{12} & … & \phi_{1I}\\
\phi_{21}& \phi_{22} & … & \phi_{2I}\\
⋮&⋮ & ⋱ & ⋮\\
\phi_{J1}& \phi_{J2} & … & \phi_{JI}
\end{bmatrix}$$

所以
$$\vec{\phi}(z) = \mathbf{P}\vec{N}(z)$$

所以
$$\mathbf{A}\mathbf{P}\frac{d\vec{N}(z)}{dz}+\mathbf{\Sigma}\mathbf{P}\vec{N}(z) = \vec{s}$$

galerkin方法

在方程两端同时右乘以基函数$[\vec{N}]^{T}(z)$，并对单元区域进行积分，得
$$\mathbf{A}\mathbf{P}_e\int_{e}\frac{\partial \vec{N}(z)}{\partial z}\vec{N}^{T}(z)dz + \int_{e}\mathbf{Σ}_e\mathbf{P}_e\vec{N}(z)\vec{N}^{T}(z)dz = \int_{e}\vec{s}_e\vec{N}^{T}(z)dz$$

利用分部积分，得
$$n\cdot \mathbf{A}\mathbf{P}_{\Gamma_n}-\mathbf{A}\mathbf{P}_e\int_{e}\vec{N}(z)\frac{\partial \vec{N}^{T}(z)}{\partial z}dz + \int_{e}\mathbf{Σ}_e\mathbf{P}_e\vec{N}(z)\vec{N}^{T}(z)dz = \int_{e}\vec{s}_e\vec{N}^{T}(z)dz$$

边界条件

令$\mathbf{A}_n = n\cdot \mathbf{A}$，将其进行特征值分解得$\mathbf{A}_n = \mathbf{A}_n^{+}+\mathbf{A}_n^{-}$

• 真空边界
  $$\mathbf{A}_n^{-} \vec{\phi}_{\Gamma_{v}} = 0$$
• 反射边界
  $$\mathbf{A}_n^{-} \vec{\phi}_{\Gamma_{r}} = \mathbf{A}_n^{-} \mathbf{R}\vec{\phi}_{\Gamma_{r}}$$
  其中矩阵 $\mathbf{R} = \int_{-1}^{1} \vec{P}_n(\mu)\vec{P}_n(-\mu)d\mu$，该矩阵是一个对角矩阵
  $$\mathbf{R} = diag[1,-1,1,-1,…,(-1)^{n}]$$

向量化

将以上矩阵方程变为向量方程为

$$-\mathbf{B}_{1z}^{T}⊗\mathbf{A}\vec{\phi} + \mathbf{B}_{11}⊗\mathbf{Σ}\vec{\phi}+\mathbf{A}_n^{+}\vec{\phi}_{Γ}+\mathbf{A}_n^{-} \mathbf{R}\vec{\phi}_{\Gamma_{r}} = \vec{b} ⊗\vec{s}$$

其中
$$\mathbf{B}_{1z} = \int_{e}\vec{N}(z)\frac{\partial \vec{N}^{T}(z)}{\partial z}dz$$
$$\mathbf{B}_{11} = \int_{e}\vec{N}(z)\vec{N}(z)^{T}dz$$
$$\vec{b} = \int_{e}\vec{N}(z)dz$$

一维线性单元

形状函数及其导数
$$\vec{N}(z) = [λ_1,λ_2]^{T}$$
$$\frac{\partial \vec{N}(z)}{\partial z} = \frac{[b_1,b_2]^{T}}{l}$$

其中
$$b_1 = -1,b_2 = 1$$

由欧拉公式

$$\int_{z}λ_1^aλ_2^bdz=\frac{a!b!}{(a+b+1)!}l$$

得
$$\int_{z}dz=l$$
$$\int_{z}λ_idz=\frac{l}{2}$$
$$\int_{z}λ_i^{2}dz=\frac{l}{3}$$
$$\int_{z}λ_iλ_jdz=\frac{l}{6}$$

所以
$$\mathbf{B}_{1z} = \int_{e}\vec{N}(z)\frac{\partial \vec{N}^{T}(z)}{\partial z}dz = \frac{1}{2}\begin{bmatrix}
-1&1\\
-1&1
\end{bmatrix}$$

$$\mathbf{B}_{11} = \int_{e}\vec{N}(z)\vec{N}(z)^{T}dz = \frac{l}{6}\begin{bmatrix}
2&1\\
1&2
\end{bmatrix}$$

$$\vec{b} = \int_{e}\vec{N}(z)dz = \frac{l}{2}\begin{bmatrix}
1\\
1
\end{bmatrix}$$

方程组装

单元矩阵

内部

$$\mathbf{K}_e = -\mathbf{B}_{1z}^{T}⊗\mathbf{A} + \mathbf{B}_{11}⊗\mathbf{Σ}$$

边界

$$\mathbf{K}_{be} = \mathbf{A}_n^{+}+\mathbf{A}_n^{-} \mathbf{R}$$

RHS

$$\vec{f}_e = \vec{b} ⊗\vec{s}$$

附录

勒让得多项式

定义式

$$P_n(\mu) = \frac{1}{2^{n}n!}\frac{d^{n}(\mu^2-1)^{n}}{d\mu^{n}}$$

正交关系式

$$\int_{-1}^{1}P_n(\mu)P_m(\mu)d\mu = \frac{2}{2n+1}δ_{nm}$$

递推关系式

$$\mu P_n(\mu) = \frac{1}{2n+1}[nP_{n-1}(\mu) + (n+1)P_{n+1}(\mu)]$$

加法定理

$$P_n(\mu_0) = P_n(\mu)P_n(\mu^{‘}) + 2\sum_{m=1}^{n}\frac{(n-m)!}{(n+m)!}P_{n}^{m}(\mu)cos(m(φ-\varphi^{‘}))$$

多项式系数

	1	$x$	$x^2$	$x^3$	$x^4$	$x^5$	$x^6$	$x^7$
1	1
2	0	1
3	$-\frac{1}{2}$		$\frac{3}{2}$
4		$-\frac{3}{2}$		$\frac{5}{2}$
5	$\frac{3}{8}$		$-\frac{30}{8}$		$\frac{35}{8}$
6		$\frac{15}{8}$		$-\frac{70}{8}$		$\frac{63}{8}$
7	$-\frac{5}{16}$		$\frac{105}{16}$		$-\frac{315}{16}$		$\frac{231}{16}$
8		$-\frac{35}{16}$		$\frac{315}{16}$		$-\frac{693}{16}$		$\frac{429}{16}$

图像

向量化

向量化算子

定义向量化算子，按列重排$V_c$，按行重排 $V_r$
$V_c(\mathbf{P}) = \begin{bmatrix}
\phi_{11}\ \phi_{21}\\⋮\ \phi_{J1}\ \phi_{12} \ ⋮\\\phi_{J2} \ ⋮\\\phi_{1I} \ ⋮\\\phi_{JI}
\end{bmatrix}$ $V_r(\mathbf{P}) = \begin{bmatrix}
\phi_{11}\ \phi_{12}\\⋮\ \phi_{1I}\ \phi_{21} \ ⋮\\\phi_{2I} \ ⋮\\\phi_{J1} \ ⋮\ \phi_{JI}
\end{bmatrix}$

矩阵方程向量化

$\mathbf{A}\mathbf{X}\mathbf{B} = \mathbf{Y} \Leftrightarrow (\mathbf{A} ⊗ \mathbf{B}^{T})V_r(\mathbf{X}) = V_r(\mathbf{Y})\Leftrightarrow (\mathbf{B}^{T}⊗\mathbf{A} )V_c(\mathbf{X}) = V_c(\mathbf{Y})$