裂縫分布對CO2-EGS熱提取的影響研究

周 慶，廖建興，徐 斌，姚 亮，周浩宇，趙 權

(貴州大學 土木工程學院，貴州 貴陽 550025)

與傳統化石能源相比，地熱能具有低碳性、可持續性、儲量豐富、耐高溫等特點[1]，被認為是一種極具潛力的清潔能源。大部分地熱能以干熱巖(Hot Dry Rock，HDR)的形式儲存在地下3～10 km深度的巖層中，其溫度通常為150～650℃[2]。增強型地熱系統(Enhanced Geothermal System，EGS)是通過注入流體與干熱巖之間發生熱交換來實現地熱提取的，被認為是目前從干熱巖中提取熱量最有效的方法[3-4]。D.W.Brown[5]首次提出了以CO2作為注入流體的新型增強型地熱系統(CO2-EGS)。K.Pruess[6]比較了水和CO2的熱物理性質，結果表明，與水相比，CO2具有更大的壓縮性與膨脹性，且CO2的低黏度與低密度更有利于傳熱。Zhong Chenghao等[7]定量分析了水與CO2作為工作流體時EGS的采熱性能，表明以CO2作為工作流體的熱提取率明顯高于水。為減少碳排放，助推我國“雙碳”目標達成，CO2-EGS的發展與應用引起了廣泛關注[8-9]。

地殼巖石中普遍存在著各種復雜裂隙，如節理、斷層等[10]。上述裂隙在不同地質作用下形成復雜的裂縫網絡[11]，裂縫網絡的空間分布通常是不均勻的[12]。P.Davy等[13]首次提出采用冪律分布函數描述裂縫網絡空間分布。C.Darcel等[14]研究了分形裂縫網絡的空間分布，表明裂縫冪律分布函數可以擴展到其他分形裂縫網絡分布。地熱儲層中裂縫網絡表現出明顯的分形特征，孫致學等[15]采用冪律分布函數刻畫裂縫的分形特征。而裂縫空間分布直接影響著裂縫網絡連通性[16]，進而影響流體流動。J.R.De Dreuzy等[17]研究了不同長度指數時裂縫中的流體分布情況，表明長度指數在1～2時，長度指數越小，流體分布越均勻。M.Gudala等[18]研究了裂縫數量分別為3、5和7條的情況下采熱性能的變化，并以儲層壽命、儲層流動阻抗和熱功率3種評價指標評估儲層熱性能，表明裂縫數量對地熱儲層有顯著影響。單丹丹[19]、王丹丹[20]等總結了裂縫數量以及裂縫長度的變化對系統采熱的影響，并以生產溫度、生產質量流率、采熱效率和熱提取率4項評價指標評估儲層熱性能，表明隨著裂縫數量的增加，儲層熱提取率越大。綜上分析，以往研究大多集中在分析具體裂縫數量或裂縫長度的變化對EGS采熱性能的影響，忽略了裂縫分布的隨機性與裂縫變形對流體通道的影響。

基于此，筆者以冪律分布的裂縫網絡為基礎，采用數值方法系統研究不同裂縫長度指數(a)、密度(β)的裂縫網絡對CO2-EGS采熱性能的影響，并以熱突破時間、EGS壽命、產熱率與總產熱能以及產熱效率5種評價指標對儲層熱性能進行詳細的評估，以期為干熱巖開發提供參考。

1 研究方法

1.1 模型構建

為探究裂縫空間分布對EGS采熱性能的影響，本文利用數值模擬軟件生成尺寸長(x)×寬(y)×高(z) 為 400 m×400 m×100 m的模型(圖1)，其底面深度為2 000 m，頂層深度為1 900 m。在模型中嵌入兩簇裂縫，其夾角(銳角)為60°。熱儲層采用一注一采的開采模式，其中，注入井(Inj.)和生產井(Pro.)的位置分別設置在x=100 m、y=100 m與x=300 m、y=300 m處，兩井間距約為282.84 m。

圖1 裂縫儲層模型的形狀與尺寸Fig.1 Geometry and dimensions of the fractured reservoir model

該模型中隨機生成裂縫的空間分布滿足冪律分布，其函數[10]可以表示為：

本研究中，D=2，對應的a取值范圍為[1，∞][10]。圖1所示模型中裂縫網絡是根據冪律分布函數所生成的，其中a=1.0，β=0.1。

圖2比較了上述模型裂縫長度模擬統計值與式(1)中裂縫長度理論統計值，顯示模擬數量與理論數量統計值基本相符，說明該模型的裂縫空間分布符合冪律函數分布。

圖2 不同裂縫長度下裂縫數量模擬結果與理論值的比較Fig.2 Comparison of the simulated and theoretical results of fracture numbers under different fracture lengths

1.2 初始與邊界條件

在該模型中，邊界不透水，且與外部沒有熱交換。模型z方向的初始壓應力為50 MPa，在水平方向上施加40 MPa的各向同性應力。

該模型采用經典的雙井布置模式，注入井深度z1=2 000 m，生產井深度z2=1 950 m，計算時間設置為20 a。計算模型的其他參數見表1。

表1 計算模型相關參數Table 1 Parameters of the fractured reservoir model

1.3 控制方程與耦合模型

本文中裂縫等效為連續多孔介質。多孔介質中流體的連續方程與質量守恒方程[22]分別為：

多孔介質中的溫度控制方程[22]為：

其中，qT=λeff?T，λeff=(1?φ)λs+φλf。

多孔介質變形的平衡方程與幾何控制方程[22]分別為：

在產熱過程中，裂縫寬度的變化取決于應力狀態與裂隙壓力。接觸狀態下裂縫寬度變化的函數[23]為：

裂縫寬度增量主要由純拉壓應力作用與剪切滑移引起的變形組成[23]。在純拉壓應力作用下，裂縫寬度增量 Δwm[23]為：

滑移引起的裂縫寬度增量 Δwd[23]為：

因此，在時間步為t+1時的裂縫寬度[23]為：

對于裂縫單元之間的流動，裂縫之間的滲透率可由以下函數[23]近似表示：

CO2的熱力學特征參數(包括焓、黏度和密度)是根據模擬軟件TOUGH2MP自帶的公式計算。若要了解詳細過程，可參考文獻[21]。

地熱儲層中的流體流動涉及復雜的熱?水?力學(THM)耦合過程。本研究中，THM耦合模型是通過耦合FLAC3D和TOUGH2MP兩個模擬軟件實現的。軟件FLAC3D解決研究過程中的力學問題，流體與熱流問題的求解則由軟件TOUGH2MP解決[22]。

1.4 模型驗證

現實中EGS儲層由復雜的裂縫網絡組成。單裂縫是裂縫網絡的基本組成部分，且單裂縫模型可通過巖石基質與裂縫之間的熱傳導充分研究裂縫熱提取的基本特征。故可建立一個模擬流體在裂縫巖石中流動的單裂縫模型來驗證裂縫網絡中的傳熱行為[1]。如圖3所示，在模型中嵌入一個寬度恒為1 mm的裂縫，x=0與入口重合。模型重要相關參數見表2。單裂縫模型內傳熱的解析解[23]如下：

表2 裂縫傳熱相關參數Table 2 Parameters related to heat transfer in fractures

圖3 2D單裂縫模型Fig.3 2D model of a single fracture

圖4a顯示了10、50和100 d后沿裂縫的空間溫度分布。并將模擬的x=10、20和50 m不同位置的溫度變化與理論值進行了比較(圖4b)。結果顯示模擬結果與解析結果吻合良好，驗證了該模型在裂隙巖中熱傳導和對流的可行性。

圖4 不同時間與不同位置的理論值和模擬值的比較Fig.4 Comparison of theoretical and simulated temperatures at different times and different locations

1.5 評價指標

連續穩定的高提取溫度以及大量的提取能量是地熱提取的必備條件[21]。因此，產熱相關因素是評估EGS性能的關鍵指標。

(1) 熱突破時間(Γb)：從運行開始至生產溫度較初始溫度降低1℃所用的總時間。

(2) EGS壽命(Γl)：從運行開始到生產溫度降至初始溫度的80%以下的總時間。EGS的開采周期受提取溫度和地熱系統壽命的影響[24]。

(3) 總產熱能W與產熱率Wh：總產熱能W[25]是指從運行開始到生產壽命期間產生的總能量，J；產熱率Wh為產熱量和注入熱量之間的差額，J；計算公式分別為：

(4) 產熱效率ηh：內部能量消耗(Wp)包括注入井消耗和生產井消耗。產熱效率可被定義為產熱率與內部能量消耗的比率[25]。內部能量消耗(Wp)和產熱效率(ηh)計算公式分別為：

2 模擬方案設置

裂縫冪律分布變化規律如圖5所示。圖中橫向表示密度β相同，裂縫長度指數a增大時裂縫變化情況。縱向為a相同，β增大時裂縫變化情況。紅色裂縫是變化后新增加的裂縫。結果顯示，當a增大，裂縫長度減小，長裂縫占比減小。而密度β增加，裂縫數量會增加，但是裂縫長短比例不會變化。

圖5 裂縫分布規律Fig.5 Distribution patterns of fractures

參考Lei Qinghua等[10]對裂縫冪律分布的研究，本文將裂縫分布情況設置了a=1.00、1.25、1.50、1.75、2.00與密度項β=0.100、0.125、0.150、0.175、0.200、0.225、0.250，正交組合共35個案例。在設置的35個案例中，密度項β=0.10、0.15、0.20、0.25的20組裂縫分布情況可清晰顯示裂縫分布隨裂縫長度指數a與密度β變化的規律(圖6)。結果顯示，隨著長度指數a的增大，裂縫長度減小，長裂縫占比減小，注采井之間形成的貫穿裂縫數量減少。隨著密度β的增大，裂縫數量增加，貫穿裂縫數量增加，裂縫連通性加強。

圖6 不同a與β組合生成的裂縫網絡Fig.6 Fracture networks generated under combinations of different a and β values

3 結果與分析

3.1 溫度分布

不同a和β條件下生產溫度(Tpro)在20 a生產過程中的演變如圖7所示。初始階段，生產溫度基本保持不變。隨著生產時間推移，注入的冷流沿裂縫逐漸流向生產井，在裂縫中與干熱巖接觸發生熱交換，從而導致生產溫度逐漸下降。a相同時，隨著密度β增大，裂縫數量增加，對流換熱面積變大，生產溫度下降越晚，降幅越小。隨著長度指數a增大，裂縫長度減小，貫穿裂縫數量減少，換熱面積減小，流體更為集中流向生產井，生產溫度下降時間越早，降幅越大。隨著長度指數a增大，運行20 a后不同密度之間的最大溫差可從a=1.00的10℃上升到a=2.00的46℃。說明a越大，裂縫密度對溫度分布的影響越明顯。

圖7 不同a和β條件下生產溫度隨時間的演變Fig.7 Time-varying production temperature under different a and β values

圖8顯示了運行20 a后不同案例的基質溫度分布情況。顯然，溫度在儲層中分布情況的變化規律與裂縫分布情況有關，而其隨密度項變化的規律可由β分別為0.10、0.15、0.20、0.25的20個案例展示。在連續注入條件下，低溫區域逐漸從注入井傳播到生產井，主要沿著連接良好的裂縫傳至生產井。與Sun Zhixue等[26]的研究結果相同，由于溫差驅動裂隙流體與圍巖之間發生對流換熱，熱量首先在靠近裂縫的區域擴散。故在裂縫連接不佳的采熱區域的溫度高于主裂縫區域。

圖8 不同a與β條件下儲層溫度分布(20 a)Fig.8 Distributions of reservoir temperature (20 a) under different a and β values

結果顯示，低溫區域隨著長度指數a增大向生產井延伸越明顯。而當a相同時，密度β越大，流體向注入井附近擴散越明顯，注入井周圍低溫區域越大。結合裂縫網絡分布圖(圖6)，與a=1.75相比，a=1.50的有效裂縫占比較大，注采井之間連通性更好，使其溫度下降更早，降幅更大。同時也使a=1.50時低溫區域向生產井延伸較a=1.75時更明顯(圖8)。

3.2 注入壓力(pinj)

圖9展示了注入井中不同a和β條件下注入壓力隨時間的演變。注入壓力是流體循環的主要驅動力，初始注入壓力為26 MPa。開始階段，注入壓力明顯上升，達到最大值后由于溫度的下降引起熱效應，加大裂縫寬度，使壓力緩慢下降。所以注入壓力的下降與溫度變化有關。

圖9 不同a和β條件下注入壓力隨時間的變化Fig.9 Time-varying injection pressure under different a and β values

結果顯示，當a=2.00時，β=0.250的最大注入壓力為28.9 MPa，高于β=0.225的28.1 MPa，結合裂縫分布情況(圖6)，β=0.250的裂縫連通性較差，未形成明顯的貫穿裂縫，流體流向生產井緩慢，使其注入壓力增加越多。其余情況下，a相同時，密度β越大，注入壓力上升越少。a=1.00時，β=0.100的最大注入壓力為28.5 MPa，而β=0.250的最大注入壓力為27.1 MPa，降低了約5%。隨著a增大，最大注入壓力也隨之增大。當β=0.100時，最大注入壓力由a=1.00的28.5 MPa上升至a=2.00時的32.8 MPa，上升了約15%。a越大，由于生產溫度的降幅變大，注入壓力的降幅也越大。在本文研究案例中觀察到的最大注入壓力32.8 MPa，這明顯高于Liao Jianxing等[23]觀察到的注入壓力，這是因為本研究采用了高注入速率。

3.3 裂縫寬度

裂縫作為注采井之間的主要通道，裂縫寬度的變化決定通道的大小，進而影響流體流動。圖10顯示了20 a時不同冪律分布的裂縫寬度情況，裂縫寬度隨裂縫密度的變化規律可由β分別為0.10、0.15、0.20、0.25的20個案例清晰展示。顯然，在注采井之間，裂縫寬度隨著與注入井距離的增加而減小。密度β相同時，a越大，大裂縫寬度向生產井延伸越明顯。而a相同時，隨著密度β的增大，裂縫寬度越小。該現象是因為裂縫寬度變化與基質溫度下降有關。由于溫度的下降，產生了熱效應，基質發生冷縮，從而使裂縫寬度增大。結合圖7與圖9比較，雖然壓力降低會使裂縫寬度變小，但溫降引起的熱應力對裂縫寬度變化的影響更顯著。溫度下降越多，產生的熱效應越明顯，使裂縫寬度增加越多。這與Zhang Xu等[27]的研究結果一致。

圖10 不同a和β條件下的裂縫寬度分布(20 a)Fig.10 Distributions of fracture width (20 a) under different a and β values

3.4 CO2與熱流動

圖11顯示了運行20 a時不同a和β條件下流體流線分布情況。顯然，流線的分布與裂縫分布有關，且流線分布隨裂縫分布變化的規律由密度項為0.10、0.15、0.20、0.25的20個案例可清晰展示。較小的裂縫密度β對應的裂縫連通性更好，為流體流動創造了主要的流經通道，并在其中觀察到明顯的竄流。且密度β越小，竄流越明顯。而裂縫密度β增加，流體流線分布越均勻。結合圖6，隨著長度指數a增大，注采井之間形成貫穿裂縫的數量越少，流體流動越集中，竄流越明顯。凸顯了裂縫分布對流體流動的影響。

圖11 不同a和β條件下CO2流體流線分布(20 a)Fig.11 Distributions of CO2 fluid streamline (20 a) under different a and β values

運行20 a時不同a和β條件下熱流的分布情況如圖12所示。熱傳導是由于溫差引起的，且對流換熱伴隨著流體流動。結果顯示，在裂縫中存在明顯的熱對流現象，a相同時，β越大，熱量流線分布越密集，向生產井延伸越緩慢；而隨著a增大，熱量流線向生產井延伸越明顯；且其規律與運行20 a時儲層溫度分布(圖8)基本一致。在熱采區邊界處，由于注入CO2與干熱巖之間的溫差顯著，熱流線由高溫區向低溫區流動。在其他區域，熱流與流體流線基本一致，其中熱對流占主導地位。這種現象與Yu Guojun等[28]的研究結果一致。因此，裂縫是CO2流體流動的主要通道。熱量首先從干熱巖中通過熱傳導傳遞給裂縫中的冷CO2，然后由CO2流以對流換熱的形式攜帶到生產井。

圖12 不同a和β條件下熱量流線分布(20 a)Fig.12 Distributions of heat-flow streamline (20 a) under different a and β values

3.5 性能評估

3.5.1 熱突破時間(Γb)與EGS壽命(Γl)

圖13a顯示了不同方案的熱突破時間。結果表明，在a相同時，密度β越大，熱突破時間越長。a=1.00時，熱突破時間由3.65 a上升至16.06 a，增加了3.4倍。a=2.00時，熱突破時間由0.10 a延長至6.99 a，增加了近70倍。密度β相同時，a越大，熱突破時間越短。β=0.100時，熱突破時間從3.65 a縮短到0.10 a，下降了3.5 a。β=0.250時，熱突破時間從16.06 a縮短到6.99 a，下降了約9 a。說明a越大，裂縫密度對熱突破時間的影響越顯著。結合圖6與圖10可知，當a=2.00、β=0.100時，貫穿裂縫數量少，裂縫寬度較大，且換熱面積較小，致使儲層溫度下降迅速，從而使其很快達到熱突破。

圖13 不同a和 β條件下的熱突破時間與儲層運行壽命Fig.13 Thermal breakthrough time and EGS life under different a and β values

EGS壽命如圖13b所示。裂縫分布不同時，流體對流換熱面積不同，達西阻力不同，有效換熱面積不同，從而導致儲層壽命不同。由此可見，裂縫分布對儲層壽命有一定影響。由于本文設置最大開采時間為20 a，即儲層最大壽命為20 a。結果表明，a=1.00和1.25時，在20 a生產時間內，最終生產溫度未下降到初始溫度的80%以下，故其壽命為20 a。其余情況下，密度β越大，EGS壽命越長。由于裂縫數量增加，對流換熱面積增加，生產溫度降幅減小，從而延長EGS壽命。當β=0.100時，EGS壽命從20 a縮短至9.87 a，下降了50.65%。顯然，當β相同時，隨著a增大，貫穿裂縫數量減少，換熱面積減小，生產溫度下降加快，使EGS壽命縮短。結合裂縫分布(圖6)與溫度分布(圖8)情況，a=1.50時的裂縫連通性較a=1.75時更好，且低溫區域向生產井延伸更明顯，溫度下降更多，所以a=1.50時熱突破時間更短，EGS壽命也更短。

3.5.2 產熱率Wh與總產熱能W

圖14顯示了不同案例在20 a生產過程中產熱率隨時間的變化。結果顯示，初始階段，產熱率基本保持不變，之后逐漸下降。根據式(15)可知，由于注入速率與注入流體的焓不變，當注采井之間形成穩定流動后，生產速率基本保持不變，故產熱率的變化與產出流體的焓有關，從而產熱率隨著生產溫度的降低而下降。

圖14 不同a和β條件下產熱率隨時間的變化Fig.14 Time-varying heat production rate under different a and β values

與圖7相比較，顯然，產熱率與生產溫度的變化規律一致，說明基質生產溫度是影響產熱率的主要因素。

在20 a的生產過程中，總產熱能隨時間線性增加(圖15)。當a=1.00和1.25時，其產熱速率(斜率)在整個運行期間基本一致，最終總產熱能也基本一致。說明當a較小時，裂縫密度的變化對總產熱能的影響可忽略不計。其余情況下，僅在開始一段時間內斜率相同，隨后，a越大，斜率相差越大。a相同時，隨著密度的增大，產熱速率越快，最終總產熱能相對越多。運行20 a后，a=1.00時，總產熱能從β=0.100的1.13×1015J增加到β=0.250的1.24×1015J，增加了9.73%。a=2.00時，總產熱能從β=0.100的1.01×1015J增加到β=0.250的1.24×1015J，增加了22.77%。說明a較大時，增加裂縫密度有利于提高總產熱能。

圖15 不同a和β條件下總產熱能隨時間的變化Fig.15 Time-varying total heat production under different a and β values

3.5.3 產熱效率ηh

20 a生產過程中不同a和β條件下產熱效率隨時間的變化如圖16所示。結果顯示，初始階段ηh發生急劇下降，后達到穩定，后期由于生產溫度的降低，使其逐漸下降。a相同時，密度β越大，前期產熱效率相對越小。β相同時，a越大，前期相對應的產熱效率越大。由于溫度下降導致流體對流換熱效率降低，所以，產熱效率后期下降的降幅與儲層溫度下降速度相關。

圖16 不同a和β條件下產熱效率隨時間的變化Fig.16 Time-varying heat production efficiency under different a and β values

4 結論

a.考慮了裂縫變形對熱生產的影響。注入壓力在初始階段明顯增加，并在注采速率穩定時注入壓力達到最大值。之后注入壓力逐漸降低，因為發生冷收縮導致的裂縫變形將提高流動效率。

b.長裂縫占比與裂縫密度(β)越大，流體流線分布越均勻，換熱面積越大，致使生產溫度下降速度變慢，增長EGS壽命。反之，竄流越明顯，生產溫度降低越快，EGS壽命越短。

c.提高長裂縫占比與裂縫密度可明顯改善儲層熱性能。在本研究設置的所有案例中，熱突破時間可從0.1 a延長至16.0 a，總產熱能從1.01×1015J提高至1.25×1015J。

d.在恒速注入的情況下，隨著裂縫長度指數a增大，生產溫度與產熱率降低越快，更早達到熱突破，從而縮短EGS壽命，降低總產熱能。當a相同時，裂縫密度β越大，生產溫度與產熱率降低越慢，熱突破時間和EGS壽命延長，總產熱能提高。熱突破時間與EGS壽命越長，生產溫度降低越少，更有利于實現商業化開采；在符合商業開采條件下，盡可能提高總產熱能，更好地實現熱提取。

符號注釋：

a為裂縫長度指數；b為擬合參數，Pa；Cd為等效比熱容，J/(kg·℃)；Cf為流體比熱容，J/(kg·℃)；Cs為巖石比熱容，J/(kg·℃)；d為最大閉合量，m；D為縫網分形維數，D取2；e為體積應變；f為與裂縫粗糙度相關摩阻系數；F為參數計算函數；g為重力加速度，9.8 m/s2；gi為加速度分量，m/s2；hinj、hpro分別為注入和產出流體的熱焓，J/kg；i=x,y,z；j=x,y,z；qinj、qpro分別為注入、產出速率，kg/s；qT為巖石與裂縫之間熱傳導，℃；QT為熱源，J/(s·m3)；k為滲透率，m2；Kn為接觸狀態下裂縫的法向剛度，Pa/m；L為研究區域的大小，m；M為流體的體積模量，Pa；n(l,L)為模型內裂縫長度l在[l,l+dl]內的裂縫數量；N為源項，L/s；p為孔隙壓力，Pa；pf為裂縫壓力，Pa；pinj、ppro分別為注入井、生產井的壓力，Pa；t為時間，s；△t為時間變化量，s；上標t、t+1為時間步；T為巖石溫度，℃；Tinj為注入流體的溫度，℃；Ts為基質的初始溫度，℃；T(x,t)為與位置和時間相關的溫度，℃；△up為剪切滑移位移，m；v為注入速率，m/s；vf為裂縫流體流速，m/s；vi為i方向上的速度，m/s；vi,j為i方向上的速度在j軸上的投影，m/s，其他同；w為裂縫寬度，m；w0為零應力對應裂縫寬度，m；△wm為純拉壓應力作用下的裂縫寬度增量，m；△wd為滑移引起的裂縫寬度增量，m；為相鄰兩裂縫單元平均寬度，m；wf為裂縫孔徑，m；z為深度，m；β為裂縫密度常數；△εij為應變張量的分量；ηp為注入泵效率，一般取80%；λf為流體導熱系數，W/(m·℃)；λeff為多孔介質的有效導熱系數，W/(m·℃)；λs為巖石導熱系數，W/(m·℃)；μ為流體黏度，Pa·s；ρ為密度，kg/m3；ρf為流體密度，kg/m3；ρm為流體平均密度，kg/m3；ρs為巖石密度，kg/m3；(ρCd)eff為多孔介質的等效體積熱容，J/(m3·℃)；σij,j為應力張量的分量，Pa；σn為裂縫面法向應力，Pa；σ'n為裂縫面法向有效應力，σ'n=σn–pf，Pa；△σ'n為有效正應力增量，Pa；φ為巖石孔隙率，%；ψ為膨脹角，(°)；▽為哈密頓算子；erf()為誤差函數；δ()為單位步長函數。