核爆炸对地冲击作用下土体运动特性等效模拟

荣吉利， 宋逸博,2， 王玺， 郭振， 项大林， 吴志培,2

(1.北京理工大学 宇航学院, 北京 100081; 2.中国运载火箭技术研究院, 北京 100076;3.北京航天发射技术研究所, 北京 100076; 4.北京宇航系统工程研究所, 北京 100076)

0 引言

核武器由于其巨大的毁伤能力及战略价值，得到了快速的发展。以往对于地面及地下结构在核爆袭击下的强度校核，主要采用冲击波超压峰值和正压作用时间或冲量值作为创伤指示参数[1]，该方法过于单一且不能完全反映出结构所遭受的威胁。随着武器技术的发展以及对核爆效应的了解日益深入，由爆炸产生的对地冲击效应引起了广泛的关注，并被认为是对地下结构最大的威胁。1968年4月26日，在美国内华达地下进行的1 200 kt当量的Boxcar核爆炸产生了强大的对地冲击效应及岩土的不可逆变形，引起最远震中距6 100 m处地表出现断层破裂，断层错动的最大垂直相对位移为1 m[2]. 美军在1964至1976年间对土体介质及岩石场地进行了一系列大当量爆炸试验以对核爆炸对地冲击效应进行论证分析[3-7]，并将研究成果拟合成经验公式，整理在TM 5-858设施抗核武器效应设计手册中的TM 5-858-2武器效应分册中[8]。Pathak等[9]利用一维弹性波传播模型，对核武器空爆引起的地运动位移进行了研究，为工程提出了一种预测地面位移的方法。我国学者也对爆炸地冲击效应进行了大量的研究[10-12]，形成了较为完备的计算方法。但由于对爆炸引起的地运动规律缺乏有效且统一的评估方法，不同计算方式对地冲击效应的计算结果存在较大差异，因此难以得到一致结论。此外大当量的对地冲击试验工程巨大，难以大范围开展获取充足的试验数据，因此对地冲击效应的认识仍存在一定的局限性。

随着计算机模拟技术的发展，对核爆炸对地冲击效应进行准确的数值仿真成为了可能，林大超等[13]提出地面随机平稳噪声的时变曲线是一条分型曲线的假设，给出了爆炸所致地面运动的仿真方法和过程；唐廷等[14]利用有限元分析软件LS-DYNA对地冲击后地下预置裂纹拱结构进行了抗爆炸分析，以研究受损的硐室抵抗二次爆炸的能力；Ma等[15]、李永胜等[16]分析了地下硐室小当量化爆下岩土响应，其结果与试验结果拟合较好，但相关工作研究的爆炸当量较小，无法应用于核爆等大当量爆炸仿真；王小盾等[17]、苗青等[18]针对天津8·12爆炸事故，考虑了爆炸空气冲击波与爆炸对地冲击波联合作用，利用有限元分析软件Abaqus对网架结构动力响应进行分析，证明了爆炸对地冲击效应对结构的影响不可忽视，但由于缺少相关实测数据，仿真结果缺乏对比验证，模型的准确性和可靠性难以保证，且没有考虑空气冲击波激发的感生冲击波的作用。

综上所述，对于大当量爆炸下土体运动准确仿真方法的建立与模型的选取亟需探讨。本文考虑了感生冲击波与直接冲击波的叠加作用，通过使用流体- 固体耦合技术，得到了由空气冲击波激发的感生冲击波，研究了二者的贡献程度，证明了爆炸引起的对地冲击效应是感生冲击波和直接冲击波联合作用的产物。提出了以Shock Hugoniot绝热关系描述的Mie-Gruneisen状态方程以及Drucker-Prager(D-P)弹塑性准则来描述黏土的行为，使用JWL状态方程描述爆轰产物的分析方法，通过有限元分析软件AUTODYN以美军代号LN302的对地冲击试验[3]为例进行数值仿真模拟。将分析结果与美军试验数据和美军手册TM5-858-2的公式计算结果进行对比，吻合良好，验证了本模型的准确性和有效性，为核爆炸对地冲击效应的仿真模拟提供了有效的方法。

1 核爆炸对地冲击试验

1968年，美军实施了LN302对地冲击试验中代号为PRAIRIE FLAT的子试验[3]，以研究黏土地质下核爆炸对地冲击效应。试验在广阔的平原区域进行，该地区由河床沉积而成，没有岩层，主要地质成分为黏土。由于核爆炸与化学爆炸主要区别在于爆炸反应发生的初期阶段经历了不同的物理过程以及有无后续辐射产物，而对于爆炸对地冲击波的产生几乎没有差别，因此美军使用高能炸药梯恩梯(TNT)作为核爆炸的等效爆炸源[12]。在该试验中，等效TNT当量为500 t，成球状堆叠置于地表。根据预测的超压峰值，在距离爆心25.603～350.520 m的距离范围内布置了地运动测量仪器，深度0.457～9.144 m. 仪器封装在铝制圆柱罐内，固定置于地下，其具体布置方案如表1所示。

表1 地运动测量仪器布置阵列[3]

核武器在地表爆炸会在空气中产生向外高速扩张的球面空气冲击波，空气冲击波掠过地表，其激波前缘会黏连并牵引地表的土体与间隙的空气产生应力波扰动，称为感生冲击波，而由于爆炸释放的能量直接向地面耦合，在地下产生从弹坑处向外扩散的应力波称为直接冲击波[11]。美军通过理论分析与试验验证，将对地冲击作用区域从爆心向外划分为超震区和超前区。距离爆心较近的区域为超震区，由于爆心附近空气冲击波速度远大于爆炸引起的直接冲击波速度，因此位于超震区内的物体先接收到空气冲击波及其激发的感生冲击波，之后接收到直接冲击波。随着距离的增加，空气冲击波速度迅速衰减，而直接冲击波在地下传播速度不变，直接冲击波逐渐追赶并超越感生冲击波，以二者同时到达的距离为界，之后进入超前区。位于超前区的物体先感受到直接冲击波，后感受到冲击波和感生冲击波。

试验测得的正向速度峰值数据如表2所示，表中垂直速度和水平速度的正方向相对于爆心分别为向上和向外。

表2 对地冲击正向速度峰值试验数据[3]

2 美军手册TM 5-858-2经验公式

美军将一系列大当量爆炸对地冲击试验的数据结果拟合出用于计算核爆炸对地冲击效应的经验公式，作为手册TM 5-858-2[8]中核爆炸对地冲击效应的计算方法。在计算对地冲击速度峰值响应之前，需要根据爆炸工况判断目标位于超震区还是超前区，其判定方法如下：

(1)

式中：pso为冲击波压力最大值与大气压力po的差值(Pa)；L为与爆心的直线距离(m)；W为有效当量(kt)；u为空爆超压冲击波速度(m/s)；tf为从爆炸发生到直接冲击波到达目标位置所需的时间(s)；s为目标位置距离爆心的直线距离(m)；t为时间(s)；Range为判据，当Range≥1时使用超震区公式计算，当Range<1时使用超前区公式计算。

根据PRAIRIE FLAT试验报告中的设置工况和黏土的实测参数，确定公式中基本参数为当量W=0.5 kt，载荷下土壤压缩波速cL=130 m/s，s为布置有速度计的测点爆心距，分别为25.6 m、42.7 m、67.1 m、100.6 m与121.9 m，经计算均处于超震区，故使用超震区的速度峰值计算公式进行计算：

(2)

(3)

式中：uV为土壤垂直速度峰值(m/s)；uH为土壤水平速度峰值(m/s)；Z为目标深度(m)。对地下3.048 m深度的速度响应进行计算，得到的结果如表3所示，表中垂直速度和水平速度的正方向相对于爆心分别为向上和向外。

表3 3.048 m深度速度峰值手册结果[8]

3 仿真分析模型

3.1 材料模型

3.1.1 空气

空气使用理想气体状态方程描述，

p=(γ-1)ρe，

(4)

式中：p为气体压力；γ为绝热指数，γ=1+R/cv，R为一般气体常数R0除以特定气体的有效分子量，cv为气体体积恒定时的比热容；ρ为密度；e为比内能。其材料参数如表4所示。

表4 空气材料参数[19]

3.1.2 TNT

TNT使用JWL状态方程来描述，

(5)

式中：pe为爆轰产物压力；v=1/ρe为比容，ρe为装药密度；C1、r1、C2、r2为常数；ω为绝热常数；ee为爆轰产物的比内能。材料参数如表5所示。

表5 TNT材料参数[19]

3.1.3 土体

由于土体在应力作用下表现出较强的流动性，必须考虑其剪切应力与剪切破坏，因此本文以冲击状态方程结合基于摩尔- 库伦准则的D-P弹塑性强度模型对土体进行建模，并通过控制点定义屈服应力随压力分段线性变化来模拟描述土体抗压不抗拉的特性，当土体受拉时，屈服应力迅速降为0 Pa. 失效使用Hydro模型，当材料压力小于定义的最大拉伸压力时判定失效，参数如表6所示。

表6 土体材料参数[20]

3.2 模型建立与网格划分

根据文献[3]中的试验布置情况(见图1)建立如图2所示有限尺寸的二维轴对称模型，在试验速度计布置处设置测点。模型需要考虑空气与大地的耦合效应，建立150 m×400 m的欧拉网格模拟空气，建立150 m×400 m的拉格朗日网格模拟土体，为与经验公式结果进行对比，对土体进行简化，不考虑河床线对土体水饱和度的影响。由于测点分布在全距离下，因此使用均匀网格划分，网格尺寸为0.5 m×0.5 m. 二者设置50 m重叠以实现空气与土体的耦合作用，在地面以上的欧拉网格填充500 t球型TNT，对球心单点起爆以模拟核弹触地爆炸。对称轴处为对称边界，将欧拉区域边界条件设置为流出边界，将拉格朗日区域除地表面外的边界均设置为透射边界以消除边界反射，防止反射波与后续冲击波叠加。定义拉格朗日域自接触以模拟弹坑附近土体回落堆积，定义欧拉/拉格朗日耦合，实现应力从空气到大地的传播。

图1 试验测点布置Fig.1 Arrangement of test measuring points

图2 模型及测点布置Fig.2 Model and arrangement of measuring points

4 结果与分析

4.1 物质运动

基于上述建模方法与材料参数建立了500 t TNT在地面爆炸的地冲击效应仿真模型。为了观察完整的地冲击传递过程，仿真时间设置为2 s，输出时间间隔为0.002 s，截取爆炸过程部分时间的物质运动状态云图如图3所示。图3中深蓝色部分为空气，浅蓝色为黏土地质，土中呈阵列分布的点状物为测点，中心绿色部位为TNT。提取0.02 s时局部地冲击应力云图及速度云图如图4、图5所示。

图3 500 t TNT触地爆炸物质运动云图Fig.3 Material motion contour of 500 t TNT contact explosion

图4 0.02 s时对地冲击应力云图Fig.4 Stress contour of ground shock at 0.02 s

图5 0.02 s时对地冲击速度云图Fig.5 Velocity contour of ground shock at 0.02 s

通过图3可以看到，爆炸初期爆轰产物迅速膨胀形成的冲击波会带走爆心附近土体并形成半球形弹坑，后续的冲击波将弹坑周围的土体压实、抛掷和冲蚀，使弹坑继续扩大并形成凸起的唇缘。周围的土体在爆炸的挤压作用下产生了向外、向上的运动。距离爆心25.6 m的测点由于位于弹坑边缘附近的破裂区和塑性区，其位置在土体的流动下被大幅抬升并高于地面，随着与爆心距离的增加，土体位移幅度逐渐衰减。

4.2 对地冲击波

根据0.02 s时对地冲击应力云图可知，在靠近爆心的位置，由爆炸直接形成并从弹坑处向外扩散的球面应力波前方出现了一段与空气冲击波前缘相连的应力波。球面应力波为直接冲击波，在土体中以纵波波速传播，前方的应力波是感生冲击波。地冲击速度云图显示了土体中粒子速度的变化，测点先接收到感生冲击波带来的地冲击速度，直接冲击波随后到达。该现象与文献[8]中结论相互印证，也说明了对于地面建筑、设备和人员来说，在冲击波来临的同时还有其激发的感生冲击波在起作用，对于结构的强度校核不应只局限于超压，还应考虑与对地冲击作用的组合工况。

对直接冲击波与空气冲击波分离情况进行分析，结果如图6所示。

图6 超震区与超前区分界图Fig.6 Boundary of superseismic region and outrunning region

通过图6可以看到，在距离爆心约158 m处直接冲击波与空气冲击波前缘分离，与报告中超震区与超前区的预测分界点(约爆心距171 m处)十分接近，误差小于10%. 随着空气冲击波速度的衰减，先前生成的感生地震波与激波前缘分离，在地下以纵波波速继续向前传播并逐渐衰减，导致170.7 m处的测点仍能先接收到较强的感生冲击信号，因此实际的分界点更接近预测分界点。

分别选取位于超震区67 m和超前区253 m处的典型位置测点，对比分析其压力变化以研究直接地冲击波和感生地冲击波的贡献，结果如图7、图8所示。由图7、图8可以看出：位于超震区小于3.048 m的深度上，感生冲击波的压力峰值超过直接冲击波压力峰值，成为地冲击效应的主导部分，随着深度的增加，直接冲击波的压力峰值逐渐增大，感生冲击波的贡献比例逐渐降低；在超前区，直接冲击波起主导作用，压力幅值在距离上的衰减程度远小于感生冲击波，且随着深度的增加压力峰值逐渐增大。产生这种现象的原因是由于岩土中接近地表的直接冲击波会被地面运动释放，而地下的压力耗散缓慢，因此随着深度的增加强度衰减缓慢，而感生冲击波的压力是由空气冲击波提供，导致其峰值在深度上无明显变化。

图7 超震区对地冲击压力变化曲线Fig.7 Ground shock pressure curves in superseismic region

图8 超前区对地冲击压力曲线Fig.8 Ground shock pressure curves in outrunning region

4.3 速度时程

由于试验TNT当量为500 t，根据物质云图显示靠近爆心的浅层测点可能存在测量失真，因此选取25.6 m远、9.144 m深测点的速度波形与美军手册指导波形进行对比。使用该测点的归一化爆距在手册TM 5-858-2中选取波形，根据手册指导通过匹配峰值与上升时间的方式进行拟合，二者对比结果如图9、图10所示。通过对比可以发现二者吻合良好，对于垂直速度波形，其先出现方向向下，幅值与向上的速度峰值相近的脉冲，随后向上运动，之后速度缓慢下降。水平速度波形向外运动，之后速度迅速下降。波形对比的吻合印证了本模型得到的速度波形正确性，侧面证明了本仿真方法的准确性。

图9 垂直速度波形对比Fig.9 Comparison of vertical velocity waveforms

图10 水平速度波形对比Fig.10 Comparison of horizontal velocity waveforms

由于试验TNT当量为500 t，根据物质云图显示靠近爆心的浅层测点可能存在测量失真，因此选取25.6 m远、9.144 m深测点的速度波形与美军手册指导波形进行对比。经计算得到该测点的归一化爆距并在手册中选取波形进行拟合，二者对比结果如图9、图10所示。通过对比可以发现二者吻合良好，对于垂直速度波形，其先出现方向向下，幅值与向上的速度峰值相近的脉冲，随后向上运动，之后速度缓慢下降。水平速度波形向外运动，之后速度迅速下降。波形对比的吻合印证了本模型得到的速度波形正确性，侧面证明了本仿真方法的准确性。

将仿真得到的各测点正向速度峰值与美军试验数据使用B样条方法进行对比，结果如图11～图16所示。

图11 地面下0.457 2 m深垂直速度Fig.11 Vertical velocity at 0.457 2 m depth

图12 地面下0.457 2 m深的水平速度Fig.12 Horizontal velocity at 0.457 2 m depth

图13 地面下1.524 m深的垂直速度Fig.13 Vertical velocity at 1.524 m depth

图14 地面下1.524 m深的水平速度Fig.14 Horizontal velocity at 1.524 m depth

图15 地面下3.048 m深的垂直速度Fig.15 Vertical velocity at 3.048 m depth

图16 地面下3.048 m深的水平速度Fig.16 Horizontal velocity at 3.048 m depth

通过曲线可以看到土体的正向速度峰值随爆心距的增加呈现出指数衰减的趋势，同一距离不同深度的速度峰值也随深度而减小。通过对比可以发现仿真结果与试验数据拟合很好，由于试验的浅层测点表面覆土由人工回填，与周围土体可能存在一定的非连续性，且靠近爆心的测点位于弹坑附近，容易在爆炸作用下被发掘出地面，结合其他测点的整体趋势判断图11中42.7 m处试验数据可能存在失真。

使用美军防核设施设计手册武器效应分册的对地冲击效应计算公式，根据试验条件及当地黏土介质特性测量结果，对地下3.048 m的测点进行计算，并与上述结果进行对比，结果如图17、图18所示。通过对比可以看出三者整体拟合情况良好，仿真结果与试验数据更接近，垂直速度拟合优于水平速度，且均呈现出随爆心距的增加指数衰减的特征。其中防核设施设计手册水平分量与仿真结果的近场计算结果相比于试验数据均存在较大误差，计算结果高于试验数据，而仿真结果低于试验数据。原因可能是由于美军的计算公式基于大量地冲击试验数据进行拟合得到，由于爆炸效应的强非线性、测量数据的离散性以及拟合方式的选取等原因，使得拟合公式在强非线性段存在一定的误差，其特征是爆炸当量越大，爆心距离越小，深度越小，误差越大。此外有限元分析基于连续介质力学理论建立，真实试验场地的岩土介质实际上是非连续的，介质之间存在孔隙、节理等使得介质流动性更大、连接更松散，对于接近爆心的表层岩土亦存在抛掷及气化蒸发的可能，导致接近爆心处的岩土实际运动幅度比有限元仿真所得到的结果更大，造成了上述误差。

图17 地面下3.048 m深的垂直速度Fig.17 Vertical velocity at 3.048 m depth

图18 地面下3.048 m深的水平速度Fig.18 Horizontal velocity at 3.048 m depth

5 结论

核爆炸引起的对地冲击效应成分复杂，难以进行准确地仿真及预测。本文提出了一种综合考虑感生冲击波与直接冲击波叠加效应的核爆炸对地冲击准确描述方法，通过流体- 固体耦合技术还原了感生冲击波，使用Shock Hugoniot绝热关系描述的Mie-Gruneisen状态方程以及D-P弹塑性准则来描述黏土在对地冲击下的流动行为，并以美军核爆炸对地冲击试验为例对核爆炸对地冲击作用下土体运动特性进行了研究。得到如下结论：

1)本文所提出的仿真分析方法准确地反映了核爆炸对地冲击作用下黏土介质的运动特性，与美军地冲击试验数据以及美军手册计算公式结果吻合良好，可用于核爆炸对地冲击效应计算。

2)仿真结果表明核爆炸引起的空气冲击波会在地面激发感生冲击波，完整的核爆炸对地冲击效应是感生冲击波和直接冲击波两部分叠加作用的结果。在超震区浅层地表范围内感生冲击波起主导作用，在深层地下及超前区范围直接冲击波起主导作用。

3)在核爆炸打击下地面物体会同时受到空气冲击波及感生冲击波的共同作用，对于地面及浅层地表的目标，除了要经受空气冲击波的冲刷，还会经历地面瞬时大范围运动，因此在核爆炸袭击下对于建筑强度的校核不应只局限于冲击波超压峰值，还应考虑其抗震性能。

4)由于大当量爆炸的强非线性、测量数据的离散性以及拟合方式的选取等原因导致在接近爆心处仿真结果、试验数据与经验公式存在较大误差，因此难以准确描述爆心附近的土体运动情况。