一起动手算原子弹的临界质量

偶然看到了 Dr. Jorge S. Diaz 关于如何计算原子弹临界质量的视频 ，觉得挺有意思，于是稍微记录一下推导过程。 下面这个公式就是整篇文章的核心。 \[ \frac{\partial N}{\partial t} = \frac{(\nu - 1)}{\tau} N + \frac{\lambda_t \lambda_f}{3\tau} \nabla^2 N \]

问题假设及目标

既然是要算原子弹的临界质量，我们要先明确问题的条件（假设）以及问题的目标。我们先看原子弹的定义：

所谓原子弹，又称裂变弹，是利用铀、钚等可裂变物质的核裂变反应释放巨大能量而制成的爆炸装置。 — Wikipedia

从定义可知，原子弹主要利用到了核裂变(fission)。所谓核裂变，以 U-235 为例，就是一个 U-235 原子在吸收一个中子(neutron)后，变成不稳定的 U-236 ，最后迅速分解为一个 Kr-92 原子、一个 Ba-141 原子， 以及释放出三个中子；而这三个中子随后会继续轰击其它 U-235 原子，继续得到更多的中子， 从而在短时间内使系统内发生裂变的原子数量呈指数极增加，在这个过程释放出巨大的能量，这个过程表现在宏观上就是“爆炸”。

在上面的分析中，我们不难看出，问题的关键在于中子的数量，准确地讲是弹体中的中子分布密度。 如果中子分布密度可以随时间指数级地递增，那么链式反应就能发生，原子弹就可以爆炸； 反之，如果中子分布密度随时间递减，则原子弹不会爆炸。因此我们将问题简化为求解关于弹体内中子分布密度函数的微分方程。

为此我们不妨做出下面假设：

1. 弹体由均匀的 U-235 组成；
2. 弹体为球形；
3. 弹体中的中子分布密度用含时函数 \(N(\vb{r}, t)\) 表示；
4. 弹体边界处的中子分布密度为0；
5. 中子在遇到 U-235 原子时可能会被吸收引发裂变，也可能发生弹性散射被弹开；

求解目标为：

• 求解关于 \(N(\vb{r}, t)\) 的微分方程，给出维持链式反应最小质量（临界质量）。 至于这个方程长什么样，我们还不得而知。

中子扩散方程的导出

取弹体中的任意一点 \(\vb{r}\)，它经中子密度为 \(N(\vb{r}, t)\) ，如果此时没有中子被原子吸收， 也没有其它的原子产生中子，则中子数守恒， \(\vb{r}\) 处的中子密度满足 流体连续性方程 ：

\begin{equation} \label{eq:continuity-eq} \pdv{N}{t} + \divergence \vb{J} = 0 \end{equation}

公式\eqref{eq:continuity-eq} 中 \(\vb{J}\) 是中子流密度，它代表了单位时间内中子流过某个点的数量， 而 \(\divergence N = \pdv{N_x}{x} + \pdv{N_y}{y} + \pdv{N_z}{z}\) 则表示求它的散度 。 这个公式的含义是在任意一点增加的中子数量等于其流出的中子数量，也即中子数守恒¹。

公式\eqref{eq:continuity-eq} 中的 \(\vb{J}\) 可以从菲克扩散定律(Fick’s law of diffusion)推出：

粒子从高浓度区向低浓度区的流动速率正比于浓度梯度。 — 菲克定律

将菲克定律应用于中子分布密度，可以得到

\begin{equation} \label{eq:fick-law} \vb{J} = - D \grad N \end{equation}

公式\eqref{eq:fick-law} 中 \(\grad N\) 表示求中子的浓度梯度； \(D\) 为扩散常数，它的量纲是 \(\frac{\text{长度}^2}{时间}\) ²，代表了中子扩散的强弱。 一般中子的扩散常数可以使用下面的式子计算

\begin{equation} D = \frac{1}{3} \lambda_t \expval{v} \end{equation}

其中 \(\lambda_t\) 表示中子输运的平均自由程（即发生两次裂变之间中子经过路程的平均值， \(t\) 表示 transport）； \(\expval{v}\) 表示中子的平均运动速率，它可以使用下面的式子近似计算

\begin{equation} \expval{v} = \frac{\lambda_f}{\tau} \end{equation}

其中 \(\lambda_f\) 表示裂变时中子的平均自由程。自由程可以使用宏观的散射截面来计算³：

\begin{equation} \begin{aligned} \lambda = \frac{1}{\Sigma} = \frac{1}{N \sigma} = \frac{1}{(\rho N_A/M) \sigma} = \frac{M}{\rho N_A \sigma} \end{aligned} \end{equation}

其中 \(\Sigma\) 是宏观散射截面，\(N\) 是单位体积内的原子核数量， \(\sigma\) 是微观散射截面， \(\rho\) 是宏观密度， \(N_A\) 是阿伏加德罗常数， \(M\) 是原子相对质量。对于裂变自由程 \(\lambda_f\) ， 其微观散射截面就是裂变散射截面 \(\sigma_f\) ；而对于输运自由程 \(\lambda_t\) ， 其微观散射截面就是裂变散射截面由裂变散射截面和弹性散射截面共同组成 (\(\sigma_t = \sigma_f + \sigma_e\)) 。

但在原子弹中，中子可以被 U-235 吸收，也可以通过裂变产生，那么中子数量并不守恒， 因此这个方程的右手边需要加入中子产生项和中子吸收项：

\begin{equation} \label{eq:continuity-eq-with-source} \pdv{N}{t} + \divergence \vb{J} = S_+ - S_- \end{equation}

公式\eqref{eq:continuity-eq-with-source} 的右手边 \(S_+\) 表示产生的中子数， \(S_-\) 表示吸收的中子数， 它们分别可以写为

\begin{align} \label{eq:neutron-source-sink-term} S_+ &={} \nu \frac{N}{\tau} \quad \text{裂变产生的中子} \\
S_- &={} \frac{N}{\tau} \quad \text{裂变消耗的中子} \end{align}

式\eqref{eq:neutron-source-sink-term} 中 \(\nu\) 表示每次裂变反应生成的中子数，而 \(\tau\) 表示平均发生两次裂变反应的间隔时间。

现在我们把公式\eqref{eq:continuity-eq-with-source} 各项都展开，有 \[ \pdv{N}{t} + \divergence (-D \grad N) = \frac{(\nu - 1)}{\tau} N \] 使用矢量运算关系式 \(\divergence (\grad N) = \laplacian N = \pdv[2]{N}{x} + \pdv[2]{N}{y} + \pdv[2]{N}{z}\) ，以及代入扩散常数 \(D\) 的表达式， 对中子扩散方程进行整理，可以得到文章开头的那个公式

\begin{equation} \label{eq:neutron-diffusion-equation} \pdv{N}{t} = \frac{(\nu - 1)}{\tau} N - \frac{\lambda_t \lambda_f}{3\tau} \laplacian N \end{equation}

这便是中子扩散方程，也是我们计算原子弹临界质量的理论基础。下面我们来解这个方程，并推导出临界质量。

中子扩散方程的求解

不妨仔细观察这个方程，我们把需要求的变量用红色标出 \[ \pdv{\color{red}N}{t} = \frac{(\nu - 1)}{\tau} {\color{red}N} - \frac{\lambda_t \lambda_f}{3\tau} \laplacian {\color{red}N} \] 不难看出这是一个关于 \(N\) 的偏微分方程，议程左边只依赖于时间 \(t\) ，右边只依赖于空间坐标 \(x,y,z\) ， 因此可以使用分离变量法来求解。我们可以将 \(N(\vb{r}, t)\) 写成两个函数的积，其中 \(f(\vb{r})\) 只与位置 \(\vb{r}\) 有关， \(g(t)\) 只与时间 \(t\) 相关：

\begin{equation} \label{eq:variable-separation} N(\vb{r}, t) = f(\vb{r}) g(t) \end{equation}

此时将公式\eqref{eq:variable-separation} 代入方程\eqref{eq:neutron-diffusion-equation} ，可得

\begin{aligned} & f \pdv{g}{t} = \frac{\nu - 1}{\tau} fg + \frac{\lambda_t \lambda_f}{3\tau} (\laplacian f)g \\
\implies & \frac{1}{g} \pdv{g}{t} = \frac{\nu - 1}{\tau} + \frac{\lambda_t \lambda_f}{3\tau} \frac{\laplacian f}{f} \equiv K = \frac{\nu'}{\tau} \end{aligned}

我们得到两个关系式

\begin{gather} \frac{1}{g} \pdv{g}{t} = \frac{\nu'}{\tau} \label{eq:g-pde} \\
\frac{\nu'}{\tau} = \frac{(\nu - 1)}{\tau} + \frac{\lambda_t \lambda_f}{3\tau} \frac{\laplacian f}{f} \label{eq:nu-prime-definition} \end{gather}

其中方程\eqref{eq:g-pde} 是只关于 \(g\) 的微分方程，且方程右边与 \(g\) 无关，故我们将等式右边重新定义为 \(\frac{\nu'}{\tau}\) ，其具体形式如公式\eqref{eq:nu-prime-definition} 所示。

时间演化方程

\(g\) 只与时间相关，所以式\eqref{eq:g-pde} 是中子分布密度的时间演化方程，它也是 一个典型的常微分方程， 我们可以轻松写出它的通解

\begin{align} \frac{1}{g} \dv{g}{t} &={} \frac{\nu'}{\tau} \notag \\
\implies \dv{\ln{g}}{t} &={} \frac{\nu'}{\tau} \notag \\
\implies g(t) &={} A e^{\nu’t/\tau} \label{eq:g-solution} \end{align}

公式\eqref{eq:g-solution} 显示 \(g\) 是一个指数函数，虽然其指前因子 \(A\) 与指数部分 \(\nu'\) 待定，但现在我们至少可以得出一个结论，即要想中子的数量短时间大量扩增， \(g\) 必须为增函数， 即指数 \(\nu'\) 必须为正，否则 \(g\) 是减函数，中子数量会逐渐衰减。 Alex Wellerstein 制作了一个网站用于演示裂变反应的具体过程 ，它可以生动地反映当 \(\nu'\) 大于零以及小于零时系统中裂变反应的剧烈程度，读者有兴趣可以尝试⁴。

空间分布方程

下面我们来求解 \(f\) ，其只与 \(x,y,z\) 有关，因此它描述的是中子的空间分布。 对公式\eqref{eq:nu-prime-definition} 进行化简，我们可以得到关于 \(f\) 的微分方程：

\begin{align} &\frac{\nu'}{\tau} = \frac{(\nu - 1)}{\tau} + \frac{\lambda_t \lambda_f}{3\tau} \frac{\laplacian f}{f} \notag \\
\implies & \frac{\lambda_t \lambda_f}{3\tau} \frac{\laplacian f}{f} + \frac{\nu - 1 - \nu'}{\tau} = 0 \notag \\
\implies & \laplacian f + \pqty{\frac{3(\nu - \nu' - 1)}{\lambda_t \lambda_f}} f = 0 \notag \\
\implies & \laplacian f + \kappa^2 f = 0 \label{eq:f-equation} \end{align}

方程\eqref{eq:f-equation} 类似振子振动的方程，其中 \(f\) 前面的系数我们用 \(\kappa^2 = \pqty{\frac{3(\nu - \nu' - 1)}{\lambda_t \lambda_f}}\) 表示。

由于弹体是球形，系统具有球对称性，我们在球坐标系下求解这个方程会更容易：

\begin{equation} f(\vb{r}) = f(r, \theta, \phi) \equiv f( r) \end{equation}

这里我们舍去角向的 \(\theta, \phi\) 部分，只考虑径向 \(r\) 的分布。拉普拉斯算符在球坐标系下可以转化为

\begin{equation} \label{eq:spherical-laplacian} \laplacian f = \frac{1}{r^2} \dv{}{r} \pqty{ r^2 \dv{f}{r} } \end{equation}

如果直接将公式\eqref{eq:spherical-laplacian} 代入方程\eqref{eq:f-equation} ， 则很难求解这个微分方程，因此这里使用一个辅助函数 \(u( r) = r f( r)\) 使得 \(f( r) = u( r) / r\) ，此时

\begin{equation} \label{eq:laplacian-auxiliary} \laplacian f = \frac{1}{r} \dv[2]{u}{r} \end{equation}

将式\eqref{eq:laplacian-auxiliary} 代入公式\eqref{eq:f-equation} ，可得

\begin{equation} \label{eq:u-function} \dv[2]{u}{r} + \kappa^2 u = 0 \end{equation}

不难解得

\begin{equation} \label{eq:u-solution} u( r) = B \cos(\kappa r) + C \sin(\kappa r) \end{equation}

将 \(u\) 代回 \(f( r) = u( r) / r\) ，可得

\begin{equation} \label{eq:f-solution} f( r) = B \frac{\cos(\kappa r)}{r} + C \frac{\sin(\kappa r)}{r} \end{equation}

此时关于 \(N\) 的两个分量都有了表达式，将式\eqref{eq:g-solution} 与式\eqref{eq:f-solution} 合并，可得 \(N\) 关于半径和时间的函数关系式

\begin{gather} N(r, t) = A e^{\nu' t/\tau} \pqty{B \frac{\cos(\kappa r)}{r} + C \frac{\sin(\kappa r)}{r}} \\
\kappa^2 = \frac{3(\nu - \nu' - 1)}{\lambda_t \lambda_f} \label{eq:kappa-squared} \end{gather}

这个式子还是有些复杂，并且还包含更多待定的因子，如 \(A, B, C\) 等，我们需要利用边界条件对它进行化简。

边界条件

1. 弹体中心处的中子密度是有限的；

   这要求 \(N(r=0,t) \ne \infty\) ， 考虑到 \(\displaystyle \lim_{r\to 0} B \frac{\cos(\kappa r)}{r} \to \infty\) ，因此 \(B\) 必须为零； 而 \(\displaystyle \lim_{r\to 0} C \frac{\sin(\kappa r)}{r} \to C \kappa\) ，因此 \(C \ne 0\) 。

2. 弹体中心处的初始中子密度为 \(N_0\) ；

   即 \(N(r=0, t=0) = N_0\) ，对应于 \(\displaystyle \lim_{\substack{r\to 0 \\ t \to 0}} N(r,t) = AC\kappa\) ，由此可以推出 \(AC = \frac{N_0}{\kappa}\) 。

3. 弹体足够大，以至边界处的中子数为零。

   即 \(N(r=R, t) = N_0 e^{\nu’t/\tau} \pqty{\frac{\sin(\kappa R)}{\kappa R}} = 0\) ， 上面的式子要求满足 \(\sin(\kappa R) = 0\) ，因此可以取 \(\kappa = \frac{\pi}{R}\) 。 结合前面我们得到的 \(\kappa^2\) 关系式，我们可以求出待定参数 \(\nu'\)

   \begin{gather} \kappa^2 = \frac{3(\nu - \nu' - 1)}{\lambda_t \lambda_f} = \frac{\pi^2}{R^2} \notag \\
   \implies \nu' = (\nu - 1) - \frac{\pi^2 \lambda_t \lambda_f}{3R^2} \notag \end{gather}

临界质量的导出

在上一节，我们得到了关于 \(\nu'\) 的表达式，而我们在时间演化方程那一节曾讨论过，原子弹爆炸的一个必要条件是 \(\nu' > 0\) ，因此我们可以得到弹体的临界半径

\begin{gather} \nu' = (\nu - 1) - \frac{\pi^2 \lambda_t \lambda_f}{3R^2} = 0 \label{eq:nu-prime-condition} \\
\implies R_c = \pi \sqrt{\frac{\lambda_t \lambda_f}{3(\nu - 1)}} % = \frac{\pi M}{\rho N_A} \pqty{\frac{1}{3\sigma_t \sigma_f (\nu-1)}}^{1/2}\\
\end{gather}

进而根据半径算出体积，以及临界质量的表达式

\begin{equation} \implies M_c = \frac{4}{3} \pi \rho R_c^3 % = \frac{4\pi^4 M^3}{3\rho^2 N_A^3} % \pqty{\frac{1}{3\sigma_t \sigma_f (\nu - 1)}}^{3/2} \end{equation}

代入实验上测得的微观数据

数据名称	符号	数值
中子寿命（平均两次裂变反应间隔时间）	\(\tau\)	16.5 ns
平均每次裂变产生的中子数	\(\nu\)	2.6
中子输运自由程	\(\lambda_t\)	3.596 cm
中子裂变自由程	\(\lambda_f\)	16.89 cm
铀-235 的密度	\(\rho\)	19.050 g/cm³

可以计算得到临界半径和临界质量分别为

\begin{gathered} R_c \approx 11 \, \text{cm} \\
M_c \approx 106 \, \text{kg} \end{gathered}

以上是一种“基础解”，它简单地要求边界处的中子数为零，这个条件过于严格。 如果我们放宽边界条件的限制，其临界质量也会有所变化。 另一种常用的边界条件是在边界处允许一半的中子流逃逸到弹体外面，即

\begin{gather} J( R) = -D \pdv{\phi}{r}\bigg|_R \ge \frac{1}{2} \phi( R), \quad \phi( r) = N( r) \expval{v} \notag \\
\implies N(r=R, t) \le - \frac{2}{3} \lambda_t \pdv{N}{r}\bigg|_{r=R} \label{eq:marshak-bc} \end{gather}

式\eqref{eq:marshak-bc} 也被称为 Marshak 边界条件⁵，此时可以得到

\begin{equation} \frac{3R}{2\lambda_t} \le 1 - \frac{\kappa R}{\tan (\kappa R)} \label{eq:marshak-critical-radius} \end{equation}

联系到我们使用公式\eqref{eq:nu-prime-condition} 求临界质量时要求 \(\nu' = 0\) ， 这里我们也使用这个条件，并将它与其它参数代入式\eqref{eq:kappa-squared} ，可以得到 \(\kappa\) 的具体值，然后我们可以使用计算机求解方程\eqref{eq:marshak-critical-radius} ，便可以得到一个新的临界半径 \(R_c\)，及其对应的临界质量 \(M_c\)： \[ R_c \approx 8.37\;\text{cm}; \quad M_c \approx 46\;\text{kg} \] 至此，我们得到了 Marshak 边界条件下的 U-235 原子弹临界质量，相比于实际临界质量(52 kg)⁶， 这个值又偏小，说明 Marshak 边界条件又过于宽松。

小结与胡扯

以上便是从中子扩散理论出发得到的两种原子弹临界质量，分别为 106 kg 和 46 kg ，尽管它们与实际值有着不小的差别，但其推导过程仍具有启发意义。在此再次感谢 Dr. Jorge S. Diaz 博主的影片， 其在视频最后还留了几道家庭作业，分别是计算钸-239的临界质量以及正方体和圆柱体弹体的临界质量。 鉴于本人对做 Homework 并不感兴趣，这几道题还是留给读者吧。 如果读者对其它推导方法感兴趣，可以参考这篇文献 ⁷。

稍微讲一些胡说八道吧，本文抓住了中子分布密度这一核心变量来建模整个原子弹， 通过求解它的含时方程与空间分布方程，我们可以得到整个原子弹的临界半径与临界质量， 这一步看似简单，但其实并不容易，需要建模人员有敏锐的直觉和丰富的经验。 另外一个点是，文中推导的中子扩散方程实际是以菲克扩散定律为基础，而后者可以从微观粒子的随机行走理论推导得到，这实际上忽略了中子与中子的相互作用而只考虑了中子与铀原子核的相互作用， 如果考虑中子与中子的相互作用，那么得到的临界质量又会是另外一个值。 最后要说的是，推导的最后，我们仅仅改变了中子扩散方程的边界条件，所得到的临界质量就大相径庭， 这说明我们在建模时不仅需要简化模型，也需要适当考虑更多的因素，这实际上相当考验建模人员的数学功底和物理直觉，至少本人自己做不到^_^。

文章结尾再送读者一篇文献⁸，如果你对海森堡算错临界质量这件事比较感兴趣可以看这篇文献， 文章系统描述了二战后海森堡在和平环境下用一周的时间内算出的结果， 他甚至还考虑了带中子反射层的边界条件，得到的临界质量与实验值“惊人地一致”， 这也侧面说明研究环境对于研究人员状态的影响有多大。

---

1. 如果把\(N\)换成质量，它也就是流体力学中的质量守恒定律； ↩︎

2. 如果单独搜索“扩散系数”，会有资料显示它的量纲是长度 (\(D = \frac{\lambda}{3}\)， \(\lambda\) 是平均自由程)，此时 \(D\) 中的速度项被 \(\vb{J} = -D \grad \vb{\phi}\) 中的 \(\vb{\phi}\) 所吸收， 详见这个链接 。 ↩︎

3. 详见 https://clive.semmens.org.uk/Energy/FreePath.html  ↩︎

4. 裂变反应模拟网站 https://blog.nuclearsecrecy.com/misc/criticality/  ↩︎

5. Noh, Taewan. “A Study on Diffusion Approximations to Neutron Transport Boundary Conditions.” Journal of the Nuclear Fuel Cycle and Waste Technology(JNFCWT) 16, no. 2 (2018): 203–9. https://doi.org/10.7733/jnfcwt.2018.16.2.203 . ↩︎

6. https://en.wikipedia.org/wiki/Critical_mass  ↩︎

7. Reed, B. Cameron. “Estimating the Critical Mass of a Fissionable Isotope.” Journal of Chemical Education 73, no. 2 (1996): 162. https://doi.org/10.1021/ed073p162 . ↩︎

8. El Naschie, M.S. “Heisenberg’s Critical Mass Calculations for an Explosive Nuclear Reaction.” Chaos, Solitons & Fractals 11, no. 6 (2000): 987–97. https://doi.org/10.1016/S0960-0779(99)00110-1 . ↩︎