趙 偉，袁 源，王 凱，徐 超，宋亞楠，馮忠凱，郭小芳，王柳懿

（1.中國(guó)礦業(yè)大學(xué)（北京）共伴生能源精準(zhǔn)開(kāi)采北京市重點(diǎn)實(shí)驗(yàn)室，北京 100083；2.中國(guó)礦業(yè)大學(xué)（北京）應(yīng)急管理與安全工程學(xué)院，北京 100083）

煤層的透氣性是影響煤層氣大規(guī)模開(kāi)采的重要因素。然而隨著埋深加大，地應(yīng)力逐步加大，進(jìn)而使得煤層滲透率顯著降低，故而如何有效率的開(kāi)采深部煤層氣資源一直是能源領(lǐng)域及安全領(lǐng)域面臨的卡脖子問(wèn)題之一[1-2]。目前我國(guó)現(xiàn)有90%以上的煤層氣井都經(jīng)過(guò)水力壓裂改造，該技術(shù)具有擾動(dòng)范圍大、增透效果好等特點(diǎn)，在煤層氣強(qiáng)化開(kāi)采中得到廣泛應(yīng)用。而水力壓裂與斷層相互作用可能會(huì)發(fā)生滑移、地震等現(xiàn)象，故經(jīng)過(guò)水力壓裂改造的煤礦對(duì)于生產(chǎn)的安全提出了嚴(yán)峻的挑戰(zhàn)。為了保證生產(chǎn)的安全，則需要人們對(duì)壓裂的機(jī)理進(jìn)行研究，而水力壓裂的過(guò)程涉及面較廣、難以觀測(cè)，因此對(duì)于水力壓裂的研究大多都基于一些假設(shè)或者在理想狀態(tài)下進(jìn)行。

數(shù)值模擬也稱為計(jì)算機(jī)模擬，即基于計(jì)算機(jī)，結(jié)合有限元與有限容積的概念，通過(guò)數(shù)值計(jì)算和圖像顯示的方法，達(dá)到對(duì)工程問(wèn)題和物理問(wèn)題乃至自然界中各類問(wèn)題研究的目的。其被廣泛應(yīng)用于涉及流體-固體相互作用的多物理場(chǎng)，是研究水力壓裂問(wèn)題的重要手段，以往已有許多學(xué)者使用有限元模擬對(duì)水力壓裂進(jìn)行了研究。Zhao 等[4]采用數(shù)值模擬分析了裂縫擴(kuò)展方向與最大主應(yīng)力和局部巖石結(jié)構(gòu)之間的關(guān)系；Lei 等[5]運(yùn)用數(shù)值模擬方法分析了流體注入裂隙巖石的流體力學(xué)行為，分析了完整巖石脆性破壞和天然裂隙摩擦滑動(dòng)引起的誘發(fā)地震活動(dòng)的時(shí)空演變；Letham 等[6]使用有限離散元方法研究了雙軸應(yīng)力下裂縫擴(kuò)展和裂縫孔隙的變化；Min 等[7]進(jìn)行了一系列數(shù)值模擬實(shí)驗(yàn)，進(jìn)而計(jì)算了不同加載條件下模擬裂隙巖體滲透率的變化；Zhang 等[8]通過(guò)數(shù)值模擬分析了裂縫的連通性對(duì)于巖體變形和滲透率的影響。

在研究壓裂造成的應(yīng)力和滲透率變化方面，Wang 等[9]對(duì)不同注水壓力和不同應(yīng)力大小邊界等壓裂條件下的煤層透氣性已有研究；Liang 等[10]研究了注水壓力對(duì)煤層應(yīng)變的影響。除上述條件外，煤層的吸附誘導(dǎo)也會(huì)導(dǎo)致煤層發(fā)生應(yīng)力應(yīng)變變化[11]。煤層在不同深度時(shí)，受到的水平應(yīng)力與垂直應(yīng)力的比例會(huì)發(fā)生變化。隨著深度的增加，水平應(yīng)力與垂直應(yīng)力的比值會(huì)由高到低，最后趨于1 個(gè)穩(wěn)定值，應(yīng)力比值變化則會(huì)對(duì)裂縫擴(kuò)展帶來(lái)影響?；诖?，同時(shí)考慮到外部壓力和煤層自身的吸附變形，構(gòu)建流—固耦合模型，通過(guò)數(shù)值模擬軟件對(duì)巖層施加不同比例的水平應(yīng)力和垂直應(yīng)力，以模擬不同深度條件下的水力壓裂，研究注壓孔附近的應(yīng)力場(chǎng)變化，對(duì)研究不同深度的裂隙擴(kuò)展情況提供一些參考。

1 幾何模型與本構(gòu)方程

1.1 幾何模型

本研究采用有限元軟件COMSOL Multiphy 模擬水力壓裂擴(kuò)展問(wèn)題，運(yùn)用流-固耦合方法來(lái)模擬井筒附近應(yīng)力場(chǎng)分布。在本次模擬中，假設(shè)目的層無(wú)限大、遠(yuǎn)場(chǎng)地應(yīng)力均勻分布、巖石為均質(zhì)的彈性體、沒(méi)有天然裂縫分布等。為達(dá)到分析應(yīng)力場(chǎng)的目的，采用二維平面模型進(jìn)行應(yīng)力分析，2D 模型由任意選擇的垂直于水平注壓孔孔軸的橫截面組成。據(jù)上述條件建立的力學(xué)模型如圖1，模型的大小為2 m×2 m，注水孔半徑為0.1 m。

1.2 本構(gòu)方程

由于水力壓裂主要是煤體變形和瓦斯運(yùn)移共同作用的結(jié)果，所以在數(shù)值模擬中通過(guò)求解煤體變形控制方程和氣體平衡方程來(lái)實(shí)現(xiàn)水力壓裂中的流-固耦合。由煤的總應(yīng)變方程、Langmiur 體積應(yīng)變方程、平衡方程，可得到煤的變形控制方程[12]：

式中：ρg為氣體密度；qg為達(dá)西速度矢量；Qs為氣體源；t 為時(shí)間；m 為氣體含量。

Zhang 等[12]考慮氣體的吸附，對(duì)式（2）進(jìn)行了擴(kuò)展，得到：

式中：pa為1 個(gè)大氣壓，101.325 kPa；ρc為煤的密度；? 為孔隙度；k 為煤的滲透率；μ 為氣體的動(dòng)力黏度；VL為朗繆爾體積常數(shù)。

進(jìn)一步地，Zhang 等提出了一般孔隙度模型[12]：

針對(duì)上式中的εs（p），Liu 和Harpalani 提出了更為簡(jiǎn)便的吸附應(yīng)變模型[13]，即：

式中：R 為氣體常數(shù)；T 為儲(chǔ)層溫度；EA為固體膨脹模量；V0是氣體摩爾體積，22.4 L/mol。

由式（4）和式（5），可得? 對(duì)t 的偏導(dǎo)數(shù)為：

將式（6）代入式（3），得到最終的耦合方程為：

式（7）中，等式左邊第1 項(xiàng)和第2 項(xiàng)分別代表孔隙度變化和孔隙壓力變化引起的瓦斯流量變化，等式右邊2 項(xiàng)分別表示氣體源以及煤體力學(xué)變形引發(fā)的瓦斯流量變化。式（1）～式（7）即為模型的控制方程。

2 仿真模擬與相關(guān)參數(shù)

COMSOL Multiphy 軟件作為有限元分析、求解和仿真模擬軟件，在物理學(xué)和工程學(xué)中有著廣泛的應(yīng)用。特別是對(duì)于多物理場(chǎng)耦合現(xiàn)象而言，其可通過(guò)對(duì)多個(gè)耦合偏微分方程組的求解來(lái)實(shí)現(xiàn)。因此使用COMSOL Multiphy 的內(nèi)置求解器來(lái)求解指定邊界和初始條件下的偏微分方程，在COMSOL Multiphy中添加邊界載荷模擬應(yīng)力與注水壓力，使用固體力學(xué)模塊處理煤的變形控制方程，使用PDE 模塊處理氣體流量控制方程，結(jié)合本構(gòu)方程進(jìn)行數(shù)值求解。

把600 m 作為深淺度的分界線，將600 m 以上的位置認(rèn)為是較淺的位置，將600 m 以下的位置認(rèn)為是較深的位置。在較淺的深度（600 m 以上）上，水平應(yīng)力往往超過(guò)垂直應(yīng)力，但是在比較深的地方（600 m 以下），水平應(yīng)力則小于垂直應(yīng)力。600 m 以下，水平應(yīng)力與垂直應(yīng)力的比值從600 m 以上的大于1 下降到趨近于0.6 的位置[14]。不同的壓裂條件會(huì)影響注水孔周圍應(yīng)力從而影響裂縫的擴(kuò)展效果，因此模擬分析了不同壓裂條件下注水孔周圍的應(yīng)力分布，結(jié)合上文中提到的不同深度應(yīng)力的比值，對(duì)上述模型進(jìn)行了6 組模擬，6 組模擬中除修改應(yīng)力大小外其余條件均不做改變。模擬中使用的物理參數(shù)見(jiàn)表1，模擬方案見(jiàn)表2。

表1 模擬所用物理參數(shù)[15]Table 1 Physical parameters used by the model

表2 模擬方案Table 2 Simulation schemes

3 應(yīng)力比值對(duì)裂紋轉(zhuǎn)向的影響

通過(guò)修改不同的邊界參數(shù)來(lái)對(duì)已建立的模型進(jìn)行數(shù)值求解。Mises 應(yīng)力分布如圖2。以下模擬中注壓孔壓力均為15 MPa。

1）水平應(yīng)力為15 MPa，垂直應(yīng)力5 MPa。此時(shí)水平應(yīng)力與垂直應(yīng)力的比值為3∶1，應(yīng)力分布如圖2（a），高應(yīng)力主要在注水孔左右兩側(cè)（水平方向）集中，觀察1/4 的區(qū)域可以發(fā)現(xiàn)在10.3°～67.6°相較其他區(qū)域有更高的應(yīng)力分布，且該區(qū)域覆蓋角度大、面積廣，注壓孔沿某一方向開(kāi)裂后更容易向高應(yīng)力分布區(qū)域延伸，且裂紋長(zhǎng)度會(huì)較長(zhǎng)。

2）水平應(yīng)力為10 MPa，垂直應(yīng)力為5 MPa。此時(shí)水平應(yīng)力與垂直應(yīng)力的比值為2∶1，應(yīng)力分布如圖2（b），高應(yīng)力同樣集中在注水孔的左右兩側(cè)（水平方向），在23.6°～60.4°之間有較其他區(qū)域更高的應(yīng)力分布，高應(yīng)力區(qū)域較2（a）更靠近水平方向，說(shuō)明當(dāng)水平應(yīng)力與垂直應(yīng)力比值下降時(shí)，裂紋的擴(kuò)展也更靠近水平方向。

3）水平應(yīng)力為7 MPa，垂直應(yīng)力為5 MPa。兩者比值在1～2 之間時(shí)，應(yīng)力分布如圖2（c），高應(yīng)力區(qū)域依然分布在注水孔的左右兩側(cè)（水平方向）。相較于上兩者，沒(méi)有明顯在周圍產(chǎn)生高應(yīng)力的分布區(qū)域，此時(shí)裂縫的主要延伸方向?yàn)樗椒较颉?/p>

4）當(dāng)水平應(yīng)力為7 MPa，垂直應(yīng)力也為7 MPa。即垂直應(yīng)力與水平應(yīng)力比值為1∶1 時(shí)，如圖2（d），此時(shí)高應(yīng)力區(qū)域均勻分布在注水孔周圍，裂紋向四周均勻擴(kuò)展。

5）當(dāng)水平應(yīng)力為7 MPa，垂直應(yīng)力為10 MPa。此時(shí)水平應(yīng)力與垂直應(yīng)力的比值在0.6～1 之間，如圖2（e），與2（a）和2（b）相反，高應(yīng)力集中在注水孔上下（垂直方向）方，在30.8°～65.5°有較其他區(qū)域更高的應(yīng)力分布，高應(yīng)力區(qū)域張開(kāi)的角度小，裂紋傾向于向該區(qū)域延伸。

6）當(dāng)水平應(yīng)力為7 MPa，垂直應(yīng)力增大到12 MPa。此時(shí)水平應(yīng)力垂直應(yīng)力之比趨近于0.6，應(yīng)力分布如圖2（f），此時(shí)高應(yīng)力區(qū)域同樣主要分布在注水孔上下（垂直方向）方，高應(yīng)力區(qū)在23.1°～69.9°之間分布，較2（e）而言，水平應(yīng)力與垂直應(yīng)力比值降低，高應(yīng)力區(qū)張開(kāi)的角度較大，裂紋延伸的角度更大，且會(huì)更加靠近垂直方向。

監(jiān)測(cè)線及線上壓力分布如圖3。為了更好研究應(yīng)力的變化規(guī)律，在模型中添加如圖3（a）所示的截線，以左下角頂點(diǎn)為坐標(biāo)原點(diǎn)，截線端點(diǎn)坐標(biāo)分別為（0.75，1）與（1，1.25），監(jiān)測(cè)截線在水平應(yīng)力/垂直應(yīng)力分別為3、2、1.4 時(shí)（即水平應(yīng)力15 MPa，垂直應(yīng)力5 MPa；水平應(yīng)力10 MPa，垂直應(yīng)力5 MPa；水平應(yīng)力7 MPa，垂直應(yīng)力5 MPa，注壓孔壓力均為15 MPa 時(shí)）的應(yīng)力分布，將結(jié)果繪制在圖3（b）上。

由截線上應(yīng)力分布可知，當(dāng)比值越大時(shí)，截線上應(yīng)力的峰值會(huì)向右上方移動(dòng)，說(shuō)明當(dāng)水平應(yīng)力與垂直應(yīng)力比值大于1 時(shí)，高應(yīng)力區(qū)會(huì)隨著比值增加逐漸向垂直方向靠攏；相對(duì)的，當(dāng)比值小于1 時(shí)，高應(yīng)力區(qū)會(huì)隨著比值降低往水平方向靠攏。

4 實(shí)驗(yàn)驗(yàn)證與討論

在Chen 等[16]的壓裂實(shí)驗(yàn)中，使用了高黏度壓裂液引發(fā)水力壓裂，分別在射孔方向0°、15°、30°、45°、60°、75°時(shí)，將壓裂液以恒定速率壓入鉆孔，直到裂縫開(kāi)始并在破裂后擴(kuò)展，然后記錄壓力下降，實(shí)驗(yàn)結(jié)束后，將每個(gè)塊體沿裂縫平面切開(kāi)觀察了裂縫的幾何形狀。選取實(shí)驗(yàn)中射孔方向0°時(shí)不同垂直應(yīng)力與水平應(yīng)力比值時(shí)的擊穿壓力，與在第3 章中模型所得的結(jié)果進(jìn)行對(duì)比，不同應(yīng)力比值下的擊穿壓力如圖4。

將2 組數(shù)據(jù)進(jìn)行描點(diǎn)擬合后可以看出，雖然因?yàn)樽⑷雺毫Φ牟煌瑢?dǎo)致模擬與實(shí)驗(yàn)擊穿壓力的大小有所區(qū)別，但是隨著垂直應(yīng)力和水平應(yīng)力比值的降低，實(shí)驗(yàn)與模擬中壓裂的擊穿壓力都近似發(fā)生線性降低，這表明模擬得出的數(shù)值符合實(shí)驗(yàn)規(guī)律，在不考慮其他條件的影響下，起裂開(kāi)始后，可以通過(guò)模型模擬得出的應(yīng)力分布推測(cè)裂縫的擴(kuò)展情況。裂縫延伸方向如圖5。

模擬結(jié)果表明，在非裂縫性儲(chǔ)層中（儲(chǔ)層中不考慮天然裂縫存在），水力裂縫的擴(kuò)展應(yīng)該為對(duì)稱型擴(kuò)展，注壓孔周圍應(yīng)力的變化會(huì)導(dǎo)致?lián)舸毫Φ淖兓?，且垂直?yīng)力與水平應(yīng)力比值的變化在裂縫的擴(kuò)展中扮演1 個(gè)重要的角色。由圖5 可知，當(dāng)位于較淺位置，水平應(yīng)力與垂直應(yīng)力比值大于1 時(shí)，高應(yīng)力主要分布在水平方向上；除此之外，深度越淺（即應(yīng)力比值越大），裂縫越容易向垂直方向發(fā)生轉(zhuǎn)向，在沒(méi)有天然裂縫等條件影響下，裂縫會(huì)向兩側(cè)對(duì)稱擴(kuò)展。隨著深度的增加，當(dāng)水平應(yīng)力與垂直應(yīng)力比值為1時(shí)，應(yīng)力均勻分布，表明裂縫此時(shí)可能會(huì)向各個(gè)均勻擴(kuò)展；當(dāng)水壓足夠大時(shí)，會(huì)在周圍同時(shí)產(chǎn)生多個(gè)裂縫。而當(dāng)比值繼續(xù)減小，即水平應(yīng)力與垂直應(yīng)力比值小于1 時(shí)，高應(yīng)力分布會(huì)由水平方向轉(zhuǎn)換為垂直方向，裂縫會(huì)主要在垂直方向上進(jìn)行擴(kuò)展；與較淺位置時(shí)相反，此時(shí)比值越小，裂縫越容易朝水平方向發(fā)生轉(zhuǎn)向。李傳亮等[17]已經(jīng)詳細(xì)闡述過(guò)壓裂井裂縫的沿伸方向，解釋了為何地層中最容易產(chǎn)生的裂縫是平行于最大水平地應(yīng)力方向的裂縫，最不容易產(chǎn)生的裂縫是垂直于最大地應(yīng)力水平方向的裂縫，也與本研究的模擬結(jié)果保持一致。

5 結(jié) 論

1）在較淺的位置，即水平應(yīng)力與垂直應(yīng)力的比值大于1 時(shí)，裂縫主要沿最大主應(yīng)力方向即水平方向上開(kāi)裂。在較深的位置，即水平應(yīng)力與垂直應(yīng)力的比值小于1 時(shí)，裂縫主要沿垂直方向開(kāi)裂。

2）在垂直應(yīng)力不變的情況下，垂直應(yīng)力與水平應(yīng)力的比值增大，裂隙擊穿壓力隨比值變化近似發(fā)生線性降低的變化。

3）在較淺位置時(shí)，水平應(yīng)力與垂直應(yīng)力比值越大，裂隙向垂直方向偏轉(zhuǎn)角度越大；深度較深時(shí)，水平應(yīng)力與垂直應(yīng)力比值越小，裂隙向水平方向偏轉(zhuǎn)角度越大。