【核反应堆物理分析】3-中子扩散理论： 一些扩散方程

⌛ 中子扩散理论

💿 1 单能中子扩散方程

🌝 1.1 中子输运的几何描述

初始具有某能量和位置和某运动方向的中子，在稍晚时运动到堆内的另一位置以另一能量和运动方向出现，这种过程称为中子输运。

中子运动的数学描述：

• 位置矢量：$\overrightarrow{r}(x,y,z)$；

• 能量：$E$，也可用运动速度$v$表示。

• 运动方向：$\overrightarrow{\mathrm{\Omega }}$，其为立体角，由中子相对坐标原点的极角$\theta$和方位角$\phi$表示。

  $$\text{d}\overrightarrow{\mathrm{\Omega }}=\frac{\text{d}S}{r^2}=\frac{r^2\sin \theta \text{d}\theta \text{d}\phi }{r^2}=\sin \theta \text{d}\theta \text{d}\phi$$

• 中子角密度：$n(\overrightarrow{r},E,\overrightarrow{\mathrm{\Omega }})$。其为在$\overrightarrow{r}$处单位体积内和能量为$E$的单位能量间隔内，运动方向为$\overrightarrow{\mathrm{\Omega }}$的单位立即角内的中子数目。

• 中子角通量密度：$\varphi (\overrightarrow{r},E,\overrightarrow{\mathrm{\Omega }})=n(\overrightarrow{r},E,\overrightarrow{\mathrm{\Omega }})v(E)$。沿$\overrightarrow{\mathrm{\Omega }}$方向在单位之间内穿过垂直于这个方向的单位面积上的中子数目。

• 中子密度（标量）：

  $$n(\overrightarrow{r},E)={\int }_{4\pi }n(\overrightarrow{r},E,\overrightarrow{\mathrm{\Omega }})\text{d}\overrightarrow{\mathrm{\Omega }}$$

• 中子通量密度（标量）：

  $$\varphi (\overrightarrow{r},E)={\int }_{4\pi }\varphi (\overrightarrow{r},E,\overrightarrow{\mathrm{\Omega }})\text{d}\overrightarrow{\mathrm{\Omega }}$$

描述中子运输过程的精确方程称为玻尔兹曼输运方程，经过近似简化得到的方程称为扩散方程。

🌾 1.2 斐克定律

中子扩散是中子与介质原子核散射碰撞的结果，中子总是从中子密度高的地方向底的地方扩散。

斐克定律的基本假设：

1. 介质是无限的、均匀的；
2. 在实验室L系中，散射是各向同性的；
3. 介质的吸收截面很小，${\mathrm{\Sigma }}_a\ll {\mathrm{\Sigma }}_s$；
4. 中子通量密度是随空间位置缓慢变化的函数。（中子通量密度函数随空间分布连续可导）

$\text{d}A$为一个垂直于$z$轴的平面，研究从某${\overrightarrow{r}}^{\prime }$附近体积$\text{d}V$散射到$\text{d}A$的中子数，以计算整个$\text{d}A$上的中子数$J_z^{+}$。每秒从$\text{d}V$内散射出来沿着$\overrightarrow{\mathrm{\Omega }}$方向未经碰撞到达$\text{d}A$上的中子数为：

$$\frac{1}{4\pi }{\mathrm{\Sigma }}_s\varphi ({\overrightarrow{r}}^{\prime })e^{-{\mathrm{\Sigma }}_a|l|}\cos \theta \text{d}A\text{d}l$$

沿$\overrightarrow{\mathrm{\Omega }}$每秒穿过$\text{d}A$的中子数为沿$l:-\mathrm{\infty }\to 0$的积分：

$$\frac{\text{d}A}{4\pi }{\int }_{-\mathrm{\infty }}^0{\mathrm{\Sigma }}_s\varphi ({\overrightarrow{r}}^{\prime })e^{-{\mathrm{\Sigma }}_s|l|}\cos \theta \text{d}l$$

用Taylor展开来替代$\varphi ({\overrightarrow{r}}^{\prime })$项，其中一阶导数项为：

$$\frac{\text{d}\varphi }{\text{d}l}=\mathrm{\Omega }\mathrm{\nabla }\varphi =({\mathrm{\Omega }}_x,{\mathrm{\Omega }}_y,{\mathrm{\Omega }}_z)\cdot (\frac{\partial \varphi }{\partial x},\frac{\partial \varphi }{\partial y},\frac{\partial \varphi }{\partial z})$$

其中向量$({\mathrm{\Omega }}_x,{\mathrm{\Omega }}_y,{\mathrm{\Omega }}_z)$是向量$\overrightarrow{\mathrm{\Omega }}$在三个轴上的投影：

$$\{\begin{array}{c}{\mathrm{\Omega }}_x=\sin \theta \cos \phi \\ {\mathrm{\Omega }}_y=\sin \theta \sin \phi \\ {\mathrm{\Omega }}_z=\cos \theta \end{array}$$

带入并计算得到：

$$J_z^{+}(\overrightarrow{r})=\frac{\varphi (\overrightarrow{r})}{4}-\frac{1}{6{\mathrm{\Sigma }}_s}\frac{\partial \varphi (\overrightarrow{r})}{\partial z}$$

同样在反方向上：

$$J_z^{-}(\overrightarrow{r})=\frac{\varphi (\overrightarrow{r})}{4}+\frac{1}{6{\mathrm{\Sigma }}_s}\frac{\partial \varphi (\overrightarrow{r})}{\partial z}$$

沿着$z$方向穿过$\text{d}A$平面上单位面积净中子数为：

$$J_z(\overrightarrow{r})=J_z^{+}(\overrightarrow{r})-J_z^{-}(\overrightarrow{r})=-\frac{{\lambda }_s}{3}\frac{\partial \varphi (\overrightarrow{r})}{\partial z}$$

同样在$x$轴和$y$轴方向，可以同样得到：

$$\begin{gathered} J_x(\overrightarrow{r})=-\frac{{\lambda }_s}{3}\frac{\partial \varphi (\overrightarrow{r})}{\partial x} \\ {J}_y(\overrightarrow{r})=-\frac{{\lambda }_s}{3}\frac{\partial \varphi (\overrightarrow{r})}{\partial y} \end{gathered}$$

中子流密度：单位时间穿过单位面积的净中子流动情况，矢量，是具有不同方向微分中子束矢量合成的量。

$$\begin{gathered} J_n=\overrightarrow{J}\cdot \overrightarrow{n} \\ \overrightarrow{n}=\cos \alpha \overrightarrow{i}+\cos \beta \overrightarrow{j}+\cos \gamma \overrightarrow{k} \\ \overrightarrow{J}=-\frac{{\lambda }_s}{3}\text{grad}\varphi \end{gathered}$$

中子流密度与中子通量密度的关系：

$$\overrightarrow{J}(\overrightarrow{r})={\int }_{4\pi }\overrightarrow{\mathrm{\Omega }}\varphi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }})\text{d}\overrightarrow{\mathrm{\Omega }}$$

当净中子流方向与法向方向一致时，$J_n>0$，$\overrightarrow{J}$取得最大值。

斐克定律：中子流密度正比于负的中子通量密度梯度：

$$\begin{gathered} \overrightarrow{J}=-D\text{grad}\varphi \\ D=\frac{{\lambda }_s}{3} \end{gathered}$$

修正的斐克定律：为了修正散射的各项异性的影响，用运输平均自由程替代散射平均自由程：

$$\begin{gathered} D=\frac{{\lambda }_{\text{tr}}}{3} \\ {\lambda }_{\text{tr}}=\frac{{\lambda }_s}{1-{\overline{\mu }}_0} \end{gathered}$$

其中${\overline{\mu }}_0$是L系中平均散射角余弦，${\mu }_0={\displaystyle \frac{2}{3A}}$。重核由于${\overline{\mu }}_0\ll 1$，${\lambda }_{\text{tr}}\approx {\lambda }_s$。

斐克定律的物理意义：中子穿过某个平面左边每秒每单位体积内发生散射碰撞的中子数比右边发生散射碰撞的中子数多，左边散射碰撞到右边的中子比右边散射碰撞到左边的多。这样就在该平面上产生了一个沿$x$方向的净中子流，平面两侧中子通量密度梯度越大，中子流就越大。

🌔 1.3 单能中子扩散方程（中子数守恒-质量守恒）

在$V$内中子数守恒的方程可写成：

$$\frac{\text{d}}{\text{d}t}{\int }_{V}n(\overrightarrow{r},t)\text{d}V=\text{Source rate}-\text{Leak rate}-\text{Absorption rate}$$

泄漏率 Leak rate：

$$L={\int }_S\overrightarrow{J}(\overrightarrow{r},t)\cdot \overrightarrow{n}\text{d}S={\int }_V\text{div}\overrightarrow{J}(\overrightarrow{r},t)\text{d}V$$

产生率 Source rate：

$$S={\int }_{V}S(\overrightarrow{r},t)\text{d}V$$

吸收率 Absorption rate：

$$A={\int }_V{\mathrm{\Sigma }}_a\varphi (\overrightarrow{r},t)\text{d}V$$

故有中子数连续方程：

$$\frac{\partial n(\overrightarrow{r},t)}{\partial t}=S(\overrightarrow{r},t)-{\mathrm{\Sigma }}_a\varphi (\overrightarrow{r},t)-\text{div}\overrightarrow{J}(\overrightarrow{r},t)$$

当斐克定律成立时，有：

$$\frac{1}{v}\frac{\partial \varphi (\overrightarrow{r},t)}{\partial t}=S(\overrightarrow{r},t)+D{\mathrm{\nabla }}^2\varphi (\overrightarrow{r},t)-{\mathrm{\Sigma }}_a\varphi (\overrightarrow{r},t)$$

上式称为单能中子扩散方程，可以用其近似确定许多情况下的中子通量密度分布。

注：散度与算子

$$\begin{gathered} div\overrightarrow{J}=\mathrm{\nabla }\cdot \overrightarrow{J}=\mathrm{\nabla }\cdot (-Dgard\varphi )=-D\mathrm{\nabla }\cdot \mathrm{\nabla }\varphi =-D{\mathrm{\nabla }}^2\varphi =-D\mathrm{\Delta }\varphi \\ \mathrm{\Delta }x(x\mathrm{的}\mathrm{增}\mathrm{量}) \\ \mathrm{\Delta }x(x\mathrm{的}Laplace\mathrm{算}\mathrm{子}\mathrm{运}\mathrm{算}) \end{gathered}$$

当中子通量密度不随着时间变化，则有稳态单能中子扩散方程：

$$D{\mathrm{\nabla }}^2\varphi (\overrightarrow{r})-{\mathrm{\Sigma }}_a\varphi (\overrightarrow{r})+S(\overrightarrow{r})=0$$

🐐 1.4 扩散方程边界条件

求解扩散方程中常用的集中边界条件：

1. 在扩散方程适用范围内，中子通量密度数值必须是正的、有限的实数。

2. 在两种不同扩散性质的介质交界面上，垂直于分界面的中子流密度相等，中子通量密度相等。

   $${\varphi }_A={\varphi }_B$$

3. 介质与真空交界的外表面上，自真空返回介质的中子流等于零：

   $$J_x^{-}{|}_{x=0}=0$$

   由于中子在空气中的平均自由程比在非气体介质中的平均自由程大得多，所以中子不可能自真空中散射会介质中来。

4. 从交界面处将中子通量密度的分布曲线，按照交界面上的斜率向真空向外直线外推，则在交界面某个位置$d$处中子通量密度为0，有：

   $$d=0.7104{\lambda }_{\text{tr}}$$

   $d$称为直线外推距离，其位置中子通量密度为零。
   
   

🐝 1.5 斐克定律和扩散理论的适用范围

斐克定律的范围：

1. 假定扩散介质是无限的。有限介质内距离表面几个自由程之外的内部区域，斐克定律是成立的，在距离真空边界三个自由程内的区域是不适用的。
2. 推导把Taylor项展开到一阶项，要求中子通量密度必须缓慢变化，或变化梯度不大。所以斐克定律在强吸收体或两种扩散性质显著不同的交界面附近的几个自由程内不适用。
3. 斐克定律只适用于吸收截面远小于散射截面的区域。
4. 由于假设不考虑源的贡献，所以在强中子源三个平均自由程的区域内，斐克定律不适用。
   
   

🐧 2 中子扩散方程的解

🐦 2.1 稳态扩散方程的通解形式

稳态情况下的扩散方程为：

$$D{\mathrm{\nabla }}^2\varphi (\overrightarrow{r})-{\mathrm{\Sigma }}_a\varphi (\overrightarrow{r})+S(\overrightarrow{r})=0$$

无外中子源情况下，扩散方程为：

$${\mathrm{\nabla }}^2\varphi (\overrightarrow{r})-\frac{\varphi (\overrightarrow{r})}{L^2}=0$$

其中：

$$L^2=\frac{D}{{\mathrm{\Sigma }}_a}$$

$L$称为中子扩散长度，因为$L^2$有长度平方的量纲。

上式(27)也称为波动方程或Helmholtz方程。集中特殊情况的Helmholtz方程的解如下表所示：

$${\mathrm{\nabla }}^2\varphi (\overrightarrow{r})\pm B^2\varphi (\overrightarrow{r})=0$$

	${\mathrm{\nabla }}^2$	$+B^2$	$-B^2$
一维平板	${\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}$	$A\sin Bx+C\cos Bx$	$A{e}^{-Bx}+C{e}^{-Bx}$或$A\sinh Bx+C\cosh Bx$
球	${\displaystyle \frac{{\text{d}}^2}{\text{d}{r}^2}+\frac{2}{r}{\displaystyle \frac{\text{d}}{\text{d}r}}}$	$A{\displaystyle \frac{\sin Br}{r}+C{\displaystyle \frac{\cos Br}{r}}}$	$A{\displaystyle \frac{e^{-Br}}{r}+C\frac{e^{-Br}}{r}}$或$A{\displaystyle \frac{\sinh Br}{r}+C{\displaystyle \frac{\cosh Br}{r}}}$
一维圆柱	${\displaystyle \frac{{\text{d}}^2}{\text{d}{r}^2}+\frac{1}{r}{\displaystyle \frac{\text{d}}{\text{d}r}}}$	$A{J}_0(Br)+C{Y}_0(Br)$	$A{I}_0(Br)+C{K}_0(Br)$

🌍 2.2 无限介质点源情况

无限介质有一各向同性点源S，用球坐标系求解，坐标原点在点源位置上：

$$\frac{{\text{d}}^2\varphi (r)}{\text{d}{r}^2}+\frac{2}{r}\frac{\text{d}\varphi (r)}{\text{d}r}-\frac{\varphi (r)}{L^2}=0,r>0$$

定解条件：

1. 有界性条件：除原点，中子通量密度在各处均为有限值；
2. 中子源条件：${\displaystyle \underset{r\to 0}{lim}4\pi r^{2}J(r)=S}$

其解为：

$$\varphi (r)=A\frac{e^{-r/L}}{r}+C\frac{e^{r/L}}{r}$$

常数$A$用中子源条件求出：

$$\begin{gathered} J(r)=-D\frac{\text{d}\varphi (r)}{\text{d}r} \\ \underset{r\to 0}{lim}4\pi r^{2}J(r)=S \end{gathered}$$

解得：

$$A=\frac{S}{4\pi D}$$

由有界性边界条件：$C=0$。

故中子通量密度为：

$$\varphi (r)=\frac{S{e}^{-r/L}}{4\pi Dr},r>0$$

🔎 2.3 无限平面源位于有限厚度介质内的情况

厚度为$a$（包括外推距离）的无限均匀平板的中心面上有一源强为$S$的平面源，有：

$$\frac{{\text{d}}^2\varphi (x)}{\text{d}{x}^2}-\frac{\varphi (x)}{L^2}=0,x\ne 0$$

定解条件：

1. 边界条件：当$x=\pm (a/2)$时，$\varphi (\pm a/2)=0$；
2. 中子源条件：${\displaystyle \underset{x\to 0}{lim}J(x)=S/2}$。

其普遍解为：

$$\varphi (x)=A{e}^{-x/L}+C{e}^{x/L}$$

由边界条件得：

$$C=-A{e}^{-a/L}$$

由中子源条件可以求出：

$$A=\frac{SL}{2D}(1+e^{-a/L}{)}^{-1}$$

其中子通量密度为：

$$\begin{gathered} \varphi (x)=\frac{SL}{2D}\frac{e^{-|x|/L}-e^{-(a-|x|)/L}}{1+e^{-a/L}} \\ =\frac{SL}{2D}\frac{\sinh [(a-2|x|)/2L]}{\cosh (a/2L)} \end{gathered}$$

对于无限介质平面源情况，令$a\to \mathrm{\infty }$：

$$\varphi (x)=\frac{SL}{2D}{e}^{-|x|/L}$$

结论1：在无限介质中没有中子泄露损失，但对有限厚平板，中子会不断从边界面泄露出去。对于较薄的平板，在边界处中子通量密度下降的会更快一些。

结论2：介质厚度等于中子扩散长度3倍时，除在边界面附近外，中子通量密度的分布与无限介质内的分布相差不多。对于单能情况，反射层厚度大于三个扩散长度时，其效果和无限厚大致相当了。

📟 2.4 包含两种不同介质的情况

一厚度为$a$，长、宽均为无限的平面中子源，源强为$S$，在平板两侧为一无限厚的另一种介质，有扩散方程：

$$\begin{gathered} \frac{{\text{d}}^2{\varphi }_1(x)}{\text{d}{x}^2}-\frac{1}{L_1^2}{\varphi }_1(x)=0,|x|<\frac{a}{2},x\ne 0 \\ \frac{{\text{d}}^2{\varphi }_2(x)}{\text{d}{x}^2}-\frac{1}{L_2^2}{\varphi }_2(x)=0,|x|>\frac{a}{2} \end{gathered}$$

定解条件：

1. 无穷边界条件：当$|x|\to \mathrm{\infty }$时，${\varphi }_2(x)$趋近于零。
2. 中子源条件：${\displaystyle \underset{x\to 0}{lim}J(x)=S/2}$。
3. 交界面第一类边界条件：${\varphi }_1(\pm a/2)={\varphi }_2(\pm a/2)$。（1函数值，2函数导数值，3导数+函数=规律）
4. 交界面第二类边界条件：${D_1{\displaystyle \frac{\text{d}{\varphi }_1}{\text{d}x}}|}_{x=\pm a/2}={D_2{\displaystyle \frac{\text{d}{\varphi }_2}{\text{d}x}}|}_{x=\pm a/2}$

其普遍解为：

$$\begin{gathered} {\varphi }_1=A_1\cosh (x/L_1)+C_1\sinh (x/L_1) \\ {\varphi }_2=A_2{e}^{-x/L_2}+C_2{e}^{x/L_2} \end{gathered}$$

由无穷边界条件得：$C_2=0$；

中子源条件可以求得：${\displaystyle C_1=-\frac{S{L}_1}{2D_1}}$；

经过交界面边界条件可以求出：

$$\begin{gathered} A_1=\frac{S{L}_1}{2D_1}\frac{D_1{L}_2\cosh (a/2L_1)+D_2{L}_1\sinh (a/2L_1)}{D_2{L}_1\cosh (a/2L_1)+D_1{L}_2\sinh (a/s{L}_1)} \\ {A}_2=\frac{S{L}_1{L}_2}{2}\frac{\exp \left(\frac{a}{2L_2}\right)}{[D_2{L}_1\cosh \left(\frac{a}{2L_1}\right)+D_1{L}_2\sinh \left(\frac{a}{2L_1}\right)]} \end{gathered}$$

🔓 3 *反照率

反照率：反射层的效率：

$$\beta =\frac{J^{-}}{J^{+}}$$

在两个区域交界面上（$A\leftrightarrow B$），$J^{-}$为介质B设入A的中子流密度，$J^{+}$为介质A射入B的中子流密度。根据斐克定律有：

$$\beta =\frac{J^{-}}{J^{+}}=\frac{{\displaystyle \frac{\varphi }{4}+\frac{D}{2}\frac{\text{d}\varphi }{\text{d}x}}}{{\displaystyle \frac{\varphi }{4}-\frac{D}{2}\frac{\text{d}\varphi }{\text{d}x}}}=\frac{1+{\displaystyle \frac{2D}{\varphi }\frac{\text{d}\varphi }{\text{d}x}}}{1-{\displaystyle \frac{2D}{\varphi }\frac{\text{d}\varphi }{\text{d}x}}}$$

对于无限平板反射层，这是反射层内中子通量密度分布为：

$$\varphi =C{e}^{-x/L}$$

因而其反射率为：

$${\beta }_{\mathrm{\infty }}=\frac{1-{\displaystyle \frac{2D}{L}}}{1+{\displaystyle \frac{2D}{L}}}$$

有反射层的内扩散方程的求解，可以得出反照率为：

$$\beta =\frac{1-{\displaystyle \frac{2D}{L}\coth (\frac{a}{L})}}{1+{\displaystyle \frac{2D}{L}\coth (\frac{a}{L})}}$$

💉 4 中子运动长度

📒 4.1 扩散长度

扩散长度：

$$L^2=\frac{D}{{\mathrm{\Sigma }}_a}=\frac{{\lambda }_s{\lambda }_{\text{tr}}}{3}=\frac{{\lambda }_s{\lambda }_{\text{tr}}}{3(1-{\overline{\mu }}_o)}=\frac{1}{3{\mathrm{\Sigma }}_a{\mathrm{\Sigma }}_{\text{tr}}(1-{\overline{\mu }}_o)}$$

由于热中子符合麦克斯韦谱，故扩散长度也符合下式：

$$L^2=L_0^2\sqrt{\frac{T_n}{293}}\frac{g_n(293)}{g_n(Tn)}$$

其中，$L_0$为中子温度在$T_n=293$K的热中子扩散长度。

扩散长度的平方等于热中子从产生地点到被吸收地点穿行直线距离均方值的六分之一：

$$L^2=\frac{1}{6}\overline{r^2}$$

对于平面源，有：$L^2={\displaystyle \frac{1}{2}\overline{x^2}}$

结论3：$L$越大，则热中子从产生到被吸收的直线平均距离越大，而热中子泄露的概率越大。

📚 4.2 慢化长度

慢化长度表征快中子从产生地点到慢化到分界能之下的穿行直线距离：

$$L_l^2=\frac{D}{{\mathrm{\Sigma }}_l}=\frac{1}{3\xi {\mathrm{\Sigma }}_s{\mathrm{\Sigma }}_{\text{tr}}}\ln \frac{E_{\text{tr}}}{E_0}$$

通常也把$L_1^2$称之为热中子年龄，用${\tau }_{\text{th}}$，$\sqrt{{\tau }_{\text{th}}}$即为慢化长度：

$$\tau (E)={\int }_E^{E_0}\frac{D(E)}{\xi {\mathrm{\Sigma }}_s(E)}\frac{\text{d}E}{E}$$

慢化长度平方$L_1^2$或热中子年龄${\tau }_{\text{th}}$等于在无限介质内中子自原点阐述出发在介质中慢化年龄时穿行的直线距离的均方值的六分之一。

$${\tau }_{\text{th}}=\frac{1}{6}\overline{r_d^2}$$

💡 4.3 徙动长度

徙动面积：表征快中子慢化至热中子直至被吸收的穿行直线距离：

$$\begin{gathered} M^2=L^2+{\tau }_{\text{th}} \\ {M}^2=\frac{1}{6}(\overline{r_s^2}+\overline{r_d^2}) \end{gathered}$$

徙动面积是中子由快中子产生至热中子被吸收的穿行的直线距离均方值的六分之一。