基于ALE 算法的連續(xù)水射流破煤特征數(shù)值模擬研究

郝從猛 ，陳學(xué)習(xí) ，劉清泉 ，金霏陽

（1.中國(guó)安全生產(chǎn)科學(xué)研究院，北京 100012；2.華北科技學(xué)院 安全監(jiān)管學(xué)院，河北 三河 101601；3.中國(guó)礦業(yè)大學(xué) 安全工程學(xué)院，江蘇 徐州 221116）

龐大的煤炭產(chǎn)能有力支撐了我國(guó)經(jīng)濟(jì)的高速發(fā)展，同時(shí)旺盛的煤炭需求使得我國(guó)礦井的開采正在以平均每年10～30 m 的速度向深部延伸，許多礦井的開采深度已超過800 m，特別是山東、安徽、黑龍江、河南等地區(qū)多對(duì)礦井已超過1 000 m[1-2]。隨著開采深度的不斷增加，煤層地應(yīng)力、瓦斯壓力和含量急劇增大，而滲透率卻大幅降低，此“三高一低”的狀況使得越來越多的非突礦井正在轉(zhuǎn)變?yōu)橥怀觥?qiáng)突出礦井，煤礦企業(yè)面臨的安全生產(chǎn)形勢(shì)愈發(fā)嚴(yán)峻[3]。作為煤炭生產(chǎn)、消費(fèi)大國(guó)，我國(guó)每年投入巨額經(jīng)費(fèi)治理煤礦瓦斯，盡管成績(jī)顯著，但瓦斯事故依然是引起煤礦重大傷亡的主要原因之一[4]。

實(shí)踐表明，瓦斯抽采是解決煤礦瓦斯事故，實(shí)現(xiàn)節(jié)能環(huán)保的最有效措施[5]。但普遍存在的低滲煤層使得瓦斯抽采一直面臨著抽采效率低、難度大、治理成本高等問題[6]；特別是進(jìn)入深部以后，高地應(yīng)力作用使得瓦斯抽采瓶頸日益顯著。為了提高煤層滲透性，廣大的科研工作者們提出了很多的增透措施和技術(shù)；其中，水力沖孔技術(shù)由于操作簡(jiǎn)單、卸壓效果顯著而得到了廣泛的應(yīng)用，相關(guān)研究成果十分豐富。劉明舉等[7]、王兆豐等[8]、陶云奇等[9]通過現(xiàn)場(chǎng)考察、實(shí)驗(yàn)研究獲得了水射流破煤條件和增透效果，證明煤體卸壓增透效果和沖孔半徑成正相關(guān)；ZHANG 等[10]結(jié)合應(yīng)用開發(fā)了高壓水射流鉆沖一體化造穴裝備，在新景礦試驗(yàn)成功后進(jìn)行了大規(guī)模推廣；葛兆龍等[11]、高亞斌等[12]通過試驗(yàn)研究了水射流沖擊性能，得到水射流的破煤細(xì)觀演化機(jī)制。梳理文獻(xiàn)可以發(fā)現(xiàn)，目前水射流破煤方面的研究主要集中在破煤機(jī)理及應(yīng)用優(yōu)化上，對(duì)于水射流破煤特征也多是對(duì)非淹沒射流條件下進(jìn)行了分析。

然而，在一些俯孔條件下進(jìn)行水射流破煤時(shí)往往出現(xiàn)孔內(nèi)積水不能及時(shí)排出，進(jìn)而形成淹沒環(huán)境下射流破煤的情況[4]。由于淹沒環(huán)境下水射流能量損耗與非淹沒環(huán)境差距較大，導(dǎo)致兩者破煤性能具有顯著的差別，而目前的研究對(duì)兩者破煤特性認(rèn)識(shí)不足，特別是對(duì)淹沒射流破煤效率及其最佳沖擊深度研究不足，這也導(dǎo)致在淹沒環(huán)境下開展水射流破煤造穴技術(shù)的應(yīng)用存在一定的盲目性。為此，采用數(shù)值模擬的方法對(duì)淹沒射流/非淹沒射流的破煤特征進(jìn)行對(duì)比分析，進(jìn)一步揭示不同環(huán)境介質(zhì)下水射流破煤特性的差異及其產(chǎn)生機(jī)制，并對(duì)淹沒環(huán)境下水射流破煤時(shí)效性及可行性展開分析。

1 基于ALE 的連續(xù)水射流破煤特性數(shù)值模型

1.1 模擬軟件及流程

多物質(zhì)流固耦合法（ALE）是一種任拉格朗日-歐拉方法，可同時(shí)實(shí)現(xiàn)單物質(zhì)和多物質(zhì)的任意耦合，可以較好地實(shí)現(xiàn)連續(xù)射流沖擊的數(shù)值計(jì)算方法。ALE 算法兼顧了Lagrange 和Euler 方法兩者的優(yōu)點(diǎn)，不僅具備Lagrange 方法在有效跟蹤物質(zhì)在結(jié)構(gòu)邊界上運(yùn)動(dòng)的優(yōu)點(diǎn)，還具有Euler 方法在內(nèi)部網(wǎng)格與物質(zhì)實(shí)體之間相互獨(dú)立的強(qiáng)項(xiàng)[13]。ALE算法不僅確保了計(jì)算過程中的網(wǎng)格單元的完整性，同時(shí)縮短了復(fù)雜問題的求解時(shí)間并提高了精度，特別適用于沖擊作用下大變形問題的處理。

LS-Dyna 模擬軟件可實(shí)現(xiàn)多種算法，并且擁有豐富的材料庫和材料本構(gòu)模型，可以提供50 多種材料間的接觸方式，適合分析高速碰撞、侵徹等沖擊問題。因此，采用LS-Dyna 進(jìn)行水射流沖擊破煤模擬，結(jié)合Ls-Prepost 軟件的前后處理功能，進(jìn)行建模、修改參數(shù)、結(jié)果顯示、結(jié)果輸出等。數(shù)值模擬流程如圖1。

1.2 幾何模型及材料參數(shù)

由于需要對(duì)非淹沒/淹沒環(huán)境下的破煤特性進(jìn)行對(duì)比分析，因此幾何模型需要建立不同的射流域，即煤體分別處于空氣域和水域中。對(duì)于對(duì)稱問題，為了提高計(jì)算效率可以進(jìn)行適當(dāng)?shù)木?jiǎn)取其一部分進(jìn)行分析。高壓水射流破煤過程就是對(duì)稱的，因此取模型的1/4 進(jìn)行模擬計(jì)算，在結(jié)果分析時(shí)按對(duì)稱面呈現(xiàn)，非淹沒/淹沒射流破煤數(shù)值計(jì)算1/4 模型如圖2。左圖為空氣域（非淹沒）模型，右圖為水域（淹沒）模型，煤體的尺寸都為100 mm×100 mm×100 mm。

圖2 非淹沒/淹沒射流破煤數(shù)值計(jì)算1/4 模型Fig.2 Numerical calculation of 1/4 model of coal breaking by unsubmerged/submerged water jet

煤體本構(gòu)模型采用高應(yīng)變率、大變形下的巖石材料的Johnaon-Holmquist-Concrete（J-H-C）模型，根據(jù)煤層特點(diǎn)并結(jié)合相關(guān)文獻(xiàn)[4,14]，煤體選取J-HC 模型，研究表明該模型適合模擬水射流沖擊下煤體發(fā)生損傷的標(biāo)準(zhǔn)，即以最大剪切應(yīng)變和最大拉伸應(yīng)變?yōu)榕袛嘁罁?jù)的破壞條件；流體（水和空氣）選取MAT-NULL 模型。

煤體材料模型（H-J-C）參數(shù)如下：①煤體密度ρc：1.4 g/cm3；②剪切模量G：1.48 GPa；③特征化黏性強(qiáng)度A：0.79；④特征化壓力硬化因子B：1.6；⑤應(yīng)變率影響系數(shù)C：0.007；⑥壓力硬化指數(shù)N：0.61；⑦無側(cè)限單軸抗壓強(qiáng)度Fc：4.8×10-3GPa；⑧抗拉強(qiáng)度T：4×10-4GPa；⑨累計(jì)塑性應(yīng)變EFmin：0.001；⑩特征化最大強(qiáng)度SFmax：7；?壓潰點(diǎn)的壓力pc：-1.6×10-4GPa；?壓潰點(diǎn)的體積應(yīng)變Uc：0.001；?壓實(shí)點(diǎn)的壓力pL：0.008 GPa；?壓實(shí)點(diǎn)的體積應(yīng)變UL：0.1；?損傷常數(shù)D1：0.04；?損傷常數(shù)D2：1.0；?壓力常數(shù)K1：0.85 GPa；?壓力常數(shù)K2：-1.71 GPa；?壓力常數(shù)K3：2.08 GPa；?無量綱系數(shù)：3.57。

流體材料模型（NULL）空氣參數(shù)如下：①空氣密度ρa(bǔ)：1.9×10-3g/cm3；②空氣壓力截止pca：-1×10-11GPa；③空氣動(dòng)力黏度系數(shù)Ca：0.033 1；④無量綱系數(shù)S1a：0；⑤無量綱系數(shù)S2a：0；⑥無量綱系數(shù)S3a：0；⑦狀態(tài)方程系數(shù)GAMA0a：0；⑧一階體積更正系數(shù)αa：0；⑨空氣初始內(nèi)能Ea：0。

流體材料模型（NULL）水參數(shù)如下：①水密度ρw：1.0 g/cm3；②水壓力截止pcw：-1×10-5GPa；③水動(dòng)力黏度系數(shù)Cw：0.148；④無量綱系數(shù)S1w：2.56；⑤無量綱系數(shù)S2w：1.986；⑥無量綱系數(shù)S3w：0.226；⑦狀態(tài)方程系數(shù)GAMA0w：0.45；⑧一階體積更正系數(shù)αw：0；⑨水初始內(nèi)能Ew：0。

1.3 控制方程

模型中流體的連續(xù)性方程描述如下[15]:

式中：ρ為流體介質(zhì)的密度，kg/m3；t為時(shí)間；vi為i方向的速度分量，m/s；wi為i方向的相對(duì)速度的分量，m/s；xi為i方向的位移分量，m。

射流的動(dòng)量方程如下：

式中：v為射流速度，m/s；σij,j為應(yīng)力張量在j方向的偏導(dǎo)數(shù)；bi為i方向的體積力，N；xj為j方向的位移分量，m。

射流的能量方程為：

式中：E為流體的內(nèi)能密度，J/m3；σi,j為在i方向和j方向體應(yīng)力矢量和；vi,j為i方向和j方向的速度矢量和；wj為j方向的相對(duì)速度的分量，m/s。

通過與狀態(tài)方程聯(lián)用而為壓力行為提供應(yīng)力組件，從而為材料提供應(yīng)力張量。對(duì)空氣域選取*EOD-LINEAR-POLYNOMIAL 狀態(tài)方程來提供壓力，將空氣域視為理想氣體后可得狀態(tài)方程為[14]：

式中：p為壓力；γ為空氣的比熱比，γ=CP/CV；CP、CV分別為定壓比熱容和定容比熱容；ρ0為初始密度。

對(duì)于水域，采用*EOS-GRUNEISEN 狀態(tài)方程來提供壓力p，表達(dá)式為：

式中：Cb為沖擊波波速與單元速度關(guān)系曲線的縱軸截距；S1、S2、S3為沖擊波波速與單元速度關(guān)系曲線的斜率系數(shù)；γ0為Gruneisen 系數(shù)；μ為泊松比；α為對(duì) γ0的一階體積修正系數(shù)。

煤的本構(gòu)模型選用H-J-C 模型，模型的規(guī)范化等效應(yīng)力描述如下：

式中：A為無量綱黏性力；B為規(guī)范化壓縮系數(shù)；p*為無量綱壓力；為無量綱應(yīng)變率；C為材料的應(yīng)變系數(shù)；N為硬化系數(shù)；D為損傷因子。

式中：Δεp為等效塑性應(yīng)變?cè)隽?，無量綱；Δμp為等效體積應(yīng)變?cè)隽?，無量綱；D1（p*+T*）D2為常壓p下材料斷裂時(shí)的塑性應(yīng)變，無量綱；T*為最大拉伸靜水壓力，MPa；D1、D2為損傷常數(shù)。

2 連續(xù)水射流破煤特性數(shù)值模擬結(jié)果

2.1 射流能量耗散分析

通過分析淹沒/非淹沒狀態(tài)下水射流的傳播速度可以初步獲得兩者的傳播特征和能量耗散情況。設(shè)定射流直徑為0.3 mm，射流速度為170 m/s，噴嘴與煤壁作用面的間距都為20 mm，水射流在非淹沒和淹沒條件下的傳播情況如圖3。

圖3 非淹沒/淹沒水射流傳播特征Fig.3 Propagation characteristics of unsubmerged/submerged water jet

從圖3 中可以看出：非淹沒和淹沒條件下水射流自噴嘴噴出后兩者的傳播過程差別巨大；首先，從水射流形態(tài)上看，非淹沒射流傳播過程中呈圓柱狀向前傳播，射流直徑基本上沒有變化，而淹沒射流的傳播過程中的前端由于水域的阻力作用逐漸向后分離，同時(shí)水射流內(nèi)部受擠壓而產(chǎn)生應(yīng)變，這是水射流周圍部分水速低于射流速度后導(dǎo)致的滯后效應(yīng)，射流受阻是導(dǎo)致淹沒條件下水射流的傳播距離較小的根本原因；其次，射流在傳播速度上的明顯差異，非淹沒射流的傳播距離和傳播時(shí)間成正比，而相同時(shí)間淹沒射流的傳播距離大幅下降；根據(jù)射流到達(dá)煤壁的時(shí)間，淹沒條件下水射流傳播10 mm 的距離所用時(shí)間是非淹沒條件下的1.6 倍，傳播20 mm 的距離所用時(shí)間是非淹沒條件下的2.18 倍，這說明隨著傳播距離的增大淹沒條件下的水射流要不斷消耗自身的質(zhì)量來抵消環(huán)境的阻力，也就是說淹沒條件下需要水射流持續(xù)不斷地噴出才能保證射流向前傳播，否則射流就會(huì)在環(huán)境水域的阻力消耗下消失。

另外，根據(jù)水射流傳播的時(shí)間變化可以推測(cè)，隨著射流噴嘴與煤體間距的增大所需射流的持續(xù)時(shí)間也會(huì)變長(zhǎng)。

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

2.2.1 水射流破煤過程分析

為了進(jìn)一步研究射流環(huán)境對(duì)破煤的影響，建立射流直徑為3 mm、靶距為20 mm 的破煤工況，沖擊計(jì)算時(shí)間1 000 μs，射流速度分別為170、190、210 m/s 時(shí)非淹沒/淹沒射流破煤過程如圖4。

圖4 非淹沒/淹沒射流破煤過程Fig.4 Coal breaking process of unsubmerged/submerged water jet

從圖4 中可以看出：非淹沒/淹沒條件下水射流破煤特性存在著明顯的區(qū)別；首先，3 組非淹沒射流從噴嘴噴出后到達(dá)煤壁的時(shí)間不同；另外，非淹沒射流模擬試驗(yàn)的煤體全部被穿透，水射流速度為170、190、210 m/s 時(shí)煤體穿透的時(shí)間分別是970、900、830 μs，而淹沒條件下3 組試驗(yàn)全部沒有穿透煤體，但隨著射流速度的增大破煤深度逐漸增加；最后，除破煤深部不同外，2 種射流的破煤基坑直徑的差別明顯，在淹沒射流模型中，隨著射流時(shí)間的延長(zhǎng)，侵蝕坑呈漏斗狀逐漸向下延伸，但是非淹沒射流的破煤侵蝕坑始終保持圓柱狀增長(zhǎng)，這可能是由于在水域環(huán)境下射流周圍產(chǎn)生強(qiáng)烈的湍流作用使得射流的動(dòng)量不斷向周圍水分子傳遞，進(jìn)而二次沖擊使得破煤范圍增大，同時(shí)，湍流作用下的返流動(dòng)量會(huì)很大程度上減弱射流自身的沖擊能力，使得射流向前沖擊煤體的能力變?nèi)酰罱K形成漏斗形的基坑。

2.2.2 煤體應(yīng)力分布特征分析

非淹沒/淹沒射流破煤時(shí)煤體等效應(yīng)力分布如圖5。

圖5 非淹沒/淹沒射流破煤時(shí)煤體等效應(yīng)力分布Fig.5 Distribution of equivalent stress of coal mass during unsubmerged/submerged water jet breaking coal

從圖5 可以看出：2 種射流模式下的最大應(yīng)力非常接近。這主要是因?yàn)閼?yīng)力的峰值主要是由于煤體自身強(qiáng)度決定的，與射流模式關(guān)系不大，但是在應(yīng)力分布上差別明顯；非淹沒條件下的應(yīng)力主要在基坑的前方，沿橫向上分布范圍較??；而淹沒條件下應(yīng)力不僅在沖擊位置分布范圍較大，在煤體的底部也出現(xiàn)了應(yīng)力的升高，這主要是淹沒射流擠壓煤體的范圍較大，而壓強(qiáng)達(dá)到煤體應(yīng)力極限后才會(huì)出現(xiàn)損傷，從而導(dǎo)致煤體整體受力較大，而底部由于是固定邊界，因此會(huì)出現(xiàn)較高的應(yīng)力集中現(xiàn)象。

2.2.3 破煤效果分析

為了獲得不同沖擊模式下水射流破煤效率，對(duì)比分析淹沒/非淹沒射流破碎煤體基坑深度隨時(shí)間的變化關(guān)系，統(tǒng)計(jì)了2 種射流破煤模式下每間隔20 μs 的侵蝕深度變化，非淹沒/淹沒射流破煤深度變化如圖6。

圖6 非淹沒/淹沒射流破煤深度變化Fig.6 Coal breaking depth variation of unsubmerged/submerged water jet

從圖6 可以看出：2 種沖擊模式下侵蝕深度都隨著沖擊時(shí)間的增加不斷增大，區(qū)別在于非淹沒射流的侵蝕深度基本與沖擊時(shí)間成正比，并且在沖擊計(jì)算結(jié)束之前穿透煤體。為了便于對(duì)比，按照規(guī)律將非淹沒射流的侵蝕深度延深至1 000 μs，此時(shí)射流速度為170、190、210 m/s 的侵蝕深度分別達(dá)到103、112、121 mm，對(duì)應(yīng)沖擊時(shí)間時(shí)下淹沒射流的侵蝕深度分別為21、25、30 mm，分別為非淹沒射流侵蝕深度的20.4%、22.35%、24.8%。另外，淹沒條件下的侵蝕深度呈現(xiàn)階梯增加，并且增加的速度遠(yuǎn)小于非淹沒條件下增加的速度。

根據(jù)淹沒條件下侵蝕深度的變化趨勢(shì)可以看出，隨著沖擊時(shí)間的延長(zhǎng)，侵蝕深度的增加速率逐漸減小，這主要是隨著射流噴嘴距離煤體逐漸增大，射流到達(dá)煤體需要克服的行程阻力逐漸增加?？梢酝茰y(cè)隨著時(shí)間的增加射流沖擊煤體最終會(huì)停止，即一定速度下的射流通過能量不斷衰減達(dá)到煤體表面后的動(dòng)能不足以破壞煤體，所形成的侵蝕坑深度就是極限破煤深度。

由于非淹沒條件下3 種射流速度的煤體均被穿透，雖然在圖6 中進(jìn)行了預(yù)測(cè)補(bǔ)充，但卻不能對(duì)2 種射流模式下的破煤效果進(jìn)行準(zhǔn)確的分析，因此以210 m/s 穿透時(shí)的沖擊時(shí)間為上限，即對(duì)比沖擊發(fā)生后至830 μs 之間的侵蝕深度和煤體損傷單元，非淹沒/淹沒射流沖擊830 μs 時(shí)的破煤效果對(duì)比如圖7。

圖7 非淹沒/淹沒射流沖擊830 μs 時(shí)的破煤效果對(duì)比Fig.7 Comparison of coal breaking effect of unsubmerged/submerged water jet impacting 830 μs

從圖7 可以看出：在沖擊了830 μs 時(shí)，非淹沒條件下射流速度為170、190、210 m/s 時(shí)的侵蝕深度分別為83、91、100 mm，單元損傷數(shù)量分別為6 640、7 360、7 996 個(gè)；而對(duì)應(yīng)速度下淹沒射流的侵蝕深度分別為16、20、25 mm，單元損傷數(shù)量分別為7 366、8 246、9 730 個(gè)，侵蝕深度分別為非淹沒射流的19.2%、21.9%、25%，單元損傷數(shù)量分別為非淹沒射流的110.9%、112%、121%。

對(duì)比2 種沖擊模式下煤體的侵蝕深度和損傷單元的數(shù)量發(fā)現(xiàn)，不同射流速度下的破煤深度與時(shí)間基本成正線性關(guān)系，根據(jù)2.1 節(jié)的分析可知，無論是淹沒環(huán)境還是非淹沒環(huán)境，水射流的傳播速度是確定的，但淹沒射流需要更多的能量克服介質(zhì)水的阻力而使其傳播速度降低。另外，2 種沖擊模式下的損傷單元雖然隨射流速度的增加而增大，但是規(guī)律性較差，特別是淹沒環(huán)境下，損傷單元由于射流的湍流作用使得破煤面積大大增加，同時(shí)也導(dǎo)致了破煤深度大幅下降，但整體損傷單元的數(shù)量在增大。

2.3 淹沒/非淹沒射流破煤深度關(guān)系時(shí)變特征

為了研究淹沒射流的破煤特性，需要先獲得淹沒條件對(duì)水射流破煤的弱化作用。以非淹沒射流的侵蝕深度為基礎(chǔ)，求出對(duì)應(yīng)沖擊時(shí)間下淹沒射流侵蝕深度和非淹沒射流侵蝕深度的關(guān)系，淹沒/非淹沒射流破煤深度時(shí)變關(guān)系如圖8。

從圖8 可以看出：在沖擊的早期淹沒射流的侵蝕坑深度與非淹沒射流的侵蝕坑深度的比值迅速增大，這是由于在淹沒射流發(fā)生沖擊破煤的時(shí)候非淹沒射流沖擊已經(jīng)持續(xù)了一段時(shí)間，因此兩者的比值是從0 開始增加；之后，隨著非淹沒射流的侵蝕深度不斷增大，兩者比值也迅速增加，但3 種射流速度下的侵蝕深度變化規(guī)律一致，都是先迅速增大，之后趨于平穩(wěn)。這里需要注意的是趨于平穩(wěn)的是淹沒射流的侵蝕深度與非淹沒射流侵蝕深度的比值，即隨著沖擊時(shí)間的延長(zhǎng)侵蝕深度都在增大，但兩者的比值卻趨于平穩(wěn)。根據(jù)圖中擬合曲線可以推測(cè)淹沒射流的破煤深度最大約等于非淹沒射流的21%，因此根據(jù)這一結(jié)果可以結(jié)合已知非淹沒射流的破煤深度求得相應(yīng)時(shí)間下的淹沒射流破煤深度。

根據(jù)研究結(jié)果可知，淹沒環(huán)境下水射流破煤效果較差，但在工程應(yīng)用中通過延長(zhǎng)沖刷時(shí)間或射流速度可以相應(yīng)的提高破煤深度。另外，由于淹沒環(huán)境下水射流破煤效果難以考察，可以基于研究結(jié)果對(duì)非淹沒環(huán)境下的效果進(jìn)行分析，初步評(píng)價(jià)淹沒射流的可行性。

3 結(jié)語

1）淹沒條件下水射流傳播過程中需要不斷消耗自身的能量來克服環(huán)境的阻力，因此需要水射流持續(xù)不斷地噴出才能保證射流向前傳播，否則射流會(huì)在環(huán)境水域的阻力下消耗消失。

2）相同條件下淹沒射流對(duì)煤體的損傷單元數(shù)略大于非淹沒射流，但侵蝕深度遠(yuǎn)小于非淹沒射流，說明淹沒射流的能量在環(huán)境介質(zhì)的阻礙下發(fā)散并使得破煤面積增大，同時(shí)導(dǎo)致破煤深度大幅下降。

3）相同沖擊時(shí)間下淹沒射流的破煤深度約等于非淹沒環(huán)境的21%，因此，在淹沒環(huán)境下開展破煤造穴時(shí)需要合理延長(zhǎng)沖擊時(shí)間或者縮小鉆孔間距。