基于嵌入式離散裂縫模型的增強型地?zé)嵯到y(tǒng)熱—流—力—化耦合分析

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

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

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

0 引言

地?zé)崮苁莾Υ嬖诘貧ぶ械囊环N穩(wěn)定可再生能源，可應(yīng)用于發(fā)電、供暖、工業(yè)和農(nóng)業(yè)生產(chǎn)等多個方面。地?zé)豳Y源可以分為淺層地?zé)?、水熱型地?zé)岷透蔁釒r等。據(jù)估計，全球范圍內(nèi)埋深3.0～10.0 km的干熱巖資源蘊藏的熱能為全球石油、天然氣和煤炭儲藏能量的30倍[1]。我國干熱巖資源量為2.52×1025J，約占世界資源總量的1/6，其中埋藏地下3.0～5.0 km深度的干熱巖資源量為我國化石能源的80倍[2-4]。由于干熱巖儲層地質(zhì)條件的差異性和復(fù)雜性，以及現(xiàn)有儲層改造技術(shù)的“不可復(fù)制”性，造成其開發(fā)難度大[5-6]。目前，我國地?zé)崮艿睦弥饕运疅嵝蜑橹?，高溫干熱巖的開發(fā)仍處于探索研究和發(fā)展階段。

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

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

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

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

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

1 模型構(gòu)建

1.1 EDFM模型

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

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

1.2 滲流模型

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

第i條裂縫質(zhì)量守恒方程：

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

1.3 換熱模型

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

基巖能量守恒方程：

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

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

1.4 力學(xué)模型

基于線彈性小變形假設(shè)，考慮孔隙壓力、熱應(yīng)力變化影響，忽略慣性項和體積力后，巖體受力本構(gòu)方程為：

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

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

圖2 XEFM 結(jié)點加強方式圖

1.5 化學(xué)模型

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

第i條裂縫的溶質(zhì)運移方程為：

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

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

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

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

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

式中?0表示儲層初始時刻的孔隙度；Ms表示礦物摩爾質(zhì)量，kg/mol?？紫抖鹊淖兓瘜B透率的影響通過Kozeny-Carman方程確定：

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

2 數(shù)值方法

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

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

圖3 THMC耦合求解流程圖

3 模型驗證

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

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

圖4 兩種方法的網(wǎng)格劃分圖

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

圖5 兩種方法的數(shù)值解對比圖

4 數(shù)值算例

在EGS長期取熱過程中，儲層中的滲流壓力、溫度、礦物離子濃度和巖體隨時間發(fā)生變化，為明晰其時空演化規(guī)律，應(yīng)用所開發(fā)和驗證的THMC四場耦合模型及程序進行模擬。本文設(shè)置算例主要用于檢驗?zāi)Ｐ秃陀嬎惴椒ǖ目尚行裕瑫簝H考慮注入流體與儲層中的礦物石英相互作用。石英在水溶液中發(fā)生水解反應(yīng)而溶解生成硅酸（H4SiO4），或者其等價形式Si(OH)4，采用下式[33]：

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

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

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

表1 模型主要參數(shù)表

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

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

所用二維裂縫型地?zé)醿訋缀文Ｐ腿鐖D8-a所示，左下角綠色圖標表示注入井，右上角紅色圖標表示生產(chǎn)井，注入井和生產(chǎn)井之間由人工裂縫貫穿，儲層共預(yù)置裂縫70條?；谒O(shè)置的模型參數(shù)對不同網(wǎng)格數(shù)量下的開采溫度隨開采年限變化進行計算，結(jié)果如圖8-b所示，可以看出當計算域網(wǎng)格數(shù)增加到24 820個單元繼續(xù)增加時，開采溫度基本不隨網(wǎng)格數(shù)量變化而變化。故本文采用該網(wǎng)格劃分方案進行計算域離散，結(jié)果如圖8-c所示，其中，基巖網(wǎng)格數(shù)為22 801，裂縫網(wǎng)格數(shù)為2 019。

圖8 裂縫型儲層模型和網(wǎng)格劃分示意圖

4.1 演化規(guī)律

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

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

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

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

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

4.2 影響因素分析

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

4.2.1 儲層巖石參數(shù)

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

圖10 儲層巖石參數(shù)對取熱溫度影響圖

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

4.2.2 注水參數(shù)

一般來講，注入流體參數(shù)與當?shù)毓こ虠l件緊密相關(guān)。為了研究注入流體參數(shù)對開采溫度的影響，本文分別模擬了注水壓力和注入流體物性對開采溫度的影響，其他參數(shù)設(shè)置不變，結(jié)果如圖11所示?？梢钥吹?，兩種情況下開采溫度均先維持約10年的穩(wěn)定階段，然后出現(xiàn)不同幅度的下降階段。當系統(tǒng)運行第40年時，開采壓力由9 MPa、12 MPa、14 MPa增大到15.5 MPa時，對應(yīng)的開采溫度分別為229.1 ℃、221.8 ℃、215.5 ℃和210.1 ℃，如圖11-a所示?？梢姡⑷雺毫υ酱?，隨著開采時間增加，開采溫度降低越快。注入壓力改變注入流體的滲流速率，影響其取熱效率。當開采壓力由9 MPa增大到15.5 MPa，系統(tǒng)運行第40年時，單位厚度開采井的凈取熱功率前者比后者低54%。

圖11 注水參數(shù)對取熱溫度影響圖

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

4.3 水—巖反應(yīng)影響

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

圖12 水—巖反應(yīng)對取熱溫度的影響圖

5 結(jié)論

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

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

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

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

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