基于ALE 算法的连续水射流破煤特征数值模拟研究

2024-02-12郝从猛陈学习刘清泉金霏阳

煤矿安全 2024年1期

郝从猛，陈学习，刘清泉，金霏阳

（1.中国安全生产科学研究院，北京 100012；2.华北科技学院安全监管学院，河北三河 101601；3.中国矿业大学安全工程学院，江苏徐州 221116）

庞大的煤炭产能有力支撑了我国经济的高速发展，同时旺盛的煤炭需求使得我国矿井的开采正在以平均每年10～30 m 的速度向深部延伸，许多矿井的开采深度已超过800 m，特别是山东、安徽、黑龙江、河南等地区多对矿井已超过1 000 m[1-2]。随着开采深度的不断增加，煤层地应力、瓦斯压力和含量急剧增大，而渗透率却大幅降低，此“三高一低”的状况使得越来越多的非突矿井正在转变为突出、强突出矿井，煤矿企业面临的安全生产形势愈发严峻[3]。作为煤炭生产、消费大国，我国每年投入巨额经费治理煤矿瓦斯，尽管成绩显著，但瓦斯事故依然是引起煤矿重大伤亡的主要原因之一[4]。

实践表明，瓦斯抽采是解决煤矿瓦斯事故，实现节能环保的最有效措施[5]。但普遍存在的低渗煤层使得瓦斯抽采一直面临着抽采效率低、难度大、治理成本高等问题[6]；特别是进入深部以后，高地应力作用使得瓦斯抽采瓶颈日益显著。为了提高煤层渗透性，广大的科研工作者们提出了很多的增透措施和技术；其中，水力冲孔技术由于操作简单、卸压效果显著而得到了广泛的应用，相关研究成果十分丰富。刘明举等[7]、王兆丰等[8]、陶云奇等[9]通过现场考察、实验研究获得了水射流破煤条件和增透效果，证明煤体卸压增透效果和冲孔半径成正相关；ZHANG 等[10]结合应用开发了高压水射流钻冲一体化造穴装备，在新景矿试验成功后进行了大规模推广；葛兆龙等[11]、高亚斌等[12]通过试验研究了水射流冲击性能，得到水射流的破煤细观演化机制。梳理文献可以发现，目前水射流破煤方面的研究主要集中在破煤机理及应用优化上，对于水射流破煤特征也多是对非淹没射流条件下进行了分析。

然而，在一些俯孔条件下进行水射流破煤时往往出现孔内积水不能及时排出，进而形成淹没环境下射流破煤的情况[4]。由于淹没环境下水射流能量损耗与非淹没环境差距较大，导致两者破煤性能具有显著的差别，而目前的研究对两者破煤特性认识不足，特别是对淹没射流破煤效率及其最佳冲击深度研究不足，这也导致在淹没环境下开展水射流破煤造穴技术的应用存在一定的盲目性。为此，采用数值模拟的方法对淹没射流/非淹没射流的破煤特征进行对比分析，进一步揭示不同环境介质下水射流破煤特性的差异及其产生机制，并对淹没环境下水射流破煤时效性及可行性展开分析。

1 基于ALE 的连续水射流破煤特性数值模型

1.1 模拟软件及流程

多物质流固耦合法（ALE）是一种任拉格朗日-欧拉方法，可同时实现单物质和多物质的任意耦合，可以较好地实现连续射流冲击的数值计算方法。ALE 算法兼顾了Lagrange 和Euler 方法两者的优点，不仅具备Lagrange 方法在有效跟踪物质在结构边界上运动的优点，还具有Euler 方法在内部网格与物质实体之间相互独立的强项[13]。ALE算法不仅确保了计算过程中的网格单元的完整性，同时缩短了复杂问题的求解时间并提高了精度，特别适用于冲击作用下大变形问题的处理。

LS-Dyna 模拟软件可实现多种算法，并且拥有丰富的材料库和材料本构模型，可以提供50 多种材料间的接触方式，适合分析高速碰撞、侵彻等冲击问题。因此，采用LS-Dyna 进行水射流冲击破煤模拟，结合Ls-Prepost 软件的前后处理功能，进行建模、修改参数、结果显示、结果输出等。数值模拟流程如图1。

1.2 几何模型及材料参数

由于需要对非淹没/淹没环境下的破煤特性进行对比分析，因此几何模型需要建立不同的射流域，即煤体分别处于空气域和水域中。对于对称问题，为了提高计算效率可以进行适当的精简取其一部分进行分析。高压水射流破煤过程就是对称的，因此取模型的1/4 进行模拟计算，在结果分析时按对称面呈现，非淹没/淹没射流破煤数值计算1/4 模型如图2。左图为空气域（非淹没）模型，右图为水域（淹没）模型，煤体的尺寸都为100 mm×100 mm×100 mm。

图2 非淹没/淹没射流破煤数值计算1/4 模型Fig.2 Numerical calculation of 1/4 model of coal breaking by unsubmerged/submerged water jet

煤体本构模型采用高应变率、大变形下的岩石材料的Johnaon-Holmquist-Concrete（J-H-C）模型，根据煤层特点并结合相关文献[4,14]，煤体选取J-HC 模型，研究表明该模型适合模拟水射流冲击下煤体发生损伤的标准，即以最大剪切应变和最大拉伸应变为判断依据的破坏条件；流体（水和空气）选取MAT-NULL 模型。

煤体材料模型（H-J-C）参数如下：①煤体密度ρc：1.4 g/cm3；②剪切模量G：1.48 GPa；③特征化黏性强度A：0.79；④特征化压力硬化因子B：1.6；⑤应变率影响系数C：0.007；⑥压力硬化指数N：0.61；⑦无侧限单轴抗压强度Fc：4.8×10-3GPa；⑧抗拉强度T：4×10-4GPa；⑨累计塑性应变EFmin：0.001；⑩特征化最大强度SFmax：7；⑪压溃点的压力pc：-1.6×10-4GPa；⑫压溃点的体积应变Uc：0.001；⑬压实点的压力pL：0.008 GPa；⑭压实点的体积应变UL：0.1；⑮损伤常数D1：0.04；⑯损伤常数D2：1.0；⑰压力常数K1：0.85 GPa；⑱压力常数K2：-1.71 GPa；⑲压力常数K3：2.08 GPa；⑳无量纲系数：3.57。

流体材料模型（NULL）空气参数如下：①空气密度ρa：1.9×10-3g/cm3；②空气压力截止pca：-1×10-11GPa；③空气动力黏度系数Ca：0.033 1；④无量纲系数S1a：0；⑤无量纲系数S2a：0；⑥无量纲系数S3a：0；⑦状态方程系数GAMA0a：0；⑧一阶体积更正系数αa：0；⑨空气初始内能Ea：0。

流体材料模型（NULL）水参数如下：①水密度ρw：1.0 g/cm3；②水压力截止pcw：-1×10-5GPa；③水动力黏度系数Cw：0.148；④无量纲系数S1w：2.56；⑤无量纲系数S2w：1.986；⑥无量纲系数S3w：0.226；⑦状态方程系数GAMA0w：0.45；⑧一阶体积更正系数αw：0；⑨水初始内能Ew：0。

1.3 控制方程

模型中流体的连续性方程描述如下[15]:

式中：ρ为流体介质的密度，kg/m3；t为时间；vi为i方向的速度分量，m/s；wi为i方向的相对速度的分量，m/s；xi为i方向的位移分量，m。

射流的动量方程如下：

式中：v为射流速度，m/s；σij,j为应力张量在j方向的偏导数；bi为i方向的体积力，N；xj为j方向的位移分量，m。

射流的能量方程为：

式中：E为流体的内能密度，J/m3；σi,j为在i方向和j方向体应力矢量和；vi,j为i方向和j方向的速度矢量和；wj为j方向的相对速度的分量，m/s。

通过与状态方程联用而为压力行为提供应力组件，从而为材料提供应力张量。对空气域选取*EOD-LINEAR-POLYNOMIAL 状态方程来提供压力，将空气域视为理想气体后可得状态方程为[14]：

式中：p为压力；γ为空气的比热比，γ=CP/CV；CP、CV分别为定压比热容和定容比热容；ρ0为初始密度。

对于水域，采用*EOS-GRUNEISEN 状态方程来提供压力p，表达式为：

式中：Cb为冲击波波速与单元速度关系曲线的纵轴截距；S1、S2、S3为冲击波波速与单元速度关系曲线的斜率系数；γ0为Gruneisen 系数；μ为泊松比；α为对 γ0的一阶体积修正系数。

煤的本构模型选用H-J-C 模型，模型的规范化等效应力描述如下：

式中：A为无量纲黏性力；B为规范化压缩系数；p*为无量纲压力；为无量纲应变率；C为材料的应变系数；N为硬化系数；D为损伤因子。

式中：Δεp为等效塑性应变增量，无量纲；Δµp为等效体积应变增量，无量纲；D1（p*+T*）D2为常压p下材料断裂时的塑性应变，无量纲；T*为最大拉伸静水压力，MPa；D1、D2为损伤常数。

2 连续水射流破煤特性数值模拟结果

2.1 射流能量耗散分析

通过分析淹没/非淹没状态下水射流的传播速度可以初步获得两者的传播特征和能量耗散情况。设定射流直径为0.3 mm，射流速度为170 m/s，喷嘴与煤壁作用面的间距都为20 mm，水射流在非淹没和淹没条件下的传播情况如图3。

图3 非淹没/淹没水射流传播特征Fig.3 Propagation characteristics of unsubmerged/submerged water jet

从图3 中可以看出：非淹没和淹没条件下水射流自喷嘴喷出后两者的传播过程差别巨大；首先，从水射流形态上看，非淹没射流传播过程中呈圆柱状向前传播，射流直径基本上没有变化，而淹没射流的传播过程中的前端由于水域的阻力作用逐渐向后分离，同时水射流内部受挤压而产生应变，这是水射流周围部分水速低于射流速度后导致的滞后效应，射流受阻是导致淹没条件下水射流的传播距离较小的根本原因；其次，射流在传播速度上的明显差异，非淹没射流的传播距离和传播时间成正比，而相同时间淹没射流的传播距离大幅下降；根据射流到达煤壁的时间，淹没条件下水射流传播10 mm 的距离所用时间是非淹没条件下的1.6 倍，传播20 mm 的距离所用时间是非淹没条件下的2.18 倍，这说明随着传播距离的增大淹没条件下的水射流要不断消耗自身的质量来抵消环境的阻力，也就是说淹没条件下需要水射流持续不断地喷出才能保证射流向前传播，否则射流就会在环境水域的阻力消耗下消失。

另外，根据水射流传播的时间变化可以推测，随着射流喷嘴与煤体间距的增大所需射流的持续时间也会变长。

2.2 淹没/非淹没射流破煤特征分析

2.2.1 水射流破煤过程分析

为了进一步研究射流环境对破煤的影响，建立射流直径为3 mm、靶距为20 mm 的破煤工况，冲击计算时间1 000 μs，射流速度分别为170、190、210 m/s 时非淹没/淹没射流破煤过程如图4。

图4 非淹没/淹没射流破煤过程Fig.4 Coal breaking process of unsubmerged/submerged water jet

从图4 中可以看出：非淹没/淹没条件下水射流破煤特性存在着明显的区别；首先，3 组非淹没射流从喷嘴喷出后到达煤壁的时间不同；另外，非淹没射流模拟试验的煤体全部被穿透，水射流速度为170、190、210 m/s 时煤体穿透的时间分别是970、900、830 μs，而淹没条件下3 组试验全部没有穿透煤体，但随着射流速度的增大破煤深度逐渐增加；最后，除破煤深部不同外，2 种射流的破煤基坑直径的差别明显，在淹没射流模型中，随着射流时间的延长，侵蚀坑呈漏斗状逐渐向下延伸，但是非淹没射流的破煤侵蚀坑始终保持圆柱状增长，这可能是由于在水域环境下射流周围产生强烈的湍流作用使得射流的动量不断向周围水分子传递，进而二次冲击使得破煤范围增大，同时，湍流作用下的返流动量会很大程度上减弱射流自身的冲击能力，使得射流向前冲击煤体的能力变弱，最终形成漏斗形的基坑。

2.2.2 煤体应力分布特征分析

非淹没/淹没射流破煤时煤体等效应力分布如图5。

图5 非淹没/淹没射流破煤时煤体等效应力分布Fig.5 Distribution of equivalent stress of coal mass during unsubmerged/submerged water jet breaking coal

从图5 可以看出：2 种射流模式下的最大应力非常接近。这主要是因为应力的峰值主要是由于煤体自身强度决定的，与射流模式关系不大，但是在应力分布上差别明显；非淹没条件下的应力主要在基坑的前方，沿横向上分布范围较小；而淹没条件下应力不仅在冲击位置分布范围较大，在煤体的底部也出现了应力的升高，这主要是淹没射流挤压煤体的范围较大，而压强达到煤体应力极限后才会出现损伤，从而导致煤体整体受力较大，而底部由于是固定边界，因此会出现较高的应力集中现象。

2.2.3 破煤效果分析

为了获得不同冲击模式下水射流破煤效率，对比分析淹没/非淹没射流破碎煤体基坑深度随时间的变化关系，统计了2 种射流破煤模式下每间隔20 μs 的侵蚀深度变化，非淹没/淹没射流破煤深度变化如图6。

图6 非淹没/淹没射流破煤深度变化Fig.6 Coal breaking depth variation of unsubmerged/submerged water jet

从图6 可以看出：2 种冲击模式下侵蚀深度都随着冲击时间的增加不断增大，区别在于非淹没射流的侵蚀深度基本与冲击时间成正比，并且在冲击计算结束之前穿透煤体。为了便于对比，按照规律将非淹没射流的侵蚀深度延深至1 000 μs，此时射流速度为170、190、210 m/s 的侵蚀深度分别达到103、112、121 mm，对应冲击时间时下淹没射流的侵蚀深度分别为21、25、30 mm，分别为非淹没射流侵蚀深度的20.4%、22.35%、24.8%。另外，淹没条件下的侵蚀深度呈现阶梯增加，并且增加的速度远小于非淹没条件下增加的速度。

根据淹没条件下侵蚀深度的变化趋势可以看出，随着冲击时间的延长，侵蚀深度的增加速率逐渐减小，这主要是随着射流喷嘴距离煤体逐渐增大，射流到达煤体需要克服的行程阻力逐渐增加。可以推测随着时间的增加射流冲击煤体最终会停止，即一定速度下的射流通过能量不断衰减达到煤体表面后的动能不足以破坏煤体，所形成的侵蚀坑深度就是极限破煤深度。

由于非淹没条件下3 种射流速度的煤体均被穿透，虽然在图6 中进行了预测补充，但却不能对2 种射流模式下的破煤效果进行准确的分析，因此以210 m/s 穿透时的冲击时间为上限，即对比冲击发生后至830 μs 之间的侵蚀深度和煤体损伤单元，非淹没/淹没射流冲击830 μs 时的破煤效果对比如图7。

图7 非淹没/淹没射流冲击830 μs 时的破煤效果对比Fig.7 Comparison of coal breaking effect of unsubmerged/submerged water jet impacting 830 μs

从图7 可以看出：在冲击了830 μs 时，非淹没条件下射流速度为170、190、210 m/s 时的侵蚀深度分别为83、91、100 mm，单元损伤数量分别为6 640、7 360、7 996 个；而对应速度下淹没射流的侵蚀深度分别为16、20、25 mm，单元损伤数量分别为7 366、8 246、9 730 个，侵蚀深度分别为非淹没射流的19.2%、21.9%、25%，单元损伤数量分别为非淹没射流的110.9%、112%、121%。

对比2 种冲击模式下煤体的侵蚀深度和损伤单元的数量发现，不同射流速度下的破煤深度与时间基本成正线性关系，根据2.1 节的分析可知，无论是淹没环境还是非淹没环境，水射流的传播速度是确定的，但淹没射流需要更多的能量克服介质水的阻力而使其传播速度降低。另外，2 种冲击模式下的损伤单元虽然随射流速度的增加而增大，但是规律性较差，特别是淹没环境下，损伤单元由于射流的湍流作用使得破煤面积大大增加，同时也导致了破煤深度大幅下降，但整体损伤单元的数量在增大。

2.3 淹没/非淹没射流破煤深度关系时变特征

为了研究淹没射流的破煤特性，需要先获得淹没条件对水射流破煤的弱化作用。以非淹没射流的侵蚀深度为基础，求出对应冲击时间下淹没射流侵蚀深度和非淹没射流侵蚀深度的关系，淹没/非淹没射流破煤深度时变关系如图8。

从图8 可以看出：在冲击的早期淹没射流的侵蚀坑深度与非淹没射流的侵蚀坑深度的比值迅速增大，这是由于在淹没射流发生冲击破煤的时候非淹没射流冲击已经持续了一段时间，因此两者的比值是从0 开始增加；之后，随着非淹没射流的侵蚀深度不断增大，两者比值也迅速增加，但3 种射流速度下的侵蚀深度变化规律一致，都是先迅速增大，之后趋于平稳。这里需要注意的是趋于平稳的是淹没射流的侵蚀深度与非淹没射流侵蚀深度的比值，即随着冲击时间的延长侵蚀深度都在增大，但两者的比值却趋于平稳。根据图中拟合曲线可以推测淹没射流的破煤深度最大约等于非淹没射流的21%，因此根据这一结果可以结合已知非淹没射流的破煤深度求得相应时间下的淹没射流破煤深度。

根据研究结果可知，淹没环境下水射流破煤效果较差，但在工程应用中通过延长冲刷时间或射流速度可以相应的提高破煤深度。另外，由于淹没环境下水射流破煤效果难以考察，可以基于研究结果对非淹没环境下的效果进行分析，初步评价淹没射流的可行性。