预加热效果对磁化套筒惯性聚变放能影响的模拟研究

赵海龙; 王强; 阚明先; 谢龙

doi:10.11858/gywlxb.20200661

预加热效果对磁化套筒惯性聚变放能影响的模拟研究

doi: 10.11858/gywlxb.20200661

中国工程物理研究院流体物理研究所，四川绵阳　621999

基金项目: 国家自然科学基金（11205145，21805262，12075226）

详细信息

作者简介:
赵海龙（1985－），男，博士研究生，副研究员，主要从事磁流体理论与数值模拟研究. E-mail：ifp.zhaohailong@qq.com

中图分类号: O539
计量
- 文章访问数: 3400
- HTML全文浏览量: 1471
- PDF下载量: 21
出版历程
- 收稿日期: 2020-12-28
- 修回日期: 2021-01-21

Simulation of the Preheating Effects on the Discharging of Magnetized Liner Inertial Fusion

Institute of Fluid Physics, CAEP, Mianyang 621999, Sichuan, China

摘要

摘要: 得益于预加热和轴向磁场的作用，理论上磁化套筒惯性聚变（MagLIF）构型能有效降低聚变实现的难度，具有极大的应用潜力。然而，当前研究过于重视对激光能量沉积效率的提升，而忽略了预加热自身参数对MagLIF过程和内爆结果的影响。为此，采用一维集成化数值模拟程序MIST，开展了MagLIF过程中预加热效果对聚变放能影响的模拟研究，基于参数扫描方法，从简单的模型着手，逐步深入探讨相关参数对内爆结果的影响。模拟结果表明：预加热是MagLIF构型能够成功的必要条件，最佳时间是套筒即将开始向内压缩燃料的时刻；燃料预加热的设计原则是让燃料获得尽可能平缓分布的高温，而中心局部加热方式对于未能达到点火条件的负载更有优势；激光预加热模式下，脉宽越短越好，对于以ZR装置驱动能力为目标的算例而言，最佳套筒高度为1.0 cm。研究结果有助于加深对MagLIF过程中预加热机制和效果的认知和理解，对于具体的负载参数设计也有较强的指导意义。
- 磁化套筒惯性聚变 /
- 预加热 /
- MIST /
- 内爆
Abstract: Benefiting from laser preheating and axial magnetization, magnetized liner inertial fusion (MagLIF) has great application potential because it can effectively reduce the difficulties to realize the controlled fusion in theory. However, much attention has been paid to the improvement of laser energy deposition efficiency in current research, while the influence of preheating parameters on the MagLIF process and implosion is ignored. For this reason, the one-dimensional integrated simulation code, MIST, is used here to study the preheating effect on the fusion discharging in MagLIF process. Based on the method of parameter scanning, starting from a simple model, the studies of the influences of relevant parameters on implosion results are gradually advanced. The simulation results show that preheating is a necessary condition for the success of MagLIF configuration, and the best preheating time is the moment when the liner is about to compress the fuel. The design principle of preheating is to allow the fuel to acquire as smoothly distributed high temperature as possible, and the central local heating mode is more advantageous when the ignition fails. If the laser is preheated, the shorter pulse width will be better. For driving ability of ZR facility, the optimal liner height is 1.0 cm. These results are not only helpful to understand laser preheating mechanism and effect during MagLIF process, but also providing useful guidance for the design of the detail load parameters of MagLIF configuration.
- Magnetized Liner Inertial Fusion /
- preheat influences /
- MIST /
- implosion

HTML全文

地铁站通常具有人流量大^[1]、结构封闭等特点，一旦发生爆炸，将会造成巨大的人员伤亡和财产损失，爆炸产生的能量还会影响地铁站结构的安全性，因此研究地铁站台及围岩爆炸响应规律对于地铁站的安全运营具有重要意义。

针对地下空间爆炸问题开展实验研究，成本高且安全风险大，因此国内外学者一般采用有限元方法进行模拟研究。例如：王勇等^[2]基于LS-DYNA3D算法，对某地铁区间隧道爆炸响应情况进行数值模拟，验证了隧道结构能抵抗10 kg TNT炸药的爆炸；Choi等^[3]对隧道内外不同位置起爆的响应情况进行了数值模拟，为坑道易损性评估提供了一种简单实用的方法；Zhang等^[4]研究了典型矿山密封在冲击载荷作用下的力学响应性能；柴永生等^[5]模拟了爆炸冲击波在地铁出入口的传播规律，为地铁及其他地下工程出入口部位的技术标准制定和防护装备设施设计提供了参考。有限元方法在模拟由爆炸引起的岩土体大变形问题时，常采用输入应力波的方式模拟爆炸应力，与气体爆炸冲击岩土体的实际情况存在较大差异，此外，岩土体大变形还会造成有限元网格畸变^[6]。物质点法是采用欧拉法与拉格朗日法双重描述的无网格粒子法，将连续体离散成一组带有质量的质点，并在背景网格节点上求解动量方程，计算空间导数。物质点法兼具拉格朗日法和欧拉法的双重优势，能够避免网格畸变和对流项处理造成的误差和计算量增大问题^[7]，对物质交界面无需特殊处理，自动满足无滑移的接触条件^[8]，能够直接模拟物体之间的冲击现象。这些特点使物质点法在处理爆炸冲击^{[6, 9]}和岩土大变形^[10-11]问题上较有限元方法更具优势。

由于地铁站结构多为混凝土结构，需要考虑爆炸作用下混凝土孔隙坍缩现象，因此在数值模拟中应选择适宜的材料模型。Holmquist、Johnson和Cook提出了用于模拟混凝土材料的HJC模型^[7]，该模型能够描述爆炸作用对地下基础设施的损伤区域大小和损伤强度，为此，本研究将HJC模型嵌入物质点法程序中，以模拟站台结构响应情况。

本研究依托上海某地铁站工程，基于张雄教授公开的开源版物质点法MPM3D-F90程序^[7]，将HJC模型嵌入程序，模拟地铁站内固体炸药爆炸冲击结构和围岩过程，探究爆炸作用下地铁站台及围岩的响应规律，以期为合理应对地铁爆炸灾害以及灾后人员救援、结构抢修和加固提供一定的指导。

1. 材料本构模型

1.1 材料本构模型基本原理

针对站台结构和大厅底板等混凝土结构，采用HJC本构模型进行模拟。该模型考虑了混凝土材料的受损情况和高压下混凝土孔隙坍缩现象，能够描述爆炸作用下地下基础设施的损伤区域大小及损伤强度^[7]。

HJC 模型的极限面被描述为损伤、应变率与静水压力的函数，其破坏失效面可以表示为

$\mathop \sigma \nolimits^ * = [A(1 - D) + B\mathop p\nolimits^{ * N} ](1 + C\ln {\dot \varepsilon ^ * })$

(1)

式中： $\mathop \sigma \nolimits^ *$ 为归一化等效强度， $\mathop \sigma \nolimits^ * = \sigma /{f_{\rm c}}$ ， $\sigma$ 为真实应力， ${f_{\rm c}}$ 为单轴抗压强度，且满足 $\mathop \sigma \nolimits^ * \leqslant {S_{\max }}$ ， ${S_{\max }}$ 为归一化最大等效屈服强度； $\mathop p\nolimits^ *$ 为归一化压力， $\mathop p\nolimits^ * = p/{f_{\rm c}}$ ，p为真实压力； ${\dot \varepsilon ^ * } = \dot \varepsilon /{\dot \varepsilon _0}$ ， ${\dot \varepsilon ^ * }$ 为等效应变率， $\dot \varepsilon$ 为真实应变率， ${\dot \varepsilon _0}$ 为参考应变率；D为损伤因子（ $0 \leqslant D \leqslant 1$ ）；A为归一化黏聚强度；B为归一化压力硬化系数；N为压力硬化指数；C为应变率系数。

损伤因子 D 由等效塑性应变和体积应变累加得到，即损伤演化方程可表示为

$D={\displaystyle \sum (\Delta {\varepsilon }_{\rm p}+\Delta {\mu }_{\rm p})}/({\varepsilon }_{\rm p}^{\rm f}+{\mu }_{\rm p}^{\rm f})$

(2)

式中： $\Delta {\varepsilon _{\rm p}}$ 为等效塑性应变增量， $\Delta {\mu _{\rm p}}$ 为塑性体积应变增量， $\varepsilon _{\rm p}^{\rm f}$ 和 $\mu _{\rm p}^{\rm f}$ 分别为常压下材料破碎时的塑性应变和塑性体积应变。多数情况下，HJC模型的损伤主要由等效塑性应变引起。

针对围岩结构，采用适用于岩石材料的Drucker-Prager（D-P）模型进行模拟，其屈服函数为

${f_{\rm s}} = \alpha \frac{{{I_1}}}{3} + \sqrt {{J_2}} - k$

(3)

式中：I₁为应力张量第一不变量，J₂为偏应力张量第二不变量，参数 $\alpha$ 和k可由岩石黏聚力c和内摩擦角 $\phi$ 确定。若取屈服面为Mohr-Coulomb准则外接圆，则 $\alpha$ 和k可由下式计算

$a = \frac{{6\sin \phi }}{{\sqrt 3 \left( {3 - \sin \phi } \right)}}, \;\;\;\;\;\; k = \frac{{6c\cos \phi }}{{\sqrt 3 \left( {3 - \sin \phi } \right)}}$

(4)

1.2 应力更新

更新后的应力 ${\sigma '_{ij}}$ 可以由应力率 ${\dot \sigma _{ij}}$ 获得，即

${\sigma '_{ij}} = {\sigma _{ij}} + \Delta t \; {\dot \sigma _{ij}}$

(5)

式中： $\Delta t$ 为时间增量。在大变形物体的本构关系中，为了确保应力率不受刚体旋转的影响，需要将其转化为在与物体固定连接的动态坐标系中的焦曼应力率 $\dot \sigma _{ij}^{\rm J}$

$\dot \sigma _{ij}^{\rm J} = \dot \sigma _{ij}^{} - {\sigma _{ik}}{\varOmega _{jk}} - {\varOmega _{ik}}{\sigma _{kj}}$

(6)

式中： ${\varOmega _{ij}}$ 为旋率张量分量。将式(6)代入式(5)，得到更新后的应力

${\sigma '_{ij}} = {\sigma _{ij}} + \Delta t\left( {\dot \sigma _{ij}^{\rm J} + {\sigma _{ik}}{\varOmega _{jk}} + {\varOmega _{ik}}{\sigma _{kj}}} \right)$

(7)

1.3 接触算法

本研究中多个物体之间的接触算法与文献[6]相同，受篇幅限制，这里不再赘述。

1.4 HJC模型验证

参照Hanchak等^[12]的实验结果，模拟子弹穿越无筋混凝土板过程，记录子弹以不同初速度穿越混凝土板后的剩余速度，并将结果与实验数据进行对比。模拟中混凝土材料采用HJC模型，子弹采用理想弹性模型，具体参数与文献[12]一致。弹体初速度为1058 m/s时的侵彻结果如图1所示，从中可以清晰地看到，在子弹的冲击作用下混凝土质点受损因子呈辐射状分布。图2显示了子弹剩余速度的对比结果，可见物质点法模拟得到的数据与实验结果吻合较好，模拟结果正确可靠。

图 1 混凝土受损因子云图

Figure 1. Nephogram of concrete damage factor

下载: 全尺寸图片幻灯片

图 2 子弹剩余速度对比

Figure 2. Comparison of residual velocity of bullet

下载: 全尺寸图片幻灯片

2. 计算模型设置

2.1 物质点法模型建立

上海某地铁站站台全长271.7 m（其中有柱区域长205.0 m，其余为无柱区域），宽19.1 m，层高6.15 m。站台有柱区域均匀分布52根方柱，站台与上部车站大厅之间的围岩厚度为5 m，方柱贯穿此围岩，将站台与车站大厅相连。站台本体结构材料为钢筋混凝土，厚度为0.9 m，标号为C40，内部结构的BIM模型如图3所示。

图 3 站台内部结构的BIM模型

Figure 3. BIM model of internal structure of the platform

下载: 全尺寸图片幻灯片

约定模型中与站台列车运行方向平行的方向为纵向，与运行方向垂直的方向为横向。同时，为模拟车站上部的建筑物和人员荷载，设置站台埋深为5 m，并在模型顶部添加90 kPa竖向荷载，模拟站台顶部建筑荷载。由于爆炸作用范围有限，为减小数值模拟的计算量，本研究按照无柱区域与有柱区域分别建立站台模型。无柱区域模型长度为99 m（5.2倍横向尺寸），如图4所示，模型整体尺寸为99 m×61 m×40 m。该地铁站台结构柱的标准间距为横向4.900 m，纵向7.125 m，距边墙6.100 m。结构柱将站台顶板和底板与大厅底板连接，据此建立站台净空内结构柱模型，如图5所示。

图 4 站台无柱区域模型

Figure 4. Model of the non-pillar area

下载: 全尺寸图片幻灯片

图 5 站台结构柱模型平面图

Figure 5. Plan view of structural column model of the platform

下载: 全尺寸图片幻灯片

2.2 模型参数选取

模型四周采用无反射边界模拟侧向无限围岩，顶部采用自由边界模拟车站大厅地面。站台底部围岩在水平方向上均匀分布，在竖直方向上，上部围岩的强度较低，而底部50 m处围岩为中低压缩性砂土，土质密实，故底部采用对称边界模拟其界面。站台本体结构和上部车站大厅采用混凝土HJC强度准则，除最大静水压力外（站台本体结构为4 MPa，车站大厅底板为2 MPa），其他材料参数均相同，如表1所示，其中： $\,\rho$ 为密度，E为弹性模量，ν为泊松比， ${f'_{\rm c}}$ 为无侧限抗压强度， $\varepsilon _{\min }^{\rm f}$ 为最小破碎塑性应变，p_crush为孔洞坍塌压力，p_lock为压实压力，D₁、D₂、K₁、K₂、K₃为材料参数，v_p为纵波波速，v_s为横波波速。围岩为淤泥质砂土，采用D-P模型描述，物理力学参数如表2所示，其中： $q_{\phi}$ 和 $K_{\phi}$ 为模型参数， $q_{\varPsi}$ 为剪胀角， $\sigma_{\rm{t}}$ 为岩石抗拉强度。内部空气采用无偏应力空气模型模拟，其状态方程采用理想气体状态方程，物理力学参数如表3所示，其中：c为波速，E₀为初始内能， $\kappa$ 为比热容比。

表 1 混凝土结构的HJC模型参数

Table 1. HJC model parameters of concrete structure

$\,\rho$ /(kg·m^–3)	E/GPa	$\nu$	${ f{'_ {\rm{c} } } }$ /MPa	S_max	A	B	N	C	$\varepsilon {_{\min }^{\rm f} }$
2439	32.5	0.2	48	7.0	0.79	1.6	0.61	0.007	0.01

p_crush/GPa	p_lock/GPa	D₁	D₂	K₁/GPa	K₂/GPa	K₃/GPa	v_p/(m·s⁻¹)	v_s/(m·s⁻¹)
0.016	0.8	0.04	1	85	−171	208	3872	2375

下载: 导出CSV

| 显示表格

表 2 围岩的物理力学参数

Table 2. Physical and mechanical parameters of surrounding rock

$\,\rho$ /(kg·m^–3)	E/GPa	$\nu$	$q{_{\phi}}$ /(°)	$K{_{\phi}}$	$q{_{\varPsi} }$ /(°)	$\sigma$ _t/kPa	v_p/(m·s⁻¹)	v_s/(m·s⁻¹)
1850	0.04	0.35	0.388	11171	0	0.1	185	89

下载: 导出CSV

| 显示表格

表 3 空气（空模型）的物理力学参数

Table 3. Physical and mechanical parameters of air (air model)

$\,\rho$ /(kg·m^–3)	c/(m·s⁻¹)	E₀/(MJ·m^–3)	$\kappa$	v_p/(m·s⁻¹)	v_s/(m·s⁻¹)
1.29	340	0	1.4	340	0

下载: 导出CSV

| 显示表格

2.3 固体爆炸模型建立

JWL状态方程已广泛应用于以TNT炸药为代表的固体炸药爆炸过程模拟^[13]，且有大量的实验参数可供参考。JWL状态方程为

$p = {A_{\rm {JWL}}}\left(1 - \frac{\omega }{{{R_1}V}}\right){{\rm e}^{ - {R_1}V}} + {B_{\rm {JWL}}}\left(1 - \frac{\omega }{{{R_2}V}}\right){{\rm e}^{ - {R_2}V}} + \frac{{\omega e}_0}{V}$

(8)

式中：e₀为单位体积炸药的初始内能， $\omega$ 、A_JWL、B_JWL、R₁和R₂为JWL状态方程参数，取值见表4，其中： $\,\rho_0$ 为初始密度，p_CJ为CJ爆轰压力， $\gamma$ 为爆轰产物的绝热指数，D_J为爆速。

表 4 固体爆炸物参数

Table 4. Calculation parameters of solid explosives

$\,\rho$ ₀/(kg·m^–3)	e₀/(GJ·m^–3)	p_CJ/GPa	$\gamma$	D_J/(m·s^–1)
1500	7.0	21	2.727	6930

A_JWL/GPa	B_JWL/GPa	R₁	R₂	$\omega$
371.2	3.23	4.15	0.95	0.30

下载: 导出CSV

| 显示表格

本研究参考李志鹏^[14]在ALE模拟中所用体积填充法的装药网格数目，选用的基础模型中爆炸物由8个质点构成，每个质点的质量为50 kg，替换相应体积的空气，炸药爆炸的TNT当量为444.4 kg^[15]。爆炸物位于站台底板处，以模拟恐怖袭击中常见的固体爆炸装置置于站台层底板的爆炸情况。

3. 爆炸作用下站台及围岩响应规律分析

3.1 爆炸作用下站台及围岩的响应压强

图6为爆炸发生后4和400 ms时站台无柱区域顶板和底板的压强分布云图。随着爆炸后时间的增加，顶板和底板的响应压强峰值均逐渐减小。顶板压强在爆炸后出现应力集中现象，这是由起爆后爆炸物粒子撞击结构导致的，撞击点偏离了起爆点在结构顶板上的投影位置，当爆炸物被继续细化时，这种偏离距离会消失。响应压强以应力波的形式从起爆点向外辐射传播，使得站台结构内部既存在正压强也出现负压强，即站台底板和顶板既存在受压区也存在受拉区。400 ms时顶板中心处出现正压强，这是由爆炸应力波向上传递到大厅底板后反射回来的反射波导致的；而底板中心处正压强与负压强分布不规则是由于在爆炸作用下底部围岩出现塌陷坑，由底部地质界面反射回来的应力波又发生多次反射，使得结构底板压强分布不规律。

图 6 爆炸作用下不同时刻结构顶板和底板的压强分布

Figure 6. Pressure distribution of the top and bottom plates of the structure at various times under the explosion loading

下载: 全尺寸图片幻灯片

图7为站台结构顶板中心质点压强随时间变化曲线。爆炸应力波在起爆后1.6 ms时到达顶板，4.5 ms时顶板质点压强出现峰值，为6.97 MPa，此后压强峰值和振幅快速减小，这是由于应力波在围岩-站台结构和站台结构-空气两个界面之间不断传播和反射，能量快速损耗，响应压强迅速减小。

图 7 站台结构顶板中心点处压强随时间变化曲线及局部放大图

Figure 7. Pressure curve at the center point of the top plate with time and its partial enlarged view

下载: 全尺寸图片幻灯片

图8为站台结构底板中心质点压强随时间变化曲线。从图8可以看出，爆炸应力波在起爆后0.2 ms时到达顶板， 7.0 ms时达到峰值200 MPa，随后迅速减小并在0 kPa附近振荡，这是由于起爆后底板质点区域下方的围岩形成塌陷坑，使该质点处于悬空状态，不受反射波影响，因此响应压强迅速减小。底板质点先于顶板质点1.1 ms受到爆炸作用，结合站台层高，可以求得爆炸应力波的传播速度为4615 m/s，小于文献[16]中描述的5800～7000 m/s，主要原因是离散的空气质点减缓了传播速度。

图 8 站台结构底板中心点处压强随时间变化曲线及局部放大图

Figure 8. Pressure curve at the center point of the bottom plate with time and its partial enlarged view

下载: 全尺寸图片幻灯片

3.2 爆炸作用下站台净空内爆炸超压分析

根据距地面3.0 m高的监测数据，采用Kriging插值法，计算得到爆炸过程中起爆点附近的超压峰值分布，如图9所示。可以看出，爆炸产生的超压峰值出现在纵向上偏离站台中轴线2.5 m处，即出现超压峰值偏移，这是由于爆炸物位于站台底部，起爆后冲击波波阵面呈半球形，与3.0 m高平面的交点与起爆点有一定的偏移距离。图9中爆炸超压峰值为6.5 MPa，且在起爆点横向上出现超压峰值突变情况，这是由于初始爆轰应力波传递至隧道边墙后产生反射，与后续爆轰反应形成的应力波叠加，致使横向超压峰值发生突变。在以往研究中，站台净空爆炸超压峰值分布计算结果未出现这样的突变区域，其原因在于他们仅考虑了结构净空内压力分布，将结构-围岩界面简化为无反射边界，忽略了应力波由结构传递至围岩时发生的折射和反射现象^[17]。

图 9 爆炸超压分布

Figure 9. Distribution of explosion overpressure

下载: 全尺寸图片幻灯片

Henrych^[18]、Brode^[19]、Sadovskyi^[20]等通过理论推导提出了不同适用范围超压峰值p_s的经验公式。Wu等^[21]基于前人的研究，提出了更加精确的超压峰值公式

${p_{\rm {s}}} = \left\{ \begin{gathered} \frac{{{1.059}}}{{{{\bar R}^{2.56}}}} - {0.051}\quad{\text{ 0}}{\text{.1 }}\leqslant\bar R \leqslant {1.0} \hfill \\ \frac{{{1.008}}}{{{{\bar R}^{2.01}}}}\quad\quad\quad\quad\;{1.0} < \bar R \leqslant {10.0} \hfill \\ \end{gathered} \right.$

(9)

式中：p_s的单位为MPa； $\bar R$ 为比例距离， $\bar R = R/{W^{1/3}}$ ，其中W为装药量（单位kg），R为距起爆点距离（单位m）。表5列出了本研究模拟得到的5和10 m范围内的超压峰值与式(9)的计算结果，可见二者相近，进一步验证了本研究提出的数值模拟方法的准确性。

表 5 超压模拟结果与经验公式计算结果对比

Table 5. Comparison of overpressure simulation results and empirical formula calculation results

Distance/m	Overpressure/MPa		Error/%
Distance/m	Theoretical formula	Numerical simulation	Error/%
5	2.35	2.25	4.44
10	0.47	0.52	−7.69

下载: 导出CSV

| 显示表格

为评价爆炸对不同范围内人员的损伤情况，基于师光达^[22]提出的爆炸对人员伤亡损失评价标准和图9，建立了人员伤亡区域云图，如图10所示。可见，人员致死区域是半径为8 m、面积为200 m²的蝴蝶形区域，实际爆炸造成的伤亡区域面积远大于此范围，这是由于爆炸过程中产生的飞散弹片以及后期救治不及时造成的^[23]。

图 10 人员伤亡区域预测分布

Figure 10. Regional distribution of casualties

下载: 全尺寸图片幻灯片

3.3 爆炸作用下站台结构及围岩位移

图11为起爆后400 ms时围岩和车站结构位移分布云图。从图11中可以看出，爆炸作用使站台结构发生整体沉降，最大沉降位移为2.57 m，起爆点下方围岩出现深度为1.64 m的凹陷深坑。站台上下覆围岩均与大厅底板或站台底板分离。图12为爆炸发生400 ms时站台结构顶板和底板竖向位移分布云图。可见，横向上顶板位移分布较一致，纵向上则随着与起爆点距离的增加而逐渐减小；底板竖向位移在横向和纵向上均随距起爆点距离的增加而减小，起爆点附近的底板位移急剧增大，并且随着时间的增加，该处质点位移继续增大，该区域发生明显破裂。以下结合爆炸后应力波的传递情况，对站台结构和围岩位移进行分析。

图 11 400 ms时横向和竖向上站台及围岩的位移分布

Figure 11. Displacement distributions of platform and surrounding rock in horizontal and vertical directions at 400 ms

下载: 全尺寸图片幻灯片

图 12 站台结构顶板和底板的竖向位移分布

Figure 12. Vertical displacement distributions of the top and bottom plates of the platform structure

下载: 全尺寸图片幻灯片

爆炸物起爆后，空气冲击波波阵面扩散到站台本体结构，此时本体结构受力发生振动并向围岩传递应力波。应力波产生的原因有两个：一是飞散的爆炸产物直接撞击本体结构，二是气流冲击波冲击至本体结构。爆炸应力波向下透过站台底板传至围岩，使其产生大量向下的位移，由于围岩之间的黏结力，站台侧向围岩也向下运动，因此站台结构发生整体沉降。与此同时，因站台底板受损，应力波向下传播较为集中，使得起爆点下方围岩产生较大变形，局部凹陷。图13显示了塌陷坑深度随时间的变化规律。0～80 ms内，受初始应力波的作用，塌陷坑迅速下沉；80～320 ms内，反射应力波传递至围岩与站台底板接触面，降低围岩的下降速度，此时塌陷坑深度进入缓慢增长阶段。另外，当站台底板与围岩整体沉降时，由于围岩的弹性模量较小，其沉降位移比站台结构大，因此围岩与站台底板发生分离。

图 13 站台底板下部深坑深度随时间变化曲线

Figure 13. Curve of pit depth under platform bottom plate with time

下载: 全尺寸图片幻灯片

站台上方围岩和大厅底板在向上传播的冲击波作用下相对于周围向上隆起，且大厅底板和围岩发生分离，这是因为冲击波传播至围岩与车站大厅底板接触面时发生反射和折射，反射冲击波将减缓围岩整体向上隆起的速度。取距离起爆点20 m处质点为参考点，分析站台纵向距起爆点5、10、15 m处围岩表面质点的相对位移随时间变化曲线，如图14所示。可以看出，各质点的相对位移随时间的增加呈现相似的平滑增加规律。首先0～350 ms内，在爆炸冲击波及其次生波的作用下，围岩加速隆起；在350～400 ms内，由于界面之间的反射和黏弹性边界的吸收，应力波能量全部耗散，围岩的位移也趋于稳定。

图 14 围岩表面质点的相对位移随时间变化曲线

Figure 14. Curves of relative displacement of surrounding rock surface particles with time

下载: 全尺寸图片幻灯片

3.4 站台本体结构损伤因子

站台结构的受损部位主要集中于靠近爆炸物的车站本体结构底板，为此，本研究着重分析结构底板受损因子的分布规律。

图15为爆炸产生的本体结构受损因子分布云图。可以看出，受损区域呈椭圆形，其中椭圆长轴平行于站台纵向，这是由于站台横向尺寸小于纵向尺寸，竖向墙壁结构增加了横向上站台的稳定性，使其受损情况较纵向上的站台轻。图15中出现了集中损伤区域，这是由于爆炸产物直接冲击本体结构造成的，而均布损伤区域是由空气冲击波冲击产生的。此外，当质点受损因子达到1时会被判定为失效，此时底板有效承重区域分布不规则，致使后续质点受损因子的计算结果出现不均质性。

图 15 站台受损因子云图

Figure 15. Nephogram of the damage factor of the platform floor

下载: 全尺寸图片幻灯片

在站台横纵截面受损因子分布云图中，横纵向受损因子均集中在固体爆炸物正下方的结构底板，且在截面上损伤因子上小下大，底板上表面的受损程度小于底板下表面的受损程度，此外，由于结构底板内应力波在上下表面内反复反射叠加，导致中部区域受损最严重。通过底板截面受损区域的梯形分布，可以看出应力波在结构内的传播趋势。

3.5 站台有柱区域与无柱区域的爆炸响应对比

表6列出了站台有柱区域和无柱区域在固体爆炸物爆炸作用下结构围岩、净空超压等各项指标的对比，其中：p_m为净空超压峰值，R_d为净空超压致死区半径，n_d为结构受损因子峰值，C_d为结构受损面积，p_r、p_f1、p_f2分别为结构顶板、结构底板和大厅底板的响应压强峰值，d_r、d_f1、d_f2分别为结构顶板、结构底板、大厅底板的位移峰值，d_p为围岩塌陷坑深度。表6中，站台有柱区域的p_m、R_d、n_d、C_d、p_f1、d_f1、d_p均较无柱区域有所下降，但p_r、d_r、p_f2、d_f2有所增加。这是由于结构柱增加了整个结构体系的刚度，使大厅底板和站台结构顶板分担了炸药对站台结构底板的冲击作用，提升了站台结构的抗爆能力。此外，结构柱的存在增加了站台结构封堵率，削弱了爆炸冲击波的传播，导致超压致死区半径减小。

表 6 站台有柱区域和无柱区域在固体爆炸作用下结构参数的对比

Table 6. Comparison of various parameters under the solid explosion in the pillared and non-pillared areas of the platform

Area	p_m/MPa	R_d/m	n_d	C_d/m²	p_r/MPa	p_f1/MPa
With volumns	6.5	8	1.00	21	5.025	218.0
Without volumns	6.0	5	0.67	13	5.667	110.5
Error/%	−7.69	−37.50	−33.00	−38.10	12.78	−49.31

Area	p_f2/MPa	d_r/m	d_f1/m	d_f2/MPa	d_p/m
With volumns	0.034	0.124	0.188	0.037	1.64
Without volumns	0.041	0.084	0.103	0.089	0.43
Error/%	20.59	−32.26	−45.21	140.54	−73.78

下载: 导出CSV

| 显示表格

通过各项数据的对比分析可以看出，站台有柱区域结构整体性优于无柱区域，抗爆能力更强。

4. 结　论

将HJC模型嵌入开源版物质点法程序，研究了爆炸荷载作用下地铁站台及围岩的压强、位移、净空超压和结构受损的响应规律，得到如下主要结论。

(1) 结构在爆炸应力波作用下既存在受拉区又存在受压区；应力波在围岩-站台结构和站台结构-空气两个界面之间不断传播和反射，能量逐渐损耗，压强逐渐减小；空间内超压分布存在峰值偏移，同时在站台边墙处应力波与反射波叠加出现超压突变区；爆炸超压致死区域半径约为8 m。

(2) 结构体系在爆炸作用下会发生整体沉降，起爆点上方围岩和大厅地板相对周围向上隆起，起爆点正下方围岩出现塌陷坑，上下覆围岩短时间内会出现与站台结构体脱离的情况。

(3) 结构受损区域主要集中在结构底板，受损区域呈椭圆形，底板中部受损程度最严重，上表面的受损程度小于下表面的受损程度，站台有柱区域的抗爆能力强于无柱区域。

图 1 ZR装置在95 kV充电电压下的驱动电流随时间演化曲线^[7]

Figure 1. Driving current from ZR facility with charging voltage 95 kV^[7]

下载: 全尺寸图片幻灯片

图 2 250 eV预加热温度与无预加热条件下计算得到的燃料内能和聚变产额随时间演化曲线

Figure 2. Demonstrations of fusion yield and fuel internal energy calculated with 250 eV and no preheat

下载: 全尺寸图片幻灯片

图 3 不同初始预加热温度下计算得到的内爆结果对比

Figure 3. Demonstrations of implosion results calculated with different preheat temperature

下载: 全尺寸图片幻灯片

图 4 余弦预加热方式计算得到的滞止时刻燃料温度和密度分布曲线

Figure 4. Distributions of fuel temperature and density at stagnation time with cosine preheat

下载: 全尺寸图片幻灯片

图 5 余弦预加热方式计算得到的聚变产额与燃料内能演化曲线

Figure 5. Demonstrations of fusion product and internal energy calculated with cosine preheat

下载: 全尺寸图片幻灯片

图 6 不同沉积半径下预加热和迟滞阶段的温度分布曲线

Figure 6. Distributions of preheat and stagnation temperatures with different preheat radius

下载: 全尺寸图片幻灯片

图 7 不同沉积半径下迟滞阶段的温度分布和聚变产额演化曲线（B_z = 5 T）

Figure 7. Distributionsof stagnation temperature and evolvement of fusion product with different preheat radii (B_z = 5 T)

下载: 全尺寸图片幻灯片

图 8 不同沉积半径下预加热时刻和迟滞阶段温度分布（B_z = 5 T, E_las = 3 kJ）

Figure 8. Distributions of preheat and stagnation temperature with different preheat radii (B_z = 5 T, E_las = 3 kJ)

下载: 全尺寸图片幻灯片

图 9 不同沉积半径下磁化强度BR随时间演化曲线(B_z = 5 T, E_las = 3 kJ)

Figure 9. Schematic of BR evolving with time with different preheat radii (B_z = 5 T, E_las = 3 kJ)

下载: 全尺寸图片幻灯片

图 10 不同脉宽时预加热功率曲线

Figure 10. Schematic of preheat power with different durations

下载: 全尺寸图片幻灯片

图 11 不同激光脉宽和功率参数下预加热时燃料中温度和密度分布

Figure 11. Distributions of fuel temperature and density with different laser power and durations at preheat time

下载: 全尺寸图片幻灯片

图 12 简化后的ZR装置的等效电路示意图^[25]

Figure 12. Schematic of simplified equivalent circuit of ZR facility^[25]

下载: 全尺寸图片幻灯片

图 13 ZR装置绝缘堆电压曲线和MIST程序计算得到的负载电流曲线

Figure 13. Voltage curve from the vacuum insulator and calculated current curve by MIST code

下载: 全尺寸图片幻灯片

表 1 不同预加热温度下计算得到的内爆结果对比

Table 1. Calculated implosion results with different preheat temperatures

Preheat temperature/ eV	Preheat energy/ kJ	Fuel temperature/ keV	Internal energy/ (kJ·cm⁻¹)	Fusion yield/ (kJ·cm⁻¹)	Q
50	3.8	4.7	310	510	1.65
100	7.2	7.5	470	1400	2.98
150	10.6	9.4	580	2000	3.45
200	14.0	10.5	650	2300	3.54
250	17.4	10.0	700	2420	3.46
300	20.8	11.6	730	2450	3.36
350	24.0	11.7	740	2390	3.23

下载: 导出CSV

表 2 不同套筒高度计算得到的内爆结果对比

Table 2. Calculated implosion results calculated with different liner heights

h/cm	Preheat temperature/eV	Peak current/MA	Internal energy/(kJ·cm⁻¹)	Fusion yield/(kJ·cm⁻¹)	Total yield/kJ
0.50	890	29.5	786	2426	1213
0.75	615	28.9	668	2133	1600
1.00	450	28.2	565	1614	1614
1.25	364	27.4	478	1172	1465

下载: 导出CSV

参考文献(28)

[1]	AYMAR R. The ITER project [J]. IEEE Transactions on Plasma Science, 1997, 25(6): 1187–1195. doi: 10.1109/27.650895
[2]	SHIMOMURA Y, SPEARS W. Review of the ITER project [J]. IEEE Transactions on Applied Superconductivity, 2004, 14(2): 1369–1375. doi: 10.1109/TASC.2004.830580
[3]	HUANG C J, LI L F. Magnetic confinement fusion: a brief review [J]. Frontiers in Energy, 2018, 12(2): 305–313. doi: 10.1007/s11708-018-0539-1
[4]	HURRICANE O A, SPRINGER P T, PATEL P K, et al. Approaching a burning plasma on the NIF [J]. Physics of Plasmas, 2019, 26(5): 052704. doi: 10.1063/1.5087256
[5]	MCCRORY R L, MEYERHOFER D D, BETTI R, et al. Progress in direct-drive inertial confinement fusion [J]. Physics of Plasmas, 2008, 15(5): 055503. doi: 10.1063/1.2837048
[6]	ROSEN M D. The physics issues that determine inertial confinement fusion target gain and driver requirements: a tutorial [J]. Physics of Plasmas, 1999, 6(4): 1690–1699. doi: 10.1063/1.873427
[7]	SLUTZ S A, HERRMANN M C, VESEY R A, et al. Pulsed-power-driven cylindrical liner implosions of laser preheated fuel magnetized with an axial field [J]. Physics of Plasmas, 2010, 17(5): 056303. doi: 10.1063/1.3333505
[8]	HARVEY-THOMPSON A J, GEISSEL M, JENNINGS C A, et al. Constraining preheat energy deposition in maglif experiments with multi-frame shadowgraphy [J]. Physics of Plasmas, 2019, 26(3): 032707. doi: 10.1063/1.5086044
[9]	PARADELA J, GARCÍA-RUBIO F, SANZ J. Alpha heating enhancement in MagLIF targets: a simple analytic model [J]. Physics of Plasmas, 2019, 26(1): 012705. doi: 10.1063/1.5079519
[10]	PERKINS L J, LOGAN B G, ZIMMERMAN G B, et al. Two-dimensional simulation of thermonuclear burn in ignition-scale inertial confinement fusion targets under compressed axial magnetic fields [J]. Physics of Plasmas, 2013, 20(7): 072708. doi: 10.1063/1.4816813
[11]	SLUTZ S A, VESEY R A. High-gain magnetized inertial fusion [J]. Physical Review Letters, 2012, 108(2): 025003. doi: 10.1103/PhysRevLett.108.025003
[12]	SEFKOW A B, SLUTZ S A, KONING J M, et al. Design of magnetized liner inertial fusion experiments using the Z facility [J]. Physics of Plasmas, 2014, 21(7): 072711. doi: 10.1063/1.4890298
[13]	SLUTZ S A. Magnetized Liner Inertial Fusion (MagLIF): the promise and challenges [C]//Proceedings of MagLIF Workshop. Albuquerque: 2012.
[14]	GOMEZ M R, SLUTZ S A, SEFKOW A B, et al. Experimental demonstration of fusion-relevant conditions in magnetized liner inertial fusion [J]. Physical Review Letters, 2014, 113(15): 155003. doi: 10.1103/PhysRevLett.113.155003
[15]	AWE T J, MCBRIDE R D, JENNINGS C A, et al. Observations of modified three-dimensional instability structure for imploding z-pinch liners that are premagnetized with an axial field [J]. Physical Review Letters, 2013, 111(23): 235005. doi: 10.1103/PhysRevLett.111.235005
[16]	GOMEZ M P, SLUTZ S A, SEFKOW A B, et al. Recent progress in magnetized liner inertial fusion (MagLIF) experiments [R]. Austin: NNSA, 2015.
[17]	SINARS D. Magnetized Liner Inertial Fusion (MagLIF) research at Sandia national laboratories [C]//The 1st Chinese Pulsed Power Society Workshop. Chengdu, 2015.
[18]	GEISSEL M, HARVEY-THOMPSON A J, AWE T J, et al. Minimizing scatter-losses during pre-heat for magneto-inertial fusion targets [J]. Physics of Plasmas, 2018, 25(2): 022706. doi: 10.1063/1.5003038
[19]	DAVIES J R, BAHR R E, BARNAK D H, et al. Laser entrance window transmission and reflection measurements for preheating in magnetized liner inertial fusion [J]. Physics of Plasmas, 2018, 25(6): 062704. doi: 10.1063/1.5030107
[20]	SLUTZ S A. On the feasibility of charged particle-beam preheat for MagLIF: SAND 2015-1515R [R]. Albuquerque, USA: Sandia National Laboratories, 2015.
[21]	赵海龙, 肖波, 王刚华, 等. 磁化套筒惯性聚变一维集成化数值模拟 [J]. 物理学报, 2020, 69(3): 035203. doi: 10.7498/aps.69.20191411 ZHAO H L, XIAO B, WANG G H, et al. One-dimensional integrated simulations of magnetized liner inertial fusion [J]. Acta Physica Sinica, 2020, 69(3): 035203. doi: 10.7498/aps.69.20191411
[22]	BASKO M M, KEMP A J, MEYER-TER-VEHN J. Ignition conditions for magnetized target fusion in cylindrical geometry [J]. Nuclear Fusion, 2000, 40(1): 59–68. doi: 10.1088/0029-5515/40/1/305
[23]	阚明先, 王刚华, 赵海龙, 等. 金属电阻率模型 [J]. 爆炸与冲击, 2013, 33(3): 282–286. doi: 10.11883/1001-1455(2013)03-0282-05 KAN M X, WANG G H, ZHAO H L, et al. Electrical resistivity model for metals [J]. Explosion and Shock Waves, 2013, 33(3): 282–286. doi: 10.11883/1001-1455(2013)03-0282-05
[24]	ZOLLWEG R J, LIEBERMANN R W. Electrical conductivity of nonideal plasmas [J]. Journal of Applied Physics, 1987, 62(9): 3621–3627. doi: 10.1063/1.339265
[25]	薛全喜, 江少恩, 王哲斌, 等. 基于神光Ⅲ原型装置开展的激光直接驱动准等熵压缩研究进展 [J]. 物理学报, 2018, 67(4): 045202. doi: 10.7498/aps.67.20172159 XUE Q C, JIANG S E, WANG Z B, et al. Progress of laser-driven quasi-isentropic compression study performed on SHENGUANG Ⅲ prototype laser facility [J]. Acta Physica Sinica, 2018, 67(4): 045202. doi: 10.7498/aps.67.20172159
[26]	JENNINGS C A, CHITTENDEN J P, CUNEO M E, et al. Circuit model for driving three-dimensional resistive MHD wire array Z-pinch calculations [J]. IEEE Transactions on Plasma Science, 2010, 38(4): 529–539. doi: 10.1109/TPS.2010.2042971
[27]	MCBRIDE R D, JENNINGS C A, VESEY R A, et al. Displacement current phenomena in the magnetically insulated transmission lines of the refurbished Z accelerator [J]. Physical Review Accelerators and Beams, 2010, 13(12): 120401. doi: 10.1103/PhysRevSTAB.13.120401
[28]	SINARS D B, SLUTZ S A, HERRMANN M C, et al. Measurements of Magneto-Rayleigh-Taylor instability growth during the implosion of initially solid al tubes driven by the 20-MA, 100-ns Z facility [J]. Physical Review Letters, 2010, 105(18): 185001. doi: 10.1103/PhysRevLett.105.185001