基于PFC2D-COMSOL 的煤層水力壓裂增透促抽瓦斯數(shù)值模擬研究

王世斌，王 剛，，陳雪暢，范酒源，遲利輝

（1.山東科技大學(xué)能源與礦業(yè)工程學(xué)院，山東青島 266590；2.山東科技大學(xué)安全與環(huán)境工程學(xué)院，山東青島 266590）

水力壓裂作為一種有效的強(qiáng)滲-增潤手段被廣泛應(yīng)用于煤層氣開采過程中，通過向煤層中打入壓裂鉆孔，注入高壓水對煤體進(jìn)行壓裂，提高煤層的滲透率，從而對煤層中的瓦斯進(jìn)行抽采[1-6]。水力壓裂及瓦斯抽采的效果往往通過煤層的裂隙擴(kuò)展、瓦斯抽采流量和抽采作用范圍等指標(biāo)進(jìn)行評價(jià)，因此，為了達(dá)到理想的壓裂和抽采效果，在實(shí)際應(yīng)用前通常會采用數(shù)值模擬軟件進(jìn)行模擬，從而確定實(shí)際應(yīng)用的工藝參數(shù)。

針對煤層水力壓裂和瓦斯抽采的數(shù)值模擬研究，國內(nèi)外研究學(xué)者開展了較多的研究工作。從?？黐7]利用FLAC3D數(shù)值模擬分析了頂板水力壓裂裂縫擴(kuò)展演化特征，討論了壓裂孔直徑和壓裂壓力對裂縫發(fā)育的影響規(guī)律；閆炎等[8]基于損傷力學(xué)與流固耦合理論，建立了水力壓裂過程中射孔段水泥環(huán)裂縫擴(kuò)展的數(shù)值模型，對水平井段壓裂工況下固井水泥環(huán)徑向裂縫沿井筒軸向的擴(kuò)展過程進(jìn)行了數(shù)值模擬；付海峰等[9]通過三維離散格子模擬方法建立了等尺度射孔水力裂縫起裂擴(kuò)展全三維數(shù)值模型，揭示了不同射孔模式下的裂縫起裂規(guī)律；張玉等[10]以多孔介質(zhì)流體滲流和圍巖應(yīng)力耦合理論為基礎(chǔ)，提出一種基于有限容積法（FVM）的水力耦合作用下射孔圍巖水力壓裂破裂數(shù)值模擬方法，該方法豐富了水力壓裂破裂機(jī)理的研究，對實(shí)際工程設(shè)計(jì)提供重要的參考；劉嘉等[11]基于多孔彈性理論和能量最小化原理建立水力耦合的相場模型，采用交錯(cuò)策略的分離式耦合方法進(jìn)行求解，并基于三維水力壓裂模擬結(jié)果進(jìn)行分析，驗(yàn)證了該方法對于模擬不同地應(yīng)力條件下的水力裂紋擴(kuò)展問題的適用性；閆曉等[12]提出了統(tǒng)一管道-界面單元法（UP-IEM），模擬可滲透性孔隙裂隙介質(zhì)中的水力壓裂，通過3 個(gè)經(jīng)典數(shù)值模型驗(yàn)證了UP-IEM 模擬流體流動(dòng)、裂紋擴(kuò)展以及水力壓裂的準(zhǔn)確性及適用性；楊晨旭等[13]針對水力壓裂過程中所涉及到的巖石變形與裂縫擴(kuò)展、縫內(nèi)流體的流動(dòng)及濾失等力學(xué)問題，基于有限元法建立了全三維全耦合水力壓裂數(shù)值模型，并采用數(shù)值模擬方法對模型進(jìn)行驗(yàn)證；郭欣等[14]以瓦斯?jié)B流理論為基礎(chǔ)，根據(jù)達(dá)西定律和質(zhì)量守恒定律，利用COMSOL Multiphysics 數(shù)值模擬軟件建立考慮煤巖變形控制方程的氣固耦合數(shù)學(xué)模型，研究相關(guān)因素對瓦斯有效抽采半徑大小的影響；郝長勝等[15]運(yùn)用FLUENT 軟件模擬分析布置于平行回風(fēng)巷且孔間距為10 m 的頂板走向長鉆孔束在不同參數(shù)條件下的瓦斯抽采效果，為寺河礦5309 工作面上隅角瓦斯超限問題的解決提供了一種方法；賈毅超等[16]采用FLAC 數(shù)值模擬軟件對高位鉆孔布置優(yōu)化進(jìn)行了研究，并在實(shí)際工程應(yīng)用中提高了瓦斯抽采濃度，減小了上隅角和回風(fēng)流中的瓦斯?jié)舛?；汪騰蛟等[17]應(yīng)用有限體積法分別離散了考慮溫度變化下的采空區(qū)瓦斯抽采多場耦合二維數(shù)學(xué)模型，基于Microsoft Visual Basic 編制了計(jì)算機(jī)解算程序，利用Tecplot軟件對求解結(jié)果進(jìn)行可視化，研究了瓦斯抽采前后的采空區(qū)瓦斯壓力、瓦斯體積分?jǐn)?shù)、氧體積分?jǐn)?shù)、固體溫度和氣體溫度分布情況及前后變化趨勢；王剛等[18]建立急傾斜煤層鉆孔周圍氣體滲流單元模型，利用數(shù)值模擬分析了急傾斜煤層瓦斯鉆孔有效抽采范圍隨時(shí)間的變化規(guī)律；王廣宏[19]利用COMSOL Multiphysics 數(shù)值模擬軟件對榆樹田煤礦的定向長鉆孔瓦斯抽采工藝技術(shù)進(jìn)行了研究，以提高瓦斯抽采效率，確保礦井安全生產(chǎn)；張育磊[20]以萬峰煤礦1201 工作面為工程背景，采用FLAC3D數(shù)值模擬軟件，研究煤層開采后對鄰近煤層的影響及合理的瓦斯抽采鉆孔布置位置；張智淵等[21]基于建立的復(fù)合儲層和單一儲層瓦斯抽采流固耦合模型，模擬研究了不同傾角的穿層鉆孔對復(fù)合儲層和單一儲層的瓦斯抽采效果，以及多鉆孔耦合下的穿層鉆孔瓦斯抽采效果。

從上述研究中可以看出，絕大多數(shù)學(xué)者僅針對煤層水力壓裂和瓦斯抽采中1 種工藝方式進(jìn)行了數(shù)值模擬，未對水力壓裂作用后的煤體開展瓦斯抽采數(shù)值模擬，也即未對2 種工藝數(shù)值模擬間的“鏈接”進(jìn)行構(gòu)建，從而無法達(dá)到先壓裂后抽采的模擬要求。針對這一問題，提出了一種基于PFC2D水力壓裂和COMSOL 瓦斯抽采結(jié)合的數(shù)值模擬方法，將PFC2D水力壓裂后裂隙擴(kuò)展的圖像導(dǎo)入到網(wǎng)格繪制的圖像處理軟件中，將圖像進(jìn)行像素級優(yōu)化，按照像素點(diǎn)將圖像中的裂隙進(jìn)行繪制，然后將繪制的圖像導(dǎo)入COMSOL 中進(jìn)行瓦斯抽采數(shù)值模擬，從而實(shí)現(xiàn)了煤層先壓裂后抽采的數(shù)值模擬研究，研究內(nèi)容對煤層氣開采具有重要的指導(dǎo)意義。

1 煤層水力壓裂PFC2D 數(shù)值模擬研究

1.1 PFC2D 數(shù)值模擬模型及參數(shù)

模型構(gòu)建采用PFC2D軟件建立200 mm × 200 mm 的二維模型，模型的中心位置設(shè)定半徑7.0 mm的水力壓裂孔，模型4 個(gè)邊界為固定邊界，即對模型邊界的顆粒進(jìn)行固定。水力壓裂模型示意圖如圖1。

PFC2D模擬時(shí)的宏觀和細(xì)觀參數(shù)如下[22-23]：

1）宏觀力學(xué)參數(shù)。單軸抗壓強(qiáng)度σmc為32.22 MPa，彈性模量E 為1.81 GPa，泊松比ν 為0.25，本體抗拉強(qiáng)度σmt為2.83 MPa。

2）顆粒力學(xué)參數(shù)。顆粒最小半徑Rmin為2.50 mm，粒徑比Rmax/Rmin為1.66，顆粒密度ρ 為1 280 kg/m3，孔隙率為9%，阻尼β 為0.50 N/（s·m）。

3）接觸模型。細(xì)觀彈性模量Ec為1.00 GPa，剛度比kn/ks為2.30，摩擦因數(shù)μ 為0.50，細(xì)觀抗拉強(qiáng)度σc為11 MPa，黏聚力c 為20 MPa，摩擦角φ 為37.50°，法向臨界阻尼比βn0.57。水力壓裂模型示意圖如圖1。

圖1 水力壓裂模型示意圖Fig.1 Schematic diagram of hydraulic fracturing model

除水力壓裂的注水壓力外，其余參數(shù)保持不變，模擬20 MPa 和40 MPa 的不同注水壓力下煤體水力壓裂，模擬時(shí)僅設(shè)置注水鉆孔邊界處的壓力為20 MPa 和40 MPa，且在注水過程中注水壓力不發(fā)生變化，不存在鉆孔注水壓力上升的過程，始終保持20 MPa 和40 MPa，且裂紋擴(kuò)展后鉆孔注水壓力也為20 MPa 和40 MPa，直到模擬時(shí)注水壓力不發(fā)生改變。計(jì)算步數(shù)為10 000，每一步的時(shí)間間隔為0.001 s，對應(yīng)總時(shí)間為10 s。

1.2 水力壓裂實(shí)現(xiàn)方式

PFC2D進(jìn)行水力壓裂時(shí)，為了定量計(jì)算流體交換時(shí)的流量，假設(shè)流體通道是由相鄰兩顆粒的接觸點(diǎn)組成的1 個(gè)平行板通道，長度為L，開度為a，厚度為單位厚度，垂直于x-y 平面[24]。則2 個(gè)孔隙之間的流量可以通過Hagen-Poiseuille 方程表示：

式中：q 為流量，m3/s；k 為滲透系數(shù)，m/（Pa·s）；△p 為2 個(gè)孔隙網(wǎng)格之間的壓力差，Pa。

在顆粒之間接觸處，假設(shè)存在1 個(gè)初始開度a0，初始開度的存在允許即使兩顆粒之間緊密接觸的情況下其構(gòu)成的流體通道也有流體流動(dòng)，從而保證了材料的基質(zhì)滲透率。流體通道的開度a 依賴于顆粒之間的接觸力，當(dāng)兩顆粒之間的法向接觸力為壓應(yīng)力時(shí)，流體通道開度a 為：

式中：F 為現(xiàn)在兩顆粒間的壓縮力，N；F0為管道開度降低為初始開度一半時(shí)的壓縮力，N。

當(dāng)兩膠結(jié)顆粒處于張拉狀態(tài)，或者兩顆粒之間的黏結(jié)已破壞，顆粒在接觸點(diǎn)斷開時(shí)，開度a 為：

式中：d 為兩顆粒之間的距離，m；R1、R2分別為兩顆粒半徑，m；λ 為1 個(gè)無量綱乘子，對于大部分的模型而言，顆粒粒徑都比實(shí)際的顆粒粒徑大得多，計(jì)算得到的開度也會偏大，因此，λ 常常取1 個(gè)小于1 的常數(shù)來獲得1 個(gè)合理的開度。

在△t 的時(shí)間步里面，由于流體流動(dòng)導(dǎo)致的孔隙流體壓力變化由流體的體積壓縮模量計(jì)算?？紤]某個(gè)孔隙，其有N 條流體通道，在△t 的時(shí)間步里面，其流體總流量為∑q，則孔隙流體壓力的變化△p 可表示為：

式中：Kf為流體的壓縮模量，Pa；Vd為“域”，即孔隙體積，m3；△Vd為孔隙體積變化，m3。

在PFC2D中流體與顆粒間的耦合作用方式主要有以下3 種：①通過接觸的張開與閉合或接觸力的變化實(shí)現(xiàn)通道孔隙的變化；②通過改變研究區(qū)域的力學(xué)特性來改變其中的壓力；③區(qū)域孔隙壓力對其內(nèi)部顆粒有推移作用。

采用顯式方法求解，在所有域內(nèi)交替應(yīng)用流動(dòng)方程和壓力方程。保證模型運(yùn)行穩(wěn)定的條件就是水流引起的壓力變化必須小于擾動(dòng)壓力，當(dāng)兩者相等時(shí)可求出臨界時(shí)間步長△t 為：

式中：N 為1 個(gè)“域”所連接的管道數(shù)；r 為1 個(gè)“域”周圍顆粒的平均半徑，m。

1.3 水力壓裂結(jié)果

不同注水壓力的壓裂范圍可通過壓裂后顆粒的位移云圖來反映，不同注水壓力的裂隙擴(kuò)展規(guī)律可通過裂隙擴(kuò)展圖表現(xiàn)。2 種注水壓力下顆粒位移云圖如圖2，2 種注水壓裂狀態(tài)下的裂隙擴(kuò)展如圖3。

圖3 2 種注水壓力下注水壓裂裂隙擴(kuò)展圖Fig.3 Fracture propagation diagrams under two kinds of water injection pressure

從圖2 中可以看出：2 種注水壓力下的壓裂范圍有所不同；20 MPa 的注水壓力時(shí)，僅有鉆孔周圍的小部分顆粒具有位移，說明20 MPa 注水壓力下，注水壓裂的范圍較??；40 MPa 的注水壓力時(shí)，鉆孔周圍的許多顆粒均產(chǎn)生了位移，且顆粒的位移量沿鉆孔徑向逐漸減小，表明40 MPa 注水壓力下，注水壓裂的范圍較大，同時(shí)說明增加注水鉆孔的注水壓力可以提高水力壓裂的壓裂效果。

圖2 2 種注水壓力下顆粒位移云圖Fig.2 Cloud diagrams of particle displacement under two kinds of water injection pressure

從圖3 中可以看出：裂隙從鉆孔附近沿徑向逐漸擴(kuò)展，不同注水壓力下生成的裂隙數(shù)量和范圍有所不同；20 MPa 注水壓力生成的裂隙數(shù)量較少，僅有8 條裂隙；而40 MPa 注水壓力生成的裂隙數(shù)量較多，有78 條裂隙；20 MPa 注水壓力下裂隙擴(kuò)展延伸的范圍較小，壓裂效果較差；40 MPa 注水壓力下裂隙擴(kuò)展延伸的范圍較大，壓裂效果較好；同時(shí)也表明增加鉆孔的注水壓力可以增加水力壓裂的裂隙擴(kuò)展延伸，提高水力壓裂的效果。

2 煤層瓦斯抽采COMSOL 數(shù)值模擬

2.1 COMSOL 模型

將PFC2D水力壓裂后裂隙擴(kuò)展的圖像導(dǎo)入到網(wǎng)格繪制的圖像處理軟件中，將圖像進(jìn)行像素級優(yōu)化，按照像素點(diǎn)對圖像中的裂隙進(jìn)行繪制，將像素長度與實(shí)際長度進(jìn)行匹配，然后將繪制的圖像導(dǎo)出為COMSOL 能夠識別的網(wǎng)格文件，最后將網(wǎng)格文件導(dǎo)入到COMSOL 數(shù)值模擬軟件中，瓦斯抽采幾何模型構(gòu)建圖如圖4。

圖4 瓦斯抽采幾何模型構(gòu)建圖Fig.4 Geometric model construction diagrams of gas extraction

2.2 控制方程

模擬采用COMSOL 軟件中的達(dá)西定律模塊，控制方程為：

式中：εp為孔隙率；ρ 為瓦斯密度，kg/m3；u 為瓦斯的達(dá)西滲流速度，m/s；Qm為源匯項(xiàng)，kg/（m3·s）；k為滲透率，m2；μ 為瓦斯動(dòng)力黏度，Pa·s；△pg為壓力梯度，Pa/m；t 為時(shí)間變量，s。

將處于滲流運(yùn)動(dòng)中的瓦斯視為理想氣體，滿足理想氣體狀態(tài)方程。該方程表示平衡狀態(tài)下，理想氣體的壓強(qiáng)、體積、物質(zhì)的量以及溫度之間的關(guān)系為：

式中：Mg為瓦斯的摩爾質(zhì)量，kg/mol；pg為瓦斯壓力，Pa；R 為氣體摩爾常數(shù)，J/（mol·K）；T 為絕對溫度，K。

2.3 邊界條件及參數(shù)

瓦斯抽采數(shù)值模擬時(shí)的邊界條件如圖5。

圖5 瓦斯抽采邊界條件示意圖Fig.5 Schematic diagram of boundary conditions of gas extraction

模型四周邊界無瓦斯流動(dòng)，煤層的初始瓦斯壓力為p0，鉆孔的瓦斯抽采壓力為pout，由于PFC2D構(gòu)建注水鉆孔時(shí)采用的是刪除顆粒的命令，當(dāng)顆粒有不存在刪除區(qū)域中的部分時(shí)，顆粒就不會被刪除，因此會使構(gòu)建的鉆孔直徑有所變化，因此在抽采時(shí)將鉆孔直徑設(shè)置為5 mm。各參數(shù)具體數(shù)值如下：①煤體孔隙率φ0為9.00%；②煤體孔隙滲透率k0為9×10-16m2；③煤體裂隙滲透率kl為9×10-13m2；④初始瓦斯壓力p0為1 MPa；⑤抽采負(fù)壓pout為0.2 MPa；⑥瓦斯動(dòng)力黏度μ 為1×10-5Pa·s；⑦瓦斯密度ρ 為0.716 kg/m3；⑧氣體摩爾常數(shù)R 為8.314 5 J/（mol·K）；⑨參考溫度T 為293.15 K；⑩瓦斯摩爾質(zhì)量Mg為16 g/mol。

2.4 瓦斯抽采結(jié)果

2 種壓裂條件下瓦斯抽采10 s 的壓力分布云圖如圖6。20 MPa 注水壓裂煤層瓦斯抽采壓力分布云圖如圖7，40 MPa 注水壓裂煤層瓦斯抽采壓力分布云圖如圖8。

圖6 抽采10 s 瓦斯壓力分布云圖Fig.6 Gas pressure distribution diagrams of 10 s

圖7 20 MPa 注水壓裂煤層瓦斯抽采壓力分布云圖Fig.7 Diagrams of gas drainage pressure distribution in 20 MPa water injection fracturing coal seam

從圖6 中可以看出：瓦斯壓力沿著鉆孔徑向向中心逐漸減小，這是因?yàn)槊簩又械耐咚贡怀椴?，靠近鉆孔周圍的煤層最先被抽采出，瓦斯壓力減小，煤層中瓦斯壓力分布不均勻，從而使得高壓瓦斯不斷向低壓處運(yùn)移。

瓦斯抽采范圍可以通過瓦斯壓力分布云圖反映出，瓦斯壓力產(chǎn)生變化的地方表明鉆孔的瓦斯抽采范圍能夠達(dá)到此處。從圖7 和圖8 可以看出：隨著抽采時(shí)間的增加，瓦斯抽采的范圍逐漸增大；當(dāng)抽采時(shí)間相同時(shí)，40 MPa 與20 MPa 水力壓裂的煤層相比，前者的瓦斯抽采范圍比后者大，表明增加水力壓裂的注水壓力不僅增加了煤層中裂隙的擴(kuò)展范圍，同時(shí)也提高了煤層的瓦斯抽采范圍。

采用COMSOL 中的線積分對鉆孔的瓦斯抽采流量進(jìn)行計(jì)算，鉆孔瓦斯抽采流量變化如圖9。由于本文COMSOL 瓦斯抽采數(shù)值模擬采用的是二維幾何模型，對于抽采流量的計(jì)算是采用線積分的方式，因此抽采流量單位為m2/s。

圖9 鉆孔瓦斯抽采流量變化Fig.9 Change of gas extraction flow rate from borehole

從圖9 中可以看出，隨著抽采時(shí)間的增加，抽采流量逐漸降低，在抽采開始后的10 s 內(nèi)，鉆孔瓦斯抽采流量下降的最快，20 MPa 水力壓裂下瓦斯抽采流量從0.006 6 m2/s 降低為0.001 6 m2/s；40 MPa 水力壓裂下瓦斯抽采流量從0.010 1 m2/s 降低為0.002 0 m2/s，這是因?yàn)榇藭r(shí)煤層的瓦斯壓力與鉆孔抽采壓力之間形成的壓力梯度最大，瓦斯運(yùn)移的速度快，從而使得鉆孔的瓦斯抽采流量大。隨著抽采的進(jìn)行，煤層中的瓦斯含量和瓦斯壓力均有所下降，瓦斯壓力與抽采鉆孔之間的壓差逐漸減小，因此會造成瓦斯流量的降低。同時(shí)，2 種壓裂條件下的瓦斯抽采流量也不相同，40 MPa 水力壓裂后的煤層瓦斯抽采流量比20 MPa 水力壓裂后的煤層大，表明增加水力壓裂的注水壓力在增加煤層裂隙擴(kuò)展的同時(shí)也增加了瓦斯抽采的流量。

文獻(xiàn)［25］的實(shí)際工程應(yīng)用中也顯示水力壓裂后鉆孔的瓦斯抽采流量有所增加，表明水力壓裂可以提高煤層的透氣性，提高煤層的瓦斯抽采效果，進(jìn)一步驗(yàn)證所提出的數(shù)值模擬方法具有可行性。

3 結(jié) 語

對煤層開展了PFC2D水力壓裂數(shù)值模擬，分析了不同注水壓力下的壓裂范圍與裂隙擴(kuò)展規(guī)律，并采用一種圖像處理方法將壓裂后的模型導(dǎo)入到COMSOL 軟件中進(jìn)行壓裂后的瓦斯抽采數(shù)值模擬研究，分析了不同壓裂條件下的瓦斯抽采范圍與瓦斯抽采流量。

1）提出的幾何模型圖像像素級處理方法成功將PFC2D水力壓裂后的裂隙幾何模型導(dǎo)入到COMSOL軟件中，完成了2 個(gè)數(shù)值模擬軟件之間的結(jié)合，該方法成功實(shí)現(xiàn)了對水力壓裂煤層進(jìn)行瓦斯抽采的數(shù)值模擬。

2）40 MPa 與20 MPa 注水壓力相比，注水壓裂的范圍較大，產(chǎn)生的裂隙數(shù)量較多，裂隙擴(kuò)展范圍較大，壓裂效果好。隨著抽采時(shí)間的增加，瓦斯抽采的范圍逐漸增大，抽采流量逐漸降低。且相同抽采時(shí)間下，40 MPa 注水壓裂的煤層瓦斯抽采范圍和抽采流量大，水力壓裂可以提高煤層的瓦斯抽采效果，與實(shí)際工程應(yīng)用效果相符，模擬效果較好，該模擬方法具有可行性，能夠?yàn)楝F(xiàn)場應(yīng)用提供指導(dǎo)。