陳 鵬， 孫可明*,2， 張 宇

(1.遼寧工程技術大學 力學與工程學院，阜新 123000； 2.青島理工大學 理學院，青島 266000)

1 引 言

我國煤層氣資源豐富，具有巨大的開發(fā)潛力，但抽采量和利用量較低，增速緩慢，主要原因是煤儲層具有低壓力、低滲透率以及低含氣飽和度的三低儲層特征[1]，滲透率一般在(0.1～0.001)×10-3μm2范圍內(nèi)[2]，導致煤層氣抽采難度大，時間長，成本高，因此煤儲層滲透率低成為商業(yè)化開采的瓶頸。超臨界CO2氣爆技術是一種新型的致裂增透技術，相比深孔化學炸藥爆破的增透方法，超臨界CO2物理氣爆增透方法具有安全、高效、環(huán)保、適用范圍廣和可控性好等優(yōu)點，特別對于高瓦斯煤層氣開發(fā)其優(yōu)勢更為明顯，超臨界CO2近似液體的密度使氣爆具有較強的沖擊效果，使煤體內(nèi)生成大量氣爆裂縫，進而提高了增透效果。

天然煤體經(jīng)過長期地質(zhì)作用，孕育了大量的割理裂隙和層理等結構弱面，這些結構弱面的存在使煤體表現(xiàn)較強的非均質(zhì)性，并影響著煤體的力學特性。宋浩然等[3]基于煤體的各向異性和非均質(zhì)性，考慮煤體應力變形場和瓦斯?jié)B流場的交叉耦合作用，分析了煤層抽采中水力割縫鉆孔周圍瓦斯壓力以及滲透率的時空演化規(guī)律；宋紅華[4]等基于波速測試、CT掃描和三維重構技術，分析了煤巖內(nèi)部原生裂隙、孔隙以及礦物夾雜分布的非均質(zhì)特征，研究了單軸受壓條件下，沿不同方向加載時，煤樣內(nèi)部結構分布的非均質(zhì)性對煤巖破壞特征的影響；賈慧敏[5]研究發(fā)現(xiàn)，煤巖孔隙結構具有分段分形特征，煤巖分形維數(shù)越大，其孔隙度和滲透率越小，孔隙分選性越差，非均質(zhì)性越強；趙瑜等[6]采用非均質(zhì)固氣耦合數(shù)值試驗的方法，研究不同瓦斯壓力作用下，非均質(zhì)煤巖抗壓強度尺寸效應的影響規(guī)律及其作用機理。孫可明等[7-10]進行了大量的超臨界CO2氣爆煤體致裂實驗研究及均質(zhì)煤體的氣爆模擬研究，且目前關于模擬非均質(zhì)煤體爆破的相關研究鮮有報道，因此有必要進一步開展超臨界CO2氣爆致裂非均質(zhì)煤體的模擬研究，對提高超臨界CO2氣爆技術的應用水平具有重要作用。

2 超臨界CO2氣爆煤體力學模型

2.1 爆轟波守恒方程

超臨界CO2氣爆瞬間產(chǎn)生的爆轟波滿足質(zhì)量、動量和能量守恒方程[11]

(1,2)

(3)

2.2 煤體變形場方程

(1) 增量型平衡微分方程

(4)

(2) 增量型幾何方程

(5)

式中Δεi j為應變的增量，Δui，j為位移對坐標偏導數(shù)的增量。

(3) 煤體損傷本構模型

由于煤體氣爆是大應變和高應變率的動態(tài)過程，煤體材料采用Johnson-Cook模型，其屈服應力可表示為[12]

(6)

3 超臨界CO2氣爆非均質(zhì)煤體有限元模型

為了建立非均質(zhì)含弱面煤體的幾何模型，基于數(shù)字圖像方法開發(fā)了對實際煤體圖像識別處理的Matlab程序，將實際煤體圖像灰度化，利用灰度值將400×400像素的圖像分隔為400×400個節(jié)點矩陣，利用此節(jié)點矩陣和像素代表的真實長度構建出非均質(zhì)煤體的幾何表征模型，將其導入Abaqus中生成非均質(zhì)煤體有限元模型，實際煤體圖像與有限元模型對比如圖1所示。

圖1 實際煤體圖像與非均質(zhì)煤體有限元模型

為減少計算工作量和便于對比分析，建立的煤體模型尺寸為400 mm×400 mm×20 mm，煤體中心處氣爆孔直徑為16 mm(與實驗模型[9]尺寸一致)，與爆孔軸線平行的四個邊界面設置固定約束，且各邊界無反射，煤體x和y方向施加的初應力如圖2所示，α為弱面傾角。煤體上的監(jiān)測點布置方式如圖3所示，x和y方向布置的監(jiān)測點到爆孔中心的距離相同，分別為10 mm，30 mm，60 mm，100 mm，140 mm及180 mm。實驗得到了煤體及弱面的物理力學參數(shù)[9]，列入表1。

圖2 氣爆煤體計算模型

圖3 監(jiān)測點布置

表1 煤體及弱面的物理力學參數(shù)

模擬中采用JWL狀態(tài)方程確定超臨界CO2氣爆過程中爆生氣體的壓力變化，其表達式為[15]

(7)

模擬過程采用了SPH與FEM聯(lián)合求解的計算方法，煤體與超臨界CO2單元類型均為C3D8R，SPH粒子直徑為1 mm，煤體模型中實體單元達到轉化閾值時轉變?yōu)镾PH粒子，轉化閾值采用應力標準，轉變?yōu)镾PH粒子的區(qū)域即為煤體破壞區(qū)；計算開始時,超臨界CO2有限單元直接轉變?yōu)镾PH粒子，粒子可以進入非有限單元區(qū)域與有限單元發(fā)生相互作用；粒子與粒子、粒子與有限單元之間接觸類型采用通用接觸來滿足質(zhì)量、動量和能量守恒。

4 含單個弱面煤體氣爆的模擬結果及分析

地下煤層往往受地應力作用，文獻[16]給出地應力公式為

σv=γH

(8)

式中σv為鉛錘應力，γ為上覆巖體的平均重力密度，H為巖體單元的深度。

本文選取地下800 m深度煤層所受地應力情況進行模擬，煤層頂板巖石的平均重力密度γ為 25 kN/m3，由式(8)可得σv為20 MPa，因此模擬中將煤體x和y方向的初應力工況設置為靜水應力狀態(tài)σx=σy=20 MPa及非靜水應力狀態(tài)σx=0 MPa，σy=20 MPa和σx=20 MPa，σy=10 MPa。

為了控制變量的統(tǒng)一，將一個含單個弱面煤體的圖像經(jīng)過圖像的拉伸和旋轉得到圖4和圖8所示的圖像，保證了這些圖像中煤體的弱面為相同弱面，弱面厚度約為10 mm。圖4弱面傾角為0°，弱面到爆孔中心距離分別為40 mm,80 mm,120 mm和160 mm；圖8弱面到爆孔中心距離為80 mm，弱面傾角分別為0°，15°，30°和45°。

圖4 實際煤體圖像

4.1 弱面到爆孔中心距離對氣爆致裂效果的影響

圖5(a)氣爆產(chǎn)生的主裂縫沿x和y方向擴展，可以看出弱面對裂縫擴展具有阻礙作用，隨著弱面到爆孔中心距離的增大，弱面對裂縫擴展的阻礙作用逐漸增大，穿過弱面擴展的裂縫密度逐漸減小；圖5(b)氣爆產(chǎn)生的主裂縫只沿y方向擴展，隨著弱面到爆孔中心距離的增大，穿過弱面擴展的裂縫尺度逐漸減小，弱面到爆孔中心距離為160 mm時，沒有裂縫穿過弱面擴展；圖5(c)氣爆產(chǎn)生的主裂縫只沿x方向擴展，主裂縫上生成的分支裂縫偏向y方向擴展。

圖5 不同初應力下弱面到爆孔中心不同距離氣爆結果

由模擬結果可知，弱面到爆孔中心距離較小時，有裂縫穿過弱面擴展，隨著距離的增大，穿過弱面擴展的裂縫尺度和密度不斷減小，當弱面距離爆孔中心足夠遠時，弱面完全阻斷了裂縫的擴展，沒有裂縫穿過弱面。

圖6和圖7給出了部分含單個弱面煤體氣爆時裂縫擴展及應力波傳播過程。在超臨界CO2氣爆初始時間段，爆孔周圍介質(zhì)在氣爆應力波作用下發(fā)生粉碎性破壞；45 μs 后爆孔周圍介質(zhì)繼續(xù)破壞，爆孔四周出現(xiàn)環(huán)向裂縫，裂縫隨著應力波傳播時間的增加不斷擴展，裂縫擴展到弱面時，弱面阻礙了裂縫的擴展，僅有少量裂縫穿過弱面繼續(xù)擴展；240 μs后主裂縫不再擴展，只有分支裂縫在擴展。初始時應力波呈現(xiàn)圓環(huán)狀向外傳播，75 μs時應力波傳播到弱面，表現(xiàn)非環(huán)狀特征，75 μs～120 μs 可以看出弱面對應力波的阻隔作用，應力波未穿透弱面，說明應力波在弱面上發(fā)生吸收和反射，135 μs后應力波傳播到弱面上側的煤體上，應力波發(fā)生衰減，165 μs后弱面下側的應力波穿過計算區(qū)域，240 μs后應力波穿過整個煤體。

圖6 裂縫擴展過程

結合應力分析可知，氣爆初始時，應力波主要表現(xiàn)為壓力波，對爆孔周圍煤體進行擠壓破壞形成粉碎區(qū)，接著爆孔周圍產(chǎn)生環(huán)向裂縫，爆生氣體進入初始裂縫驅(qū)動其繼續(xù)擴展形成裂隙區(qū)；應力波傳播到弱面時，應力值出現(xiàn)不連續(xù)，入射應力波分化為反射應力波和透射應力波，反射應力波表現(xiàn)為拉應力，使弱面內(nèi)側產(chǎn)生向爆孔方向擴展的反射裂縫，在入射壓應力波和反射拉應力波的作用下，造成煤體拉壓累積損傷，反射裂縫與爆孔處產(chǎn)生的徑向裂縫交匯貫通，使弱面和爆孔之間形成裂縫密集區(qū)；透射壓應力波繼續(xù)在煤體內(nèi)傳播，由于弱面對裂縫擴展的阻礙，弱面外側煤體上只生成少量裂縫。

圖7 應力波傳播過程

圖8 實際煤體圖像

4.2 弱面傾角對氣爆致裂效果的影響

由圖9(a，b)可以看出，隨著弱面傾角的增大，弱面對裂縫擴展的阻礙作用逐漸增大，穿過弱面擴展的裂縫尺度和密度逐漸減??；由圖9(b，c)可以看出，45°傾角弱面對裂縫擴展的阻礙作用最大，幾乎沒有裂縫穿過弱面。

圖9 不同初應力下弱面不同傾角氣爆結果

對比圖5和圖9的3種初應力工況下的模擬結果可得，氣爆主裂縫沿最大初應力方向擴展，初應力對裂縫的萌生和擴展具有導向作用。

根據(jù)圖3監(jiān)測不同弱面傾角下y方向各個監(jiān)測點處的振動速度峰值，監(jiān)測數(shù)據(jù)列入表2。

表2 質(zhì)點振動速度峰值

圖10 質(zhì)點振動速度峰值曲線

由圖10可以看出，監(jiān)測點1和監(jiān)測點2與弱面距離較大，受弱面的影響較小，不同傾角時同一監(jiān)測點的振動速度相差不大；弱面傾角為0°和15°時，監(jiān)測點3接近弱面，傾角為30°和45°時，監(jiān)測點4接近弱面，由曲線可看出，接近弱面處的監(jiān)測點振動速度發(fā)生了突變，因為應力波傳播到弱面時，在弱面上產(chǎn)生能量聚集，造成應力集中，導致弱面附近的質(zhì)點振動速度大幅度增大；弱面后的質(zhì)點振動速度峰值隨弱面傾角的增大而減小。

應力波通過弱面后會發(fā)生衰減，弱面不同傾角時，監(jiān)測點6均位于弱面之后，監(jiān)測得到該點的Mises應力峰值列入表3。數(shù)據(jù)表明，弱面后Mises應力峰值隨其傾角的增大而減小，由圖11可知應力波通過弱面時，弱面傾角越大，應力波衰減越多。

表3 監(jiān)測點6的Mises應力峰值

圖11 監(jiān)測點6的Mises應力峰值曲線

文獻[17]給出爆炸應力波透射率為透射應力與入射應力的比值，即

T=σt/σi

(9)

式中T為應力波透射率,σi為入射應力,σt為透射應力。

按照應力波向弱面入射和透射的時間及其入射和透射方向，監(jiān)測弱面上及接近弱面兩端的煤體上的應力，得到應力波由煤體傳入弱面及由弱面?zhèn)魅朊后w時的入射應力和透射應力，根據(jù)式(8)計算出弱面不同傾角時的應力波透射率(表4)，并得到 圖12 所示的曲線。

表4 應力波透射率

結果表明，應力波由煤體傳入弱面時的透射率小于1，應力波由弱面?zhèn)魅朊后w時的透射率大于1；隨著弱面傾角的增大，應力波透射率不斷減小，應力波透射率衰減度不斷增大。

5 含多弱面煤體氣爆的模擬結果及分析

圖1給出了含多個弱面煤體的實際圖像及其經(jīng)過圖像識別處理得到的有限元模型，給煤體x和y方向均施加20 MPa初應力，得到氣爆裂縫擴展過程及應力波傳播過程如圖13和圖14所示。

對比含單個弱面的煤體氣爆結果可知，多個弱面的存在改變了煤體氣爆的粉碎區(qū)和裂隙區(qū)形貌，應力波波形變得不規(guī)則；煤體中弱面越多，應力波衰減越多，弱面處生成向爆孔擴展的反射裂縫越多，爆孔與弱面之間煤體上的裂縫密度越大。

圖12 應力波透射率曲線

圖14 應力波傳播過程

圖15 x方向應力時程曲線

對比同一監(jiān)測點上的速度和應力時程曲線可知，質(zhì)點在應力波作用下產(chǎn)生振動，質(zhì)點起振時刻與應力波到達時刻一致。x方向監(jiān)測點7處介質(zhì)在氣爆沖擊載荷作用下瞬間達到抗壓強度發(fā)生壓縮破壞，承載應力得以釋放，應力變?yōu)?恒定不變；監(jiān)測點8處介質(zhì)受到較大沖擊應力的作用并發(fā)生屈服，之后回彈，應力逐漸減?。粦ΣS時間的推移向外傳播，傳到弱面時發(fā)生反射和吸收，反射應力波為拉力波，煤體發(fā)生拉壓累積破壞，應力波通過弱面后發(fā)生衰減，造成監(jiān)測點處的速度和應力大幅度減小，曲線波動較小。y方向在爆孔邊緣存在一個弱面，受弱面影響，監(jiān)測點1的應力和速度發(fā)生突變，其后監(jiān)測點的速度和應力時程曲線變化規(guī)律與x方向的規(guī)律比較接近。

圖16 y方向應力時程曲線

圖17 x方向速度時程曲線

圖18 y方向速度時程曲線

6 結 論

本文基于數(shù)字圖像方法開發(fā)了對實際煤體圖像識別處理的Matlab程序，建立了非均質(zhì)煤體有限元模型，模擬了不同初應力條件下超臨界CO2氣爆非均質(zhì)煤體的致裂過程，分析了初應力及煤體的弱面對超臨界CO2氣爆裂隙演化和應力波傳播的影響規(guī)律，得到以下主要結論。

(1) 弱面傾角相同時，弱面到爆孔中心距離距離越大，弱面對裂縫擴展的阻礙作用越大，穿過弱面的裂縫尺度和密度越小。

(2) 弱面到爆孔中心距離相同時，弱面傾角越大，弱面對裂縫擴展的阻礙作用越大(弱面傾角為45°時阻礙作用最大)，穿過弱面的裂縫尺度和密度越小，應力波透射率越小，應力波透射率衰減度越大。

(3) 煤體受初應力作用時，主裂縫會沿最大初應力方向擴展，初應力對氣爆裂縫的萌生和擴展有導向作用。

(4) 天然煤體中的弱結構面對超臨界CO2氣爆致裂效果影響較大，弱面數(shù)量越多、形態(tài)越復雜，煤體上生成的氣爆裂縫密度越大，可較好地提高煤體的致裂增透效果。