非對稱應力作用下干熱巖水力壓裂裂縫擴展機理研究

曹 涵， 朱珍德， 周露明

（1.河海大學巖土力學與堤壩工程教育部重點實驗室，南京 210098；2.河海大學江蘇省巖土工程技術工程研究中心，南京 210098）

近年來能源安全問題已經逐漸引起人們的重視［1］. 地熱能是一種可再生性的天然清潔能源，因為其具有低污染、低排放的特性，所以對于地熱資源的開發利用引起了人們的高度關注. 現階段，人類能夠利用的地熱能大多儲存于高溫高壓、低孔低滲的干熱巖中. 人們利用增強型地熱系統（Enhanced Geothermal Systems，EGS）將這類干熱巖通過人工手段改造為高滲透率的巖體，以此加熱低溫液體開采地熱能. 由于干熱巖高溫高壓、低孔低滲的特性，壓裂液主要在水力壓裂形成的裂縫網絡（簡稱縫網）中流通，而建立增強型地熱系統就是為了人為形成裂縫網絡，以便更高效地提取干熱巖中的地熱能［2-5］.

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

為研究非對稱應力作用下干熱巖雙井筒模型的水力壓裂裂縫擴展機理，本研究依據線彈性力學、能量守恒定律、質量守恒定律和彈-脆性損傷力學等理論建立干熱巖水壓致裂的熱流固耦合損傷模型，然后使用有限元方法對非對稱地應力和熱應力影響下的雙井筒模型的水力壓裂裂縫擴展機理進行數值模擬研究.

1 熱流固耦合損傷模型

1.1 應力場控制方程

根據彈性損傷本構關系，考慮有效應力、位移和熱應力時線彈性體的靜力平衡方程［9-10］如公式（1）所示.

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

1.2 溫度場控制方程

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

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

1.3 滲流場控制方程

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

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

1.4 損傷分析

1.4.1 破壞準則

當巖體的應力狀態達到最大拉應力準則(F1≥0)時，發生拉伸損傷；當巖體的應力狀態達到摩爾-庫倫準則(F2≥0)時，發生剪切損傷，F1和F2的表達式分別如公式（4）和公式（5）所示［12-13］.

式中：F1和F2分別為兩個表示巖體應力狀態的函數；fc和ft分別為巖體的單軸抗壓強度和單軸抗拉強度，Pa.

選用彈-脆性損傷模型對巖體進行損傷判斷，當巖體的應力狀態滿足F2＜0 且F1≥0 時，巖體發生拉伸損傷，損傷變量D可用公式（6）表示［14］.

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

當巖體的應力狀態滿足F2≥0 且F1＜0 時，巖體發生剪切損傷，損傷變量D可用公式（7）表示.

式中：εc0為巖體發生剪切損傷時的壓應變.

1.4.2 損傷對模型參數的影響

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

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

巖體的應力狀態影響著巖體的孔隙度，其關系式［17］可用公式（9）表示.

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

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

當巖石發生損傷后，其孔隙度和滲透性的演化規律較為復雜，可根據公式（11）描述其關系.

式中：αD為損傷對滲透率的影響系數，αD=5.0；k0為初始滲透率，m2；k為受到應力作用后的滲透率，m2.

雖然目前研究學者還沒有在損傷變量對巖體熱物性參數的影響方面形成統一的認識，且相關的理論和試驗研究也比較少，但可以肯定的是，損傷必然會引起巖體熱傳導系數的升高. 本研究假定損傷是以指數關系的形式影響巖體熱傳導系數的［11］，其關系式如公式（12）所示.

式中：αλ為損傷對巖體熱傳導系數的影響系數；λs(T,D)為損傷發生后巖體的熱傳導系數；λs(T)為損傷發生前巖體的熱傳導系數.

2 數值模擬程序

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

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

3 非對稱應力下水力壓裂的數值模擬研究

3.1 數值模型的建立

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

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

表1 數值模擬中采用的相關參數Tab.1 Relevant parameters used in numerical simulation

3.2 數值模擬的結果分析

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

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

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

為了模擬注水壓力對干熱巖雙井筒模型水力壓裂裂縫擴展的影響，在工況2的條件下，分別模擬了注水壓力為3.5、20、43 MPa（對應的迭代步數分別為7、40、86步）時干熱巖雙井筒模型的水壓致裂過程，其最終的損傷區域分布如圖3所示. 由圖3可知，工況2下，隨著注水壓力的增大，巖體首先沿著最大地應力方向發生破壞，隨著裂縫的不斷發展，雙井筒內側區域中產生的裂縫發生交匯，裂縫擴展后期，雙井筒外側的裂縫網絡逐漸趨于穩定，但雙井筒內側縫網交匯區的裂縫繼續發展，最終導致巖體在雙井筒內側形成的損傷區域明顯大于其在雙井筒外側形成的損傷區域，由此可見，雙井筒內側裂縫網絡的交匯對后續裂縫的發展有明顯的促進作用.

圖3 工況2注水壓力對干熱巖水力壓裂損傷區域分布的影響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 結論

本研究建立了描述巖體細觀結構的熱流固耦合損傷模型，并對非對稱應力和熱應力作用下干熱巖雙井筒模型的水力壓裂裂縫擴展過程進行了數值模擬研究. 得到如下結論：

1）當最大地應力的方向垂直于雙井筒連線方向且基巖溫度較低時，裂縫首先沿最大地應力方向發展，且此時兩個井筒周圍產生的裂縫網絡互不干擾；隨著基巖溫度的升高，巖體破裂所需的注入壓力降低，巖體的損傷區域面積則隨之增大，當基巖溫度升高到一定程度時，兩個井筒周圍的縫網發生交匯，在交匯處裂縫會繼續擴展，且裂縫在交匯處的擴展速度快于其在無交匯損傷區域的擴展速度.

2）當最大地應力的方向平行于雙井筒連線方向且注水壓力較低時，裂縫同樣首先沿著最大地應力方向發展，隨著注水壓力的增大，裂縫會在雙井筒內側區域發生交匯，內側交匯區的裂縫擴展范圍明顯大于雙井筒外側.

3）地應力場對雙井筒模型的圍巖破裂形態影響顯著，且最大地應力的方向直接決定了裂縫發展的方向和大致范圍.