1.国核电力规划设计研究院有限公司 北京 100095；2.山东电力工程咨询院有限公司 山东济南 250000
摘要：本次计算的目的是为了模拟岩体在爆炸冲击荷载工况下的损伤断裂，最终的结果是得到裂纹的发展情况。采用的工具是国际通用的大型有限元分析软件ABAQUS。ABAQUS在非线性、爆炸模拟方面具有强大的功能[1，2]，并且在ABAQUS的计算实例中给出了一些相关的爆炸模拟例子[3，4]，所以采用ABAQUS作为这次运算的工具是足够而充分的，最终的数值结果将与ABAQUS自带的爆炸计算实例进行对比分析来确定本次计算的合理性。
关键词：爆炸冲击；ABAQUS；裂纹；数值模拟
1 引言
由于数值计算在爆炸冲击下岩石材料的裂纹产生必须具有普遍性，所以首先要解决的问题便是几何模型，由于几何模型必须具有一定的实际意义，所以在参考了前人实验数据后，决定采用如下的几何模型。
这一几何模型是一般爆炸试验的标准尺寸，并且利用这一模型所做的爆炸试验结果完全可以用于对照本次计算，来检验ABAQUS模拟情况的合理性。
需要注意的是：虽然实际的爆炸试验具有一定的参考价值，但本文参考的文献是混凝土爆炸的模型尺寸，而我们所需要作的是的岩石的爆炸模型，在本质上二者具有很大的差别。
2．ABAQUS算例概况
算例中存在2个同心的圆环管，在2个环形管的间隔中塞满高速炸药（HE）。内钢管的半径为10mm，外钢管的半径为20mm，两个钢管的厚度均为2mm。两个钢管都采用6个元素来模拟发射方向，但在HE区域则采用24个元素模拟爆炸辐射方向。钢管材料为理想弹塑性的，
它的杨氏模量为221.6GPa，泊松比为0.279，屈服强度为430MPa，密度为7846kg/m3。爆炸时，HE 采用JWL 状态方程来计算应力及能量参数。爆炸时爆炸波的速率为7596m/s，参数A=520.6GPa，B=5.3GPa，R1=4.1，R2=1.2，ω=0.35，HE密度=1900kg/m3，初始具体能量=3.63J/kg。控制管内液体拉力为认为是0。
在这一分析的过程中可以采用两种方法解决问题，第一种采用平面形式，利用4节点简缩积分二次单元方法。第二种采用立体形式，利用8节点简缩积分二次单元方法[5]。
具体每一步骤的参数设置此处不做过多赘述
- Property部分：
这一部分需要对于刚性套管柱与HE（高速炸药）的材料性质进行定义。由于材料的内部与外部均为钢，所以这一部分的定义采用前述所介绍的钢管参数，同样HE也采用前面介绍的炸药参数。
需要补充的是，在利用ABAQUS计算爆炸是需要在钢管的内外周边设立4个爆炸点，而我们所选取的部分为整个圆环的八分之一，在我们这一部分上只存在一个爆炸点。在材料分配上我们将两种性质的材料分配给钢管与HE。
我们在采用平面计算时，爆炸点坐标为（1.41421，1.41421，0），即在为钢管与HE接触的最边缘位置。
- Step 部分：
在这一次的计算中我们共将动力学分析部（ABAQUS /Explicit）分为两部分，第一部分的时间长度为6μs，第二部分的时间长度为1.5μs。
- Interaction部分：
在立方体的8节点减缩积分单元计算实例中需要定义接触面，在interaction property的定义中选为“hardcontact”，而在平面问题中则不需要定义接触。
- Load 部分：
定义两个边界条件，两个边界条件分别加在扇形的两条直边上，约束类型选择YSYMM。
2.1 算例一
2.1.1 几何模型与材料参数
整个几何模型在空间上沿X、Y轴均对称，故选取整体试件（500×500×500mm）的1/4作为研究对象。在一角的位置预留半径为20mm的1/4圆弧作为爆炸的填放炸药位置。利用part中的切割功能建立好三维可变形的实体模型。材料性质的定义尤为重要，因为爆炸时的受力情况比较特殊，材料在load 步骤中并没有为载荷的作用，爆炸产生的应力是由炸药点火后才产生的，而定义爆炸材料的参数是在property下完成的。根据前文的ABAQUS计算实例，输入参数，确定爆炸点，在这一计算中爆炸点共有两个，分别在炸药的中间部位与试件接触的位置。
试件的材质定义存在着这样的问题，为模拟爆炸作用下的试件开裂情况，在定义试件材料时，需给出给损伤与断裂参数，而爆炸问题又属于动力学问题，在使用ABAQUS中explicit模块下是不能计算断裂的，这就限制了本次计算。同时这一问题还用等效的方法进行模拟，就是在试件材料定义时加入粘弹性参数，并利用ABAQUS中“删除破坏单元”的显示方法来实现模拟断裂的目的。
2.1.2 边界条件及接触定义
定义接触与边界条件决定了计算模拟的准确性，由于在爆炸数值计算参考文献：较少，无法找到较为合理的接触方式，所以决定采用最为保守的方式定义接触[6]，即接触方式为“Hard Contact”。定义岩石的边界面为接触的主平面，爆炸面为副平面，两者间的接触没有摩擦，并且在发生作用后两个面不能够分离，而边界条件可根据实际模型简化为四组，分别为：第一组是作用在如图5，图6所示的两个平面上的YSYMM方式约束；第二组是作用在如图5、图6所示的两个表面上的XSYMM方式约束；第三组是作用在整个模型底面上的完全约束，这是为了模拟真实情况中底端固定的作用；第四组是作用在如图5、图6所示的两个爆炸上下表面竖向约束。
2.1.3 网格划分及结果
网格的划分见图1、图2（Mises等效应力分布图），从图2中可以看到应力沿放射线传播，在这里解释一下爆炸时应力为何会以波的形式传播。炸药爆炸时，爆炸物质几乎是瞬时转变为高温，高压的爆炸气体。爆炸气体膨胀很快，并挤压周围空气，且占其容积。这样一层层压缩空气在爆炸气体前沿形成和发展，这就是爆炸波。爆炸初期介质受到冲击荷载，后期介质受到静压作用。在介质中引起的应力、应变和位移主要以波动的形式传播。爆炸气体中的压力逐步减小到等于大气压，然后爆炸波不再放出爆炸气体支持而继续独立传播。出于惯性，爆炸气体质点继续运动，它的压力下降到低于大气压力时，又处于周围气体的压力高，爆炸气体逐步停止并往回运动。知道它的压力又逐渐增加到稍大于大气压力，以后又膨胀，如此往复。这说明了爆炸应力是以波的形式传播的。
图2 Mises等效应力示意图
在临近炸药实体的位置应力为0，因为在发生单元破裂或失效时，单元自动删除，所以我们看到的炸药周围应力为0的部位便是裂纹出现的位置。爆炸属于大位移变形，应该分析对数应变的分布情况。图3-5以下为对数应变各分量的云图。
图5 对数第一主应变云图
2.2 算例二
这一次的计算采用了模拟整个岩石的空间尺寸的数值计算，模型仍为500×500×500mm。在其中间部位预留直径为40mm的装药孔，另外我们在整个的炸药内设立了4个爆炸点，分别将它们的爆炸时间作一定的延迟。
在材料方面我们加入了粘弹性参数，在动力学问题中不能够加入损伤断裂，所以要通过加粘弹性的方法来模拟岩石质材料的力学性能。
边界条件方面，加上了两个方向轴压。在真实情况下模拟一个岩石的装药爆炸过程并不具有实际的工程意义，因为岩石的爆炸实在实际的矿山工程上应用的，岩石的周围必定具有多向的轴压，所以我们在利用ABAQUS计算爆炸的过程中，人为地加上了两个侧向的轴压，这样计算的岩石爆炸问题便可以作为实际工程的参考。模型边界条件如图12-图14所示：
不采用三轴的模拟压实作用，是因为三轴压缩时，炸药在爆炸后，我们很难得到一个收敛的结果。而缩短时间步骤后计算虽然收敛，可位移场，应力场的变化都较小，无法满足出现裂纹的要求。
在爆炸刚刚开始的几百毫秒时间内，在于爆炸接触的岩石部位及利状况较为均匀，这是因为爆炸在点火的瞬间以后炸药爆炸，爆炸轻体以及应力波向四周传播，这时岩石的整体应力状态比较均匀，如图15所示。其后一段过程，对于岩石来讲是应力增加的过程，我们在计算时加上了删除破坏单元的命令，所以在岩石单元的盈利达到我们给定的参数时，此单元体的应力状态显示为0，而它的位移变量仍然可以显示。
由于采用了这样的程序，我们在利用ABAQUS分析爆炸的过程中我们发现了一个较为特殊的情况，那就是在爆炸结束后，岩石的应力状况达到破坏后，岩石的裂纹仍然在开展。
在爆炸理论中可以解释为：岩体内的最初裂纹是由于炸药在瞬间的的爆炸，产生了高温高压，这一过程只是在几毫秒甚至几十微秒的时间内完成，而在此后爆炸气体渗入裂纹并在静压的作用下，使原始形成的裂缝进一步扩展。所以可以用这一理论解释我们数字模拟过程中产生的这一现象。另外在计算的最后阶段，大位移变形从岩石的边缘出现，最后于岩石的中间大位移变形连为一体，形成开裂。
在位移场的变形情况中可以看到，模型的上表面存在着几条位移较大的位移带，同时再比较模型的应变场，发现在这些部位我们划分的单元失效了，所以可以肯定在这些位移带上出现了裂纹，要强调的是由于应力场上基本所有区域都为失效部位，我们很难找出哪一位置或哪一条带是裂纹出现的部位，仍需在模型的建立分析上下功夫，使得能够在应力场中可以清楚地看到应力失效单元，而在位移场上又可以看到在应力场上失效的部位而位移值又是最大的，这样可以得到裂纹带一定是发生在该区域。
2.3 算例三
最后为了得到裂纹的开展，经过几百次的试算后决定在约束，单元变形控制上增加一些参数。
首先是约束，这一次是采用了较灵活的Displacement/Rotation方式约束，只限制法向的平动，这样可以使岩石单元发生活动余地更大，计算时的结果可以很好地满足收敛要求。在step分析步骤中，我们考虑到爆炸在热力学上属于绝热的过程，在本步中勾选了adiabatic选项。而本次计算的最大改进点在于采用了单元控制法，我们将如下属性赋到炸药与岩石体中：
*section Controls，name=EC-1，DISTORTION CONTROL，hourglass=relax stiffness，
second order accuracy=YES
通过这些命令来限制定义单元的变形和刚度矩阵。计算得到的应力云如下图18、图19所示：
图19 位移云图
3．结论
本次计算过程中，划分单元节点数量受到了一定限制，但仍可以定性地得到裂纹产生的合理部位。裂纹主要是以一条或几条主裂纹的形式开始扩展，然后贯穿整个试件，与实验的结果较为吻合。本次计算还加入了失效单元的删除命令，来模拟岩石材料的失效，这一点上克服了在explicit模块中不能够添加断裂参数的缺陷，为以后的计算提供了一条新的途径。
参考文献：
[1]王勖成著.有限单元法.
[2]庄茁. ABAQUS有限元软件6.4 版入门指南. 北京：清华大学出版社，2004-09.
[3]庄茁，张帆，岑松. ABAQUS非线性有限元分析与实例. 北京：科学出版社，2005-03.
[4]石亦平，周玉蓉. ABAQUS有限元分析实例详解. 北京：机械工业出版社，2006-07.
[5]庄茁. ABAQUS/Standard 有限元软件入门指南. 北京：清华大学出版社，1999.
[6]杨曼娟. Abaqus 用户材料子程序开发及应用.
论文作者:段修斌1,朱磊磊2
论文发表刊物:《防护工程》2019年15期
论文发表时间:2019/11/29
标签：应力论文; 岩石论文; 裂纹论文; 炸药论文; 位移论文; 这一论文; 单元论文; 《防护工程》2019年15期论文;