工質密度和儲層滲透率演化對EGS系統傳熱效能的影響

黃涵釵,孟凡震,李沐子,修占國,張樹翠,李致遠

(青島理工大學,山東 青島 266000)

在所有清潔能源中,地熱資源以其儲量大、范圍廣、綠色清潔的特點備受青睞,被視為化石能源枯竭后最具潛力的戰略接替能源[1]。其中人類能利用的地熱資源90%都來自于干熱巖[2]。目前干熱巖的熱能主要通過增強型地熱系統(Enhanced Geothermal System,下簡稱EGS)開發,進行熱能開采時,先從注入井泵入高壓冷水,高壓冷水流過熱儲干熱巖中的裂縫,吸收巖石的熱量,直到被迫以高壓熱蒸汽或汽水混合物的方式從生產井出來,實現采熱。EGS的壽命以冷卻區到達出口為止,此時稱熱貫通或熱突破[3-4]。維持流體采熱及延緩熱突破是EGS系統的關鍵,直接影響EGS系統的產能和壽命。

流體進入地熱儲層后的流動和采熱過程包含了復雜的熱流固耦合效應,而熱流固耦合效應必然導致流體密度和滲透率的演化,進而影響流體的采熱和儲層的熱突破。針對這一現象,目前國內外學者基于數值模擬開展了大量耦合研究,其中大多數是基于商業軟件實現的,其熱流固耦合考慮的耦合項不完整[5-8],或者僅為建立在預先設定的耦合變量上的間接耦合,即不同物理場的求解系統是獨立的[9-10],這一定程度上影響了耦合結果的準確性。少數EGS系統的熱流固耦合模擬是基于自編程程序實現的。與商業軟件基于正壓流體模型建立的耦合程序不同,這些都是基于斜壓流體密度模型建立的THM完全耦合自編程程序。Salimzadeh等[11]基于CSMP++設計了含單裂隙的熱儲層THM直接耦合程序,并基于立方定律探討了裂隙滲透率演化對裂隙溫度的影響。Cui Xin等[12]開發了等效多孔介質模型的THM雙向耦合程序,并探討不同耦合模式對收斂性的影響,結果表明了HM-T的雙向耦合在保證耦合結果準確性的同時,可以比THM的直接耦合有更好的收斂。不過,對于方程中還殘余的流體實時密度與參考密度之比,為了在計算流固場的耦合系數矩陣時,利用對稱性來減少計算量,則忽略了實時流體密度的變化,這實際上是對流體密度演化的不完全實現。斜壓流體密度模型在正壓流體模型基礎上進一步考慮了溫度對于流體密度的影響,實現了求解系統中THM不同主要變量之間的直接耦合,更加符合實際情形。總體來說,基于斜壓流體密度模型的THM完全耦合目前應用不多,基于該模型的傳熱效能分析也鮮有報道。

綜上,為了提高EGS儲層傳熱效能評價的準確性,本文建立包含傳熱效能指標評價的HM-T的雙向耦合自編程程序,探究THM完全耦合作用下流體密度演化、滲透率演化對EGS系統和傳熱性能指標的影響。

1 THM耦合的數學模型

1.1 控制方程

對于變形場,基于飽和多孔彈性介質的應力平衡,建立了熱彈性力學模型。在小變形假設下,考慮溫度應力和孔隙壓力的準靜態條件下的應力平衡模型由式(1)給出:

(1)

式中:D為彈性張量;βv為多孔介質的體積熱膨脹系數;Ks為多孔介質的體積模量;f為體力;未知量u為變形張量,T為溫度,p為孔壓;α表示Biot液力耦合系數;二維時,I=[1 1 0]T。

利用關于溫度和壓力的一階泰勒展開式與參考流體密度ρf 0可以近似表征斜壓流體的實際密度ρf。其中ρp和ρT分別為流體的體壓縮系數和熱膨脹系數。

ρf 0+ρf 0βp(p-p0)-ρf 0βT(T-T0)

(2)

則斜壓流體在多孔介質內達西流動時的控制方程由式(3)給出[12]：

(3)

對于傳熱場,假設巖石基質中的流體速度足夠慢,使得巖石基質中固體顆粒和流體始終處于局部熱平衡。此時,用統一的變量表征多孔介質流體和固體的溫度時,其傳熱的物理模型如式(4)[12]：

(4)

其中,

(ρc)m=φρfcf+(1-φ)ρscs

(5)

λm=φλf+(1-φ)λs

(6)

(7)

考慮滲透率演化時,使用李盼[13]基于室內實驗得到的花崗巖滲透率演化模型對EGS系統進行模擬。

k=k0e0.23pm-0.15pc-0.0023Ts-0.21Tin

(8)

式中:pm為孔隙壓力,MPa;pc為圍壓,MPa;Ts為巖樣溫度,K;Tin為注入溫度,K。k0對應pm=0 MPa,pc= 0 MPa,Ts=0 K,Tin=0 K時的滲透率。不過這里不需要進行換算,實時滲透率k與k0之間是指數關系,使得k0可以認為是任意狀態時的滲透率,k計算時僅需對pm、pc、Ts、Tin代入兩狀態間的增量即可。

1.2 EGS熱開采的傳熱效能指標

生產井出口水溫是評估EGS出力和壽命的重要指標之一[4,5,14]。出口平均水溫定義為式(9):

(9)

式中:v為法向流速;T為溫度;Γ為出口邊界。

還有一些文獻[13,15,16]定義凈熱提取率來討論和評估不同輸入參數下熱提取效果。凈熱提取率是流經出口和入口的單位時間能量之差,通過式(10)得到:

Enet=Qouthout-Qinhin=cf(ρfoutvoutTout-ρfinvinTin)

(10)

式中:Qout、Qin分別為流出端和流入端工質的法向實時質量流率;hout、hin分別為流出端和流入端工質的實時比焓;ρfout、ρfin分別為流出端和流入端工質的實時流體密度;cf為流體工質的比熱;Tout、Tin分別為流出端和流入端工質的溫度。兩種方法中,后者要實現更精確的計算還應當在THM耦合模擬中考慮比焓、相變等的變化,這里忽略相變,僅用比熱和溫度的乘積代替比焓。

1.3 THM耦合的程序實現

通過對控制方程組進行弱形式的推導,可以建立斜壓流體THM完全耦合有限元模型,并進一步建立有限元程序。參考Cui和Wong[12],本文建立了基于等效多孔介質模型的HM-T的雙向耦合求解有限元程序,如圖1所示。

圖1 HM-T的雙向耦合求解有限元流程圖

2 算例分析與討論

2.1 計算模型及參數設置

本文使用的EGS模型基于Cui和Wong[12]的地熱開采模型,見圖2,原模型從五點注采模型中取了1/4進行了三維模擬。由于本文的有限元程序目前僅針對二維單元進行了開發,還未拓展到三維,本章在Cui等地熱模型的基礎上取A-A截面進行模擬。假設儲層已經通過水力激勵充分激發裂隙縫網,可以視作均勻等效連續介質處理。注入井流體出口位于儲層左下角,而生產井流體入口位于儲層右上角。

圖2 五點注采模型

這里對EGS的模擬分兩步進行,第一步僅對變形場和滲流場進行穩態計算,第二步以UP-T的雙向耦合求解方式對儲層模型進行持續40年的瞬態模擬計算。兩步模擬的邊界條件設置如圖3。其中上邊界σy=-25 MPa,對應地層深度約1 km,右邊界σx=-37.5 MPa,即1.5倍豎直載荷。對于滲流場和溫度場,未特殊說明的邊界分別為不透水邊界和絕熱邊界。假設以水為傳熱工質,但僅以20℃、0 MPa下的設置水流密度為1 000 kg/m3,其他輸入參數見表1。單元設置為矩形單元,總計80(寬)×100(高)個。本文參考Cui和Wong[12],將瞬態模擬時間步長設計成序列形式,前100步的步長為0.01年,后390步的步長為0.1年。

圖3 EGS儲層模擬分兩步進行的初邊值條件

表1 地熱開采模型中蓋層和儲層的輸入參數

2.2 主要變量的時空演化特征

本文同時進行了考慮斜壓流體密度演化和不考慮密度演化的數值模擬。不考慮流體密度演化時,認為流體密度為常數1 000 kg/m3。兩種情況下主要變量的等值線形態相似。為避免贅余,在此僅繪制考慮斜壓流體密度演化的模擬數據(0.01、1、10、20、30、40年后)進行觀察,如圖4所示。根據圖4可以分析得到孔隙壓力、溫度、豎直位移、流體密度的時空演化特點如下:

圖4 考慮流體密度演化時,40年內EGS模型的豎直位移v、孔壓p、溫度T、流體密度ρ演化

(1) 對于孔隙壓力,在儲層出口處的低壓使得流體流出,區域內孔隙壓力下降。周圍流體向出口附近補充,來自滲透率較大的儲層流體補充速度要比滲透率低的上蓋層快得多;隨著EGS運行,上蓋層也有發生向出口的少量流動,使得孔隙壓力逐漸演變,直至獲得能夠協調上蓋層和儲層不同滲透率的孔隙壓力梯度分布。而在儲層入口處,流體注入使得孔隙壓力增大。總體來說,除了進出口處孔隙壓力等值線稍微扭曲外,孔隙壓力在空間分布上基本表現為沿豎直方向由下到上遞減。

(2) 溫度分布體現了對流占優型問題的明顯特點,即存在完全冷卻區且完全冷卻區逐漸向出口擴張。完全冷卻區逐步向出口擴散,其中上下蓋層的滲透率較低,限制了從儲層向上下蓋層的流動,因此完全冷卻區的擴散方向從入口處從水平向右,逐漸轉為斜向上,最終指向出口。

(3) 對于豎向位移而言,儲層和上蓋層的豎向位移都是負值,即熱開采導致了儲層和上蓋層的下沉。域內豎向位移的演化以完全冷卻區的邊界為界,完全冷卻區兩側距離冷卻區邊界較遠處,其局部變形幾乎不受以溫度梯度表征的熱應力影響。完全冷卻區的邊界上,豎直位移等值線沿垂直于邊界的方向扭曲。從地表上看,注入井和生產井的沉降量通常比其他區域多,開采初期,生產井由于需要提供低壓而下沉顯著;開采過程中,注入冷水導致注入井顯著下沉,隨開采進行,沉降逐漸從注入井擴散;從完全冷卻區到達出口處(約第30年)開始,生產井的沉降量反而更顯著。

(4) 對于流體密度而言,盡管依據斜壓流體密度模型,溫度和孔隙壓力都會引起流體密度的變化,但是溫度對流體密度的影響更大。具體體現在流體密度的等值線圖中溫度等值線圖形狀相同、趨勢相反。在完全冷卻區,流體溫度較低,密度較大;在完全冷卻區外,流體在流動過程中隨著從巖層吸收熱量其溫度的升高,其密度也隨之減小。總體來說,在密度演化模型下,實時流體密度比給定的常數密度(1 000 kg/m3)小。

2.3 流體密度演化對EGS系統的影響

在斜壓流體模型中,流體密度的變化對熱流固三個物理場的變量產生了影響。為探究這一影響規律,這里進一步繪制觀測點B(見圖2(c))處幾個主要變量的變化如圖5。基于孔隙壓力的演化特點,可以將演化過程分為四個階段。

(1) 第一階段,即開始注入后0.01年內左右,由于此時出口處的孔隙壓力條件比域內其他位置都小得多,較低的出口壓力驅動流體克服重力向出口流動,域內的孔隙壓力因此迅速下降。而此時完全冷卻區還遠在流體入口附近,變形場中有效應力僅隨孔隙壓力下降而增大,變形量增大,表現為水平和豎向位移向y軸負方向下降。考慮流體密度演化時,流體的孔隙壓力下降的幅度比不考慮時大。這是因為考慮流體密度演化使得從B所在高度到出口處的流體實時密度比常數流體密度1 000 kg/m3小,維持從B到出口處的流動所需的壓差比不考慮流體密度演化時小。

(2) 第二階段,即完全冷卻區到達B之前的階段,從圖5(c)和圖5(d)可得這一時間大約在開始注入后的第0.1年至第12年。完全冷卻區從流體入口處逐漸擴散導致儲層溫度降低。冷卻區邊界上的熱應力作用使儲層向流體入口處收縮,表現為水平位移和豎直位移的進一步減小。考慮流體密度演化與否對此時儲層收縮的影響不大,且流體密度演化對完全冷卻區的到達時間也沒有顯著影響。

圖5 考慮流體密度演化與否對模擬結果的影響

(3) 第三階段,即完全冷卻區從B擴張到出口處的階段,結合圖4和圖5(d)可知這一時間大約在開始注入后的第12年至第30年。對于孔隙壓力而言,考慮流體密度演化時,從在這一階段開始,B所在高度以上的流體開始逐漸冷卻而密度增大,從而維持到出口處的流動所需的壓差增大,B處的孔隙壓力開始增大。這一階段還包含了B的冷卻。從圖5(d)可知B的冷卻時間大約在始注入后的第12年至第19年。考慮流體密度演化時,B的冷卻速度比不考慮密度演化時的慢一些。

(4) 第四階段,即完全冷卻區越過出口處所在高度后的擴張階段,這一時間開始于注入30年之后,此時B至出口區域的流體密度基本穩定,故孔隙壓力已不再受冷卻區擴張的影響。

考慮密度演化與否對傳熱性能指標的影響見圖6。如圖6(a),盡管兩者差距不太明顯,但考慮密度變化比不考慮時延緩了出口溫度的冷卻。如果儲層管理時設計某一標準如200℃作為合格的出口溫度,考慮流體密度變化時儲層壽命的模擬評估將比不考慮時延長近1年。對于凈熱提取率而言,如圖6(b),開采的第一階段EGS由于出口流速大而具有較高的凈熱提取率,進入穩定開采期后凈熱提取率演化主要受溫度影響,凈熱提取率考慮密度演化時模擬評估凈熱提取率比不考慮時提高4.09%。

圖6 考慮密度演化與否對傳熱效能的影響

2.4 流體密度和儲層滲透率耦合演化對EGS的影響

這里初始狀態對應第一步穩態計算后的狀態,認為此時具有初始滲透率k0=2.5×10-13m2,此后儲層滲透率因熱流固耦合效應發生演化。由于本身就比較小,上下蓋層的滲透率演化可以忽略。繪制40年內熱儲層的滲透率k演化如圖7,其他主要變量的演化見圖8。從圖7和圖8可以看出,滲透率k與溫度的等值線圖形狀比較相似,即在熱開采過程中,滲透率k受溫度的影響最大。以完全冷卻區為界,區內滲透率高,區外滲透率低。但在流體入口和出口處,由于局部壓力梯度大,孔隙壓力對滲透率的影響也比較明顯,表現為入口處滲透率高而出口處滲透率低。但總體來說,由于前0.01年域內孔隙壓力迅速降低,導致滲透率普遍減小,因此滲透率在演化條件下時,其大小比給定的常數滲透率小。

圖7 EGS開采40年內熱儲層的滲透率k演化

考慮滲透率與否對傳熱性能指標的影響見圖9。如圖9(a),考慮滲透率演化時,出口水溫的開始冷卻時間提前1年左右。滲透率演化加快了出口平均溫度的冷卻速度。如圖9(b),進入穩定開采階段后,考慮滲透率演化時模擬所得的凈熱提取率比不考慮時提高了24.82%。

圖9 考慮滲透率演化與否對傳熱效能的影響

3 結 論

(1) 隨著EGS開采的進行,THM主要變量的演化如下:地層發生沉降,注入井和生產井的沉降量通常比其他區域多;孔隙壓力的演化相比其他變量較快達到穩定,除了進出口處孔隙壓力等值線稍微扭曲外,孔隙壓力在空間分布上基本表現為沿豎直方向由下到上遞減;溫度分布體現了對流占優型問題的明顯特點,即存在完全冷卻區且完全冷卻區逐漸向出口擴張。

(2) 考慮流體密度演化時,溫度對流體密度演化的影響比孔隙壓力的影響更大。在對比考慮流體密度演化與否對系統平均出口溫度和凈熱提取率的影響時發現,考慮流體密度演化比不考慮時儲層出口溫度的開始冷卻時間差別不大,但延緩了出口溫度的冷卻速度,而凈熱提取率比不考慮流體密度演化時提高約4.09%。

(3) 考慮流體密度和滲透率耦合演化時,溫度對滲透率演化的影響比孔隙壓力、圍壓等其他因素大。在對比考慮滲透率演化與否對系統平均出口溫度和凈熱提取率的影響時發現,考慮滲透率演化時出口水溫的開始冷卻時間提前1年左右,加速了出口溫度的冷卻,而凈熱提取率比不考慮滲透率演化時提高約24.82%。