摘要:磁化套筒惯性聚变(magnetized liner inertial fusion, MagLIF)结合了传统磁约束聚变和惯性约束聚变的优点, 理论上可以显著地降低聚变实现的难度, 具有极大的应用潜力. 以研究MagLIF中的关键问题为目标, 建立能够综合考虑磁化、预加热、套筒内爆、聚变反应、端面效应、磁通压缩等多种复杂机制在内的集成化物理模型, 特别是通过引入流体喷射模型, 使得可以在一维计算条件下考虑具有二维特性的端面损失情况, 并额外考虑Nernst扩散项对磁通损失的影响. 在此基础上编写实现一维集成化MagLIF数值模拟程序MIST (magnetic implosion simulation tools), 与FP-1装置(2 MA, 7.2 μs)上铝套筒内爆实验结果的对比验证了程序磁流体模块的正确性; 将聚变模块纳入后与国外同类程序LASNEX和HYDRA计算结果进行整体比较, 所得数值结果总体接近, 主要差异体现在燃料温度的计算上, 对可能影响的原因进行了简要分析. 所建立的集成化模型与程序将为未来开展MagLIF聚变实验研究提供坚实的理论基础和重要工具.
关键词: 磁化套筒惯性聚变/
集成化数值模拟/
端面损失/
Nernst效应

English Abstract

全文HTML

--> --> -->

能源是人类社会和文明发展的基础. 进入21世纪后, 核能在人类社会消耗能源中所占比例^[1]越来越高, 而可控核聚变被认为是解决未来能源问题的重要手段之一. 国际可控核聚变领域的研究多集中于磁约束聚变(magnetic confinement fusion, MCF)^[2-4]和惯性约束聚变(inertial confinement fusion, ICF)^[5-7]这两大主流方向上, 它们主要区别在于实现聚变条件的具体方式不同. 与此同时, 国内聚变研究领域也有较快的发展, 其中有代表性的如: 中国科学院等离子体物理研究所的东方超环(EAST)装置, 主要用于MCF物理过程^[8,9]研究; 中国工程物理研究院建立的神光系列装置和“聚龙一号”装置, 广泛应用于激光^[10,11]和Z箍缩^[12,13]驱动ICF研究. 然而, 传统的聚变方式都面临着不同的工程和技术问题, 距离点火的目标仍有一定的差距, 因此有必要探索新的聚变构型, 以尽量降低实现聚变的难度.
2010年由美国圣地亚实验室提出^[14]的磁化套筒惯性聚变(magnetized liner inertial fusion, MagLIF), 是一种在ICF中引入轴向磁场, 并结合传统MCF与ICF优势的新聚变构型, 该聚变整个过程包括3个主要阶段: 燃料磁化(magneti-zation)、激光预加热(laser pre-heating)和套筒压缩(compression), 如图1所示. 通过简单的计算可以看出MagLIF构型对氘氚(DT)反应产生的α粒子能量利用率很高, 因此在理论上能够显著降低聚变实现的难度, 有极大的应用潜力. 初始时刻MagLIF中引入了较强(～30 T)的预置轴向磁场, 在外套筒的内爆压缩作用下, 磁场与DT燃料一起被压缩, 在压缩最紧密时刻磁场可达到数千特斯拉. 一般条件下DT反应产生的α粒子的平均自由程^[6]为${\lambda _{\rm{\alpha }}}\! \approx 0.5 \times {10^{12}}T_{\rm{e}}^{3/2}\!/(n\ln \varLambda )$$ (T_{\rm e} < 5 \times $ $10^8~{\rm K})$. 若燃料温度为7 keV, 粒子数密度10²²/cm³, 取库仑对数lnΛ = 9时λ_α大约是5 cm, 远大于MagLIF实验负载尺寸(～1 cm). 引入轴向磁场(假设峰值为2000 T)后, 根据公式${r_{\rm{L}}} = 2.71 \times {10^5}/B$, 可计算得到α粒子径向回旋半径为0.135 mm, 此时燃料内部的α粒子能量沉积效率显著提高.
图 1 MagLIF过程示意图(包含3个主要阶段)
Figure1. Schematic of MagLIF process, including three main stages.

MagLIF构型自提出后在理论模拟与实验上均取得快速的进展^[15-20], 圣地亚实验室取得的成就引发全世界科学家对该领域的关注, 未来将朝着实现点火的目标继续发展. 在本文作者前期开展的零维聚变模型的研究中^[21], 氘氚燃料区被简化为质点(即密度、压强、磁场等各种物理量在燃料区都被认为是均匀分布, 燃料能量方程由解析表达式描述), 模型相对粗糙. 本文工作在前期基础上有了较大的进步, 在建模中考虑了α粒子沉积、输运以及Nernst效应等多种物理模型, 特别是引入流体喷射模型使得程序可以在一维计算中考虑二维端面损失效应, 并将燃料区和金属套筒统一纳入磁流体力学控制方程的描述中, 构建了完整的一维集成化数值模拟程序.
本文结构安排如下: 首先介绍基本物理模型, 包括对于整个物理过程的认知、磁流体动力学(magneto-hydro-dynamic, MHD)方程组、α粒子能量方程、聚变模型以及磁场作用下的径向热扩散等; 随后引入了流体喷射模型, 用于在一维计算条件下考虑端面效应的影响, 并在磁扩散方程中增加了由于温度梯度导致的Nernst扩散项; 然后程序校验部分通过与FP-1装置上铝套筒内爆实验结果的对比验证了磁流体模块的正确性, 并将聚变模块纳入后与国外同类程序LASNEX和HYDRA计算结果进行了整体比较, 以检验物理建模与程序代码的合理性; 最后是总结与展望.

对于MagLIF而言, 发生聚变时刻状态燃料处于高温(5—10 keV)、高密度(n_e > 10²⁰/cm³), 可视为全电离、理想导电流体, 驱动源来自外部电流, 同时有轴向磁场的存在, 整体上使用MHD方程组描述是合适的. 但MagLIF包含非常复杂的多个物理过程, 为便于一体化的描述, 必须适当简化:
1)认为程序计算时, DT燃料为完全电离的等离子体, 燃料光性薄, 金属套筒光性厚, 燃料中的轫致辐射完全沉积在套筒内壁上.
2)估算结果表明, MagLIF典型参数设置条件下, 因聚变反应DT燃料损失质量不超过1%, 因此程序中不考虑燃料质量损失.
3)根据对激光预加热$r\text{-}\theta $平面二维数值模拟结果^[16], 在预加热后很短时间(～20 ns)内电子和离子温度快速取得平衡, Z装置上首批氘氘(DD)实验的结果^[18]也证实迟滞时刻电子与离子温度差别很小, 因此程序使用单温MHD模型.
2

-->

2.1.MHD方程组

基于上述考虑, 使用拉氏描述方式, 建立描述MagLIF过程的单温、单流体、多介质MHD方程组:

$\frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} + \rho \nabla \cdot {{u}} = 0,$

$\rho \frac{{{\rm{d}}{{u}}}}{{{\rm{d}}t}} = {{J}} \times {{B}} - \nabla p,$

$\rho \frac{{{\rm{d}}e}}{{{\rm{d}}t}} + p\nabla \cdot {{u}} + \nabla \cdot {{q}} = \eta {J^2} + \rho w, $

$\;\;\frac{{{\rm{d}}{{B}}}}{{{\rm{d}}t}} = - \nabla \times \left(\frac{\eta }{{{\mu _0}}}\nabla \times {{B}}\right) - {{B}}(\nabla \cdot {{u}}) + ({{B}} \cdot \nabla ){{u}}, $

式中, ρ为密度, u为速度矢量, p为压强, e为比内能, q为热流密度矢量, w为外部能量沉积率, ${\mu _0}$为真空磁导率, $\eta $为电阻率, B为磁场强度, $ {{J}} \equiv $ $ \dfrac{1}{{{\mu _0}}}\nabla \times {{B}}$为电流密度.
套筒和DT燃料区可统一采用方程组(1)—(4)描述, 二者的区别在于能量方程(3)中的外部能量沉积项$\rho w$. 对于金属套筒$\rho w = 0$. 对于DT燃料, 由于聚变反应, 需要考虑轫致辐射和α粒子能量沉积效应, 则能量沉积项为

$ \rho w = {G_{\text{α}}}{E_{\text{α}}} - {Q_{{\rm{rad}}}}, $

式中, $E_{\text{α}}$为单位体积的α能量; $G_{\text{α}}$为α能量耗散系数; ${Q_{{\rm{rad}}}}$为单位体积的轫致辐射损失项.
2 -->

2.2.α粒子能量方程

α粒子能量$E_{\text{α}}$单独满足一个演化方程, 根据文献[22]该方程可写为

$\frac{{{\rm{d}}{E_{\text{α}}}}}{{{\rm{d}}t}} = \nabla \cdot ({D_{\text{α}} }\nabla {E_{\text{α}} }) - {G_{\text{α}} }{E_{\text{α}} } + \dot n{E_{\text{α} 0}}, $

式中, $D_{\text{α}}$为方程扩散系数, $G_{\text{α}}$为方程耗散系数, $\dot n$为DT聚变反应率, ${E_{\text{α}0}}$ = 3.5 MeV为聚变反应产生的α粒子初始能量.
(6)式中扩散系数${D_{\text{α}}}$与磁场强度相关, 具体表达式(高斯制)为

${D_{\text{α}}} = \frac{{{V_{{\text{α}}}}{l_{\text{α}}}}}{{18 + 2{{\left( {\dfrac{{{\omega _{\text{α}}}}}{{{\nu _{\text{α}}}}}} \right)}^2}}}\,,$

式中, $ {l_{\text{α}}} = 0.107\dfrac{{T_{{\rm{keV}}}^{3/2}}}{{\rho {L_{\text{α}}}}}\;{\rm{(cm}}),\; {\omega _{\text{α}}} = \dfrac{{{2e_{\rm{g}}}{B_{\rm{g}}}}}{{{m_{\text{α}}}{c_0}}},\;{\nu _{\text{α}}} = \dfrac{{{V_{{\text{α}}0}}}}{{{l_{\text{α}}}}} $, 其中c₀为真空中的光速; ${L_{\text{α}}}$为α粒子库仑对数, 聚变条件下一般^[22]取7; ${m_{\text{α}}}$为α粒子质量; e_g为高斯制下电子电荷; B_g为高斯制下磁场强度.
能量方程耗散系数为

${G_{\text{α}}} = 2{\nu _{\text{α}}}. $

2 -->

2.3.DT聚变模型

为便于计算, 本文仅考虑DT聚变初级反应:

${\rm{D}} + {\rm{T}} \to { {\text{α} }}\;({\rm{3}}{\rm{.5~MeV}}) + {\rm{n}}\;({\rm{14}}{\rm{.1~MeV}}). $

相应的聚变反应率为

$\dot n = {n_{\rm{D}}}{n_{\rm{T}}}{\left\langle {\sigma v} \right\rangle _{{\rm{DT}}}} = \frac{1}{4}n_{\rm{e}}^{\rm{2}}{\left\langle {\sigma v} \right\rangle _{{\rm{DT}}}}, $

式中${\left\langle {\sigma v} \right\rangle _{{\rm{DT}}}}$为DT聚变反应截面, 根据文献[23]可近似写作:

${\left\langle {\sigma v} \right\rangle _{{\rm{DT}}}} = {C_1}{A^{ - 5/6}}{B^2}\exp ( - 3{A^{1/3}}B),$

式中, $A = 1 - \dfrac{{{C_2}T + {C_4}{T^2} + {C_6}{T^3}}}{{1 + {C_3}T + {C_5}{T^2} + {C_7}{T^3}}}, B = \dfrac{{{C_0}}}{{{T^{1/3}}}}$, T为燃料温度, 系数C₀到C₇为拟合实验曲线得到的数值, 具体取值如表1所列.

C0 /keV1/3	C1/cm3·s–1	C2/keV–1	C3/keV–1
6.661	643.41×10–16	15.136×10–3	75.189×10–3
C4/keV–2	C5/keV–2	C6/keV–3	C7/keV–3
4.6064×10–3	13.5×10–3	–0.10675×10–3	0.01366×10–3

表1系数C₀—C₇的取值
Table1.Values of coefficient C₀–C₇

2 -->

2.4.磁场影响下的径向热扩散

轴向强磁场的引入还会影响燃料中电子和离子沿径向的热扩散, 根据文献[22], 磁场影响下径向热扩散系数(高斯制)如下:

$\begin{split}& {K_{{\rm{ce}}}} = 3.16{P_{\rm{e}}} \cdot f({x_{\rm{e}}}),{K_{{\rm{ci}}}} = 3.9{P_{\rm{i}}} \cdot f({x_{\rm{i}}}),\\& \quad {P_{\rm{e}}} = \dfrac{{nT{\tau _{\rm{e}}}}}{{{m_{\rm{e}}}}}, ~{P_{\rm{i}}} = \dfrac{{nT{\tau _{\rm{i}}}}}{{{m_{\rm{i}}}}};\\&f({x_{\rm{e}}}) = \frac{{1 + 0.39x_{\rm{e}}^{\rm{2}}}}{{1 + 3.9x_{\rm{e}}^2 + 0.26x_{\rm{e}}^{\rm{4}}}}, \\& \qquad{x_{\rm{e}}} = {\omega _{\rm{e}}}{\tau _{\rm{e}}}, ~{\omega _{\rm{e}}} = \dfrac{{{e_{\rm{g}}}{B_{\rm{g}}}}}{{{m_{\rm{e}}}{c_0}}}, \\& \qquad {\tau _e} = \dfrac{{3{m_{\rm{e}}}^{1/2}{T^{3/2}}}}{{4\sqrt {2{\text{π}}} {e_{\rm{g}}}^4 n{L_{\rm{e}}}}},~ {L_{\rm{e}}} = 7; \\& f({x_{\rm{i}}}) = \frac{{1 + 0.756x_{\rm{i}}^{\rm{2}}}}{{1 + 3.99x_{\rm{i}}^{\rm{2}} + 1.48x_{\rm{i}}^{\rm{4}}}}, \\& \qquad {x_{\rm{i}}} = {\omega _{\rm{i}}}{\tau _{\rm{i}}}, {\omega _{\rm{i}}} = \dfrac{{{e_{\rm{g}}}{B_{\rm{g}}}}}{{{m_{\rm{i}}}{c_0}}}, \\& \qquad {\tau _{\rm{i}}} = \dfrac{{3{m_{\rm{i}}}^{1/2}{T^{3/2}}}}{{4\sqrt {\text{π}} {e_{\rm{g}}}^4 n{L_{\rm{i}}}}},~ {L_{\rm{i}}} = 9.\end{split}$

这里, c₀为真空中的光速; m_e为电子质量; m_i为1∶1混合DT燃料等效离子质量, m_i = 2.5m_p; L_e和L_i分别为电子与离子的库仑对数, 聚变条件下DT燃料密度约为1 g/cm³, 温度约为10 keV, 此时L_e和L_i分别取值^[22]为7和9; e_g为高斯制下电子电荷; B_g为高斯制下磁场强度.

为了实现对MagLIF中的氘氚燃料的激光预加热, 需要在聚变靶顶部开孔(用聚合物薄膜进行覆盖密封), 称为激光注入孔(laser entrance hole, LEH). 外部激光器所产生的激光束通过此孔进入燃料内部, 将燃料快速(～2 ns)加热至等离子体态, 满足预加热所需的温度要求(50—300 eV). 但是这个开孔的存在会导致在套筒压缩阶段内部的燃料向外喷射, 造成燃料质量和内能的损失, 这就是所谓的端面损失效应, 如图2所示. 端面效应的准确模拟需要至少在二维情形下进行, 本文中为了在一维中考虑端面效应, 提出近似的简化模型并在一维方程中建立描述端面效应的质量、能量损失项. 下面予以详细讨论.
图 2 端面效应简化模型示意图
Figure2. Schematic of simplified model describing end loss effect.

LEH仅占整个燃料区间的一部分, 半径用r_LEH来表示, 程序中仅当燃料网格位置r小于r_LEH时, 才考虑端面效应, 其他情况下忽略. 对于r_LEH范围内的燃料, 忽略流体黏性后, 可采用喷射模型来描述. 查询物理手册可知, 管道内相对于外界超压为P的流体, 其喷射速度${C_{{\rm{gs}}}}$为

${C_{{\rm{gs}}}} = \sqrt {2P/\rho } .$

在程序计算的单个时间步dt内, 对于r < r_LEH内任一层网格, 如果单独提取出来并考虑轴向运动情况, 如图2中右侧柱体所示. 认为每个时间步内顶部${C_{{\rm{gs}}}} \times {\rm{d}}t$高度内的燃料质量dM和内能dE全部向外喷射丢失, 剩余的质量和内能在高度h内重新均匀分配并计算密度. 由此可以写出单个时间步dt内喷射质量和内能的损失为

${\rm{d}}M = {2\text{π}}r \cdot {\rm{d}}r \cdot {C_{{\rm{gs}}}} \cdot {\rm{d}}t \cdot \rho , $

${\rm{d}}E = - P \cdot S \cdot {C_{{\rm{gs}}}} \cdot {\rm{d}}t,\;\;\;\;E = e \cdot M, $

式中, r为网格起始沿径向所在位置, $\rho $为燃料密度, E为网格内燃料的总内能, e为网格内燃料的比内能, M为网格内燃料的总质量, S为网格的底面积.
受端面效应影响, 原有MHD方程组中质量和能量方程应增加一项损失项, 由(14)和(15)式推导后可得考虑端面效应后新的控制方程如下:

$\frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} + \rho \nabla \cdot {{u}} = - \frac{{{C_{{\rm{gs}}}}}}{h}\rho, $

$\rho \frac{{{\rm{d}}e}}{{{\rm{d}}t}} + P\nabla \cdot {{u}} + \nabla \cdot {{q}} = \rho w - \frac{{{C_{{\rm{gs}}}}}}{h}P, $

式中h为燃料高度.
端面效应的流体喷射模型未考虑沿轴向流体分布和压力做功细节, 从物理模型来看, 可能高估了质量和内能损失, 损失率上应乘以系数Φ₁; 若考虑流体摩擦(黏性), 应乘以速度指数Φ₂ (对于水, Φ₂约为0.97); 若考虑喷嘴几何构型收缩, 应再乘以收缩系数Φ₃ (对于锐喷射口, Φ₃约为0.61).

为了将α粒子能量约束在氘氚燃料区的范围以内, 需要数千特斯拉的超强磁场, 在MagLIF中这是通过压缩冻结在燃料中的轴向磁场来实现的. 磁冻结效应由(4)式等号右端的第2项和第3项描述. 轴向磁场还会发生扩散损失(如图3所示), 这包含了两种机制, 一种是氘氚燃料和金属套筒的电阻率引起的普通磁扩散, 另一种是径向温度梯度导致的Nernst磁扩散效应. 普通磁扩散由(4)式等号右端的第1项描述.
图 3 磁通压缩与扩散过程示意图
Figure3. Schematic of magnetic flux compression and diffusion process.

Nernst效应的基本原理是: 在完全电离的等离子体中, 磁场被冻结在承担传导功能的电子上, 如果垂直于磁场方向存在温度梯度, 电子将沿热流方向运动, 并导致磁场沿热流方向发生额外的扩散. 该效应与温度梯度成正比, 借鉴文献[14]的写法并考虑磁通量守恒的要求后, 给出的表达式为

$\frac{{{\rm{d}}{{B}}}}{{{\rm{d}}t}} = \nabla \times \left[ {\frac{{{B}}}{{e\left| {{B}} \right|{\omega _{{\rm{ce}}}}{\tau _{\rm{e}}}}} \times \nabla ({k_{\rm{B}}}T)} \right], $

式中, B为轴向磁场强度, ${\omega _{{\rm{ce}}}}{\tau _{\rm{e}}}$为电子霍尔参量, k_B为玻尔兹曼常数.
将(18)式代入MHD方程组中磁扩散方程内, 可以写出MagLIF过程中燃料内部轴向磁场扩散方程的整体表达形式:

$\begin{split} \frac{{{\rm{d}}{{B}}}}{{{\rm{d}}t}} =\, & - \nabla \times \left(\frac{\eta }{{{\mu _0}}}\nabla \times {{B}}\right) - {{B}}(\nabla \cdot {{u}}) + ({{B}} \cdot \nabla ){{u}}\\ & + \nabla \times \frac{{{B}}}{{e\left| {{B}} \right|{\omega _{{\rm{ce}}}}{{{\tau}} _{\rm{e}}}}} \times \nabla ({k_{\rm{B}}}T).\\[-15pt] \end{split}$

可以进一步推导出(19)式在一维柱坐标下的表达形式:

$\begin{split} \frac{{{\rm{d}}B}}{{{\rm{d}}t}} =\, & \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{\eta }{{{\mu _0}}}\frac{{\partial B}}{{\partial r}}} \right) - \frac{1}{r}\frac{\partial }{{\partial r}}(ru)B \\ & + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{{k_{\rm{B}}}}}{{e{\omega _{{\rm{ce}}}}{\tau _{\rm{e}}}}}\frac{{\partial T}}{{\partial r}}} \right). \end{split}$

综合上述讨论, 得到了由(2)式、(6)式、(16)式、(17)式、(19)式组成的描述MagLIF演化过程的集成化物理模型. 基于该模型, 自主编写实现一维集成化数值模拟程序MIST. 为了检验建模和程序的正确性, 本节首先对程序中的MHD模块进行了单独校验; 然后与圣地亚实验室公开发表MagLIF模拟结果进行整体比对, 以检验聚变模块的合理性.
MHD模块的校验使用磁驱动铝套筒内爆实验的测量结果, 分别进行了单层套筒(0607发次)和双层套筒(052发次)实验结果的比对, 该实验在中国工程物理研究院脉冲功率驱动装置FP-1上进行, 驱动电流峰值约2 MA, 上升时间约7.2 μs, 两个发次的实验电流曲线如图4所示.
图 4 FP-1装置0523与0607发次实验驱动电流测量曲线
Figure4. Experimental current curves of shot 0523 & 0607.

0607发次实验负载采用单层铝套筒, 外半径为15.5 mm, 厚度为0.5 mm; 0523发次实验采用双层铝套筒, 外层套筒外半径为15.5 mm, 厚度为0.5 mm, 内层套筒外半径为11.2 mm, 厚度为0.7 mm. 使用上述参数作为MIST程序的输入条件, 计算得到套筒自由面速度与测量结果的比较如图5所示.
图 5 MIST计算得到的自由面速度曲线与实验测量结果的比较　(a) 0607发次; (b) 0523发次
Figure5. Comparison of inner surface velocity curves between the calculations and measurements: (a) 0607 shot; (b) 0523 shot.

由于实验中测量电流使用的磁场线圈安装位置距离套筒负载仍有一定距离, 且高功率脉冲驱动器汇流区存在额外的电流损失机制^[24-25], 因此为表征测量结果与真实流经负载的电流差异, MHD计算中通常对实验电流进行一个系数修正, 上述计算中该修正系数分别是0.82 (0607发次)和0.76 (0523发次). 又考虑到实验中测量设备之间触发信号的抖动可能导致所测量结果的零点不一定完全对准, 因此计算结果对速度曲线从整体上在时间轴向右进行了微小的平移, 对0607发次右移0.38 μs, 占加载时间比例为1.9%, 对0523发次右移0.23 μs, 占加载时间比例为1.9%. 图5的对比结果显示, MIST程序计算与实验测量结果吻合一致, 平均误差小于1%, 验证了MHD模块的正确性.
将聚变模块纳入后, 进行了MagLIF集成化测试算例的计算, 参数设置如下: 金属铝套筒, 初始半径为3 mm, 位型比(定义为套筒半径与厚度的比值)AR = 10, 高度为5 mm, DT燃料初始密度为3 mg/cm³, 初始轴向磁场为30 T, 预加热温度为250 eV, 加载电流为30 MA, 上升时间为120 ns, 网格数为60. 在不考虑端面效应和Nernst效应的情况下, 计算结果(如图6所示)表明迟滞时刻燃料温度为8.4 keV, 聚变产额约为550 kJ, 能量增益为1.6, 压缩比CR约为14, 峰值压力为3.6 Gbar (1 bar = 10⁵ Pa), 密度为0.6 g/cm³, 压缩磁场峰值为2500 T, 峰值内爆速度为75 km/s.
图 6 MIST程序计算得到的(a)监测点、(b)套筒内爆速度、(c)燃料温度, 以及(d)聚变产额等随时间的演化曲线
Figure6. Calculated results of (a) grid position, (b) implosion velocity, (c) fuel temperature, and (d) fusion yield evolving with time.

端面效应与Nernst效应会分别导致燃料内部的质量、内能和磁通发生额外的损失, 考虑两种效应后所计算得到的燃料质量和套筒内磁通随时间的变化如图7所示. 从30 ns左右随着激光预加热的加载, 端面效应开始发挥作用, 套筒内燃料质量不断从顶部喷射损失, 到迟滞时刻125 ns左右, 仅有58.5%的初始质量残留; 此时受Nernst效应与普通磁扩散的综合影响, 套筒内相应的磁通保有量锐减至64%. 计算结果表明端面效应与Nernst效应对于MagLIF内爆结果影响很大, 即使在一维模拟中也应尽量给予考虑, 迟滞时刻套筒内燃料质量与磁通的残余量也与文献[18]中给出的结果相符合.
图 7 端面效应与Nernst效应影响下, 套筒内(a)燃料质量与(b)磁通随时间的演化
Figure7. (a) Fuel mass and (b) magnetic flux evolving with time with consideration of end loss and Nernst effect.

为进一步校验程序, 我们选择圣地亚实验室报告^[26,27]中所使用的一维计算结果, 采用同样的初始参数和驱动电流开展计算, 与LASNEX和HYDRA程序的计算结果进行比对. 国外程序所使用的计算参数为: 金属铍套筒, LASNEX程序中初始半径为3.24 mm、HYDRA程序中初始半径为3.1 mm, 位型比AR = 6, 高度为5 mm, DT燃料初始密度为3 mg/cm³, 初始轴向磁场为30 T, 预加热温度为250 eV, 加载电流为27 MA, 上升时间约120 ns.
MIST程序与LASNEX和HYDRA程序的一维计算结果对比如表2所列, 可以发现MIST程序计算数值与国外程序计算结果总体接近. 计算结果的对比中, 有直接联系的燃料温度和聚变产额等都普遍大于国外程序计算结果, 而压缩比与磁场强度也一致地小于国外程序计算结果, 这说明我们程序的计算结果在物理上是自洽的. 导致计算结果差异的主要原因就在于聚变发生后燃料温度计算值偏高, 造成这种情况的原因可能有两点: 1)程序未考虑因聚变反应本身造成的质量损失, 聚变反应率的计算有所高估; 2)程序所使用的材料状态方程库在聚变温度范围内数据点比较稀疏, 调用时采用线性差值方式得到的燃料温度值可能偏高, 导致聚变反应率相应增大, 而聚变放能反过来又会继续推升温度, 使得温度计算的差异被放大.

程序名称	燃料密度/g·cm–3	燃料温度/keV	磁场强度/103 T	压缩比	峰值压力/Gbar	聚变产额/kJ
LASNEX	0.5	8	6—13	23	3	500
MIST	0.47	8.5	7	16	2.7	620
HYDRA	0.8—1.0	6—8	8—22	22	5	565
MIST	0.56	9.5	8	17	3.3	725

表2MIST与LASNEX和HYDRA程序一维计算结果的对比
Table2.Comparison of calculated results between MIST and LASNEX, HYDRA.

MagLIF具备科学、技术和工程3个方面的应用前景, 不仅有望为未来商业能源提供低成本的聚变实现方案, 其本身也是多物理场、多学科交叉的集成化创新典范, 且能够拓展现有大型脉冲功率驱动装置应用范畴. 然而, 国内对于MagLIF领域的相关研究尚处在起步阶段, 集成化数值模拟能力亟需进一步发展和提升.
本文从基本磁流体力学方程组出发, 建立了能够综合考虑磁化、预加热、套筒内爆、聚变反应、端面效应、磁通压缩等多种复杂机制在内的集成化物理模型, 编写实现了一维集成化MagLIF数值模拟程序MIST. 物理模型中对流体喷射模型和Nernst扩散项的引入, 使得程序能够计算得到燃料质量与磁通的损失过程. 与磁驱动铝套筒内爆实验结果的对比, 验证了程序中磁流体模块的正确性, 而聚变模块与国外同类程序LASNEX和HYDRA计算结果的对比说明程序对MagLIF的模拟总体合理, 主要差异体现在燃料温度的计算上, 并对可能的原因进行了简要分析.
本文的工作为开展MagLIF关键问题数值模拟研究搭好了框架, 但是其中的部分模型细节尚未十分完善, 例如端面效应中的流体喷射模型, 并未考虑流体性质、LEH的几何形状甚至是套筒对燃料压缩过程等因素的影响, 仍需进一步进行修正和检验. 此外, 对影响MagLIF聚变过程的关键物理问题如质量和内能损失、Nernst效应的影响、预加热机制、α粒子输运等也需进一步深入的物理分析, 这些方向将是我们下一步工作的主要目标.