基于嵌入式離散裂縫模型的增強型地熱系統熱—流—力—化耦合分析

韓東旭 張煒韜 焦開拓 宇 波 李庭宇 鞏 亮 王樹榮

1. 北京石油化工學院機械工程學院 2. 中國石油大學（華東）新能源學院 3. 西安交通大學動力工程多相流國家重點實驗室

4. 中國空氣動力研究與發展中心低速空氣動力研究所 5. 浙江大學能源高效清潔利用全國重點實驗室

0 引言

地熱能是儲存在地殼中的一種穩定可再生能源，可應用于發電、供暖、工業和農業生產等多個方面。地熱資源可以分為淺層地熱、水熱型地熱和干熱巖等。據估計，全球范圍內埋深3.0～10.0 km的干熱巖資源蘊藏的熱能為全球石油、天然氣和煤炭儲藏能量的30倍[1]。我國干熱巖資源量為2.52×1025J，約占世界資源總量的1/6，其中埋藏地下3.0～5.0 km深度的干熱巖資源量為我國化石能源的80倍[2-4]。由于干熱巖儲層地質條件的差異性和復雜性，以及現有儲層改造技術的“不可復制”性，造成其開發難度大[5-6]。目前，我國地熱能的利用主要以水熱型為主，高溫干熱巖的開發仍處于探索研究和發展階段。

在對深部干熱巖熱能開發中，需采用人工壓裂的方式提高儲層內的導流能力，即在低滲透性儲層建造高滲透性的人工裂縫，形成增強型地熱系統（Enhanced Geothermal System，EGS）。不同于常規/非常規油氣藏的“取物”過程，地熱能開采為“取熱”過程。在人工壓裂過程中，EGS注采循環井會出現工質“注不進，采不出”的問題，實質是裂縫導流能力不足或裂縫網絡溝通半徑不夠，導流能力不足是流、熱、力和化等因素造成的裂縫開度變化的綜合結果，裂縫網絡溝通半徑不足是熱及形變誘導應力對巖體的剪切破裂和裂縫的起裂、擴展的綜合影響[7]。而在長期取熱過程中，冷流體從注入井注入后，工質在裂隙巖體中以滲流的形式流向采出井，并將儲存在巖體中的熱量提取出來。在此過程中，因溫度變化產生的熱應力會導致巖體變形，同時，在滲透壓力和溫度的共同作用下，巖體中礦物溶解/沉淀會緩慢改變流體通道，即EGS儲層是一個典型受熱—流—力—化綜合作用的系統。所以有必要針對EGS開展熱力、水力、巖體變形以及水—巖反應（Thermal-Hydraulic-Mechanical-Chemical，THMC）耦合機理的綜合研究。

目前，國內外學者已經做了大量相關實驗研究：在多場耦合作用方面，針對固—熱耦合作用下花崗巖的力學、滲流特性進行實驗，指出巖體力學性質變化是高溫熱破裂和石英晶格相變導致，如水力壓裂、熱壓裂、化學刺激等儲層增產的方法，實際是在改變儲層的滲透率[8]。Kamali等[9]通過對熱—流—力—化多場耦合作用下的裂縫滲透率和孔徑變化進行測試，發現卸載時裂縫滲透率和孔徑不完全恢復，升高注入流體的溫度會提高裂縫恢復百分比。Yang等[10]對花崗巖進行冷熱處理后，測試了熱—流—力耦合參數，如導熱系數、滲透率、Biot系數等，建立了參數預測模型，預測結果與實驗吻合很好。由于實驗室使用的巖體試樣多為地表或淺層地下采集，即使取自深部巖體，也存在擾動應力的影響，且實驗室多為高溫冷卻后的性質測試，很難將深部地層高溫高壓環境實時結合[8]。趙陽升[11]認為巖體力學特性在實驗室尺度與實際工程特征尺度之間存在尺度效應和時間效應，實驗室力學參數應用于巖體變形分析時，折減系數取值差異會導致工程變形預計存在很大偏差。

針對實驗研究上的上述缺陷，數值模擬作為EGS研究的有效手段之一，與實驗研究相互補充，對于干熱巖儲層的開發、方案制定、運行優化等有重要意義[12]。EGS儲層由基巖和裂縫共同組成，是一種人工地熱系統。針對裂隙巖體，常用的裂縫表征模型有等效介質模型、雙重孔隙介質模型和離散裂縫模型（Discrete Fracture Model，DFM）。相比于前兩者，DFM能夠更加準確地描述裂縫的分布和走向，該模型適合以裂縫為主要滲流通道的裂隙巖體。Song等[13]通過COMSOL軟件，基于DFM對其所提出的多分支井的EGS開發方式進行了熱—流兩場的數值模擬；Sun等[14]針對二維EGS儲層，基于DFM建立了熱—流—力三場耦合模型；Yao等[15]在考慮“局部非平衡傳熱”的基礎上，基于DFM構建了EGS開發過程中三維熱—流—力耦合模型，并對Desert Peak EGS 項目進行了模擬。

雖然DFM描述裂隙巖體具有較高的精度，但是當裂縫較復雜時，該模型需要大量的計算網格，很難在計算精度和計算成本之間做到平衡。為了解決這一問題，研究人員提出并發展了嵌入式離散裂縫模型（Embedded Discrete Fracture Model，EDFM）[16]。在該模型中，基巖和裂縫網格獨立生成，裂縫網格降維嵌入到基巖網格系統中，裂縫和基巖之間的信息傳遞通過交換系數與傳遞系數進行。同DFM相比，EDFM使裂隙巖體的網格劃分更加靈活，同時還兼顧了計算精度和計算效率。因此，近年來EDFM被越來越多的應用于裂隙巖體的數值模擬中。Wang等[17]提出了一種將EDFM和擴展有限元（The Extended Finite Element Method，XFEM）相結合的多場耦合策略，用于描述裂縫和多孔介質中的流動、變形和裂縫擴展。Sangnimnuan等[18]基于EDFM開發了非常規油藏中兩相流—力耦合模型。Tran等[19]基于EDFM和改進的Bandis模型，求解了流動—溶質運移—誘導應力—斷裂力學的耦合問題。Li等[20]基于EDFM和XFEM框架，建立了EGS的熱—流—力耦合模型。此外，Xu等[21-23]也基于EDFM或改進的EDFM對復雜裂縫型儲層中的流動、傳熱和傳質進行了數值模擬。

盡管近幾年，學者們基于EDFM在裂隙巖體流—力、流—力—化、熱—流—力的耦合建模方面，以及針對復雜裂縫型儲層的模擬計算方面取得了進展，但是還缺少基于EDFM的熱—流—力—化四場耦合模型研究。因此，本文基于EDFM框架，建立了增強型地熱系統的THMC耦合模型，采用有限體積法（Finite Volume Method，FVM）和XFEM相結合的方式進行求解，并對二維增強型地熱系統儲層中的流場、溫度場、離子濃度場和位移場的時空演化過程及不同參數影響進行模擬與分析。

1 模型構建

1.1 EDFM模型

Li[24]首先提出了EDFM模型，隨后Moinfar等[25]和Tene等[26]對其進行了發展。在該模型中，裂縫降維嵌入到基巖網格中，基巖和裂縫網格在幾何上幾乎完全獨立，如圖1所示。在實際物理問題中，為了實現基巖和裂縫、裂縫和裂縫之間的物理性交互，包括質量、能量、濃度等，研究人員引入了兩個關鍵參數：交換系數和傳遞系數，根據基巖和裂縫的連接關系對其進行計算[27]。基巖和裂縫之間存在4種連接關系[16]：①相鄰基巖網格單元（M1、M2、M3、M4）之間的連接；②同一條裂縫相鄰裂縫網格單元之間的連接（F1、F2、F3、F4、F5）；③裂縫網格單元與所在基巖網格單元之間的連接；④交叉裂縫段網格單元之間的連接。該模型能夠實現對裂縫型儲層采用結構化網格進行剖分，避免生成比較復雜的非結構化網格，能夠降低計算成本。因此，本文基于該方法的思想，構建了描述EGS取熱過程的流動、傳熱、變形和溶質運移的THMC耦合模型。

圖1 二維嵌入式離散裂縫模型示意圖

1.2 滲流模型

EGS的巖性大部分是滲透率極低的基地花崗巖，流體在儲層中流動滿足達西定律，其流動過程采用質量守恒方程和達西方程描述，二維平面上不考慮重力影響，基于EDFM的基本思想，基巖質量守恒方程為[16]：

第i條裂縫質量守恒方程：

式中上標m和fr分別表示基巖和裂縫系統；pm、pfr分別表示基巖和裂縫中的壓力，Pa；ρf表示流體密度，kg/m3；?m、?fr分別表示基巖和裂縫的孔隙度；ct表示總壓縮系數，Pa－1；μf表示流體動力黏度，Pa·s；k表示滲透率張量，m2；Qm-fri、Qfri-m、Qfri-frj分別表示裂縫到基巖的單位體積流量、基巖到裂縫的單位體積流量、裂縫j到裂縫i的單位體積流量，s－1；為確保守恒，基巖和裂縫間的流量交換滿足V表示裂縫或基巖的網格單元體積，m3；Qw表示井源項，s－1。

1.3 換熱模型

同樣基于EDFM思想，采用兩套能量守恒方程分別描述基巖系統、裂縫系統以及它們之間的能量交換，并認為流體和基巖之間處于局部熱平衡狀態[29]。

基巖能量守恒方程：

第i條裂縫能量守恒方程：

式中T表示溫度，℃；v表示滲流速度，m/s；cpf表示流體比熱容，J/(kg·℃)；(ρcp)eff表示有效物性參數；λeff表示有效導熱系數，W/(m·℃)；Efri-m表示裂縫i到基巖的能量交換，W/m3；Efri-frj表示裂縫j到裂縫i的能量交換，W/m3，Ew表示井源項，W/m3。

1.4 力學模型

基于線彈性小變形假設，考慮孔隙壓力、熱應力變化影響，忽略慣性項和體積力后，巖體受力本構方程為：

式中σ表示柯西應力張量，Pa；λ、G分別表示第一、第二拉梅系數，Pa；α表示Biot系數；βT表示熱膨脹因子，Pa/℃；I表示單位張量；p0表示參考壓力，Pa；T0表示參考溫度，℃；ε表示柯西應變張量，；u表示位移矢量，m。對裂隙巖體的變形問題進行計算時，裂縫兩側的位移不連續，屬于強不連續問題，而XFEM能夠有效地求解強不連續問題，位移求解采用下式[30]：

式中上標 FE、H、tip 和 J 分別表示常規單元、貫穿單元、裂尖單元及交叉單元（圖2）。

圖2 XEFM 結點加強方式圖

1.5 化學模型

化學模型包括兩部分，一部分描述礦物離子在儲層中的運移和擴散，另一部分描述儲層中的水—巖反應。礦物離子在儲層中運移和擴散通過溶質運移方程描述，基巖的溶質運移方程為：

第i條裂縫的溶質運移方程為：

式中C表示溶質濃度，mol/m3；D表示擴散系數，m2/s；ψfri-m（或ψm-fri）、ψfri-frj分別表示基巖與裂縫間的離子濃度交換、裂縫與裂縫間的離子濃度交換，mol/(m3·s)；為確保守恒，基巖與裂縫間的離子濃度交換滿足Mw表示井源項，mol/(m3·s)；Qr表示化學反應源項，Qr=?rnρs，mol/(m3·s)；ρs表示礦物密度，kg/m3；rn表示化學反應速率，mol/(kgw·s)。

EGS中流體與巖石的化學反應速率極低，一般遇到的反應類型為表面控制的化學反應，同TOUGHREACT模擬軟件一樣，本文應用的是Lasaga提出的化學反應速率方程[31]：

式中下標n表示礦物索引；kn表示礦物反應速率常數，mol/(m2·s)；An表示礦物反應比表面積，m2/kgw；Qn表示礦物離子活度；Keq平衡常數；θ、η由實驗決定。礦物反應速率常數與溫度相關，采用Arrhenius的形式描述：

式中k25表示25 ℃下的速率常數，mol/(m2·s)；Ea表示活化能，J/mol；R表示氣體常數，J/(mol·K)；T表示絕對溫度，K。

由水—巖反應造成的儲層中基巖和裂縫的孔隙度變化通過下式確定：

式中?0表示儲層初始時刻的孔隙度；Ms表示礦物摩爾質量，kg/mol。孔隙度的變化對滲透率的影響通過Kozeny-Carman方程確定：

式中k0表示儲層初始時刻的滲透率，m2。

2 數值方法

上述數學模型中包括p、T、u、C等4個未知變量，涉及流場、溫度場、力場、濃度場的求解，各場求解過程中相互耦合，相互作用，耦合過程設置如下：①在流場求解中，耦合基巖和裂縫的能量守恒方程、溶質運移方程和水—巖反應方程，依賴于溫度變化對流體密度、黏度的影響以及溶質運移和水—巖反應對儲層滲透率、孔隙度的改變；②在溫度場求解中，耦合基巖和裂縫的質量守恒方程、溶質運移方程和水—巖反應方程，依賴于壓力場變化對滲流速率的影響以及溶質運移和水—巖反應對儲層滲透率、孔隙度的影響；③在化學場求解中，耦合基巖和裂縫的質量守恒方程、基巖和裂縫的能量守恒方程，依賴于壓力場變化對滲流速率的影響以及溫度變化對流體密度、黏度的影響；④在力場求解中，考慮滲流壓力和溫度變化對儲層變形的影響，因此其求解是在熱、流、化三場求解收斂基礎上進行。

數值求解時，首先采用FVM求解基巖和裂縫的滲流場，獲得壓力場的解（pm，pfr），再求解基巖和裂縫的達西方程，獲得速度場的解（vm，vfr）；然后，采用同樣方法求解基巖、裂縫的能量守恒方程，獲得溫度場的解（Tm，Tfr）；將得到的流場的解和溫度場的解帶入溶質運移方程和水—巖反應中，求解得到礦物離子濃度場的解（Cm，Cfr）；最后采用XFEM求解力平衡方程，獲得位移場的解（u）。每個時層采用順序迭代耦合策略，求解流程和求解方法如圖3所示。

圖3 THMC耦合求解流程圖

3 模型驗證

本文建立THMC耦合模型是在本文參考文獻[20]所建立THM耦合模型的基礎上添加了化學模塊，關于THM耦合模型及數值方法驗證已在文獻中進行了詳細描述，在此不再敘述。由于化學反應及礦物離子的運移受儲層溫度、流體流動影響較大，因此本文將獲得的壓力場、溫度場、離子濃度場的解與MRST[32]中Fully Resolved Solution Method（簡 稱FRSM）在細尺度網格上獲得的解進行對比。區別于EDFM的裂縫處理方法，FRSM以裂縫寬度作為網格單元的離散尺寸對整個計算域和裂縫進行剖分，然后分別對基巖和裂縫所在的網格單元定義其物理屬性。

所用物理模型尺寸和裂縫位置如圖4-a所示，模型尺寸設定為10 m×5 m，裂縫寬度設定為0.04 m。兩種方法得到的計算域網格劃分結果如圖4-b和圖4-c所示，采用FRSM劃分計算域獲得的基巖和裂縫細尺度網格數為31 250，采用EDFM劃分網格獲得的計算域基巖網格數為20 000，裂縫網格數為209。

圖4 兩種方法的網格劃分圖

模型左邊界模擬注入井，純水注入，注入壓力為1.5 MPa，注入溫度為20 ℃，右邊界模擬生產井，開采壓力為0.1 MPa，其余邊界為非滲透、絕熱邊界。初始壓力為0.1 MPa，初始溫度為180 ℃。假設儲層礦物離子在運移擴散之前處于平衡狀態，同時儲層中礦物離子以二氧化硅為主，擴散系數為1.0×10－9m2/s，礦物離子的初始濃度為3.8 mol/m3，基巖和裂縫的滲透率、孔隙度、密度等參數與本文參考文獻[20]保持一致。時間步長取0.1 d，總模擬時間365 d，儲層中溫度變化和溶質運移求解是在儲層壓力場分布穩定的基礎上進行，計算結果如圖5所示。可以看到，本文基于EDFM計算的壓力場、溫度場、濃度場與采用FRSM在細尺度網格上獲得的計算結果基本一致。恒定注采壓力下采用EDFM和采用FRSM計算的生產井溫度最大相對誤差為0.71%，平均為0.70%，溶質濃度相對誤差最大為0.4%，平均為0.12%，采用EDFM計算用時為2 652 s，采用FRSM計算用時為3 650 s，前者用時較后者減少了27.3%，證明本文采用的模型和數值方法具有可行性。

圖5 兩種方法的數值解對比圖

4 數值算例

在EGS長期取熱過程中，儲層中的滲流壓力、溫度、礦物離子濃度和巖體隨時間發生變化，為明晰其時空演化規律，應用所開發和驗證的THMC四場耦合模型及程序進行模擬。本文設置算例主要用于檢驗模型和計算方法的可行性，暫僅考慮注入流體與儲層中的礦物石英相互作用。石英在水溶液中發生水解反應而溶解生成硅酸（H4SiO4），或者其等價形式Si(OH)4，采用下式[33]：

式中，正反應過程表示礦物溶解，逆反應過程表示礦物沉淀。由式（9）可知，礦物化學反應中除了確定反應速率常數外，還需要確定礦物離子活度與給定溫度下平衡常數的比值。本文石英水解反應速率常數k25、反應活化能Ea參考文獻[34]選取，不同溫度平衡狀態下石英的溶解度采用下式[35]：

式（14）可以確定溫度在0～300 ℃時飽和水壓下的石英溶解度。

采用五點布井方式，即一個注入井周圍分布4個開采井，反之亦然。考慮所有井的流場呈對稱分布，筆者取五點開采區域的1/4進行研究。注入井壓力為15.5 MPa，注入溫度為60 ℃，生產井壓力為4.77 MPa，四周邊界為非滲透、絕熱、位移約束條件，礦物離子濃度梯度為0。儲層初始壓力為5.0 MPa，初始溫度為236 ℃，石英占儲層總礦物體積分數為66%。對于儲層中孔隙度和滲透率的非均勻性，采用高斯隨機分布和Carman-Kozeny 經驗式生成[36-37]（圖6）。儲層初始孔隙度分布為3.5%～4.5%，平均值為3.99%。初始滲透率分布為0.17～0.38 mD，平均值為0.26 mD。儲層巖石物性參考恰卜恰地熱田設置[38]，部分參數有改動，具體見表1。在取熱過程中，注入水的密度、黏度隨溫度、壓力不同而發生變化，圖7給出了水的密度、黏度隨溫度、壓力的變化情況[39]。其中，圖7-a兩條黑色實線與邊界圍成的區域表示壓力在4.77～15.5 MPa、溫度在20～236 ℃時，流體密度的變化情況；圖7-b兩條黑色實線與邊界圍成的區域表示在相同壓力和溫度范圍下的流體黏度的變化情況。

表1 模型主要參數表

圖6 儲層孔隙度和滲透率非均勻分布圖

圖7 水的密度、黏度隨溫度、壓力變化圖

所用二維裂縫型地熱儲層幾何模型如圖8-a所示，左下角綠色圖標表示注入井，右上角紅色圖標表示生產井，注入井和生產井之間由人工裂縫貫穿，儲層共預置裂縫70條。基于所設置的模型參數對不同網格數量下的開采溫度隨開采年限變化進行計算，結果如圖8-b所示，可以看出當計算域網格數增加到24 820個單元繼續增加時，開采溫度基本不隨網格數量變化而變化。故本文采用該網格劃分方案進行計算域離散，結果如圖8-c所示，其中，基巖網格數為22 801，裂縫網格數為2 019。

圖8 裂縫型儲層模型和網格劃分示意圖

4.1 演化規律

注入壓力可改變儲層中流體的滲流速率和孔隙壓力，而滲流速率會影響流體與儲層高溫巖體的對流換熱、礦物離子的運移，因此有必要對恒定注入壓力下儲層中壓力場的時空演化進行分析。圖9-a展示了系統運行40年過程中3種不同年份儲層中滲流壓力空間分布。可以看到，裂縫作為儲層中流體滲流的主要通道，使儲層中的壓力分布呈現非均勻性。因為裂縫介質滲透率大，滲流速度較快，基巖系統滲透率小，滲流速度滯后，因此裂縫介質中的壓降較快。在裂縫介質與基巖系統水力交換影響下，裂縫介質附近基巖系統的壓降也較快。當注入井和生產井之間存在貫穿裂縫時，這種影響更加明顯；同時，還可以發現隨著冷流體不斷注入，注入井附近的溫度下降，流體的黏滯性增強，流動阻力增大，表現為隨開采年限增加，注入井附近的壓力梯度逐漸增大。

圖9 儲層中壓力場、溫度場、濃度場和位移場隨時間演化規律圖

儲層中的溫度不僅對流體的密度、黏度產生影響，而且還對水—巖反應速率產生影響，同時溫度變化產生的熱應力變化也會導致儲層中的巖體變形，產生位移變化。圖9-b展示了系統運行40年過程中3種不同年份溫度場的空間分布。可以看到，在恒定注入壓力下，冷流體不斷被注入，開采初始階段受導熱和對流換熱作用，注入井附近的低溫水與高溫基巖發生熱量交換，水溫升高，而高溫基巖的溫度下降，隨著開采年限增加低溫區域逐漸擴大，此時開采井附近巖體溫度變化緩慢，在熱突破形成之前，系統可維持較長時間的高熱量輸出；由于裂縫中流體滲流速度高于周圍基巖，裂縫介質中熱對流明顯，沿著滲流方向在裂縫附近巖體的溫度變化較快，尤其是注采井之間存在貫穿裂縫時，該現象更加明顯；此外，需要指出的是雖然儲層部分區域裂縫密集，但是由于這些裂縫并未完全貫穿，因此在這些裂縫附近流體與高溫巖體的對流換熱作用體現并不明顯。

在長期取熱過程中，儲層中水—巖反應產生的礦物離子會隨著流體滲流而運移、擴散，儲層中礦物離子運移、擴散會打破儲層內礦物離子固有的平衡狀態，導致在儲層不同位置將發生礦物離子的溶解和沉淀反應，使儲層中的孔隙度隨之改變，反過來影響流場和溫度場。圖9-c展示了系統運行40年過程中3種不同年份儲層礦物離子濃度空間分布。可以看到，儲層中礦物離子濃度分布趨勢與溫度變化趨勢比較類似，兩者都受流體對流作用影響，同時溫度分布還會直接影響水—巖反應強弱，儲層溫度高的區域，礦物溶解速率大，礦物離子濃度高；但是不同于溫度場的變化，受溶質擴散系數影響，在裂縫及附近區域礦物離子濃度分布明顯低于周圍基巖系統，裂縫對儲層中礦物離子運移、擴散作用明顯。

儲層中滲流壓力、溫度、化學溶蝕長期作用會導致儲層巖體發生變形、產生位移變化，高溫巖體變形會造成儲層結構顆粒脫落、裂縫滲透性改變等。圖9-d展示了系統運行40年過程中3種不同年份儲層巖體產生的x方向位移和y方向位移空間分布。可以看到，隨著開采時間的增加，x方向位移和y方向位移的區域范圍和最大值都在增加。當儲層中不斷注入流體時，孔隙水壓力會改變，孔隙水壓力作用于裂縫面從而使巖體有效應力降低，導致裂縫介質的法向張開度發生改變，同時高溫巖體也會因溫度下降而發生收縮，裂縫面在熱應力作用下進一步張開，儲層中x方向位移和y方向位移主要集中在裂縫面附近。此外，在溫度、注入壓力發生較大梯度變化的區域，如注入井附近，越易產生較大的位移。

4.2 影響因素分析

本節主要從儲層巖石參數、注水參數兩個方面，分析各參數變化對EGS取熱的影響，其中儲層巖石參數包括基巖滲透率、裂縫滲透率和裂縫開度；注水參數包括注水壓力和流體物性。

4.2.1 儲層巖石參數

圖10-a展示的是不同基巖滲透率對應的開采溫度變化，模擬時其他參數保持不變，選用非均勻分布滲透率的最大值、最小值和平均值進行對照。可以看到，在開采前10年中，四種不同基巖滲透率對應的開采溫度都能保持穩定的高溫度輸出，隨后逐年下降，并且基巖滲透率越大開采溫度下降越快，當系統運行第40年時，基巖滲透率為0.17 mD時對應的開采溫度下降了15℃，單位厚度開采井的凈取熱功率由1.86 kW降為0.85 kW，而基巖滲透率為0.38 mD時對應的開采井溫度下降幅度增大了10.6%，單位厚度開采井的凈取熱功率由3.5 kW降為1.26 kW。可見，儲層低滲透率雖然能夠維持較長年限較高的開采溫度，但這并不表明儲層滲透率越低越好，因為低滲透率會造成流體流動緩慢，儲層中的開采井質量流率減小，凈采熱功率降低，同時還可能造成流體通道堵塞，影響系統運行。

圖10 儲層巖石參數對取熱溫度影響圖

同樣，保持其他參數設置不變，模擬改變裂縫滲透率、裂縫開度對開采溫度的影響，結果如圖10-b和圖10-c所示。可以看到，4種不同裂縫滲透率和裂縫開度都可使系統維持一定年限穩定的高溫度輸出，但是這種穩定的高溫度輸出時間比改變基巖滲透率獲得的時間要短。可見，改變裂縫參數比改變基巖參數對開采溫度的影響強烈。而在開采溫度下降階段，當系統運行第40年時，裂縫滲透率由10－9m2繼續增大到10－8m2時，前者相比于后者，開采溫度和單位厚度開采井的凈取熱功率分別僅高了2 ℃和0.01 kW；裂縫開度由0.000 5 m增大到0.001 m時，前者相比于后者，開采溫度僅高了4 ℃，單位厚度開采井的凈取熱功率方面僅低了0.05 kW。可見，當裂縫滲透率在較大范圍內變化或者裂縫開度在較小寬度范圍內變化時，其對開采溫度和凈取熱功率的影響變弱。

4.2.2 注水參數

一般來講，注入流體參數與當地工程條件緊密相關。為了研究注入流體參數對開采溫度的影響，本文分別模擬了注水壓力和注入流體物性對開采溫度的影響，其他參數設置不變，結果如圖11所示。可以看到，兩種情況下開采溫度均先維持約10年的穩定階段，然后出現不同幅度的下降階段。當系統運行第40年時，開采壓力由9 MPa、12 MPa、14 MPa增大到15.5 MPa時，對應的開采溫度分別為229.1 ℃、221.8 ℃、215.5 ℃和210.1 ℃，如圖11-a所示。可見，注入壓力越大，隨著開采時間增加，開采溫度降低越快。注入壓力改變注入流體的滲流速率，影響其取熱效率。當開采壓力由9 MPa增大到15.5 MPa，系統運行第40年時，單位厚度開采井的凈取熱功率前者比后者低54%。

圖11 注水參數對取熱溫度影響圖

模擬過程中，假定注入流體的密度、黏度不隨儲層的溫度、壓力變化而變化，與初始值保持一致，即選用常物性流體注入時，在溫度下降階段，開采溫度隨年限增加下降幅度較小，如圖11-b所示。當系統運行到第40年時，采用常物性流體注入獲得的開采溫度為231.7 ℃，與采用變物性流體注入獲得的開采溫度相比，兩者偏差22 ℃。分析可知，在給定注入流體物性時，變物性流體會因為儲層的溫度、壓力分布不均而發生變化，結合圖7可知，當壓力恒定時，注入水的密度、黏度會隨溫度增大而降低。初始階段受注入壓力、溫度和儲層溫度場的分布影響，注入流體物性變化對開采溫度的影響較少，兩者可維持一定年限穩定高溫度輸出。隨著開采年限的增加，注入井附近的壓力梯度逐漸增加，見圖9-a。開采井附近仍舊維持較高的溫度場分布，造成變物性流體在開采井附近黏度較低而滲流速度偏大，單位時間內輸出的熱量比常物性流體大，儲層的溫度下降快，使變物性流體計算的開采溫度較常物性流體偏低。

4.3 水—巖反應影響

為探究EGS系統在取熱過程中受化學反應的影響程度，利用熱—流—力—化（THMC）耦合模型和熱—流—力（THM）耦合模型模擬儲層熱開采過程中水—巖反應對開采溫度的影響。水—巖反應中石英溶解速率參考文獻[40]選取，其中，實驗條件25 ℃、pH值為7時，對應的溶解速率為－10.4，其他參數設置保持不變，兩種耦合模型模擬的開采溫度隨時間變化如圖12所示。可以看到，兩種耦合下熱儲層均可維持約10年穩定的高溫度輸出，但是在溫度下降階段，隨著開采時間的增加，開采溫度變化差別較大。當系統運行第40年時，THMC耦合下的開采溫度為195 ℃，與THM耦合相比，兩者偏差為15 ℃，偏差幅度為7%。水—巖反應對EGS儲層中的孔隙度產生影響，由式（12）可知，孔隙度改變會導致儲層固有滲透率發生變化，進而影響流體的滲流速率以及與高溫巖石的對流換熱效率。設置算例中儲層孔隙度為非均質分布，平均孔隙度由初始的3.99%變為4.49%，孔隙度增大了12.5%。可見，在取熱過程中，礦物石英的水—巖反應主要表現為溶解過程。此外，需要說明的是EGS人工熱儲中的花崗巖主要成分除了石英外，還包括堿性長石、斜長石、角閃石及云母等，這些礦物部分水—巖反應強度比石英更大，長期取熱過程中對儲層孔隙結構改變不可忽略。因此通過數值模擬研究EGS系統取熱性能時，應該綜合考慮儲層各不同礦物化學反應造成的影響。

圖12 水—巖反應對取熱溫度的影響圖

5 結論

1）裂縫作為儲層中流體的主要滲流通道，在恒定注采條件下，裂縫介質及其附近基巖系統的滲流壓力、溫度、離子濃度變化較快，當注采井之間存在貫穿裂縫時，該現象更加明顯，而位移變化主要集中在裂縫面附近，在滲流壓力、熱應力梯度大的區域，如注入井附近，易產生較大位移。

2）基巖滲透率較低時，開采溫度下降緩慢，EGS可維持較長年限的較高開采溫度，但是滲透率較低時會造成儲層中流體流動緩慢，開采井質量流率減小，進而影響凈取熱功率。裂縫滲透率或者裂縫開度發生改變時，EGS均可維持一定年限的較高溫度穩定輸出，隨后出現不同幅度的下降，但是當裂縫滲透率增大到一定值或者裂縫開度減小到一定值時，改變裂縫參數對開采溫度的影響減弱。

3）考慮注入流體物性隨滲流壓力、儲層溫度變化時，開采溫度變化明顯，模擬EGS運行40年時，采用變物性流體比采用常物性流體計算的開采溫度偏低22 ℃。因此，對于實際儲層取熱分析時，需要考慮注入流體物性參數的影響。

4）EGS長期運行過程中水—巖反應會造成儲層孔隙度改變，進而對開采溫度產生影響，模擬EGS運行40年時，考慮水—巖反應比不考慮水—巖反應計算的開采溫度偏低為15 ℃，儲層平均孔隙度值增大12.5%。可見，長期運行的EGS水—巖反應的影響不可忽略。

本文通過自主編程，將所建立的THMC耦合模型和數值方法在二維EGS儲層問題中進行了實施，相關方法和結論可為EGS多場耦合研究提供借鑒和參考，未來將考慮儲層其他礦物組分的化學反應、力學變形對流熱化三場影響的雙向耦合并推廣至三維問題中，為EGS儲層的工程開發和利用提供一些準確、可靠的預測。