触地核爆炸冲击波效应的数值模拟研究★

王英杰，闫凤国，陈建营，薛海恩，梁基冠

(1.陆军工程大学,江苏 南京 210007; 2.中山大学,广东 珠海 519082)

1 概述

当核武器在地面发生爆炸时,核爆炸产生的能量分为两种形式传播,一是爆炸产生的能量直接耦合入土地面介质中并形成直接的冲击,是地下工程所要考虑的主要荷载;另外一种是以空气冲击波的形式在空气中传播,并会在地表面区域的介质中耦合产生感生冲击波[1-2]。同时触地爆炸时还需要考虑地面的自由表面作用,而自由表面效应和介质强度的综合效应,会导致接近自由表面处的峰值应力和峰值速度低于下方轴线上等径向距离处的峰值应力与峰值速度。因此,不同于在空气中爆炸的球对称传播模式,核武器触地爆炸时冲击波传播较为复杂。本文首先通过数值模拟结果与经验公式对比的方式,用TrueGrid软件和LS-DYNA软件建立更加复杂的工况并进行模拟分析,分析地面材料为多种介质时,在触地爆炸荷载作用下介质中应力波的传播过程。

2 触地爆效应计算的经验公式

目前对触地爆的研究相对较少,但在少量的文献中,依然可以找到相关的经验公式,并且触地爆的经验公式总是与地下封闭爆炸存在一定的联系,一般情况下地下爆炸的经验公式是炸药的质量乘上某一折减系数。文献[3]中就对触地爆的经验公式做了较为详细的阐述,如沃克夫、彻尔内、斯米尔诺夫、廖哈夫、克拉维茨、莱普逊和其他的学者给定了关于土在爆炸荷载作用下的相关参数。在本文中主要采用廖哈夫和克拉维茨提出的经验公式进行对比,分别表示为式(1)，式(2):

(1)

(2)

其中,Pi为爆炸冲击波在土中形成的超压,kg/cm2,约等于0.1 MPa;R为炸药中心到测试点的距离,m;W为炸药的质量,kg。

另外,克拉维茨还提出了触地爆时,土中最大质点速度的经验公式,如式(3)所示：

(3)

3 材料本构模型

3.1 空气材料模型

空气的成分复杂,在计算时我们需要根据空气的特性将其简化。在大当量的爆炸情况下,我们可以将空气理想化为非黏性理想气体,选用9号材料*MAT-NULL[4]来作为空气的本构模型,为了描述冲击波压力与空气初始内能密度的关系,我们引入状态方程*EOS-LINEAR-POLYNOMIAL,如式(4)，式(5)所示:

Pa=C0+C1u+C2u2+C3u3+(C4+C5u+C6u2)Ea0

(4)

(5)

其中，Pa为空气冲击波压力;Ea0为空气初始单位体积内能;ρ为扩散密度;ρ0为标准密度;C0,C1,C2,C3,C4,C5,C6均为计算常数,当该状态方程描述空气时,C4和C5可由式(6)计算:

C4=C5=λ-1

(6)

其中，λ为空气的绝热指数;空气的各参数取值见表1。

表1 空气的主要计算参数

3.2 炸药材料模型

由于重核裂变和轻核聚变等核武器的链式反应的复杂性,目前国内外关于冲击波效应,特别是核武器产生的冲击波效应大多采用等效TNT当量来近似替代[5]。因此我们用8号材料*MAT-HIGH-EX·PLOSIVE-BURN作为等效TNT的本构模型,并通过引入状态方程*EOS-JWL用于表示TNT爆轰波阵面压力与初始相对体积V0、单位体积内能E0之间的关系,如式(7)所示:

(7)

其中，P为爆轰压力;V0为初始相对体积;E0为初始单位体积内能;A，ω，R1，B，R2均为材料常数。炸药的各参数取值如表2所示。

表2 炸药的主要计算参数

3.3 岩石材料模型

国内外对于岩石在大当量爆炸理论的研究较少,为了方便计算,经过综合考虑，我们将岩石视为弹塑性材料,因为3号模型*MAT-PLASTIC-KINEMATIC能充分考虑岩石的应变率影响,适合用于描述岩石材料模型的各向同性硬化和随动硬化塑性[6]。因此我们在模拟中选用3号模型*MAT-PLASTIC-KINEMATIC作为岩石的材料模型。在爆炸荷载的作用下材料的应变速率的影响不可忽略,材料模型的屈服强度可分为初始屈服强度和硬化部分强度,屈服强度可按式(8)计算:

(8)

(9)

表3 岩石的主要计算参数

3.4 土壤材料模型

土壤的成分复杂,在受到冲击后会被迅速压缩,为了能较好的模拟出土壤在受到冲击波后的压缩性能,我们一般都用5号材料模型*MAT-SOIL-AND-FOAM来作为土壤的本构模型。并可通过定义一些参数以考虑压力与体积应变关系中热效应对其所产生的影响。该材料模型的屈服函数如式(10)所示:

φ=J2-(a0+a1p+a2p2)

(10)

其中，a0，a1，a2均为材料常数,a0考虑了土壤摩擦角的影响,a1考虑了土壤黏聚力的影响,a2考虑土壤受动荷载效应的影响;p为静态压力。该屈服函数中引入第二个不变量J2来描述:

(11)

Sij=σij-δijp

(12)

(13)

其中，Sij为偏应力张量;σij为应力张量;δij为克罗内克函数。土壤的各参数取值如表4所示[7]。

表4 土壤的主要计算参数

4 触地爆的有限元模型

为了方便建模及提高计算效率,触地爆采用1/4对称模型,炸药采用等效TNT方形炸药,尺寸为85 m×85 m×85 m,所等效的TNT当量约为100万t。建模过程结合TrueGrid软件和LS-DYNA软件,分为两个部分,分别是空气和炸药、土壤与岩石。两个部分均采用TrueGrid软件建模,再将两个部分分别导入LS-DYNA软件中形成完整模型,如图1所示。图1(a)是炸药与空气模型,整体尺寸为2 600 m×1 300 m,图1(b)是土、岩混合介质模型,整体尺寸为1 200 m×1 200 m。该模型均采用Solid164单元,土与岩石接触面采用自动接触。在建模过程中,炸药与空气的单元网格均采用共节点的方式建模[8-9]。

5 纯土介质中的触地爆效应

为了验证数值模拟方式的准确性,首先开展纯土介质中核武器触地爆的数值模拟,用于与经验公式给出的触地爆炸的结果进行对比。由于与经验公式进行对比的数值模拟的纯土介质的有限元模型与章节4中所述模型基本相同,仅将土、岩石全部由土代替,因此不再赘述模型的建立过程,直接对数值模拟结果进行对比分析。

5.1 土中应力波的传播过程分析

图2呈现了土中应力波传递的大致过程。在爆炸后的初始阶段(0.06 s之前),应力波近似以半球形对称的方式向四周传播,沿地表面方向和沿径向深度方向的应力波峰值差别不大,但从应力波传播区域来看,沿纵向深度方向的应力波传播明显大于沿地表面传播方向,可以发现在触地爆炸过程中更多的能量是直接耦合入土介质中;随着冲击波的传播(0.06 s～0.14 s),冲击波波阵面以近似椭球型的方式继续向四周传播,沿着纵向深度方向的冲击波传播较为规律,冲击波超压峰值不断衰减;但由于地表面的自由表面效应和介质衰减的综合效应,地表面的冲击波传播过程较为复杂,会出现冲击波的间断和不规律传播;在0.14 s以后,冲击波传播到所建立的模型之外的区域,在目标土介质区域,应力波主要传播至较深深度处,另外,地表面处由于空气冲击波影响而产生的感生地冲击使得地表面处的应力峰值相对较大,并且土层表面在爆炸荷载的作用下,会形成一个圆形的爆坑,爆坑的直径随着时间的增长逐渐增大,在0.4 s后逐渐趋向于稳定。

5.2 土中应力波的衰减过程分析

为了更加清楚地观测应力波在土层中的衰减过程,并且对比数值模拟结果的准确性,利用origin软件对土层中不同比例距离处的入射超压变化曲线进行曲线拟合。图3(a)中表示了土层中入射超压在不同比例距离的情况下随时间变化的过程,我们可以发现,入射超压在土层中传递时,呈快速递减到缓慢递减的一个过程。结合图3(b),当比例距离为0.2 m/kg1/3～0.4 m/kg1/3时,超压峰值相较于其他段递减趋势更加明显,这就更加反映了在近距离情况下爆炸荷载作用衰减较快,且其过程较为复杂。图3(b)还对比了经验公式计算得到的不同比例距离处的入射超压峰值与数值模拟结果:从整体趋势看,数值模拟计算得到的结果与廖哈夫和克劳维茨经验公式的整体趋势均较为吻合,但当比例距离小于0.4 m/kg1/3时,经验公式计算得到的入射超压峰值均稍小于数值模拟得到的结果,这可能是由于我们采用的土层材料模型的密度与经验公式的材料模型的密度有所差别。

5.3 质点运动速度分析

在对地表面中爆炸引起的动力响应进行研究时,质点的运动速度也是一个需要关注的重要变量,因此除了对土层中的入射超压峰值进行验证,我们还对比分析了不同比例距离深度处测点的运动速度。图4(a)给出了选定的不同深度土层中质点运动速度变化的情况,以某一深度土层质点来看:其质点运动速度表现为快速上升并缓慢下降的过程,这是因为爆炸冲击波到达土层中某一点时,会引起该点处的速度瞬间发生变化,当冲击波离开以后,该点处速度逐渐减小,而由于土层中距离炸药越近的部分产生的塑性变形越大,因此土层中某一点处的速度减小时,会受到前面土层的阻碍,消耗的时间也相对增长。结合图4(a)、图4(b)也可以发现,随着比例距离深度的增加,其最大质点速度也迅速衰减,当比例距离达到0.6 m/kg1/3之后,其最大质点速度逐渐保持稳定。由于现有的触地爆研究资料较少,关于触地爆在土中引起的最大质点速度成熟的经验公式也较少,因此在图4(b)中仅用克拉维茨提出的最大质点速度经验公式作为对比参考。从图4中可以看出,数值模拟得到的最大质点速度与经验公式的计算结果吻合较好。

综合土介质中爆炸应力波的传播过程和对应力波超压峰值、质点最大速度的分析来看,数值模拟得到的应力波传播过程和经验研究较为符合,且计算得到的应力比超压峰值和质点最大速度与经验公式拟合较好,验证了所采用的数值模拟方法的合理性和数值模型参数的准确性。

6 土、岩混合介质中的触地爆效应

从上述纯土介质中的结果对比分析来看,此数值模拟方法对触地爆的研究分析的结果是合理的。所以,我们进一步将该数值模拟方法推广至土、岩混合介质的情形,并分析爆炸荷载作用下土、岩混合介质的动力响应。土层与岩层之间的接触均采用自动面面接触。

6.1 土、岩混合介质中应力波的传播过程分析

图5描述了爆炸冲击波在多种介质中的传递情况。从传播过程可以看出,触地爆炸产生的应力波整体传播规律和在纯土介质中的传播规律差别不大,主要呈现出以下几个特点:爆炸产生的冲击波同时沿着深度方向和地表面方向传播,且主要沿深度方向传播;在沿地表面深度方向,依然表现出入射超压递减的趋势和沿地表面传播方向,会受到自由表面效应的影响;但和纯土介质不同的是,由于不同介质性质的差异,在应力波传播过程中在土层与岩层的接触面处会存在应力积累,这是因为土层与岩层在爆炸荷载作用下产生的塑性变形不连续,造成了该现象的产生,使得混合介质中应力波的分布结果相对于纯土介质较为复杂。

图6显示了混合地质的质点速度随时间变化的云图。可以看出,起爆后,和炸药接触的土壤层的质点速度迅速上升,在比例爆距0.1 m/kg1/3可达到接近100 m/s,随着离地越远,质点速度逐渐下降,当碰到地下岩层的时候,质点速度开始出现急剧下降,这是因为岩层的弹性模量较大,其相对于土层坚硬了许多,因此也导致了冲击波经过岩层时质点速度发生锐减。随着距离地面深度的增加,在冲击波传至第二层深埋土时,其最大质点速度已经衰减至1 m/s以下。

6.2 土、岩混合介质中应力波的衰减过程分析

图7中给出了应力波在土层与岩层混合介质地层中传播的整体过程,并将数值模拟得到的入射超压峰值与纯土介质的计算结果进行对比。从整体上看,在多种介质混合的情况下,入射超压在传播过程中会发生较大的变化,但总体的应力衰减较慢,特别是在比例距离为0.202 m/kg1/3时,该深度处的应力水平会长时间保持较高的水平,这是因为该深度恰好在土层与岩层的接触位置附近。从图7(b)中我们可以发现,在土层、岩层混合的地质中,不同比例距离处的质点的入射超压峰值均大于纯土地质中的相应结果,并且这一现象在比例距离较小时更为显著,结合压力传播云图来看,这也是由于不同岩层接近位置处的应力积累造成的。综合以上对入射超压的分析,可以得出,混合地质中的入射超压比单一地质衰减要慢,其超压峰值相对较高,且在不同岩层的接触位置会出现应力积累的现象。

6.3 土、岩混合介质中质点速度和介质变形分析

图8给出了土岩混合介质工况中不同深度的质点速度及与单一介质计算结果的对比,可以明显看出,混合地质中质点速度的衰减速率要远大于在纯土地质中质点速度的衰减。由于选点多集中于岩层中,并结合对入射超压的描述,可以得出以下结论:应力波在岩石中的传播速度要大于在土介质中的传播速率,并且在岩石中质点速度衰减较快,这一现象在比例距离较小时差别较为显著,而随着比例距离的增加,二者差别逐渐减小。

图9是混合地质的质点位移随时间变化的云图。从中可以看出,其变化趋势和质点速度变化趋势大致相同,起爆后,炸药产生的压缩波迅速向地层传播,与炸药表面接触的土层位移较大,当冲击波传递至岩层时,质点位移出现明显的下降。从质点速度和质点位移的云图来看,到最终应力波破坏效果达到稳定后,土介质整体位移较大,地面形成爆坑,而岩石介质中开始爆炸导致的质点位移有部分逐渐恢复,综合来看,土岩混合岩层的不同介质差异对爆炸冲击波会产生的明显削弱作用[10]。

图10是混合地质的有效塑性应变随时间变化的云图。可以看出,土壤和岩层的有效塑性应变变化情况截然不同,与炸药接触的土层,仅在与炸药接触的较小范围内有明显的塑性应变,而与爆炸中心水平距离较远处的塑性应变不大,这是由于触地爆的冲击波主要是沿着地层垂直向下传播的原因;在岩层中,也仅与第一层土层接触的200 m范围内有较明显的塑性应变,随着深度的增加,冲击波在岩层中所产生的有效塑性应变发生骤减。整体来看,仅在触地爆较近的位置处有一定的等效塑性应变,而当距离较远时,在地介质中几乎无塑性变形产生。

7 结语

1)通过LS-DYNA软件,对核武器触地爆的过程进行了数值模拟,并将数值模拟得到的结果与前苏联廖哈夫和克拉维茨等人根据实验数据拟合提出的经验公式进行对比,数值模拟得到的结果和经验公式得到的结果整体规律是一致的,数据结果差别不大,验证了数值模拟方法的可靠性。

2)开展了纯土介质中核武器触地爆的数值模拟,分析了土介质中爆炸应力波的传播过程和对应力波超压峰值、质点最大速度等,数值模拟得到的应力波传播过程和经验研究较为符合,且计算得到的应力比超压峰值和质点最大速度与经验公式拟合较好,验证了所采用的数值模拟方法的合理性和数值模型参数的准确性。

3)通过对触地爆工况下所得实验结果的验证,建立了土岩介质混合的模型并对不同的工况结果进行了对比分析,通过对应力波传播过程的分析以及对爆炸冲击波所产生的峰值超压、质点最大速度、最大位移、有效塑性应变分布等特征参数的分析比较来看,在地下设置不同特性的岩层,会在岩层交界处存在应力积累线性,但其最大质点速度和质点位移均会衰减较快,因此发现不同岩层混合能够大幅降低爆炸产生的冲击波的动态破坏效应。