王超申，康道辉，游世辉，谢纯凯，曽宪任

(1.湘潭大学 土木工程与力学学院,湘潭 411105; 2.枣庄学院 机电工程学院,枣庄 277160)

1 引 言

近年来，如何提高地下防护工程在防御战时抵抗各种武器杀伤破坏的能力是国内外学者一直研究的问题[1]。李砚召等[2]通过平面装药加载试验,证明了分配层分层结构对核爆炸防护效果优于单层结构。文献[3，4]都以空壳复合材料为分配层进行集团装药化爆试验，证明这类材料既对爆炸波有衰减弥散作用，也具备承受多次打击的能力。赵跃堂等[5]模拟发现，分配层厚度越大，对冲击波的削弱作用越明显，并且支撑结构顶部荷载越均匀。叶中豹[6]提出新型三维Kelvin-Maxwell双粘性系数的等效本构模型，并用于模拟分配层中的空壳复合材料，通过试验和模拟结果对比，证明了这种本构模型的可靠性。然而，以上研究的各层材料之间的界面是水平的，波的入射角是固定的，但根据波动理论，应力波的几何扩散效应会影响波的衰减弥散[7]，利用这种效应，可改变分配层的几何形状来扩散冲击荷载。其次，分配层厚度在一定的情况下，各层材料的顺序还有待优化。

光滑粒子流体动力学法SPH(Smoothed Particle Hydrodynamics)可用于解决爆炸冲击的问题，但存在精度与效率低和边界难以处理的缺陷。而FEM法具有计算效率高和边界条件设置方便的优点，可以很好地弥补SPH法的不足[8]。SPH-FEM耦合法是在LS-DYNA的框架内建立FEM单元和SPH粒子并耦合的数值分析方法，SPH法中不同材料粒子的相对运动可形成界面滑移，光滑连续地处理物质交界面，能很好地模拟爆炸冲击引起的成坑和物质飞散现象。

典型的成层式防护工程的结构主要有伪装土层、遮弹层、分配层和支撑结构4部分。本文采取有限元网格与光滑粒子相耦合的方法，重点研究不同排列方式的马蹄形分配层结构对爆炸波衰减弥散的影响，旨在寻找更优良的分配层结构来提高地下防护工程的抗爆防震性能。

2 SPH-FEM耦合模拟方法和模型建立

2.1 SPH-FEM耦合数值模拟方法

SPH法的核心是核函数近似法。对于任意函数f(x),都可以采用积分形式

(1)

式中Ω是包含x积分域，δ是狄拉克函数。采用光滑函数W(x-x′,h)取代狄拉克δ(x-x′)函数，f(x)的积分表达式为

(2)

式中W(x-x′,h)称为光滑函数或者核函数kernel function，h定义为核函数的光滑长度。用粒子体积ΔV对dx′进行替换，粒子的质量mj为

(3)

粒子在i处的表达式为

(4)

SPH-FEM耦合算法多采用固定耦合算法，设置点面固连接触，根据线动量守恒和角动量守恒调整粒子与单元节点的速度和位置，保证两者间的位移协调，并且将粒子区的力学信息传递到有限元网格。耦合界面处的粒子与单元节点一一对应，粒子与其对应的单元节点具有相同的加速度，其大小是由粒子与单元节点的质量以及作用于两者的力共同确定[9,10]。在本文的模拟中，近场SPH区域和远场FEM区域是相同的岩石介质，粒子与其对应的单元节点将在整个计算过程中保持完全一致，如图1所示。

2.2 材料模型

本文采用的减震吸能的材料有砂土和泡沫混凝土。砂土是含水细砂，含水率控制在3%左右。泡沫混凝土模型使用可压扁材料泡沫模型*MAT_CRUSHABLE_FOAM。泡沫混凝土相关参数如下[11]，密度为0.799 g/cm3，泊松比为0.1，杨氏模量为342.2 MPa，抗拉强度为0.2 MPa，阻尼系数为 0.1。岩石材料和混凝土材料均采用适应高压、高应变率的模型*MAT_JOHNSON_HOLMQUIST_CONCRETE(HJC),HJC材料模型在爆炸荷载下能够很好地反映岩石等脆性材料损伤效应的动态响应，其等效屈服强度是压力、应变率和损伤的函数，岩石材料[12]、钢筋混凝土[13]和混凝土[14]的相关参数列入表1。本文使用Wang[15]推荐的细砂，其本构模型采用*SOIL_AND_FOAM，主要的参数值如下，密度为1800 kg/m3，剪切模量为 63.8 MPa，初始体积模量为134.0 MPa,内摩擦角为35°，表观粘聚力近似取 5 kPa，含水率为3.4%。另外，为增加分配层的抗爆性和增大材料间的不匹配度，在分配层中加入钢铁。

图1 SPH与FEM的耦合

Fig.1 Coupling of SPH to FEM

2.3 几何模型

在本文数值模拟中，考虑到岩体触地爆炸的问题具有对称特点，故采用1/4模型。采用SPH-FEM耦合法建模，先在LS-DYNA中建立TNT炸药-岩体-锚杆-衬砌的有限元模型，数值模型尺寸为470 cm×400 cm×840 cm，地下硐室采用马蹄形结构，为典型的成层式防护工程。地下巷道的最外层是100 mm的喷射混凝土层，向内是 300 mm 的分配层，再将分配层分为6层，如图2所示，这 6层 分配层每层的厚度均为50 mm。耦合模型中钢筋材料采用BEAM161四节点梁单元，其余材料均为SOLID164六面体单元。整个模型共划分为496960个有限单元和30762个SPH粒子，其中炸药粒子为216个，爆炸近区岩体粒子为30546个。图3为Lagrange单元和SPH粒子的分布情况。

SPH粒子与有限单元面采用自动点面接触，对称面施加全局约束以及SPH对称边界，为描述半无限空间，其余面定义无反射边界条件，并设置吸收剪切波和膨胀波，防止爆炸波在边界处发生反射，并对时间步长和人工体积粘性系数作相应调整。为了更好地观测钢筋混凝土破坏情况，设置了*MAT_ADD_EROSION失效准则,以最大拉应力为破坏准则。

图2 分配层结构

Fig.2 Distribution layer structure

图3 Lagrange单元SPH粒子分布情况

Fig.3 Lagrange unit and SPH particle distribution

3 模拟结果对比分析

对分配层结构的Lagrange单元进行监测，研究11个方案下各单元的位移、应力和应变。根据应力波传播的特点，距爆源中心越近，应力波峰值越大，因此选取模型对称轴线上的 7个 测点(A= H371761，B= H371600，C= H437961，D= H349280，E= H394561，F= H21440，G= H511280)作为研究对象，每个测点的单元在厚度上长度为 1.7 cm，表2是各测点应力峰值。测点分布如图4所示。

3.1 含不同材料的分配层结构对比分析

当分配层采用砂土时，等效应力呈现折线增大的趋势，图5(a)中，应力波在大约750 μs时开始从0增加，1500 μs时幅值开始降低，2250 μs时才达到峰值。而图5(b)应力波大约在600 μs时开始从0增加，1750 μs时幅值开始降低，3500 μs时到达峰值。从图6可以看出，5000 μs时，钢筋混凝土支撑结构出现失效单元(浅灰色单元)，表明砂土能更大程度地削弱到达钢筋混凝土支撑结构的应力波强度。从表2中测点F的应力峰值可知，方案〈1〉和〈2〉的爆炸波应力峰值分别衰减了99.2%和 96.4%。由于钢铁与砂土的波阻抗比值要大于钢铁与泡沫混凝土的波阻抗比值，分配层厚度一定的情况下，含水砂土对爆炸应力波的衰减明显好于泡沫混凝土。

表2 测点应力峰值

Tab.2 Stress peak at the measuring point

方案应力峰值/MPa遮弹层ABCDEFG〈1〉ggssss180.0180.00180.001.501.351.611.4220.58〈2〉gghhhh180.0180.00180.006.886.066.156.4722.57〈3〉gghhss180.0180.00180.006.845.990.890.7816.70〈4〉gghshs180.0180.00180.006.966.926.861.0120.89〈5〉ggshsh180.0180.00180.001.196.571.546.8120.12〈6〉gshgsh180.0180.002.116.1315.361.304.376.03〈7〉gshshg180.0180.001.386.827.011.1415.5322.12〈8〉ghsghs180.0180.006.971.8023.526.861.0816.68〈9〉shshsh180.02.166.900.776.521.176.128.46〈10〉hshshs180.06.855.356.725.245.700.787.25〈11〉shgshg180.02.216.9621.620.443.5620.194.77注:将钢铁材料缩写为g,泡沫混凝土为h,砂土为s。

图4 测点分布

Fig.4 Survey point distribution map

3.2 不同排列方式分配层结构的对比分析

爆炸波在不同介质中能量损耗和衰减速度相差很大，分配层中相邻介质材料的波阻抗、应变率影响因子黏度系数和材料的屈服强度等相差越大，对爆炸波的吸能弥散作用越明显[16]，并且爆炸波从硬材料入射到软材料中对波的削弱效果更强。另外，分配层材料采用周期性布置时，弥散效应会得到增强，本文的数值模拟也采取这种布置方式。

3.2.1 测点应力对比分析

图7是方案〈6〉～〈11〉测点G的等效应力曲线。可见，冲击波大约在800 μs时到达钢筋混凝土支撑结构，随着时间的推移，冲击波在分配层结构内不断发生反射和透射，使得钢筋混凝土支撑结构的应力一直在增大。当分层材料相同并且各种材料厚度相同时，方案〈9〉和〈10〉都只将砂土和泡沫混凝土采用周期性布置的方式，方案〈9〉中每个周期内材料波阻抗值是递增的，而方案〈10〉每个周期内材料波阻抗值是递减的，应力波峰值分别衰减了95.4%和96.0%。

图5 砂土和泡沫混凝土测点F有效应力时程曲线

Fig.5F-point effective stress time-course curve of sand and foam concrete

方案〈9〉和方案〈11〉的泡沫混凝土和砂土排列顺序一致，方案〈9〉中无钢铁，其应力波峰值衰减了95.3%；方案〈11〉分配层从上往下，第三层和第六层为钢铁，应力波衰减了97.4%，削波能力较好。分配层中加入钢铁后，更大程度削弱了应力波的传播，这是由于钢铁受反射波和大变形的作用，使得其也随多孔材料一起变形，吸收了一部分能量，因此局部的破坏作用加大，而钢筋混凝土支撑结构受到的破坏减小。

图6 砂土与泡沫混凝土分配层等效应力云图

Fig.6 Equivalent stress cloud diagram of sand and foam concrete distribution layer

图7 方案〈6〉～〈11〉测点G的等效应力曲线

Fig.7 Scheme 〈6〉～〈11〉 equivalent stress curve of pointG

方案〈6〉和方案〈8〉分配层从上往下，第一层和第三层为钢铁，其他四层分别按波阻抗值递增和递减来排列，应力波峰值分别衰减了96.7%和90.9%。分层结构中材料按波阻抗比递增的顺序排列时，会增强应力波强度，因此，本节模拟只考虑按波阻抗比顺序排列减小的情况。方案〈8〉中一个周期的材料顺序是按波阻抗值递减和波阻抗比减小的方式来布置的，应力波衰减了 90.9%；方案〈11〉一个周期的材料顺序是采用相邻界面材料波阻抗值递增和波阻抗比减小的方式来布置，应力波衰减97.4%,略优于方案〈8〉。很明显，分配结构中各层材料顺序按波阻抗值递增和波阻抗比减小、周期性布置和加入钢铁是较优异的方案。

3.2.2 能量对比分析

图8是六种方案中钢筋混凝土支撑结构的总能量变化曲线。可以看出，方案〈6〉和〈9〉～〈11〉在5000 μs时钢筋混凝土支撑结构的总能量大约是方案〈7〉和〈8〉的1/10。其中，方案〈11〉的支撑结构总能量最低，爆炸开始后，动能上升速度加快，内能上升速度较缓，5000 μs时，支撑结构动能占总能量的75%；而方案〈7〉的支撑结构总能量最高，5000 μs时，其动能占总能量的60%。对比两种方案可知，动能是总能量最重要的组成部分，如图9所示。方案〈7〉的支撑结构动能为39636 J，内能为28455 J；方案〈11〉的支撑结构动能为5578 J，较方案〈7〉约低86%，内能为1838 J，较方案〈7〉约低94%。方案〈8〉每一周期的材料是按波阻抗递减的方式来布置的，总能量约为70450 J，对比分析可知，此方案吸能效果较差。从图8能量曲线的斜率来看，5000 μs后，方案〈7〉和〈8〉中钢筋混凝土结构的总能量增加的速率要比方案〈11〉大，而方案〈11〉能量的增长明显变缓。

图8 方案〈6〉～〈11〉支撑结构总能量曲线

Fig.8 Scheme 〈6〉～〈11〉 total energy curve of the support structure

3.3 拱顶位移对比分析

表3是11种方案测点G的Y方向位移。在爆炸载荷作用下，钢筋混凝土支撑结构离爆心最近的部位竖向位移会最大，巷道拱顶的竖向位移与支撑结构顶部的位移很接近，因此，可以通过控制点G的竖向位移来减少爆炸波对支撑结构的破坏。11种方案中，除了方案〈7〉和〈11〉位移峰值偏大，方案〈6〉和〈8〉测点G的位移峰值最小以外，其他方案的位移峰值比较接近，约为最大位移的60%，可以说明这两个方案中的支撑结构受到来自分配层的能量较小，而方案〈6〉对支撑结构的保护更佳。分配层的变形不断调低应力峰值，减缓并减弱了爆炸波对支撑结构的破坏，对于遭受爆炸荷载的防护工程来讲，巷道拱顶的挠度值(沉降量)是判定防护结构能否抵抗二次爆炸冲击的关键指标之一。因此，提高钢筋混凝土支撑结构的抗压弯破坏能力的同时,还需科学地优化分配层结构。

图9 方案〈7〉和方案〈11〉支撑结构动能和内能曲线

Fig.9 Scheme 〈7〉 and 〈11〉 kinetic energy and internal energy curve of support structure

表3 测点G竖向位移峰值

Tab.3 Peak displacement ofGpoint vertical displacement

方案位移峰值/cmABCDEFG〈1〉ggssss2.212.191.220.840.580.210.03〈2〉gghhhh1.941.911.611.040.530.440.15〈3〉gghhss1.951.921.631.071.200.450.05〈4〉gghshs1.931.911.610.680.620.440.14〈5〉ggshsh1.911.901.191.150.560.390.14〈6〉gshgsh2.681.400.780.570.230.120.04〈7〉gshshg3.511.971.710.580.060.060.06〈8〉ghsghs2.441.690.600.690.480.340.10〈9〉shshsh2.191.390.820.470.150.100.04〈10〉hshshs1.460.870.590.230.160.120.04〈11〉shgshg1.900.860.180.060.040.030.03

3.4 误差分析

在强烈爆炸作用下，爆炸近区的动态响应非常复杂，现有的连续介质力学和线弹性力学等理论基础不足以准确描述爆炸近区的岩体变形特点，我国对于常规武器在掩体的爆炸效应的研究大多是在前苏联的试验成果的基础上进行的。如谢苗基(Щемякин Е И)和麦德维捷夫(Медведев Н С),通过大量试验得出冲击波压力和距离的关系[17]。

δm=p2/γα

(6)

(7)

测点依然是从11种方案的地下硐室中心线上选择，应力峰值误差分析选择距爆心 1.30 m,1.45 m,1.55 m,1.65 m,1.75 m,1.85 m和1.95 m处的测点。根据模拟的数据显示，距爆心1.80 m处，应力峰值变小的速率加快，1.95 m时急剧变小，与经验值的误差变大,列入表5。这是由于材料具有一定的粘性，随着应力波传播距离的增加，不断消耗能量的同时还随着波面面积的增加造成能量的扩散而衰减，冲击波逐渐变为压缩波。经计算及对比分析，在爆炸近区，测点的数值模拟结果与计算值之间误差在2.1%～3.5%之间，考虑到参数设置以及网格划分不是很准确，两者计算结果较为吻合。

表4 经验公式参数

Tab.4 Empirical formula parameters

装药初始密度/g·cm-3爆速/m·s-1岩石初始密度/g·cm-3声传播速度/m·s-1应力衰减指数1.60175302.83223501.95

表5 数值模拟结果与经验公式对比

Tab.5 Comparison of numerical simulation results with empirical formula

项目冲击波压力/MPa1.30m1.45m1.55m1.65m1.75m1.85m1.95m数值模拟161.53129.60107.66100.8885.5066.6634.56经验公式157.07126.94111.4798.6787.9878.9471.24误差/%2.82.13.52.22.915.651.5

4 结 论

(1)从含不同材料的分配层结构对比分析结果发现，低含水砂土对波的衰减和弥散作用要优于泡沫混凝土，可见砂土的波阻抗值较小。

(2)采用有利于削弱应力波峰值的周期性布置分配层材料方式模拟。从方案〈9〉和〈10〉模拟的结果来看，每一周期材料按波阻抗值增加顺序排列要优于其按波阻抗值减小削波能力。而对比方案〈6〉和〈8〉，泡沫材料和砂土组合顺序不变，加入钢铁后，每一周期材料按波阻抗值递增顺序排列要更优于其按波阻抗值递减的削波能力。可见，分配层结构的软材料中加入硬材料后，更能削弱应力波的强度。

(3)大波阻抗的钢铁与低波阻抗值的砂土紧密贴合极大增强了分配层的削波能力，没加钢铁与加钢铁的方案即〈9〉和〈11〉对比，方案〈11〉的钢筋混凝土支撑结构G处的应力峰值比方案〈9〉约小44%，表明加入钢铁的分配层吸能效果更佳。与方案〈8〉相比，方案〈11〉采用周期材料按相邻界面材料波阻抗值递增和波阻抗比减小的方式来布置更优。

(4)通过对11种方案钢筋混凝土支撑结构的总能量和竖向位移进行对比分析，方案〈11〉中支撑结构的总能量和位移响应最小，并且方案〈11〉的分配层结构主要是削弱了应力波的动能。因此，方案〈11〉的分配层结构抗震吸能效果最佳，对支撑结构的保护作用最大。

(5)将数值模拟得到的不同测点的应力峰值与前苏联经验公式的计算值对比发现，在距爆心较近位置的应力峰值计算结果和模拟结果比较吻合。