非對(duì)稱(chēng)應(yīng)力作用下干熱巖水力壓裂裂縫擴(kuò)展機(jī)理研究

曹 涵， 朱珍德， 周露明

（1.河海大學(xué)巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室，南京 210098；2.河海大學(xué)江蘇省巖土工程技術(shù)工程研究中心，南京 210098）

近年來(lái)能源安全問(wèn)題已經(jīng)逐漸引起人們的重視［1］. 地?zé)崮苁且环N可再生性的天然清潔能源，因?yàn)槠渚哂械臀廴?、低排放的特性，所以?duì)于地?zé)豳Y源的開(kāi)發(fā)利用引起了人們的高度關(guān)注. 現(xiàn)階段，人類(lèi)能夠利用的地?zé)崮艽蠖鄡?chǔ)存于高溫高壓、低孔低滲的干熱巖中. 人們利用增強(qiáng)型地?zé)嵯到y(tǒng)（Enhanced Geothermal Systems，EGS）將這類(lèi)干熱巖通過(guò)人工手段改造為高滲透率的巖體，以此加熱低溫液體開(kāi)采地?zé)崮? 由于干熱巖高溫高壓、低孔低滲的特性，壓裂液主要在水力壓裂形成的裂縫網(wǎng)絡(luò)（簡(jiǎn)稱(chēng)縫網(wǎng)）中流通，而建立增強(qiáng)型地?zé)嵯到y(tǒng)就是為了人為形成裂縫網(wǎng)絡(luò)，以便更高效地提取干熱巖中的地?zé)崮埽?-5］.

增強(qiáng)型地?zé)嵯到y(tǒng)水壓致裂過(guò)程中的熱破裂是其所面臨的重要問(wèn)題，研究干熱巖開(kāi)發(fā)過(guò)程中的熱流固耦合損傷（Thermal-Hydraulic-Mechanical-Damage，THM-D）問(wèn)題有重要的科研價(jià)值. 李士斌等［6］分別分析了地層在拉伸破壞和剪切破壞下的裂縫擴(kuò)展規(guī)律，就不同地應(yīng)力形態(tài)對(duì)裂縫擴(kuò)展的影響進(jìn)行了初探. 袁志剛［7］通過(guò)研究實(shí)現(xiàn)了流固耦合損傷模型的有限元求解，同時(shí)對(duì)煤巖體水壓致裂及裂縫擴(kuò)展的影響因素進(jìn)行了分析. Homand-Etienne和Houpert［8］對(duì)高溫處理后的巖體試樣進(jìn)行了電鏡掃描分析，同時(shí)分析了巖體微觀結(jié)構(gòu)的損傷狀況對(duì)巖體物理力學(xué)性能的影響. 綜合來(lái)看，現(xiàn)有的大量研究采用的都是熱流固耦合損傷模型，且均為單井筒模型，同時(shí)關(guān)于熱應(yīng)力在水力壓裂裂縫擴(kuò)展中的作用的研究相對(duì)較少.

為研究非對(duì)稱(chēng)應(yīng)力作用下干熱巖雙井筒模型的水力壓裂裂縫擴(kuò)展機(jī)理，本研究依據(jù)線(xiàn)彈性力學(xué)、能量守恒定律、質(zhì)量守恒定律和彈-脆性損傷力學(xué)等理論建立干熱巖水壓致裂的熱流固耦合損傷模型，然后使用有限元方法對(duì)非對(duì)稱(chēng)地應(yīng)力和熱應(yīng)力影響下的雙井筒模型的水力壓裂裂縫擴(kuò)展機(jī)理進(jìn)行數(shù)值模擬研究.

1 熱流固耦合損傷模型

1.1 應(yīng)力場(chǎng)控制方程

根據(jù)彈性損傷本構(gòu)關(guān)系，考慮有效應(yīng)力、位移和熱應(yīng)力時(shí)線(xiàn)彈性體的靜力平衡方程［9-10］如公式（1）所示.

式中：G為巖體剪切模量，Pa；v為巖體的排水泊松比；ui和Fi(i=x,y,z)分別為位移和體積力在i方向的分量；αp,i為水壓力作用項(xiàng)；K′αTT,i為熱應(yīng)力項(xiàng)；K′為介質(zhì)的排水體積模量，K′=2G(1+v)/3(1-2v)，Pa；αT為巖體的熱膨脹系數(shù)，℃-1.

1.2 溫度場(chǎng)控制方程

假設(shè)巖體與孔隙流體之間處于熱平衡狀態(tài)，則考慮熱傳導(dǎo)作用、能量守恒和體積變化的溫度場(chǎng)控制方程［11］如公式（2）所示.

式中：λM=λs(1-φ)+λ1φ，λs和λ1分別為巖體和孔隙流體的熱傳導(dǎo)系數(shù)，W/（m·K）；(ρC)M為含流體的多孔介質(zhì)的熱容，(ρC)M=ρsCs(1-φ)+ρ1C1φ，kJ·m3/K；C1為孔隙流體的比熱系數(shù)，kJ/（kg·K）；Cs為巖體的比熱系數(shù)，kJ/（kg·K）；T0為零應(yīng)力狀態(tài)下的參考溫度，K；ρ1為孔隙流體的密度，kg/m3；k為巖體的滲透率，m2；t為時(shí)間，s.

1.3 滲流場(chǎng)控制方程

基于流體的質(zhì)量守恒方程以及狀態(tài)方程，假設(shè)巖體與孔隙流體之間處于熱平衡狀態(tài)，則考慮變形和溫度的滲流場(chǎng)控制方程如公式（3）所示.

公式（1）～（3）共同構(gòu)成了熱流固全耦合的控制方程組.

1.4 損傷分析

1.4.1 破壞準(zhǔn)則

當(dāng)巖體的應(yīng)力狀態(tài)達(dá)到最大拉應(yīng)力準(zhǔn)則(F1≥0)時(shí)，發(fā)生拉伸損傷；當(dāng)巖體的應(yīng)力狀態(tài)達(dá)到摩爾-庫(kù)倫準(zhǔn)則(F2≥0)時(shí)，發(fā)生剪切損傷，F(xiàn)1和F2的表達(dá)式分別如公式（4）和公式（5）所示［12-13］.

式中：F1和F2分別為兩個(gè)表示巖體應(yīng)力狀態(tài)的函數(shù)；fc和ft分別為巖體的單軸抗壓強(qiáng)度和單軸抗拉強(qiáng)度，Pa.

選用彈-脆性損傷模型對(duì)巖體進(jìn)行損傷判斷，當(dāng)巖體的應(yīng)力狀態(tài)滿(mǎn)足F2＜0 且F1≥0 時(shí)，巖體發(fā)生拉伸損傷，損傷變量D可用公式（6）表示［14］.

式中：λ為抗拉強(qiáng)度的殘余系數(shù)，λ=ftr/ft0；ft0為巖體的單軸抗拉強(qiáng)度，Pa；ftr為殘余強(qiáng)度，Pa；εt0為巖體發(fā)生拉伸損傷時(shí)的拉應(yīng)變；εtu為極限拉應(yīng)變，εtu=ηεt0；在復(fù)雜應(yīng)力條件下，ε=ε3，即為第三主應(yīng)變.

當(dāng)巖體的應(yīng)力狀態(tài)滿(mǎn)足F2≥0 且F1＜0 時(shí)，巖體發(fā)生剪切損傷，損傷變量D可用公式（7）表示.

式中：εc0為巖體發(fā)生剪切損傷時(shí)的壓應(yīng)變.

1.4.2 損傷對(duì)模型參數(shù)的影響

根據(jù)彈性損傷理論和應(yīng)變等效原理，巖體修正后的彈性模量［15-16］可用公式（8）表示.

式中：E0和E分別為巖體發(fā)生損傷前和發(fā)生損傷后的彈性模量，Pa.

巖體的應(yīng)力狀態(tài)影響著巖體的孔隙度，其關(guān)系式［17］可用公式（9）表示.

式中：φ0為巖體的初始孔隙度，干熱巖可取值為0.01；φ為巖體當(dāng)前的孔隙度；φr為高壓縮狀態(tài)下孔隙度的極限值；αφ為孔隙度的應(yīng)力敏感系數(shù)，αφ=5.0×10-8Pa-1；σˉv為平均有效應(yīng)力，其計(jì)算公式［17］如公式（10）所示.

式中：σ1、σ2、σ3分別為最大主應(yīng)力、中主應(yīng)力和最小主應(yīng)力，Pa；p為孔隙水壓力，Pa；α為Biot系數(shù).

當(dāng)巖石發(fā)生損傷后，其孔隙度和滲透性的演化規(guī)律較為復(fù)雜，可根據(jù)公式（11）描述其關(guān)系.

式中：αD為損傷對(duì)滲透率的影響系數(shù)，αD=5.0；k0為初始滲透率，m2；k為受到應(yīng)力作用后的滲透率，m2.

雖然目前研究學(xué)者還沒(méi)有在損傷變量對(duì)巖體熱物性參數(shù)的影響方面形成統(tǒng)一的認(rèn)識(shí)，且相關(guān)的理論和試驗(yàn)研究也比較少，但可以肯定的是，損傷必然會(huì)引起巖體熱傳導(dǎo)系數(shù)的升高. 本研究假定損傷是以指數(shù)關(guān)系的形式影響巖體熱傳導(dǎo)系數(shù)的［11］，其關(guān)系式如公式（12）所示.

式中：αλ為損傷對(duì)巖體熱傳導(dǎo)系數(shù)的影響系數(shù)；λs(T,D)為損傷發(fā)生后巖體的熱傳導(dǎo)系數(shù)；λs(T)為損傷發(fā)生前巖體的熱傳導(dǎo)系數(shù).

2 數(shù)值模擬程序

由于熱流固三場(chǎng)和損傷之間的耦合作用非常復(fù)雜，且考慮到熱流固三場(chǎng)耦合控制方程組的高度非線(xiàn)性以及巖體的多個(gè)物性參數(shù)都處于動(dòng)態(tài)的變化過(guò)程中，因此一般會(huì)使用有限元方法對(duì)熱流固三場(chǎng)耦合控制方程組進(jìn)行求解，其中較常用的有限元數(shù)值模擬軟件為COMSOL，其易于進(jìn)行耦合方程的建立，且具有較高的計(jì)算精度，能高效快速地求解非線(xiàn)性熱流固三場(chǎng)耦合控制方程組.

本研究首先采用COMSOL軟件建立并求解熱流固三場(chǎng)耦合控制方程組，然后根據(jù)損傷對(duì)巖體物性參數(shù)的影響對(duì)各參數(shù)進(jìn)行修正，之后再對(duì)熱流固三場(chǎng)耦合控制方程組進(jìn)行求解，如此重復(fù)，直到不再發(fā)生損傷為止.

3 非對(duì)稱(chēng)應(yīng)力下水力壓裂的數(shù)值模擬研究

3.1 數(shù)值模型的建立

為了研究非對(duì)稱(chēng)應(yīng)力作用下干熱巖雙井筒模型的水力壓裂裂縫擴(kuò)展機(jī)理，建立了干熱巖雙井筒數(shù)值模型，如圖1 所示. 模型尺寸為1 m×1 m 的正方形，雙井筒的直徑均為0.1 m，且為非中心對(duì)稱(chēng). 模型的左邊界和下邊界固定位移，右邊界和上邊界施加水平地應(yīng)力，雙井筒筒壁施加水壓力. 考慮到巖體的非均質(zhì)性，可假定巖體的彈性模量服從Weibull 分布［18-20］. 因?yàn)楸狙芯坎捎玫碾p井筒模型為非中心對(duì)稱(chēng)的，所以地應(yīng)力的分布情況將極大地影響裂縫擴(kuò)展形態(tài).故本研究分別在兩種工況條件下對(duì)干熱巖雙井筒模型的水力壓裂過(guò)程進(jìn)行模擬，以得出非對(duì)稱(chēng)應(yīng)力作用下干熱巖雙井筒模型的水力壓裂裂縫擴(kuò)展機(jī)理. 兩種工況分別是：最大地應(yīng)力方向垂直于雙井筒連線(xiàn)方向（記為工況1）；最大地應(yīng)力方向平行于雙井筒連線(xiàn)方向（記為工況2）. 在進(jìn)行數(shù)值模擬計(jì)算時(shí)，均固定初始地應(yīng)力大小，為確定不同情況下水力壓裂的起裂壓力，以每步0.5 MPa 的速度增大孔壁的注水壓力，直至模型發(fā)生破壞. 表1為數(shù)值模擬中采用的相關(guān)參數(shù).

圖1 干熱巖雙井筒數(shù)值模型示意圖Fig.1 Schematic diagram of numerical model of hot dry rock double wellbore

表1 數(shù)值模擬中采用的相關(guān)參數(shù)Tab.1 Relevant parameters used in numerical simulation

3.2 數(shù)值模擬的結(jié)果分析

溫度差是引起熱應(yīng)力的最主要因素，為了模擬熱應(yīng)力對(duì)干熱巖雙井筒模型水力壓裂裂縫擴(kuò)展的影響，在工況1的條件下，分別模擬了基巖溫度為100、150、200 ℃時(shí)干熱巖雙井筒模型的水力壓裂過(guò)程，其最終的損傷區(qū)域分布如圖2所示.D=0表示巖體未發(fā)生損傷，0＜D≤1表示巖體發(fā)生了拉伸損傷，-1≤D＜0表示巖體發(fā)生了剪切損傷.

研究表明，當(dāng)基巖溫度為100、150、200 ℃時(shí)，巖體的破裂壓力分別為35、32、30 MPa，說(shuō)明隨著巖體和注入壓裂液之間溫差的上升，巖體的破裂壓力逐步降低. 由圖2可知，工況1下，基巖溫度為100 ℃時(shí)，裂縫首先沿著最大地應(yīng)力方向發(fā)展，隨后低溫壓裂液由裂縫進(jìn)入巖體內(nèi)部，低溫壓裂液與高溫的裂縫壁面接觸，在兩者的溫差作用下產(chǎn)生了熱應(yīng)力，熱應(yīng)力的存在加劇了后續(xù)的裂縫擴(kuò)展，最終形成裂縫網(wǎng)絡(luò)；基巖溫度為150 ℃和200 ℃時(shí)，巖體在裂縫壁面產(chǎn)生了更多的分支裂縫，這使得整個(gè)裂縫網(wǎng)絡(luò)更為復(fù)雜. 從圖2還可以看出，當(dāng)巖體和注入壓裂液之間的溫差升高到一定程度時(shí)，兩個(gè)井筒壓裂產(chǎn)生的裂縫網(wǎng)絡(luò)會(huì)在損傷發(fā)展后期發(fā)生交匯，并在縫網(wǎng)交匯處迅速產(chǎn)生新的裂縫，這說(shuō)明溫度更高的巖體會(huì)產(chǎn)生范圍更大、更廣的裂縫網(wǎng)絡(luò). 綜上可知，當(dāng)最大地應(yīng)力方向垂直于雙井筒連線(xiàn)方向時(shí)，隨著基巖溫度的增加，巖體和注入壓裂液之間的溫差會(huì)逐漸增加，從而使熱應(yīng)力也隨之增加，而熱應(yīng)力越大，巖體越容易發(fā)生損傷.

圖2 工況1基巖溫度對(duì)干熱巖水力壓裂損傷區(qū)域分布的影響Fig.2 Influence of different bedrock temperatures on the distribution of hydraulic fracturing damage area of hot dry rock under working condition 1

為了模擬注水壓力對(duì)干熱巖雙井筒模型水力壓裂裂縫擴(kuò)展的影響，在工況2的條件下，分別模擬了注水壓力為3.5、20、43 MPa（對(duì)應(yīng)的迭代步數(shù)分別為7、40、86步）時(shí)干熱巖雙井筒模型的水壓致裂過(guò)程，其最終的損傷區(qū)域分布如圖3所示. 由圖3可知，工況2下，隨著注水壓力的增大，巖體首先沿著最大地應(yīng)力方向發(fā)生破壞，隨著裂縫的不斷發(fā)展，雙井筒內(nèi)側(cè)區(qū)域中產(chǎn)生的裂縫發(fā)生交匯，裂縫擴(kuò)展后期，雙井筒外側(cè)的裂縫網(wǎng)絡(luò)逐漸趨于穩(wěn)定，但雙井筒內(nèi)側(cè)縫網(wǎng)交匯區(qū)的裂縫繼續(xù)發(fā)展，最終導(dǎo)致巖體在雙井筒內(nèi)側(cè)形成的損傷區(qū)域明顯大于其在雙井筒外側(cè)形成的損傷區(qū)域，由此可見(jiàn)，雙井筒內(nèi)側(cè)裂縫網(wǎng)絡(luò)的交匯對(duì)后續(xù)裂縫的發(fā)展有明顯的促進(jìn)作用.

圖3 工況2注水壓力對(duì)干熱巖水力壓裂損傷區(qū)域分布的影響Fig.3 Influence of different water injection pressures on the distribution of hydraulic fracturing damage area of hot dry rock under working condition 2

4 結(jié)論

本研究建立了描述巖體細(xì)觀結(jié)構(gòu)的熱流固耦合損傷模型，并對(duì)非對(duì)稱(chēng)應(yīng)力和熱應(yīng)力作用下干熱巖雙井筒模型的水力壓裂裂縫擴(kuò)展過(guò)程進(jìn)行了數(shù)值模擬研究. 得到如下結(jié)論：

1）當(dāng)最大地應(yīng)力的方向垂直于雙井筒連線(xiàn)方向且基巖溫度較低時(shí)，裂縫首先沿最大地應(yīng)力方向發(fā)展，且此時(shí)兩個(gè)井筒周?chē)a(chǎn)生的裂縫網(wǎng)絡(luò)互不干擾；隨著基巖溫度的升高，巖體破裂所需的注入壓力降低，巖體的損傷區(qū)域面積則隨之增大，當(dāng)基巖溫度升高到一定程度時(shí)，兩個(gè)井筒周?chē)目p網(wǎng)發(fā)生交匯，在交匯處裂縫會(huì)繼續(xù)擴(kuò)展，且裂縫在交匯處的擴(kuò)展速度快于其在無(wú)交匯損傷區(qū)域的擴(kuò)展速度.

2）當(dāng)最大地應(yīng)力的方向平行于雙井筒連線(xiàn)方向且注水壓力較低時(shí)，裂縫同樣首先沿著最大地應(yīng)力方向發(fā)展，隨著注水壓力的增大，裂縫會(huì)在雙井筒內(nèi)側(cè)區(qū)域發(fā)生交匯，內(nèi)側(cè)交匯區(qū)的裂縫擴(kuò)展范圍明顯大于雙井筒外側(cè).

3）地應(yīng)力場(chǎng)對(duì)雙井筒模型的圍巖破裂形態(tài)影響顯著，且最大地應(yīng)力的方向直接決定了裂縫發(fā)展的方向和大致范圍.