1. 引言
海洋关系着中华民族崛起的安全和发展利益，我国的核电工程建设、航道疏通、水坝拆除、桥梁施工、港口码头建设、水下炸礁、油气田勘探等建设项目也日益增多。考虑到水下钻孔爆破具有爆炸能量利用率好、岩石破碎均匀、工效高等优点，因此上述工程的前期基础开挖时大多采用水下钻孔爆破的方法进行施工。但水下钻孔爆破作业伴随产生的安全问题也较为突出，如水中爆破冲击波、爆破振动、涌浪、飞石、噪声等次生危害容易给爆破区域周边的生态环境及水生物等造成不同程度的安全威胁。尤其是水中爆破冲击波具有冲击作用特性且幅值大、传播远，是影响水下爆破安全的主要危害效应。水中大量渔业生物在受到水中爆破冲击波作用时不仅导致形态学改变、器官损伤，甚至可能导致大面积的死亡。因此，进行水下爆破工程下水中冲击波的传播规律研究，对指导水下爆破施工设计及安全防护有重要意义。
在水下爆破理论研究领域，S. Temkin [1] 通过对小药包水中爆炸冲击波衰减特性的研究。P. Cole [2] 深入研究了水中爆炸的物理现象与基本规律，推导了水中爆破冲击波压力峰值、比冲量等基本参数的计算公式。在水下爆破数值模拟研究领域，M. Kamegai等[3] 采用CARE程序对水下爆破冲击波产生的水面非规则稀疏反射现象进行了数值模拟分析。M.B. Liu等[4] 采用率先采用SPH方法对炸药爆轰、产物与水介质的相互作用等问题进行了数值分析。
本文以水下深孔爆破工程为例，利用ANSYS/LS-DYNA有限元程序对水下不耦合深孔装药台阶爆破的冲击波传播过程进行数值模拟，分析研究了水下深孔爆破工况下，冲击波在岩石介质与水介质中的传播规律，为爆破设计和研究提供了理论和技术上的支持。
2. LS-DYNA程序简介
LS-DYNA程序是一个通用显式非线性动力分析程序，可以求解各种二维、三维非线性结构的高速碰撞、爆炸和金属成型等接触非线性、冲击荷载非线性和材料非线性问题。DYNA程序系列最初是1976年在美国Lawrence Livermore National Lab由J.O. Hallquist 主持开发完成的。LS-DYNA具有几何非线性、材料非线性、摩擦和接触分离状态非线性等程序，以Lagrange算法为主，兼有ALE、Euler和SPH算法，主要进行结构非线性动力分析，同时可以进行热分析、静力分析。是目前使用最为广泛的有限元显示非线性动力分析软件[5] 。
3. 水下深孔爆破的数值计算模型
3.1. 材料模型
LS-DYNA材料库中提供了多种材料模型可对炸药、金属、土壤、流体、岩石、混凝土等材料进行描述。炸药及爆轰产物的材料模型采用关键字为MAT_HIGH_EXPLOSIVE_BURN的高能燃烧模型，状态方程采用JWL方程[6] ，其形式为：
(1)
式中，A、B、R1、R2和ω为材料参数，p为压力，V为爆轰产物的相对体积，E0为爆轰产物的初始化内能。本文选取2号岩石炸药，其参数如表1。
岩石介质采用Holmquist-Johnson-Cook (HJC)本构模型，在LS-DYNA中的关键字为MAT_ JOHNSON_HOLMQUIST_CONCRETE。HJC模型是Holmquist、Johnson和Cook提出的可以描述工程和地质材料在大应变、高变率下的非线性变形及断裂特性的损伤本构模型。其强度模型以特征化等效应力描述，如图1，表达式为[7] ：
(2)
式中：
为特征化等效应力(
为材料所能达到的最大特征化等效应力)，σ为实际等效应力，fc为准静态单轴抗压强度；D为损伤变量。
为特征化压力，其中p为单元内的静水压力，
为特征化应变率，
为实际应变率，
为参考应变率。C为应变率影响参数，A、B、N和Smax统称为极限面参数，其中A为特征化黏性强度系数，B为特征化压力硬化系数，N为压力硬化系数。本文选取岩石参数如表2。
水是流体介质，关键字为MAT_NULL的材料模型，状态方程采用GRUNEISEN状态方程，其形式如下[8] ：
(3)
式中：p为压力，C为剪切—压缩波速曲线截距，S1、S2、S3为剪切—压缩波速曲线斜率系数，γ0为GRUNEISEN常数，E0为初始化内能，ρ、ρ0为材料密度和材料初始密度。
3.2. 计算模型
为了考察水下深孔爆破的炸药冲击波在岩石介质与水介质中的传播规律，拟建立二维模型参数如下：装药为不耦合装药，炮孔直径14 cm，装药直径12 cm，孔深1200 cm，填塞长度200 cm，台阶高度1000 cm，
Table 1. Parameters for 2# rock explosive
表1. 2号岩石炸药的材料参数
表2. 岩石材料HJC模型的材料参数
超深200 cm，最小抵抗线250 cm，起爆方式为中心单点起爆；水面距离台阶顶层1000 cm，距台阶底层2000 cm；二维模型的长高分别为5000 cm、3000 cm，如图2。
计算模型利用SOLID164六面体单元对装药、空气、水、岩石进行网格划分。计算方法采用ALE算法，装药、水、空气为ALE网格，共节点连接。岩石介质为Lagrange网格，与ALE网格通过关键字CONSTRAINED_LAGRANGE_IN_SOLID进行耦合。由于计算模型采用三维六面体单元，建立单层网格模型，通过在模型的前后面施加约束实现对二维模型的计算效果。
4. 水下深孔爆破的数值模拟结果
4.1. 数值模拟结果
装药起爆方式为中心起爆，起爆时刻为800 us。随着装药的起爆与爆轰波的传播，爆炸能量传入岩石介质，岩石中冲击波也呈近似椭圆形状从中心向两端传播。在起爆点水平方向上，冲击波率先抵达岩
石介质与水介质交界面，爆炸能量开始传入水介质中，形成水中冲击波，并在岩石介质中形成反射波。图3为水下深孔爆破数值模拟不同时刻的应力云图。
4.2. 水中冲击波传播规律分析
装药起爆后，冲击波在岩石介质中呈近椭圆面传播，冲击波波震面依次达到岩石与水交界面，并由
此开始向水中传播。图4为以起爆点为中心，沿角度等距离取测点示意图。图5为各测点上的压力时程曲线。
由图5可知，冲击波抵达各测点的时间大致相同，其中C点的峰值压力最大。这是由于装药起爆后
Figure 4. The schematic of shockwave measuring points
图4. 冲击波测点示意图
冲击波在岩石介质中传播，波阵面与岩石水介质交界面呈一定角度传入水中，交界面质点依次受到岩石冲击波作用发生振动，冲击波以此为中心向水中传入，沿交界面一定角度的方向上，先传入水中的冲击波与后传入水中的冲击波叠加使峰值压力增加。
深孔爆破采用200 cm水介质填塞，为分析装药爆炸后的能量从炮孔处溢出多少，在炮孔及距离炮孔100 cm处选取测点如图6，图7为炮孔处各测点的压力时程曲线，图8为距离炮孔100 cm处各测点的压力时程曲线。
Figure 6. The schematic of shockwave measuring points
图6. 冲击波测点示意图
(a) 左侧测点 
(b) 中间测点 
(c) 右侧测点
Figure 7. Time history curve of pressure on blasting hole
图7. 炮孔处测点的压力时程曲线
(a) 左侧测点 
(b) 中间测点 
(c) 右侧测点
Figure 8. Time history curve of pressure 100 cm from blasting hole
图8. 距炮孔100 cm处测点的压力时程曲线
由图7、图8分析可知，在炮孔处，靠近临空面一侧的峰值压力较大，炮孔中心处及远离临空面处的峰值压力较小；距离炮孔100 cm处的三个测点的冲击波参数基本一致。这说明200 cm的水介质炮孔填塞起到了充分的填塞作用，装药爆炸产生的冲击波主要是在岩石介质内传播，且爆炸能量主要传递给靠临空面一侧的岩石介质，并从该侧岩石介质传入水介质。
5. 结语
利用ANSYS/LS-DYNA有限元程序水下爆破进行数值模拟是非常有效的。本文通过对水下深孔爆破的数值模拟，预测了在实施爆破过程中冲击波在岩石和水介质中的传播规律。在今后的水下爆破实施之前，可以先通过ANSYS/LS-DYNA有限元程序数值模拟预计爆破效果和冲击波效应是否符合预先的爆破安全设计从而进步改进实施方案，对工程实践有重要的指导意义。
基金项目
国家自然科学基金资助项目(51304218, 51339006, 51208508)。