劉先珊，曾南豆，李濤，張立君

(1.山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點(diǎn)實驗室(重慶大學(xué))，重慶，400045；2.重慶大學(xué) 土木工程學(xué)院，重慶，400045；3.湖北省漳河工程管理局，湖北 荊門，448156)

頁巖氣兼具高效、清潔、儲量豐富的優(yōu)點(diǎn)，其作為一種新型非常規(guī)油氣資源的地位將日益突出[1]。全球頁巖氣資源量達(dá)到456×1012m3，其中可采資源量僅有187×1012m3。我國頁巖氣資源豐富，氣儲量達(dá)到134×1012m3，其中技術(shù)可采儲量將近25×1012m3。然而，儲層頁巖層理明顯且天然裂縫發(fā)育，賦存環(huán)境復(fù)雜，具有物性差、孔隙度低、滲透率低等特點(diǎn)，水力壓裂作為頁巖氣商業(yè)開發(fā)的核心技術(shù)，其強(qiáng)擾動誘導(dǎo)水力裂縫的擴(kuò)展路徑復(fù)雜多變，加上缺乏準(zhǔn)確有效的現(xiàn)場監(jiān)測手段，裂縫網(wǎng)絡(luò)的形成規(guī)律難以準(zhǔn)確把握。探明儲層頁巖水力壓裂裂縫的起裂、擴(kuò)展和空間展布特征，可為水力壓裂設(shè)計、儲層改造及采收率提升提供科學(xué)依據(jù)。

近年來，一些學(xué)者針對儲層頁巖的水力壓裂過程開展了物模試驗和理論分析，取得了較豐富的研究成果。李芷等[2-3]基于頁巖真三軸壓裂試驗，探討了層理對壓裂縫網(wǎng)形成的重要影響；LIN等[4]研究了不同層理角度的頁巖起裂力，闡明層理角度對裂縫擴(kuò)展模式的影響顯著；XU 等[5]基于相似原理制作巖樣，研究層理引起的各向異性對縫網(wǎng)形成的影響；ZHAO等[6]基于頁巖水力壓裂試驗得到密集裂縫網(wǎng)絡(luò)出現(xiàn)的泵入速率極限值，揭示復(fù)雜裂縫網(wǎng)絡(luò)形成的主控因素；LI等[7]開展脈沖水力壓裂試驗，深入研究脈沖頻率對裂縫形態(tài)的影響；HU 等[8]基于水力壓裂試驗研究天然裂縫與水力裂縫的相互作用，闡明巖體剪切滑移程度與天然裂縫傾角和摩擦角密切關(guān)聯(lián)；LIU等[9]基于真三軸水力壓裂試驗研究預(yù)制裂縫參數(shù)、壓裂液性質(zhì)和巖石力學(xué)特性對定向水力壓裂過程的影響；GUNAYDIN 等[10]開展多水力裂縫壓裂試驗，闡明裂縫間距對水力裂縫擴(kuò)展的抑制作用。在此基礎(chǔ)上，一些學(xué)者提出半解析模型定量研究頁巖壓裂裂縫的擴(kuò)展過程，如CHENG等[11]發(fā)展了一種基于PKN 的非平面裂縫擴(kuò)展模型，但只涉及單一裂縫的擴(kuò)展，不能準(zhǔn)確描述復(fù)雜儲層壓裂縫網(wǎng)的形成；孫可明等[12]利用復(fù)變函數(shù)和保角變換方法，分析裂紋沿各方向擴(kuò)展的水壓力，研究水力裂紋與層理斜交的擴(kuò)展判據(jù)及擴(kuò)展規(guī)律；SUN 等[13]建立頁巖儲層水力壓裂三維計算模型，研究不同層理方向的裂縫延伸規(guī)律；ZOU 等[14]提出基于三維離散元的新裂縫網(wǎng)絡(luò)模型，研究層理對縫網(wǎng)形成的影響；王志榮等[15]根據(jù)最大周向應(yīng)力準(zhǔn)則及斷裂力學(xué)原理建立起始裂縫轉(zhuǎn)向模型，得出裂縫轉(zhuǎn)向與射孔角度、地應(yīng)力差的相關(guān)性。實際測量監(jiān)測、試樣尺寸及離散性、試驗儀器的局限性、試驗過程的重復(fù)性以及理論模型的簡化，不易獲得巖樣內(nèi)部多尺度缺陷的改變，難以準(zhǔn)確解釋上述水力壓裂儲層頁巖的多尺度缺陷演化、流-固耦合機(jī)制及裂縫形成規(guī)律。

由此，較多學(xué)者借助數(shù)值模擬技術(shù)開展儲層頁巖的水力壓裂過程研究，如擴(kuò)展有限元方法XFEM 得到了廣泛應(yīng)用，DAHI-TALEGHANI 等[16]模擬了均質(zhì)地層水力裂縫與天然裂縫的相交，運(yùn)用Picard 方法迭代解耦巖石位移場和縫內(nèi)流場；GORDELIY等[17]針對單條水力裂縫的非平面擴(kuò)展，引入隱式水平集方法確定裂縫尖端位置，推進(jìn)了XFEM 在水力壓裂模擬中的應(yīng)用；盛茂等[18]基于XFEM考慮滲流應(yīng)力耦合效應(yīng)，模擬了壓裂裂縫的擴(kuò)展演化；WANG 等[19]研究了水力裂縫間距及壓裂順序?qū)p網(wǎng)形成的影響，并將XFEM 拓展到復(fù)雜縫網(wǎng)研究；張進(jìn)科等[20]基于XFEM 建立全耦合水力壓裂模型，研究壓裂順序和裂縫間距的裂縫擴(kuò)展規(guī)律，闡明順序壓裂能夠增大裂縫寬度并有效提高采收率。可見，XFEM在水力裂縫擴(kuò)展模擬中優(yōu)勢明顯，但裂紋追蹤仍依靠水平集函數(shù)，該函數(shù)在很大程度上影響裂縫擴(kuò)展路徑的準(zhǔn)確性。為了解決如上問題，顆粒流軟件(particle flow code,PFC)能避免裂縫尖端的奇異性，有助于細(xì)致描述壓裂過程中的巖體組構(gòu)改變、流固耦合效應(yīng)及裂縫擴(kuò)展過程，其中細(xì)觀尺度的多場耦合機(jī)理研究是關(guān)鍵。THALLAK等[21]為了更準(zhǔn)確地再現(xiàn)流固耦合過程，開發(fā)了一種流固耦合顆粒流離散元模型，研究水力壓裂無黏結(jié)巖石的損傷破裂過程；BRUNO[22]改進(jìn)上述多場耦合算法，為網(wǎng)絡(luò)單元添加了應(yīng)力相關(guān)的幾何結(jié)構(gòu)和傳輸特性；ALBUSAIDI[23]基于上述模型模擬花崗巖水力壓裂裂縫擴(kuò)展過程，闡明水力裂縫的擴(kuò)展以拉張破壞為主；WANG 等[24]基于PFC 流固耦合模型模擬水力壓裂裂縫的擴(kuò)展過程，對比現(xiàn)場實測數(shù)據(jù)，驗證了顆粒流離散元流固耦合模型的可行性；倪小東等[25]深入研究PFC-CFD求解法，進(jìn)行裂紋擴(kuò)展的動態(tài)分析，驗證該方法的可行性；ZENG等[26]基于新的PFC2D 管域模型，從細(xì)觀角度研究水力壓裂裂縫形態(tài)的萌生與擴(kuò)展；FATAHI等[27]基于PFC軟件模擬不同尺寸及滲透性巖體的水力壓裂過程，對比室內(nèi)試驗以驗證模型的可行性；DUAN 等[28]基于PFC 軟件研究了注水參數(shù)與原位應(yīng)力條件對砂巖水力壓裂破裂壓力及裂紋擴(kuò)展的影響；YOON等[29-30]建立考慮流固耦合效應(yīng)的PFC2D 模型，模擬多水力裂縫干擾的擴(kuò)展過程。

為了充分描述巖體裂隙與孔隙流動的差異，本文作者從細(xì)觀尺度描述巖體壓裂過程中的流固耦合效應(yīng)，提出改進(jìn)的流固耦合算法，對比層理頁巖的應(yīng)力-應(yīng)變及滲透試驗結(jié)果，標(biāo)定裂紋優(yōu)先導(dǎo)流模型的合理性，闡明改進(jìn)耦合算法的可靠性。進(jìn)一步，基于改進(jìn)流固耦合方法建立儲層頁巖水力壓裂模型，從細(xì)觀層面研究不同主控因素影響的儲層頁巖水壓裂過程，描述水力壓裂裂縫的擴(kuò)展特性，其研究成果對優(yōu)化水力壓裂參數(shù)及提高頁巖氣開采效率具有重要的工程價值。

1 壓裂誘發(fā)頁巖破壞的改進(jìn)流固耦合算法

1.1 基于PFC2D的流固耦合算法

基于PFC2D 描述頁巖內(nèi)流體與顆粒間的耦合作用效應(yīng)，通過循環(huán)函數(shù)在整個模型中建立空隙與空隙周圍顆粒的關(guān)系，利用額外屬性的設(shè)定將顆粒與對應(yīng)圍成的空隙進(jìn)行綁定，在顆粒閉合鏈的中央形成的空隙便成了流體域(如圖1所示)。由于綁定過程中，2個顆粒無法構(gòu)成閉合鏈，即無法形成流體域，需要刪除上述類型的顆粒。然后，在流體域之間需要設(shè)定通道來實現(xiàn)流動，通道又是通過顆粒接觸來實現(xiàn)，同樣采用上述的額外函數(shù)功能，將接觸添加管道的額外屬性，由此定義了“流動管道”，并將流動通道和流體域關(guān)聯(lián)起來建立整個滲流網(wǎng)絡(luò)。流動通道相當(dāng)于一個平行板通道，可基于立方定律獲得單位時間內(nèi)該流動通道的流量q：

圖1 PFC2D 中流體域和流動通道示意圖Fig.1 Sketch map of fluid domain and flow channel in PFC2D

式中：μ為流體動力黏滯系數(shù)；a為流動管道的開度；L為兩相關(guān)顆粒半徑和，作為管道長度時L=r1+r2，r1和r2分別為兩相關(guān)顆粒的半徑；P1和P2分別為兩空隙區(qū)域內(nèi)的流體壓力；t為數(shù)值模型試樣的厚度，二維數(shù)值模擬取t=1；q為單位時間內(nèi)兩流體域間傳遞的流量。假定流動管道的初始開度為a0，表示模型完好且無荷載作用時的開度，當(dāng)流動通道承受正向壓力F作用后，其流動管道的開度a為

式中：a0為接觸法向作用力為零時的殘余開度；F0為流動通道開度減少為殘余開度a0的一半時，流動管道接觸位置法向方向的力。F為0，則開度a=a0；當(dāng)法向力F=F0時，開度a=a0/2，若F無限大，則開度a無限接近零。

當(dāng)流體域處的顆粒受拉時，流動管道的開度a為

式中：g為對應(yīng)兩顆粒表面間的距離；m為模型顆粒與實際巖石微粒之間的轉(zhuǎn)換系數(shù)。

如圖2所示流體流動結(jié)構(gòu)中的流體域，可得相連通的流動管道在單位時間內(nèi)流入的總流量為∑q，則流體域中顆粒流體壓力的改變量為

式中：Kf為流體體積模量；Vd為流體域表觀體積；為受流量改變的影響項，為受體積改變的影響項。

基于上述PFC2D模型獲得流體流動結(jié)構(gòu)網(wǎng)絡(luò)，進(jìn)行相應(yīng)的流固耦合計算。首先，流域內(nèi)的水壓會在每次模型計算循環(huán)時對顆粒產(chǎn)生力作用，同時每次循環(huán)會對流體域內(nèi)的體積進(jìn)行更新，讓其壓力產(chǎn)生改變，且在荷載作用下的顆粒間接觸力實時改變，引起了流動管道開度的不斷改變，即對流體運(yùn)動特性產(chǎn)生影響。實際計算過程中，水壓看作均勻分布，流體域看作周圍顆粒接觸點(diǎn)的連線封閉區(qū)，則顆粒j所受該區(qū)域的力為

式中：P為流體域內(nèi)流體壓力；nj為顆粒j兩接觸點(diǎn)連線的單位法向量；s為該連線的長度(如圖2所示)。

圖2 流體流動結(jié)構(gòu)示意圖Fig.2 Sketch map of fluid flow structure

計算中，為了更合理確定計算時間步，首先基于式(1)，結(jié)合域內(nèi)的壓力改變量計算流速q：

式中：N為該流體域所連通的管道數(shù)量；k為滲透率；R為所屬流體域顆粒的平均半徑。根據(jù)式(4)，可以根據(jù)壓力波動ΔPP得到一定時間Δt內(nèi)產(chǎn)生的壓力響應(yīng)ΔPR：

當(dāng)壓力響應(yīng)ΔPR小于壓力波動ΔPP時，提前穩(wěn)定計算系統(tǒng)。令ΔPR=ΔPP，求得臨界時間Δt：

通過式(8)可獲得多個臨界時間Δt，計算過程中選取的計算時步比臨界時間要小。

1.2 改進(jìn)的流固耦分析方法

上述流固耦合模型均是基于顆粒與顆粒間的接觸建立頁巖內(nèi)部流體流動管道網(wǎng)絡(luò)，該模型由顆粒接觸定義了流體管道，當(dāng)顆粒間黏結(jié)接觸破壞之前，各流體的區(qū)域與管道間的力學(xué)計算是互通，但當(dāng)顆粒間接觸不存在時，之前所定義的管道將不復(fù)存在，接觸破壞處的流固耦合計算將出錯。一些學(xué)者提出相應(yīng)的處理措施，ZHAO等[31]假設(shè)在生成裂隙后(顆粒與顆粒之間的黏結(jié)遭到了破壞)，原流動通道所連接的2 個流體域內(nèi)的流體交換是瞬間完成的，原兩者流體域內(nèi)部的水壓力進(jìn)行中和，取值某水壓力Pw，即為顆粒黏結(jié)遭到破壞以前原兩流體域內(nèi)部水壓的平均值；ZENG等[26]對初始流固耦合算法進(jìn)行了改進(jìn)，基于顆粒流離散元法進(jìn)行流固耦合計算時，接觸破壞生成裂縫的過程會被裂縫生成函數(shù)記錄下來，其中包括兩顆粒的位置、指針等信息，很好地區(qū)分接觸破壞與未破壞的管道區(qū)域，避免在裂縫產(chǎn)生后由于顆粒間接觸破壞引起流動網(wǎng)絡(luò)的破壞，改進(jìn)算法能很好地解決裂隙產(chǎn)生后流動管道網(wǎng)絡(luò)的工作問題。因此，當(dāng)裂隙產(chǎn)生時，除了將記錄裂隙處位置顆粒信息記錄于裂隙系統(tǒng)外，還需將以前該處的信息從黏結(jié)系統(tǒng)中刪除。為了更方便地運(yùn)用上述流體流動網(wǎng)絡(luò)，本文進(jìn)一步完善上述流固耦合算法，定義流體流動管道并不使用接觸，而僅僅使用接觸遍歷搜索周圍孔隙來定義流體域，則有：

1) 用經(jīng)典流固耦合的算法遍歷顆粒間的接觸搜索空隙區(qū)域定義流體域，并依據(jù)接觸的位置定義一個geometry集合體；

2) geometry 集合體的基礎(chǔ)上定義流體通道，并設(shè)置額外屬性crack_active；

3) 過生成裂隙的函數(shù)記錄顆粒信息，并給該處的geometry 的額外屬性crack_active(取0 或1，1代表黏結(jié)破壞，0代表黏結(jié)未破壞)賦值；

4) 據(jù)額外屬性crack_active 的值進(jìn)行黏結(jié)破壞后的流體網(wǎng)絡(luò)運(yùn)行。

通過geometry 的方式建立流體流動管道后，顆粒間黏結(jié)破壞與否將與流體流動管道的存在與否無聯(lián)系，即任何情況下的流動網(wǎng)絡(luò)都保持完整。因此，由一個新系統(tǒng)代替之前的流體流動管道系統(tǒng)，當(dāng)模型產(chǎn)生新裂縫時，新系統(tǒng)會根據(jù)有無裂縫產(chǎn)生對該位置賦值不同參數(shù)，以此判斷模型中該位置是孔隙流還是裂隙流，對裂縫產(chǎn)生處的管道開度乘一個放大系數(shù)，以此表征裂隙流，由此完成流固耦合過程中不同類型滲流過程的表征。具體的計算流程如圖3所示。

圖3 遍歷獲取流動通道和流體域信息流程圖Fig.3 Flow chart of getting flow channel and fluid domain

1.3 模型驗證

對現(xiàn)場采集的石柱縣漆遼海相志留系統(tǒng)龍馬溪組頁巖加工形成圓柱形試樣，其直徑為50 mm，高為100 mm，開展常規(guī)三軸壓縮-滲流試驗，試驗數(shù)據(jù)用于數(shù)值計算模型的驗證?？紤]流固耦合計算工作量大，選擇顆粒粒徑范圍為0.64～1.02 mm，顆粒密度為2 650 kg/m3，首先導(dǎo)入線性接觸模型，基于contact model命令將顆粒與顆粒間的接觸變?yōu)槠叫叙そY(jié)接觸，再通過geometry 設(shè)置一定間距的DFN(discrete fracture network，離散裂隙網(wǎng)絡(luò))模型，完成整個儲層頁巖標(biāo)準(zhǔn)模型的構(gòu)建，如圖4所示，共生成1 938個顆粒。

圖4 基于PFC的頁巖物理模型Fig.4 Physical model of shale rock based on PFC

進(jìn)一步建立儲層頁巖滲透模型，如圖5所示。數(shù)值模型的頂部位置5 mm 區(qū)間內(nèi)施加進(jìn)口壓力P1，底部5 mm 區(qū)間內(nèi)施加出口壓力P2，用fish 函數(shù)遍歷模型中央10 mm 區(qū)域內(nèi)的流體域以及流體管道，分別記錄該區(qū)域的平均孔壓pj和總流量qj(規(guī)定從上往下為正向)。跟蹤每個時間的平均孔壓pj和總流量qj，當(dāng)pj和qj都保持恒定值時，說明達(dá)到了穩(wěn)態(tài)滲流狀態(tài)。依據(jù)達(dá)西定律，得到模型滲透系數(shù)K：

圖5 頁巖穩(wěn)態(tài)法滲透率測試模型Fig.5 Steady-state model of shale permeability

式中：Lg為模型試樣流體進(jìn)口端與流體出口端間的模型試樣高度；A為模型試樣的橫截面積，二維模型中取厚度為1，即橫截面積的值即為試樣寬度；γ為流體的重度。

根據(jù)式(9)，得到模型試樣的滲透率k：

根據(jù)三軸力學(xué)試驗得到的偏應(yīng)力-應(yīng)變曲線，采用試錯法進(jìn)行參數(shù)標(biāo)定。表1所示為細(xì)觀力學(xué)參數(shù)。建立層理傾角0°的頁巖數(shù)值模型，運(yùn)用改進(jìn)的流固耦合方法進(jìn)行滲透試驗?zāi)M，對比數(shù)值模擬結(jié)果與試驗結(jié)果，得到表2所示的水力學(xué)參數(shù)，進(jìn)一步驗證計算參數(shù)的合理性。圖6所示為圍壓10 MPa 組合滲透壓5 MPa 的層理角度0°頁巖滲透試驗與數(shù)值模擬對比。從圖6可知：滲透率曲線初期緩慢降低，流體通過頁巖基質(zhì)系統(tǒng)管道流動，隨著顆粒間的接觸壓應(yīng)力增大，流體流動管道開度降低；繼續(xù)增大顆粒間應(yīng)力達(dá)到峰值應(yīng)力，試樣內(nèi)部開始出現(xiàn)大裂縫，此時裂縫處的管道計算會進(jìn)入裂縫系統(tǒng)。圖6中顯示數(shù)值模擬曲線與試驗曲線吻合較好，驗證了改進(jìn)流固耦合算法的可靠性和合理性。

圖6 圍壓10 MPa組合滲透壓5 MPa的層理角度0°頁巖滲透試驗與數(shù)值模擬對比Fig.6 Comparison of tests and simulations of shale with bedding angle 0° considering combination of confining pressure valued 10 MPa and water pressure of 5 MPa

表1 頁巖數(shù)值模型的細(xì)觀參數(shù)Table 1 Micro parameters of numerical simulation of shale

表2 頁巖滲透模型水力學(xué)參數(shù)Table 2 Numerical hydraulic parameters of shale model

為了進(jìn)一步驗證該數(shù)值模型的合理性和細(xì)觀參數(shù)的準(zhǔn)確性，模擬圍壓5 MPa組合滲透壓4 MPa下層理傾角45°頁巖的滲透演化過程，并與試驗結(jié)果進(jìn)行對比，結(jié)果如圖7所示。從圖7可知：低圍壓與低滲透壓下，儲層頁巖滲透率演化曲線與圖7顯示的曲線一致。實際水力壓裂誘發(fā)頁巖損傷破裂成縫，其間的流體運(yùn)動由孔隙流動到裂隙流動，本文提出的改進(jìn)流固耦合算法可用來分析不同流動結(jié)構(gòu)的流動過程，能更準(zhǔn)確地描述考慮流固耦合效應(yīng)的裂縫形成過程。

圖7 圍壓5 MPa組合滲透壓4 MPa的層理角度45°頁巖滲透試驗與數(shù)值模擬對比Fig.7 Comparison of tests and simulations considering combination of confining pressure valued 5 MPa and water pressure of 4 MPa

2 頁巖水力壓裂模型構(gòu)建及驗證

2.1 水力壓裂裂縫開度確定

根據(jù)本文提出的改進(jìn)流固耦合算法，通過設(shè)定額外屬性crack_active(取1 或0，1 為裂隙系統(tǒng)，0 為黏結(jié)系統(tǒng))來區(qū)別管道流和裂隙流，并根據(jù)該管道所處的狀態(tài)選擇不同的開度計算公式。當(dāng)crack_active為0時，該管道的開度a通過接觸上的法向力p來計算：

其中：a0B和a0C分別為頁巖基質(zhì)和裂縫的殘余開度；σt為顆粒黏結(jié)的抗拉強(qiáng)度；pclB為頁巖基質(zhì)閉合強(qiáng)度；CB為頁巖基質(zhì)的應(yīng)力敏感系數(shù)。

圖8所示為流動通道的開度-法向應(yīng)力作用示意圖。圖9所示為顆粒黏結(jié)沒有破壞時的流動管道開度與顆粒間法向力的關(guān)系。從圖9可知：當(dāng)顆粒間接觸受壓時，壓力增大導(dǎo)致對應(yīng)流動管道的開度變小，由應(yīng)力敏感系數(shù)CB控制減小速率，當(dāng)CB變小時，開度的減小速度變慢，反之變快。圖9所示的藍(lán)色線段表明若顆粒間接觸壓力較大，管道已處于閉合區(qū)段。顆粒間接觸拉力跟開度呈線性正相關(guān)，當(dāng)顆粒間接觸拉應(yīng)力增大到近黏結(jié)破壞時，則開度無限接近于破壞后的裂縫殘余開度。圖10所示為3 種黏結(jié)破壞后的情況。當(dāng)顆粒黏結(jié)遭到破壞，該處產(chǎn)生裂縫DFN，賦予geometry 額外屬性crack_active的值為1，管道通道由接觸破壞前的黏結(jié)系統(tǒng)轉(zhuǎn)變?yōu)轲そY(jié)破壞后的裂縫系統(tǒng)。在此處分為顆粒間有無壓力狀態(tài)2種情況，當(dāng)顆粒間存在壓力時，要判斷壓力與裂縫基質(zhì)強(qiáng)度的大小，若壓力并未達(dá)到裂縫的基質(zhì)強(qiáng)度，則管道的開度a依然與壓力存在一定的關(guān)系，即隨著壓力的變化而變化；當(dāng)壓力超過裂縫基質(zhì)強(qiáng)度，即管道會完全封閉，開度取為0；當(dāng)顆粒間無壓力存在時，此時管道的“通透性”將會非常大，流體的交換規(guī)模會變大，也就是說此時的管道開度a會很大，且與顆粒間距離成正比。

圖8 流動通道的開度-法向應(yīng)力作用示意圖Fig.8 Sketch map of aperture and normal stress for flow channels

圖9 流動通道的開度-法向應(yīng)力關(guān)系曲線Fig.9 Curves of aperture and normal stress for flow channels

圖10 裂縫系統(tǒng)中的流動通道的開度確定示意圖Fig.10 Determination of aperture in flow channels for cracks

當(dāng)crack_active 為1 時，組成管道通道的顆粒受壓力作用時，法向應(yīng)力p影響開度a；當(dāng)管道通道兩邊的顆粒不存在力作用時，通道兩邊顆粒間的表面距離g決定對應(yīng)管道的開度。黏結(jié)破壞后開度為

其中：a0C為裂縫殘余開度；pclC為裂縫閉合強(qiáng)度；CC為裂縫的應(yīng)力敏感系數(shù)；m為尺度放大系數(shù)。

圖11所示為顆粒黏結(jié)破壞后流動管道的開度隨顆粒接觸應(yīng)力的變化規(guī)律。從圖11可知：顆粒間的應(yīng)力為壓應(yīng)力時，壓應(yīng)力增大，開度減小，減小速率由裂縫應(yīng)力敏感系數(shù)CC控制，CC的增大導(dǎo)致開度減小的速率增大，而CC的減小反過來會引起開度減小的速率變小。圖11中藍(lán)色線段表示顆粒間接觸壓應(yīng)力過大使其進(jìn)入閉合階段。圖11(b)所示為當(dāng)顆粒間沒有擠壓作用時，此處管道開度依據(jù)顆粒表面距離g來控制。

圖11 裂縫系統(tǒng)中流動通道的開度確定曲線Fig.11 Curves of aperture in flow channels for cracks

2.2 水力壓裂模型驗證

建立如圖12所示的水力壓裂模型，刪除中心半徑為4 mm區(qū)域內(nèi)顆粒，作為施加水壓或流速的注水孔，x和y方向分別施加恒定地應(yīng)力σx和σy，規(guī)定最大主應(yīng)力方向為x向。首先施加邊界應(yīng)力σx和σy，當(dāng)應(yīng)力保持恒定并達(dá)到平衡時，注水區(qū)以恒定的流量注水，當(dāng)試件開始出現(xiàn)裂縫時，不再以恒定流量注入，以恒定水壓(此壓力記錄為起裂應(yīng)力)施加并繼續(xù)計算。

圖12 水力壓裂模型Fig.12 Hydraulic fracturing model

調(diào)用PFC2D 伺服系統(tǒng)保持恒定圍壓σx=σy=1 MPa，注水孔以恒定流量8 m3/h施加水壓，直至裂縫擴(kuò)展，如圖13所示(紅色代表壓裂裂縫，灰色代表顆粒，藍(lán)色代表水壓)。圖14所示為水力壓裂過程中的孔內(nèi)壓力變化曲線(圍壓為1 MPa)。從圖14可知：中心孔內(nèi)水壓從初始水壓力逐漸增加到峰值，其峰值即為起裂壓力，隨后突然下降，主要在于水力壓裂誘導(dǎo)裂縫出現(xiàn)后，中心孔內(nèi)的流體會沿著裂縫流入其他管域內(nèi)，使得中心孔內(nèi)的水壓力驟降，其驟降程度與殘余裂縫開度密切相關(guān)。當(dāng)中心孔內(nèi)持續(xù)注入，流體持續(xù)進(jìn)入裂縫中，中心孔內(nèi)壓力逐漸減小，主要在于注水孔所在區(qū)域的水壓力增高，對注水孔周的顆粒應(yīng)力加大，孔周出現(xiàn)細(xì)小裂紋，隨著中心孔壓不斷增大，壓裂裂紋不斷向模型邊界處擴(kuò)展，并與內(nèi)部顆粒接觸較薄弱的區(qū)域形成貫通，最后延伸到邊界，形成圖13所示的3 條主干裂縫相交于注水孔中心，其中2 條裂縫向模型兩底角方向延伸至邊緣，3 條裂縫間的夾角都近似為120°。另外，計算得到水力壓裂過程中的顆粒速度分布，如圖15所示。從圖15可知：裂縫擴(kuò)展處相鄰顆粒的顏色差異顯著，說明速度相差大，結(jié)合圖16所示的顆粒位移圖，可知顆粒位移量與位移方向差異顯著的位置易出現(xiàn)裂縫，闡明顆粒相對位移增大是導(dǎo)致顆粒間接觸黏結(jié)破壞并形成裂縫的主要原因。

圖13 水力壓裂裂縫圖Fig.13 Cracks caused by hydraulic fracturing

圖14 水力壓裂過程中的孔內(nèi)壓力變化曲線Fig.14 Variation of water pressure in hole caused by hydraulic fracturing

圖15 顆粒速度圖Fig.15 Diagram of particle velocity

圖16 顆粒位移圖Fig.16 Diagram of particle displacement

進(jìn)一步考慮不同地應(yīng)力影響的水力壓裂裂縫形成過程，考慮4 種水平應(yīng)力σx及豎直應(yīng)力σy組合：σx=5 MPa，σy=10 MPa；σx=5 MPa，σy=15 MPa；σx=10 MPa，σy=5 MPa；σx=15 MPa，σy=5 MPa，得到不同工況下的壓裂裂縫，如圖17所示。從圖17可知：壓裂裂縫沿著較大地應(yīng)力的方向擴(kuò)展延伸，主要在于大主應(yīng)力作為應(yīng)力主控方向，使得非大主應(yīng)力方向壓縮變形增大，滲流通道變窄，滲流速度減小，裂縫擴(kuò)展減弱，壓裂難度增大。另外，對比圖形內(nèi)左右側(cè)的壓裂裂縫形態(tài)，可知右側(cè)壓裂裂縫延伸長度明顯長于左側(cè)壓裂裂縫的延伸長度，說明地應(yīng)力差對裂縫擴(kuò)展長度的影響顯著，地應(yīng)力差越大，水力壓裂裂縫的延伸長度越大。

圖17 不同地應(yīng)力組合下的裂縫擴(kuò)展形態(tài)Fig.17 Crack propagation considering different stress combinations

3 層理頁巖的水力壓裂裂縫擴(kuò)展特性

3.1 層理頁巖的水力壓裂模型構(gòu)建

建立層理頁巖數(shù)值模型，采用節(jié)理模型描述層理結(jié)構(gòu)，如圖18所示，模型長×寬為100 mm×100 mm，基于DFN建立光滑節(jié)理模型，設(shè)置n行n列均布節(jié)理組，注水孔處不設(shè)節(jié)理，節(jié)理長度均設(shè)為20 mm，各節(jié)理與水平方向的夾角為α。

圖18 層理頁巖的水力壓裂數(shù)值模型(黑色為DFN模型，綠色為光滑節(jié)理接觸)Fig.18 Numerical model of hydraulic fracturing for bedding shale (black is DFN model,green is smooth joint contact)

3.2 不同圍壓下的頁巖壓裂裂縫擴(kuò)展研究

設(shè)置3行3列節(jié)理組的層理頁巖模型，研究不同圍壓下的頁巖水力壓裂裂縫擴(kuò)展特性。圖19所示為各種圍壓下的頁巖起裂注水壓力。從圖19可知：同一圍壓下不同層理傾角頁巖的起裂注水壓力差異較小，主要在于數(shù)值計算中的節(jié)理為光滑節(jié)理模型，標(biāo)定后的節(jié)理剛度取值相對較大，當(dāng)水從中心點(diǎn)注入時，光滑節(jié)理的滑移變形相對較小，闡明層理傾角的控制較弱；此時控制層理頁巖破壞的主要因素為圍壓，則圍壓一定時，層理傾角的改變對起裂壓力的影響較小，起裂壓力主要隨著圍壓增大而增大，起裂壓力呈現(xiàn)左高右低的變化，說明起裂注水壓力與圍壓呈正相關(guān)。同時，計算同一圍壓下層理頁巖平均起裂注水壓力。圖20所示為平均起裂壓力-圍壓曲線。從圖20可知：曲線呈線性增長，斜率為2.090 75，與經(jīng)典水壓致裂破壞準(zhǔn)則對應(yīng)的起裂力Pini=σt+2σ3中的斜率2接近，再次驗證了基于改進(jìn)流固耦合算法的水力壓裂模型是可行的。

圖19 不同圍壓下的注水起裂壓力變化Fig.19 Variation of initiation pressure for shales with different angles of bedding planes considering different confining pressures

圖20 平均起裂壓力與圍壓的關(guān)系Fig.20 Curves of average initiation pressure with confining pressure

圖21所示為σx/σy=1.0 組合不同層理傾角的頁巖水力壓裂裂縫形態(tài)(紅色代表水力壓裂裂縫，黑色代表層理，灰色代表顆粒)。從圖21可知：當(dāng)?shù)蛧鷫簽? MPa且無主應(yīng)力差時，水力壓裂誘導(dǎo)裂縫易被層理捕獲，當(dāng)壓裂產(chǎn)生的裂隙遇到天然裂隙或?qū)永頃r會順結(jié)構(gòu)面產(chǎn)生開裂。當(dāng)圍壓大于或等于5 MPa 時，層理傾角為60°，75°和90°時，水力壓裂主裂縫遇到層理后偏離一小段距離，然后離開層理繼續(xù)朝原方向發(fā)展。對次裂縫而言，層理存在將影響其擴(kuò)展，層理角度為0°～30°時，次裂縫延伸方向與層理方向一致，次裂縫沿著層理發(fā)展。當(dāng)層理傾角為60°，75°和90°時，次裂縫順層理延伸一段后沿原方向延伸，并最終延伸到模型邊界，少數(shù)情況下，次裂縫以90°逼近層理時，直接穿過層理并以原方向延伸。從圖21還可知：在圍壓作用下不同層理傾角導(dǎo)致水力壓裂裂縫以不同的逼近角逼近層理，引起水力壓裂裂縫的擴(kuò)展方式不同，水力壓裂裂縫形態(tài)差異顯著，可見層理傾角對水壓裂縫擴(kuò)展過程及裂縫形成狀態(tài)起著主控作用。當(dāng)圍壓大于等于5 MPa 時，同一條件下層理傾角為60°，75°和90°的最終裂縫分布比較接近?？梢姴煌瑢永韮A角對水壓裂縫擴(kuò)展的影響在低圍壓時更加顯著，且水力裂縫的最終形態(tài)也受圍壓影響大。

圖21 σx/σy=1.0時的水壓裂縫演化圖Fig.21 Evolution of hydraulic fractures with σx/σy=1.0

3.3 不同注水壓力下的層理頁巖壓裂裂縫擴(kuò)展研究

在儲層頁巖壓裂過程中，不同注水壓力對裂縫擴(kuò)展產(chǎn)生顯著影響，考慮2種圍壓組合不同的中心注水壓力(圍壓為1 MPa和圍壓為20 MPa的一倍起裂壓力分別為9.357 MPa 和49.543 MPa)。圖22所示為各層理傾角頁巖的水力壓裂裂縫擴(kuò)展形態(tài)(紅色為壓裂裂縫，裂隙外，黑色圓盤為水壓，綠色為管道流速，黑色短線為層理)。從圖22可知：當(dāng)圍壓為20 MPa、組合注水壓力為49.543 MPa時，層理傾角為0°～30°的頁巖模型只生成2條主要裂縫，其他模型最終都形成了3條主要裂縫。當(dāng)圍壓為1 MPa、組合注水壓力為9.357 MPa 時，裂縫擴(kuò)展到邊界時，管道流量只達(dá)到了距離層理前一定距離的位置，并沒有跨過層理流到邊界外，闡明裂縫擴(kuò)展主要靠基質(zhì)應(yīng)力的傳遞，主要在于低圍壓下，頁巖強(qiáng)度較低，裂縫起裂擴(kuò)展的速度很快，裂縫擴(kuò)展相對裂縫水流傳導(dǎo)優(yōu)先。當(dāng)圍壓為20 MPa、組合注水壓力為49.543 MPa 時，管道流量也滯后于裂縫的擴(kuò)展，但裂縫系統(tǒng)中的管道流動方向大多靠近頁巖層理，裂縫未延伸到層理位置時，裂縫的擴(kuò)展主要依靠管域內(nèi)孔壓的增大，主要在于圍壓較大導(dǎo)致裂縫擴(kuò)展較慢，管道流量的傳遞緊跟裂縫擴(kuò)展，當(dāng)裂縫擴(kuò)展到層理位置時，由于基質(zhì)的應(yīng)力傳遞及層理強(qiáng)度較弱，裂縫通過基質(zhì)應(yīng)力傳遞并繼續(xù)延伸，此時的管道流體遇到開度未擴(kuò)大的基質(zhì)管道系統(tǒng)，流量傳遞受到阻礙。當(dāng)組合注水壓力為二倍起裂壓力時，管道流量的分布較裂縫擴(kuò)展速度滯后，此時中心孔注水壓力突增引起基質(zhì)裂縫快速萌生并迅速擴(kuò)展，隨后流體流入萌生擴(kuò)展的裂縫中。另外，由圖22還可知：裂縫起裂方向基本一致，但不同層理傾角頁巖的最終裂縫擴(kuò)展方向并不完全一致，說明不同注水壓力下層理傾角對頁巖壓裂裂縫擴(kuò)展形成具有控制作用。

針對不同圍壓及注水壓力下的層理頁巖壓裂裂縫擴(kuò)展機(jī)理進(jìn)行了細(xì)致的研究，其結(jié)果表明：當(dāng)圍壓一定時，層理傾角的主控作用較弱，此時的起裂壓力變化較小，當(dāng)層理傾角一定時，圍壓對起裂壓力的控制作用顯著，圍壓增大則起裂壓力大，闡明了增大的圍壓抑制了層理頁巖破壞，且不同層理傾角頁巖在不同圍壓下的起裂擴(kuò)展形態(tài)表明層理傾角對最終裂縫形態(tài)的主控作用；當(dāng)注水壓力增大時，中心孔處的起裂加快且裂縫迅速擴(kuò)展，不同層理傾角頁巖的擴(kuò)展路徑不同，進(jìn)一步闡明層理傾角對最終壓裂裂縫形態(tài)的顯著影響。可見，為了形成復(fù)雜的裂縫網(wǎng)絡(luò)以獲得更優(yōu)的采收率，考慮地應(yīng)力及注入模式的層理頁巖壓裂設(shè)計是關(guān)鍵。

3.4 不同層理密度頁巖的水力壓裂裂縫擴(kuò)展研究

選用層理密度為4行4列且長度保持為20 mm，建立層理頁巖水力壓裂數(shù)值模型，力學(xué)參數(shù)與3行3列模型一致。圖23所示為起裂注水壓力隨層理傾角的變化曲線。從圖23可知：高密度層理頁巖的起裂注水壓力明顯低于低密度層理頁巖的起裂注水壓力，主要在于層理存在弱化了頁巖強(qiáng)度，水力壓裂過程中，較大層理密度頁巖內(nèi)注的水壓力一旦達(dá)到起裂壓力，則頁巖破損形成裂縫，且層理密度較高，更快形成水力壓裂裂縫網(wǎng)絡(luò)。當(dāng)層理傾角為30°～60°區(qū)間時，較高層理密度頁巖的起裂注水壓力波動較大，特別是當(dāng)層理傾角為45°時，由于層理靠近注水孔附近，導(dǎo)水能力增強(qiáng)，此時注水壓力波動變大。

圖23 不同層理密度頁巖的起裂壓力Fig.23 Variation of initiation pressure for shales with different densities of bedding planes

圖24所示為含4行4列層理頁巖在σx/σy=1.0的壓裂裂縫形狀圖。從圖24可知：層理密度增大后，與低層理密度頁巖的最終壓裂裂縫網(wǎng)絡(luò)差異顯著，層理傾角為0°和90°時，水力壓裂裂縫完全被靠近注水孔附近的天然層理捕獲，導(dǎo)致中心位置的水力壓裂裂縫擴(kuò)展受限，使水力壓裂裂縫沿層理擴(kuò)展到邊界，再沿初始方向擴(kuò)展。當(dāng)圍壓大于等于5 MPa時，45°，60°和75°巖樣的水力壓裂裂縫由3條主裂縫變?yōu)? 條以中心為交點(diǎn)的主裂縫，2 條都沿層理擴(kuò)展一段后再沿初始方向擴(kuò)展，特別是圍壓為10 MPa 和20 MPa 時，頁巖層理傾角為45°時存在1 條主裂縫，以逼近角90°直接穿過裂縫。由此說明，在同一壓裂參數(shù)工況下，高密度層理影響的頁巖水力壓裂裂縫更易出現(xiàn)，且形成的裂縫網(wǎng)絡(luò)更復(fù)雜。

圖24 σx/σy=1.0水壓裂縫演化圖Fig.24 Evolution of hydraulic fractures with σx/σy=1.0

4 結(jié)論

1) 針對水力壓裂裂縫擴(kuò)展中的管道滲流與裂縫滲流特性，分析原顆粒流的流固耦合方法的不足，提出改進(jìn)的流固耦合算法，模擬荷載作用下的層理頁巖力學(xué)響應(yīng)及滲透特性，對比試驗結(jié)果，以驗證改進(jìn)算法的合理性。

2) 基于改進(jìn)的流固耦合模型，通過設(shè)定額外屬性區(qū)分管道流和裂隙流，提出水力壓裂過程中裂縫開度的計算公式。

3) 不同地應(yīng)力下的水力壓裂裂縫擴(kuò)展模式不同，但大致沿最大主應(yīng)力方向延伸，且兩方向上的地應(yīng)力差越大，效果越明顯，驗證水力壓裂模型的可行性。

4) 不同地應(yīng)力組合的壓裂裂縫擴(kuò)展模式不同，兩地應(yīng)力相差較大時，裂縫依然沿著最大主應(yīng)力方向延伸，當(dāng)層理傾角不同時，裂縫以較大角度逼近層理面，且較易穿過層理面，當(dāng)裂縫逼近角較小時，易被層理捕獲并延伸；當(dāng)注水壓力增大時，基質(zhì)裂縫迅速擴(kuò)展萌生并擴(kuò)展，流體流入萌生擴(kuò)展的裂縫中，裂縫延伸模式差異顯著；層理密度越小，起裂注水壓力越大，但層理密度較大時，對裂縫延伸的阻礙效應(yīng)更顯著。