吴志军1,2，张鹏林1,2，刘泉声1,2，李万峰3，江维中3
(1．武汉大学土木建筑工程学院，武汉 430072；2．岩土与结构工程安全湖北省重点实验室，武汉 430072；3．淮南矿业(集团)有限责任公司，安徽 淮南232001)
摘 要：军火库或危险品仓库存在着偶然爆炸的威胁，而钢筋混凝土是这些建筑物的主要构成材料，因此研究钢筋混凝土结构在爆炸荷载下的破坏过程具有重要意义。该文基于LS-DYNA动力有限元程序，利用任意拉格朗日–欧拉(ALE)方法，以及多物质流固耦合方法对混凝土结构在爆炸荷载作用下的动态破坏过程进行研究。为了更好分析混凝土结构在爆炸荷载作用下的动力响应，采用了考虑应变率影响的钢筋和混凝土材料本构模型，并引入零厚度粘聚力单元来模拟混凝土的动态破坏过程，克服基于侵蚀算法单元删除带来的质量损失问题。该文首先介绍零厚度粘聚力单元模型的生成过程并对比试验结果，验证所建立的零厚度粘聚力单元模型的合理性。其次，对比不同爆炸荷载下基于侵蚀算法以及零厚度粘聚力单元两种不同模拟方法的模拟结果，验证基于零厚度粘聚力单元模拟的优越性。最后基于零厚度粘聚力单元模型，分析不同爆炸荷载对混凝土结构动态破坏过程以及碎片抛射的影响。
关键词：爆炸荷载；钢筋混凝土；零厚度粘聚力单元；数值模拟；LS-DYNA
军火库或危险品仓库等存在着偶然爆炸的威胁，这些建筑物一旦发生爆炸，其产生的碎片将对周围建筑物以及人群产生极大的危险。由于钢筋混凝土是这些建筑物的主要构成材料，因此研究其在爆炸荷载作用下的动态力学响应及破坏过程，对于提高钢筋混凝土结构的抗爆安全性具有重要的指导意义。由于爆炸试验危险性高，且试验结果难以量测，因此对于爆炸破坏分析大多采用数值模拟进行。龚顺风等[1]采用侵蚀算法研究了不同炸药量作用下钢筋混凝土板的损伤机理和破坏特征；同样基于侵蚀算法，Tai等[2]探讨了不同炸药量对于混凝土板的动力响应以及破坏方式的影响，并分析了网格尺寸对于爆炸波传播的影响；张舵等[3]采用约束自由度的方法分析了爆炸荷载下钢筋与混凝土的相互作用，并对比了加载函数法与流固耦合法的模拟结果；李天华[4]基于单自由度体系动力分析和数值模拟方法对不同条件钢筋混凝土板抗爆性能进行了系统研究。目前，对于混凝土板在爆炸荷载作用下的破裂破坏分析，通常采用侵蚀单元模拟[1－3,5－9]，或者通过等效塑性应力或应变来实现[4,10－11]，亦或通过将处于同一位置的节点进行分离来模拟[12]。侵蚀单元方法为了计算结果精确，有限元网格必须足够小，导致计算非常耗时，更为严重的是，在强爆炸荷载下，尽管单元体积很小，仍然会有不少单元被删除，导致模型质量损失，从而使模拟结果偏离实际；等效塑性应力或应变方法则是一种近似的表示方法，无法真实反映裂隙特征，裂隙所在位置依旧被单元连接，造成模拟结果偏差，无法正真形成碎片；分离节点法虽然能够通过将处于同一位置的节点分离而轻易实现断裂破坏模拟，但其模拟过程的物理意义不够明确。
本文基于 LS-DYNA显示动力有限元分析软件，利用任意拉格朗日-欧拉(ALE)方法，以及多物质流固耦合方法对混凝土结构在爆炸荷载作用下的动态破坏过程进行研究。考虑到钢筋混凝土结构在爆炸荷载作用下的最大应变率可达10 s-1～100 s-1，为了更好分析混凝土结构在爆炸荷载作用下的动力响应，采用考虑应变率影响的钢筋材料本构和混凝土材料本构模型，并引入零厚度粘聚力单元来模拟钢筋混凝土结构的动态破坏过程，克服基于侵蚀单元算法带来的质量损失问题。文章首先介绍零厚度粘聚力单元模型的生成过程并对比试验结果，验证所建立的梁厚度粘聚力单元模型的合理性。其次，对比不同爆炸荷载下基于侵蚀算法以及零厚度粘聚力单元两种不同模拟方法的模拟结果，验证基于零厚度粘聚力单元模拟的优越性。最后基于零厚度粘聚力单元模型，分析探讨不同当量爆炸荷载对混凝土结构动态破坏过程以及碎片抛射的影响。
在 LS-DYNA中一般用线性多项式状态方程(*EOS_LINEAR_POLYNOMIAL)来描述爆炸后空气的性能[13]：
式中：P为空气压力；E为单位体积空气的初始能量；Ci(i=1～6)是系数；
。对于理想气体 C0=C1=C2=C3=C6=0，C4=C5=γ-1，可将式(1)简化为Gamma准则：
式中：
为相对密度； ρ 为当前空气密度；ρ0为初始时刻密度；γ为比热比；E为初始时刻单位体积的内能。对于理想气体ρ0=1.29 kg/m3，γ=1.4，E=2.5×105J/m3。
炸药在爆轰过程中压力与内能的关系用 JWL状态方程(*EOS_JWL)表示[13－15]：
式中：A、B为线性爆炸系数； ω、R1、R2为非线性爆炸系数；p为爆炸产生的压力；V为当压力为p时的体积与初始体积的比值； ψ为炸药的比内能。对于 TNT 炸药，A=3.712×1011Pa，B=3.231×109Pa，ω=0.30，R1=4.15，R2=0.95，eρ=1630 kg/m3，ψ=4.29×106J/kg。
钢筋采用塑性随动模型(*MAT_PLASTIC_KINEMATIC)[4,13,16]，该模型能同时模拟各向同性硬化、随动硬化以及混合硬化混合模型，通过在0(仅随动硬化)和1(仅各向同性硬化)之间调整参数β来实现；同时，该模型还可考虑材料以及应变率对材料强度的影响，通过设定有效塑性应变值FS以及采用Cowper-Symonds模型实现，如下式所示：
式中：σY为考虑应变率影响后的屈服强度；σ0为静载屈服强度；
为应变率；C 和 Pr为 Cowper-Symonds模型参数；
为有效塑性应变；Ep为塑性硬化模量。模型参数取自长安大学李天华[4]的博士论文，并根据孙文彬[17]所做试验进行了相应的修改，具体参数如表1所示。
表1 钢筋材料参数
Table 1 The material parameters for the reinforced bar
对于混凝土，受到爆炸冲击荷载作用时需要考虑大应变、高应变率和高围压下材料损伤实效的动态响应，而 JHC[13,18－20]模型是一种适用于高应变率、大变形下混凝土与岩石的材料模型，它与金属材料中应用广泛的Johnson-Cook材料模型相类似，其等效屈服强度是压力、应变率及损伤的函数，损伤量则是塑性体应变、等效塑性应变和压力的函数。JHC模型的等效屈服强度(标准化等效应力)如下式所示：
式中：
为标准化等效应力；σ为等效应力；
为单轴抗压强度；D 为损伤量(0≤D≤1)；
为标准化压力；
为无量纲应变率；
为参考应变率；
˙为等效应变率；A、B、N和C均是由试验确定得出的常数。
模型的损伤量是等效塑性应变、压力以及塑性体积应变的函数，如下式所示：
式中：ΔμP和ΔεP为等效塑性体积应变增量和等效塑性应变增量；
为常压下材料断裂时的塑性应变；
为标准化情况下的最大抗静水压力，Tc为最大抗静水压力；D1和D2为通过试验测得的损伤常数。模型参数取自张舵等[3]的论文，具体参数如表2所示。
表2 混凝土材料参数
Table 2 The material parameters for the concrete model
吴智敏等[21]研究了混凝土Ⅰ-Ⅱ复合型裂缝扩展准则，通过试验及数值模拟得到混凝土破坏过程中荷载-裂纹张开位移曲线，如图1所示。为使建立的零厚度粘聚力单元模拟能够反映混凝土材料的破坏特征，采用的粘聚力本构模型应符合混凝土破坏过程的荷载-位移曲线特征。考虑到LS-DYAN材料库中 186号材料模型*MAT_COHESIVE_GENERAL[13,22－23]的应力-位移曲线可以任意定义，同时粘聚力单元主要反映混凝土破坏过程的响应而对于破坏前的响应主要由实体单元模型控制，因此本文定义如图2所示的粘聚力单元本构模型来模拟混凝土的破坏特征。该模型在峰值应力前采用非常大的刚度使得实体单元之间的应力能够完全传递，而对于峰值应力之后的响应采用图1所示的混凝土峰后荷载-位移曲线进行模拟。对于采用的广义三维粘聚力本构模型其涉及的参数有：垂直于粘结面的拉或压应力tn、粘结面内两个互相垂直方向的剪应力ts1、ts2、垂直于粘结面的法相位移δ3以及粘结面内两个互相垂直的切向位移δ1、δ2；简单起见，图2中坐标做了归一化处理，即横坐标缩小δF倍，纵坐标缩小tmax倍，使其成为横纵坐标极限都为1的曲线。根据图2，可以得到曲线与坐标轴围成的面积、峰值应力、极限位移与断裂能的关系：
式中：G为材料断裂能，为Ⅰ型断裂能
或Ⅱ型断裂能
；tmax为峰值法向应力tnmax或峰值切向应力tsmax；δF为法向极限位移δnmax或切向极限位移δsmax；
，ATSLC为曲线与坐标轴围成面积；由此可得：
图1 荷载-裂纹张开位移曲线(摘自文献[18])
Fig.1 Load-crack separation curve (From reference [18])
图2 T-S曲线
Fig.2 Traction-separation curve
由图2可知在应力达到极限应力tmax之前，应力随位移的增大快速增大，当粘聚力达到临界值时，则开始进入应变软化阶段，零厚度粘聚力单元开始损伤，随着外力的进一步增加，单元相对位移增加，粘聚力则开始降低。牵引力-位移曲线具有不可逆性，即有加载/卸载发生时，卸载曲线为卸载点到原点的直线。
引入一个无纲量的有效分离参数λ，以此来表示切向相对位移(δ1，δ2)与法向相对位移(δ3)之间的关系：
式中：
和
分别为临界法向和切向位移分量，McCauley括号用来区分混凝土受压(δ3＜0)与受拉(δ3≥0)，由于混凝土抗压强度远远大于其抗拉强度，当受压时不考虑混凝土的应变软化，即：
本文采用的起裂准则的表达式为λ＞1。
本文的模拟方法是将零厚度粘聚力单元嵌入各实体单元边界上，通过零厚度粘聚力单元的破坏来实现混凝土材料裂隙萌生、扩展、汇合破坏过程模拟。为了实现这一过程，需要对初始有限元网格进行处理，图3为嵌入零厚度粘聚力单元前后的有限元网格以及裂缝单元的示意图，需要说明的是，插入的零厚度粘聚力单元在几何上厚度为零，为了示意粘聚力位置，在图中显示为带一定厚度的单元。具体插入过程如下[24]：
实体单元离散：将整个模型的单元信息与节点信息读入程序(图3(a))，通过对单元节点重新编排，使得每个单元都由不与别的单元共享的节点组成，如图3(b)所示。
粘聚力单元节点排序：将每个单元与其相邻单元通过坐标匹配办法搜索出两个重合的单元面，并将两个面上的节点按照LS-DYNA对于零厚度粘聚力单元节点排列顺序的要求进行排列，对于八节点零厚度粘聚力单元，要求八个节点的前四个属于一个面，后四个属于与之对应的另一个面，且前四个节点的排列顺序符合按右手螺旋法则判定的大拇指方向指向另一个面的要求，如图3(c)所示。
零厚度粘聚力单元生成：将重合面上的八个节点按Ⅱ中的排列顺序输出，并赋予新的单元编号以及单元类型，这样就完成了一个零厚度粘聚力单元的插入，如图3(d)所示。
图4为在实体边界插入零厚度粘聚力单元的效果图。
图3 零厚度粘聚力单元生成流程图
Fig.3 Flow chart of inserting zero thickness cohesive element
图4 插入零厚度粘聚力单元模型图
Fig.4 Sketch of numerical model with zero thickness cohesive elements
为了验证插入粘聚力单元模型的有效性，参考文献[25]进行一个简单的数值试验。如图5所示，该正方形板边长为 2w=0.8 m，厚度 D=0.01 m，在其中心位置有一条长度为2a=0.2 m的预制裂纹，在平行于裂纹的两端施加有均布的拉应力荷载σf=0.10125 MPa。
图5 数值试验模型
Fig.5 Numerical test model
采用零厚度粘聚力单元模型进行分析，模型的网格尺寸dp=0.01 m，在板中心裂纹位置不插入粘聚力单元，其余位置插入零厚度粘聚力单元进行模拟。图6给出了模拟结果的拉应力的分布云图，由图可知，模型很好的扑捉到了裂纹尖端的拉应力集中现象。为了进一步定量的反映所模拟结果的精确性，我们提取距裂纹尖端最近的五个粘聚力单元上的拉力，将它同基于线弹性力学的理解解进行对比，如图7所示。
图6 拉应力云图数值试验结果
Fig.6 Tensile stress contour of numerical test model
图7 距裂纹尖端不同位置处单元所受拉力
Fig.7 The tensile force at different locations from the crack tip
根据线弹性断裂力学(LEFM)，距裂纹尖端很近位置(r＜＜a)，其拉应力为：
式中：a为1/2裂纹长度；w为1/2板的边长；r为距离裂纹尖端长度；σn为距裂纹尖端r处的拉应力；C(a/w)为模型的几何函数：
取r=(2m-1)dp/2，dp为单元宽度，m为距离裂纹尖端的粘聚力单元个数，即每个粘聚力单元中心，则每个粘聚力单元所受到的拉力为：
由图7可得，模拟得到裂纹尖端附近的拉力基本与基于线弹性理论结果吻合，这说明所采用的零厚度粘聚力单元模型模拟断裂问题的有效性。
为了验证所插入粘聚单元模型计算的准确性，本文根据孙文彬[17]所做的实际试验模型(图8)，用LS-DYNA软件进行分析对比，模型的尺寸为1.3 m×1.0 m×0.12 m，钢筋采用双层双向布置，受力钢筋直径 10 mm，中心间距 100 mm，分布钢筋直径也为 10mm，中心间距 200mm，炸药量为2.0kg CompB炸药，相当于TNT当量为2.09 kg[17]。由于模型完全对称，为了节省计算时间，取1/2模型进行建模分析，如图9所示。文献[2]中指出网格质量对爆炸荷载时程曲线有很大的影响，且流体与固体的网格尺寸对流固耦合作用有较大影响，通过试验性模拟分析，最后选定炸药与空气的网格大小为4 mm，分别生成9828个和9213672个单元，都采用任意拉格朗日-欧拉(ALE)算法，ALE算法是先执行一个或几个拉格朗日时步计算，此时单元网格随流体流动产生变形，开始执行欧拉时步计算，既网格固定不动，物质在网格中流动，这种方法可以很好处理流体流动及大变形问题[26－27]；混凝土沿着荷载的方向划为12层，单元最小边长为10 mm，最大边长为 31 mm，共 62880个实体单元，采用Lagrange算法；LS-DYAN可以通过关键字*CONSTRAINED_LAGRANGE_IN_SOLID将ALE网格与Lagrange网格方便的耦合[13]，可以很好的模拟爆炸气体与混凝土结构的相互作用。通过2.1节建立的程序在混凝土实体单元边界插入零厚度粘聚力单元共 181792个，插入粘聚力单元的模型如图10所示；钢筋采用梁单元模拟，单元大小为10 mm，共生成2302个单元。为了验证基于粘聚力单元模型的优越性，对插入零厚度粘聚力单元之前的模型采用侵蚀算法进行对比分析。根据实际模型试验的边界条件，钢筋混凝土板为简支约束，采用将混凝土板自由搁置在两端的刚体单元上进行模拟。由于爆炸荷载的作用时间非常短暂，通常假定钢筋与混凝土的粘结是完好的，因此本文采用加速度和速度约束的方法[3]将钢筋单元与混凝土单元进行耦合，从而实现一种近似于固结的钢筋与混凝土的连接方式，避免了分离式钢筋混凝土建模时，钢筋与混凝土单元节点必须一致的要求，大大简化了建模过程。空气边界可能会造成爆炸冲击波的折射和反射，从而对冲击波产生增强或者减弱效果，为了避免这种现象产生，将空气边界定义为*BOUNDARY_PRESSURE_OUTFLOW[13]，使得压力在边界上流出，避免产生折射与反射。
图8 试验模型
Fig.8 Test model
图9 LS-DYNA模型
Fig.9 LS-DYNA model
图10 零厚度粘聚力单元模型
Fig.10 Zero thickness cohesive element model
图11、图12(a)和图12(e)分别为钢筋混凝土板在2.09 kg TNT当量爆炸荷载下的试验、侵蚀算法以及零厚度粘聚力单元模型模拟得到板正反表面破坏结果，图中深灰色(图12(a)～图12(d))以及黑色(图12(f)～图12(h))为已删除单元，可视为裂纹，蓝色为没有被删除的单元。通过比较发现，在2.09 kg TNT当量爆炸荷载下，两种方法都能给出基本与试验结果吻合的破坏结果。即在试件正面中心受爆炸荷载冲击，在中心四周产生同心环状的拉伸裂纹，试件正面中心产生波阵面为半球形的强压缩波，压缩波到达背面经背面反射后形成拉伸波，在试件背面形成辐射状的径向裂纹。两种方法的数值模拟结果其正面环形裂纹的直径都比试验结果稍大，这可能是由于试验环境与数值模拟环境不同所造成，试验在室内进行，冲击波会在室内墙壁上发生反射然后继续作用在混凝土板上，可能会对冲击波起到增强或减弱的作用，而数值模拟中则设置成完全自由边界的空气爆炸，不会出现冲击波的反射与折射。
图11 2.09 kg TNT当量爆炸荷载试验结果
Fig.11 Test results of 2.09 kg TNT equivalent explosive load
图12 不同当量爆炸荷载下混凝土板表面破坏模拟结果
Fig.12 Simulation results of concrete slab surface failure under different equivalent explosive loads
为了探索两种算法在高当量爆炸荷载作用下的破坏模拟效果，对上述混凝土板分别采用5.0 kg、7.5 kg、10.0 kg TNT爆炸荷载进行分析。侵蚀算法结果如图12(b)～图12(d)所示，由图12(b)发现当炸药的质量增加至5 kg时，板正面环状裂纹最小直径明显减小，裂纹密度增加，由背面扩展至正面的径向裂纹增多，且在中心出现直径约为0.13 m圆形完全破碎区，板背面中心裂纹密度增大，径向裂纹大部分并未扩展至板边缘，支座处产生密集裂纹；通过图12(c)发现，当炸药当量增至7.5 kg时，其正面最小环状裂纹的直径与 TNT=5.0 kg时接近，但其密度增大，且位于板中心的完全破碎区域直径增至约 0.5 m，板背面裂纹形成一个较破碎的直径约为1 m的圆形区域。通过图12(d)发现，当炸药当量增至 10.0 kg时，板正面表面基本完全破坏，背面也产生密集裂纹。由于侵蚀算法通过删除实体单元来模拟破坏，在低爆炸荷载下单元破坏不多时其也能获得较好的效果。但随着爆炸荷载的增加，有大量实体被单元删除导致气体逃逸，从而影响模拟结果的准确性。图13为统计得到的在不同当量TNT炸药荷载作用下基于侵蚀算法所造成的质量损失占总质量百分比随时间的变化曲线，由图可以看出四条曲线在初始阶段以较大斜率上升，并且很快达到拐点，此后基本为水平状，这说明钢筋混凝土在爆炸荷载的作用下，混凝土的破坏在短时间内完成，由连接碎片的钢筋约束，这与实际情况相符[5]。由图14可以更明显的发现，随着爆炸荷载的增加，混凝土板的质量损失也在增加，当炸药量达到10 kg时，删除单元的质量占总质量的比例高达 46%，这说明在高爆炸荷载下采用侵蚀算法模拟混凝土在爆炸荷载作用下的破坏并不适用。与其对应的是基于零厚度粘聚力单元的模拟结果，如图12(f)～图12(h)所示，其裂纹分布规律与侵蚀算法结果在爆炸荷载较低时相近，但在高爆炸荷载下，由于引入零厚度粘聚力单元算法通过粘聚力单元的删除来实现混凝土材料在爆炸荷载作用下的破坏，且零厚度粘聚力单元本身没有厚度，因此不会如侵蚀算法一样造成整个板面破坏，其模拟结果裂纹更加真实。同时由于零厚度粘聚力单元没有质量，不存在质量损失，这使得零厚度粘聚力单元模型不仅能够很好的模拟混凝土材料在低爆炸荷载下的破坏过程，同时也能很好地避免侵蚀算法在高爆炸荷载作用时的质量损失问题。
图13 不同当量爆炸荷载侵蚀算法质量损失百分比随时间变化图
Fig.13 Time history of mass loss percentage for different equivalent explosive load
图14 不同当量爆炸荷载侵蚀算法质量损失百分比
Fig.14 Percentage of mass loss of different equivalent explosive load by erosion algorithm
龚顺风等[5]采用侵蚀算法模拟了密闭空间中爆炸荷载对爆炸碎片的抛射速度的影响，发现由于侵蚀算法造成的质量减小导致其预测的碎片速度值大于理论计算值。本文为了探索零厚度粘聚力单元模型对混凝土爆炸碎片抛射的模拟效果，首先对自由空气爆炸荷载下混凝土板碎片的抛射速度进行理论推导。对于自由空气爆炸，爆炸冲击波传递至钢筋混凝土板的能量转换为钢筋混凝土的断裂能以及碎片的动能，如下式所示：
式中：Ws为冲击波传递到板的能量；Uk为钢筋混凝土碎块动能；Uc为混凝土断裂所消耗能量。由此可得：
式中：v为碎片的平均抛射速度；M为板的总质量。Ws根据Xu[28]等提出的计算公式计算：
式中：I为板的动量；Ps为超压峰值； τ 为正压持续时间；A为板的迎爆面面积。其中 Ps根据Henrych[29]给出的经验公式计算：
式中：
为折合距离；r炸药中心距结构的实际距离；m 为炸药质量。τ根据 Wu[30]所得到经验公式计算：
混凝土断裂的能量则由 Xu[28]等给出的公式计算：
式中：t为混凝土板厚度；fyd为钢筋动态屈服强度；
y d为钢筋动态屈服应变；Rr钢筋体积占板总体积分数；fcd为钢筋动态屈服强度；
c d为混凝土动态屈服应变；Rc混凝土体积占板总体积分数。
根据上述平均速度计算方法，最终计算得到理论值。将混凝土板在不同爆炸荷载作用下的碎片平均抛射速度采用零厚度粘聚力单元模型进行了预测，并对比了相应的基于侵蚀算法的预测结果与理论计算结果，如图15所示。由图15可以发现采用在低爆炸荷载作用下，基于零厚度粘聚力单元预测的碎片抛射速度与基于侵蚀算法得到的抛射速度都与理论值接近，但是随着爆炸荷载的增加，基于粘聚力单元预测的速度较基于侵蚀算法得到的抛射速度更接近理论值。由此可以推断，基于零厚度粘聚力单元模型模拟预测的碎片抛射速度更加精确。
图15 不同炸药当量爆炸荷载碎片平均抛射速度对比
Fig.15 Comparison of debris average launching velocity of RC slab under different equivalent explosive load
为了进一步探索不同爆炸荷载对上述混凝土板的破坏过程以及破坏模式的影响，采用零厚度粘聚力单元模型对上述混凝土板在不同爆炸荷载下的破坏过程以及破坏模式进行了分析，并与采用侵蚀算法模拟结果比较，如图16和图17所示。
图16 不同当量爆炸荷载下混凝土板破坏过程模拟结果
Fig.16 Failure process of the concrete plate under different equivalent explosive load
图17 不同当量爆炸荷载下混凝土板最终破坏模拟结果
Fig.17 Failure results of the concrete plate under different equivalent explosive load
由图16(a)可知，当爆炸荷载为 TNT=2.09 kg时，在时间=0.2 ms时，在板内部靠近背面产生平行于板的连续的圆面裂纹，直径约为 0.2 m，在时间= 0.3 ms时，产生多条环状裂纹，最内部的直径约为 0.7 m，最后在背面产生径向裂纹，随着裂纹的扩展，最终混凝土被切成许多独立的碎块，局部位于板背面中心处碎块脱落，但大部分碎块仍被钢筋连接成一个整体，钢筋混凝土板向下弯曲，钢筋并未发生断裂，如图17(a)所示，最终统计共有19403个粘聚力单元删除，根据文献中对于爆炸破坏的分类，起初产生层状裂纹具有明显的震塌破坏特征[6]，但由于钢筋的约束，最终转变为整体弯曲破坏，这也与试验结果吻合[17]；通过图16(b)发现，当爆炸荷载增至TNT=5.0 kg时，在时间=0.2 ms时钢筋混凝土板首先产生直径约为0.7 m的环状裂纹，在时间= 0.3 ms时，环状裂纹增多且向板背面扩展，并且产生平行于板的层状裂纹，最终在板内产生密集裂纹，最终统计共有53422个粘聚力单元删除，由17(b)发现在支座处有较明显的剪切破坏；通过图16(c)发现，当炸药当量继续增加至TNT=7.5 kg时，在时间=0.2 ms时，混凝土板迎爆面就已产生比较密集的辐射状的裂纹，随后在板内产生密集裂纹，最终统计共有66630个粘聚力单元删除，由图17(c)可知最终在支座处发生彻底剪切破坏，且整个板在荷载方向有较大位移；通过图16(d)发现，当TNT=10.0 kg时，在时间=0.2 ms时，混凝土板迎爆面就已形成一个直径约为1m的完全破碎区，随后在板内产生密集裂纹，最终统计共有 102231个粘聚力单元删除，由图16(d)可知，最终混凝土板在中间与支座处发生明显的冲切破坏与剪切破坏。
由图16可知，在2.09 kg及5.0 kg TNT当量爆炸荷载下，两种算法结果的最初破坏位置以及最终破坏状态都比较接近，且由图17发现板的最终位移也相差不大。但在7.5 kg及10.0 kg TNT当量爆炸荷载下，侵蚀算法最初在混凝土板正面及背面产生环状裂纹，而零厚度粘聚力模型则在正面产生环状及辐射状裂纹，并且基于侵蚀算法位移分析明显大于零厚度粘聚力单元模型分析结果，这是由于在高爆炸荷载下，侵蚀算法质量损失过大造成的。
本文基于LS-DYNA动力有限元程序，利用任意拉格朗日–欧拉(ALE)方法，以及多物质流固耦合方法，建立爆炸气体与结构耦合的三维数值模型，提出了建立零厚度粘聚力单元模型的算法并验证了该模型模拟钢筋混凝土结构在爆炸荷载作用下破裂破坏的可行性，比较了采用侵蚀算法和零厚度粘聚力单元模型模拟钢筋混凝土在不同爆炸荷载作用下的破坏过程以及碎片抛射速度，并基于零厚度粘聚力单元模型模拟结果研究了不同当量爆炸荷载对钢筋混凝土板破坏过程以及破坏模式的影响，结果表明：
(1)在低爆炸荷载下，无论是侵蚀算法还是零厚度粘聚力单元模型都能很好的模拟钢筋混凝土在爆炸荷载作用下的破坏过程，但是随着爆炸荷载的增加，侵蚀算法删除实体单元所造成的质量损失逐渐增大，从而影响计算结果的准确性，而零厚度粘聚力单元本身没有厚度与质量，因此能很好地避免侵蚀算法在高爆炸荷载作用时的质量损失问题。
(2)通过将两种模拟方法预测的不同爆炸荷载作用下钢筋混凝土板的碎片平均抛射速度与理论计算结果比较，发现基于零厚度粘聚力单元预测的碎片抛射速度更接近理论值，由此推断基于零厚度粘聚力单元模型预测的碎片抛射速度更加合理。
(3)随着炸药当量的增加，不仅钢筋混凝土板的破坏更加迅速，其破坏方式也在发生改变；随着爆炸荷载的增加，混凝土板的破坏模式由低爆炸荷载时的整体弯曲破坏，逐渐转变为支座处的剪切破坏，再到支座处剪切破坏与中心处冲切破坏混合的破坏模式。
参考文献：
[1]龚顺风, 邓欢, 朱升波, 等. 近爆作用下钢筋混凝土板动态破坏的数值模拟研究[J]. 振动与冲击, 2012,31(2): 20―24.Gong Shunfeng, Deng Huan, Zhu Shengbo, et al.Numerical simulation for dynamic failure of a reinforced concrete slab under close-in explosion [J]. Journal of Vibration & Shock, 2012, 31(2): 20―24. (in Chinese)
[2]Tai Y S, Chu T L, Hu H T, et al. Dynamic response of a reinforced concrete slab subjected to air blast load [J].Theoretical & Applied Fracture Mechanics, 2011, 56(3):140―147.
[3]张舵, 卢芳云, 王瑞峰. 钢筋混凝土板在爆炸作用下的破坏研究[J]. 弹道学报, 2008, 20(2): 13―16.Zhang Duo, Lu Fangyun, Wang Ruifeng. Study on destruction of reinforced concrete slab under blast loading [J]. Journal of Ballistics, 2008, 20(2): 13―16.(in Chinese)
[4]李天华. 爆炸荷载下钢筋混凝土板的动态响应及损伤评估[D]. 西安: 长安大学, 2012.Li Tianhua. Dynamic response and damage assessment of reinforced concrete slabs subjected to blast loading[D]. Xi’an: Chang’an University, 2012. (in Chinese)
[5]龚顺风, 金伟良. 内部爆炸荷载作用下钢筋混凝土板碎片抛射速度的预测[J]. 工程力学, 2009, 26(9):225―230.Gong Shunfeng, Jin Weiliang. Prediction of debris launch velocity of reinforced concrete slab subjected to internal blast loading [J]. Engineering Mechanics, 2009, 26(9):225―230. (in Chinese)
[6]张想柏, 杨秀敏, 陈肇元, 等. 接触爆炸钢筋混凝土板的震塌效应[J]. 清华大学学报: 自然科学版, 2006,46(6): 765―768.Zhang Xiangbai, Yang Xiumin, Chen Zhaoyuan, et al.Explosion spalling of reinforced concrete slabs with contact detonations [J]. Journal of Tsinghua University,2006, 46(6):765―768. (in Chinese)
[7]武海军, 黄风雷, 付跃升, 等. 钢筋混凝土中爆炸破坏效应数值模拟分析[J]. 北京理工大学学报, 2007,27(3): 200―204.Wu Haijun, Huang Fenglei, Fu Yuesheng, et al.Numerical simulation of reinforced concrete breakage under internal blast loading [J]. Beijing Ligong Daxue Xuebao/ transaction of Beijing Institute of Technology,2007, 27(3): 200―204. (in Chinese)
[8]赖建中, 朱耀勇, 谭剑敏. 超高性能混凝土在埋置炸药下的抗爆试验及数值模拟[J]. 工程力学, 2016,33(5): 193―199.Lai Jianzhong, Zhu Yaoyong, Tan Jianmin. Experiment and simulation of ultra-high performance concrete subjected to blast by embedded explosive [J].Engineering Mechanics, 2016: 33(5): 193―199. (in Chinese)
[9]Yao S, Zhang D, Chen X, et al. Experimental and numerical study on the dynamic response of RC slabs under blast loading [J]. Engineering Failure Analysis,2016, 66: 120―129.
[10]阎石, 张亮, 王丹. 钢筋混凝土板在爆炸荷载作用下的破坏模式分析[J]. 沈阳建筑大学学报自然科学版,2005, 21(3): 177―180.Yan Shi, Zhang Liang, Wang Dan. Failure mode analysis for rc slab under explosive loads [J]. Journal of Shenyang Architectural & Civil Engineering Institute, 2005, 21(3):177―180. (in Chinese)
[11]Zhao C F, Chen J Y. Damage mechanism and mode of square reinforced concrete slab subjected to blast loading[J]. Theoretical and Applied Fracture Mechanics, 2013,63/64(1): 54―62.
[12]张晓天, 贾光辉, 黄海. 基于节点分离 Lagrange有限元方法的超高速碰撞碎片云数值模拟[J]. 爆炸与冲击,2010, 30(5): 499―504.Zhang Xiaotian, Jia Guanghui, Huang Hai. Simulation of hypervelocity-impact debris clouds using a Lagrange FEM with node separation [J]. Explosion & Shock Waves, 2010, 30(5): 499―504. (in Chinese)
[13]LS-DYNA keyword user’s manual, version 970 [M].Livermore, California: Livermore Software Technology Corporation, April, 2007.
[14]Lee E L, Hornig H C, Kury J W. Adiabatic expansion of high explosive detonation products [R]. Livermore,USA: Livermore Lawrence Radiation Laboratory, 1968.
[15]Baker W E. Explosions in air [M]. Austin, TX, US:University of Texas Press, 1973.
[16]Krieg R D, Key S W. Implementation of a time independent plasticity theory into structural computer programs [J]. Constitutive equations in viscoplasticity:Computational and Engineering Aspects, 1976:125―137.
[17]孙文彬. 钢筋混凝土板的爆炸荷载试验研究[J]. 辽宁工程技术大学学报, 2009, 28(2): 217―220.Sun Wenbin. Experimental studies on reinforced concrete(RC) slabs subjected to blast loads [J]. Journal of Liaoning Technical University, 2009, 28(2): 217―220.(in Chinese)
[18]Holqmuist T J. A computational constitutive model for concrete subjected to large strains, high s train rate, and high pressures [C]. 14th International Symposium on Ballistics. Canada: BIC, 1993: 591―600.
[19]孙其然, 李芮宇, 赵亚运, 等. HJC模型模拟钢筋混凝土侵彻实验的参数研究[J]. 工程力学, 2016, 33(8):248―256.Sun Qiran, LI Ruiyu, Zhao Yayun, et al. Investigation on parameters of HJC model applied to simulate perforation experiments of reinforced concrete [J]. Engineering Mechanics, 2016, 33(8): 248―256. (in Chinese)
[20]方秦, 孔祥振, 吴昊, 等. 岩石 Holmquist-Johnson-Cook 模型参数的确定方法[J]. 工程力学, 2014, 31(3):197―204.Fang Qin, Kong Xiangzhen, Wu Hao, et al.Determination of the Holmquist-Johnson-Cook consitiutive model parameters of rock [J]. Engineering Mechanics, 2014, 31(3): 197―204. (in Chinese)
[21]吴智敏, 董伟, 许青. 混凝土Ⅰ-Ⅱ复合型裂缝扩展准则及扩展全过程的数值模拟[J]. 水利学报, 2009,40(2): 1453―1459.Wu Zhimin, Dong Wei, Xu Qing. Propagation criterion for mixed mode I-II crack propagation and numerical simulation of the whole propagation process [J]. Journal of Hydraulic Engineering, 2009, 40(2): 1453―1459. (in Chinese)
[22]Tvergaard V, Hutchinson J W. On the toughness of ductile adhesive joints [J]. Journal of the Mechanics and Physics of Solids, 1996, 44(5): 789―800.
[23]Tvergaard V, Hutchinson J W. Effect of strain-dependent cohesive zone model on predictions of crack growth resistance [J]. International Journal of Solids and Structures, 1996, 33(20/21/22): 3297―3308.
[24]徐海滨, 杜修力, 杨贞军. 基于预插黏性界面单元的Koyna重力坝强震破坏过程分析[J]. 振动与冲击,2014, 33(17): 74―79.Xu Haibin, Du Xiuli, Yang Zhenjun. Seismic failure analysis of Koyna gravity dam using cohesive interface elements [J]. Journal of Vibration & Shock, 2014,33(17): 74―79. (in Chinese)
[25]Kazerani T. Micromechanical study of rock fracture and fragmentation under dynamic loads using discrete element method [D]. Switzerland: École Polytechnique Fédérale De Lausanne, 2011.
[26]年鑫哲, 张耀, 孙传怀, 等. 空气冲击波作用于柔性防爆墙的透射和绕射效应分析[J]. 工程力学, 2015,32(3): 241―248.Nian Xinzhe, Zhang Yao, Sun Chuanhuai, et al. Analysis of transmission and diffraction effects of air shock waves upon flexible explosion-proof walls [J]. Engineering Mechanics, 2015: 32(3): 241―248. (in Chinese)
[27]Basting S, Quaini A, Glowinski R. Extended ALE method for fluid-structure interaction problems with large structural displacements [J]. Journal of Computational Physics, 2017, 331(C): 312―336.
[28]Xu K, Lu Y. Debris velocity of concrete structures subjected to explosive loading [J]. International Journal for Numerical and Analytical Methods in Geomechanics,2006, 30(9): 917―926.
[29]Henrych J, Abrahamson G R. The dynamics of explosion and its use [M]. New York, Elsevier Scientific Publishing Company, 1979.
[30]Wu C, Hao H. Modeling of simultaneous ground shock and airblast pressure on nearby structures from surface explosions [J]. International Journal of Impact Engineering, 2005, 31(6): 699―717.
DYNAMIC FAILURE ANALYSIS OF REINFORCED CONCRETE SLAB BASED ON COHESIVE ELEMENT UNDER EXPLOSIVE LOAD
WU Zhi-jun1,2, ZHANG Peng-lin1,2, LIU Quan-sheng1,2, LI Wan-feng3, JIANG Wei-zhong3
(1. School of Civil Engineering, Wuhan University, Wuhan 430072, China;2. The Key Laboratory of Safety for Geotechnical and Structural Engineering of Hubei Province, Wuhan 430072, China;3. Huainan Mining Industry (Group) CO, LTD. Anhui, Huainan 232001, China)
Abstract:For arsenal or dangerous goods warehouse, accidental explosion may happen. Since the warehouse is mainly made of reinforced concrete materials, it is therefore essential to study the failure process of reinforced concrete under explosive loading. In this study, based on the LS-DYNA software, the paper studies the dynamic failure process of reinforced concrete under explosive loading, which was investigated by using the Arbitrary Lagrange Euler (ALE) method and the multi-materials Fluid-Solid coupling method. In order to analyze the dynamic response of concrete structure under explosive load, strain rate dependent constitutive models were adopted for the rebar and concrete. In addition, zero thickness cohesive element was introduced to simulate the dynamic breakup process of concrete structure, which can overcome the mass loss problem caused by the erosion algorithm. In this paper, the generation process of the cohesive element model was firstly introduced. The developed model was then verified by comparing the numerical results with the experimental results. To demonstrate the superior capability of the proposed cohesive element model, the simulation results from the proposed model were compared with those obtained from erosion methods. Finally, based on the cohesive element models, the influence of different explosive loads on the dynamic failure process of concrete structure and the debris launching were analyzed.
Key words:explosive load; reinforced concrete; zero thickness cohesive element; numerical simulation;LS-DYNA
中图分类号：TU352; O383
文献标志码：A
doi:10.6052/j.issn.1000-4750.2017.04.0270
文章编号：1000-4750(2018)08-0079-12
收稿日期：2017-04-06；修改日期：2017-08-14
基金项目：国家973项目(2014CB046900)；国家自然科学基金青年基金项目(41502283)
通讯作者：刘泉声(1962―)，男，江苏人，教授，博士，博导，主要从事深部工程围岩稳定性控制研究(E-mail: liuqs@whu.edu.cn).
作者简介：
吴志军(1983―)，男，浙江人，教授，博士，博导，混凝土结构爆炸灾害和岩石破坏数值方法研究(E-mail: wuzh0013@e.ntu.edu.sg);
张鹏林(1994―)，男，甘肃人，博士生，主要从事混凝土结构爆炸灾害分析研究(E-mail:zhang_penglin@outlook.com);
李万峰(1962―)，男，河南人，高工，学士，从事采矿工程研究(E-mail: liwanfeng@sohu.com);
江维中(1973―)，男，安徽人，工程师，从事采矿工程研究(E-mail: 642278380@qq.com).