熱補償條件下雙井EGS產(chǎn)能和壽命預(yù)測方法研究

張 霄，孫玉學(xué)，張慶松，李相輝，尹占超，李 壯，周 政

(山東大學(xué) 巖土與結(jié)構(gòu)工程研究中心，山東 濟南 250061)

近年來，化石資源臨近枯竭，環(huán)境問題日益突出[1]。干熱巖資源因儲量巨大、清潔環(huán)保，決定了其必然會成為未來重要能源之一[2]。增強型地熱系統(tǒng)（EGS）作為開發(fā)干熱巖資源的有效手段，受到了廣泛關(guān)注[3-5]。

研究EGS取熱過程時，若利用建造干熱巖電站的方法獲得實際工程的地質(zhì)資料、工程經(jīng)驗和數(shù)據(jù)，對研究干熱巖儲層熱提取極為有利，但目前的干熱巖儲層改造技術(shù)尚不完善，投資巨大、周期長。因此，現(xiàn)階段研究大多基于數(shù)值模擬方法建立數(shù)學(xué)模型，利用多場耦合過程對EGS取熱過程進行評估，指導(dǎo)EGS的設(shè)計與優(yōu)化[6-8]。在熱儲改造和EGS運行過程中，熱儲巖體的物理化學(xué)性質(zhì)都會隨時間和空間不斷變化[9-10]，改造后熱儲尺寸多以km為量級，內(nèi)部裂隙通常以mm為量級，在利用精細化網(wǎng)格直接模擬取熱過程時計算量巨大[11]，需要對熱儲內(nèi)的裂隙進行簡化。目前，常用的簡化模型有單孔隙模型[12]、雙孔隙模型[13-14]和多孔隙模型。其中，雙孔隙模型和多孔隙模型因精度要求較高，需要大量試驗數(shù)據(jù)支撐，在實際中應(yīng)用困難[15]。

干熱巖熱量資源預(yù)測方面，目前主要有兩種預(yù)測方式：一種是利用目標熱儲體積及熱儲溫度和某一確定溫度（如地表溫度）的差值計算，按百分比的可開采量預(yù)測[16-19]，主要有體積法、蒙特卡羅法及地表流量法等[20-23]，但這些方法在應(yīng)用于EGS產(chǎn)能預(yù)測時，熱儲開采前后溫差的設(shè)定值偏大，導(dǎo)致預(yù)測不準確；另一種是數(shù)值法[24]，根據(jù)詳細的熱儲參數(shù)計算和評價勘查程度較高、具有一定開采歷史的地熱田[23]，此種方法較為繁瑣，計算周期長，需要詳細的熱儲參數(shù)，不適合未開采的熱儲產(chǎn)能的預(yù)測。在EGS運行系統(tǒng)可控因素的研究方面，目前的研究多集中于以不同的熱媒介質(zhì)注入溫度、不同的熱媒介質(zhì)注入速度及不同井間距等人為可控因素控制EGS的運行[8,25]。以注水速度控制時，隨著EGS的運行，注入井的壓力不斷產(chǎn)生變化，大大增加了EGS運行時的不可控性[26]。EGS取熱效果方面，大多數(shù)研究將熱儲視為單孔隙模型，針對熱儲改造后的孔隙比、滲透率、注入井溫度、布井方式等多種因素分析EGS取熱效果[27-28]。如：王昌龍[29]考慮了11種因素對EGS取熱效果的影響。在這些研究中，大多局限于單因素的定性研究，較少考慮多因素組合的定量研究，且鮮有將大地熱流的熱補償作為初始條件加入到模型之中[30-33]。有研究發(fā)現(xiàn)，熱儲周圍的大地熱流對熱儲層長期積累的熱補償熱量會對EGS運行壽命產(chǎn)生不可忽略的影響[34]。

以吸放熱公式[35]和Dupuit公式為基礎(chǔ)，建立EGS產(chǎn)能和壽命控制方程；依據(jù)產(chǎn)能和熱儲參數(shù)，通過調(diào)節(jié)熱媒介質(zhì)注入溫度和注采壓差兩個可控因素取值預(yù)測EGS壽命，增加EGS運行時的可控性。同時，根據(jù)地下水流動方程和熱傳導(dǎo)方程分析出熱儲初始溫度、熱媒介質(zhì)注入溫度、EGS壽命、熱儲體積和熱儲比表面積5個因素對發(fā)生熱突破（生產(chǎn)井產(chǎn)出溫度低于423.16 K[36]）時對控制方程中熱儲減少熱量、大地熱補償熱量、生產(chǎn)井平均產(chǎn)出溫度和熱儲形狀系數(shù)4個未知參數(shù)的影響。在考慮大地熱補償?shù)那闆r下，利用雙井EGS的3維數(shù)值模型，對不同工況的EGS在發(fā)生熱突破時的4個未知參數(shù)的變化規(guī)律進行研究，得到4個未知參數(shù)隨5個影響因素變化的具體預(yù)測公式，實現(xiàn)EGS產(chǎn)能預(yù)測和多因素定量化分析。

1 流熱耦合下EGS換熱過程分析

如圖1所示，假設(shè)熱儲周圍圍巖滲透率為0，熱媒介質(zhì)在流動過程中無熱量損失，人工儲層原有熱量Qm和大地補償?shù)臒崃縌t完全由熱媒介質(zhì)攜帶并提取，其中，熱媒介質(zhì)攜帶熱量由Qw表示。

圖1 EGS換熱過程的流熱耦合機制Fig.1 Flow-heat coupling mechanism of EGS heat transfer process

根據(jù)熱量平衡，由吸放熱公式（式（1））結(jié)合Dupuit公式（式（2））可得熱媒介質(zhì)攜帶熱量（EGS產(chǎn)出熱量），如式（3）所示：

式（1）～（3）中：t為EGS壽命；k為熱儲滲透率；mw為熱媒介質(zhì)質(zhì)量； ρw為熱媒介質(zhì)密度；L為熱媒介質(zhì)路徑；A1為熱媒介質(zhì)在熱儲中流過的截面面積；cw為熱媒介質(zhì)的熱容；vw為熱媒介質(zhì)平均流速；P為注采壓差；ΔT1為注采溫差，ΔT1=Tout-Tin。

若熱儲在tn時刻發(fā)生熱突破，tm時刻的產(chǎn)出溫度Tout(tm)可作為Tout(t)在t0至tn時刻的平均值，即：

整理式（3）和（4），可得EGS產(chǎn)能和壽命控制方程：

利用達西定律和多孔介質(zhì)傳熱方程描述流熱耦合過程，如式（6）～（9）所示：

連續(xù)性方程：

流體動量方程：

流體能量方程：

周圍圍巖熱補償能量方程：

牛頓冷卻公式：

式（5）～（10）中，t為熱媒介質(zhì)流動時間，xi為熱媒介質(zhì)在i方向上的位移，ρw、ρm分別為流體、圍巖的熱儲密度，vi為熱儲中流體的流速，F(xiàn)i為體積力， ?2為拉普拉斯算子，cw、cm為流體和圍巖比熱容，A2為熱儲與圍巖接觸面積，η為動力黏度，T為熱儲溫度場，T∞為圍巖初始溫度，λw、λp分別為流體和圍巖的導(dǎo)熱系數(shù)，h為表面?zhèn)鳠嵯禂?shù)，Φ為熱流量，A3為流體與固體接觸面積，Ts、Tw分別為固體和流體溫度。

利用式（6）～（9），結(jié)合牛頓冷卻公式（10），對式（5）中的4個參數(shù)進行分析，可確定出影響4個參數(shù)的主要因素，如表1所示。

表1 4個未知變量的影響因素匯總Tab. 1 Summary of influencing factors of four unknown parameters

Qm由熱儲終溫Tml和熱儲初始溫度Tm控制，所以在預(yù)測時僅需考慮Tml即可，整理后如式（11）所示。

若利用式（5）達到討論Tin和P與t之間關(guān)系的目的，可利用數(shù)值模擬預(yù)測4個參數(shù)的具體表達式，并代入式（5）中得到Tin、P、t三者關(guān)系。

2 EGS取熱模型概述

2.1 地質(zhì)模型構(gòu)建

如圖2所示，EGS地下部分由人工熱儲、熱儲周圍的圍巖、注入井和生產(chǎn)井組成。為了更準確地模擬EGS地下熱交換過程，在模型中作出主要假設(shè)：1）熱儲裂隙小而均勻，可看做單孔隙率的多孔介質(zhì)；2）熱媒介質(zhì)流動為單相流；3）熱儲中熱媒介質(zhì)流動是飽和層流流動（0.003 m/s以下），初始時刻熱儲層中充滿了與周圍巖石溫度相同的熱媒介質(zhì)；4）周圍圍巖滲透率為0，熱媒介質(zhì)在流動過程中無損失；5）熱儲及周圍圍巖不會產(chǎn)生形變和各種化學(xué)反應(yīng)。

圖2 EGS地下部分概念圖Fig.2 Concept picture of EGS underground part

圖3 雙井EGS模型3維圖Fig.3 Three-dimensional pictures of the double well EGS model

如圖3所示，熱儲周圍包覆200 m厚的花崗巖圍巖，在外表面施加以3 K/100 m隨深度線性增加的溫度邊界，模擬大地熱流給予熱儲的熱量補給；分別建立不同邊長的正六面體，研究體積對4個參數(shù)的影響，如圖3（a）所示；分別建立同體積下不同形狀的熱儲，研究比表面積對4個參數(shù)的影響，如圖3（b）所示；其余模型邊界條件如表2所示。模型采用“下注上采”的開采方式，最大模擬時間50 a，采用瞬態(tài)求解方式，計算結(jié)果每年采集一次，且最終結(jié)果取符合EGS運行25 a，生產(chǎn)井出水速率達到80 kg/s以上，保證EGS商業(yè)化利用[35]。

表2 模型邊界條件取值Tab. 2 Values of model boundary conditions

2.2 模型熱物性參數(shù)取值

模型其他熱物性參數(shù)如表3所示[30-31]。

表3 模型熱物性參數(shù)取值Tab. 3 Values of model thermal property parameters

2.3 網(wǎng)格獨立性驗證

本文主要對計算結(jié)果中的生產(chǎn)井平均產(chǎn)出溫度和平均速度、熱儲平均溫度和熱儲中流體平均速度進行獨立性檢驗。

進行網(wǎng)格無關(guān)性檢驗時，按時間步長0.01 s，對不同自由四面體的網(wǎng)格數(shù)量進行無關(guān)性檢驗。采用熱儲邊長為500 m正六面體，熱儲初始溫度-熱媒介質(zhì)注入溫度-注采壓差為565.66 K-353.16 K-20 MPa的工況，分別選取第15年和第30年的計算結(jié)果進行比較（圖4）。由圖4可知：網(wǎng)格數(shù)量對計算結(jié)果有著較大的影響，當網(wǎng)格數(shù)量小于25萬個時，計算結(jié)果存在較大差異；當網(wǎng)格數(shù)量大于30萬個時，模擬結(jié)果變化不大。

圖4 網(wǎng)格無關(guān)性檢驗結(jié)果Fig.4 Results of the grid independence verification

進行時間獨立性檢驗時，以第12種（36.76萬個）網(wǎng)格數(shù)量為網(wǎng)格模型，時間步長為0.007 5、0.010 0、0.012 5和0.015 0 s，其余條件與網(wǎng)格獨立性檢驗時相同，結(jié)果如表4所示。由表4可知，當網(wǎng)格數(shù)量達到30萬以后，時間步長在0.012 5 s以下時對模擬結(jié)果影響不大。

由圖4和表4可知，在計算時，最大網(wǎng)格取40 m，最小網(wǎng)格取0.3 m，網(wǎng)格增長率取2.4，單元格數(shù)量41.21萬個，時間步長選擇0.01 s較為合理。

表4 時間獨立性檢驗結(jié)果Tab. 4 Results of the time independence verification

3 結(jié)果分析

圖5 Tm、Tin、V、SSa、t對Tml、Qt、、 α等參數(shù)的影響結(jié)果Fig.5 Influences of five factors (Tm、Tin、V、SSa、t) to four parameters (Tml、Qt、 Tout、 α)

由圖5（a）、（c）、（d）、（e）可知：當Tm和Tin增加時，Tml逐漸上升，且Tm對Tml的影響要低于Tin對Tml的影響；隨著SSa的增加，Tml也會逐漸升高；隨著V的增大，Tml表現(xiàn)出逐漸減小的趨勢；以Tm=565.66 K、Tin=293.16 K為例，EGS不同壽命所對應(yīng)的Tml相差不大，近似服從高斯分布，與時間項無關(guān)。由于熱突破條件并非與熱媒介質(zhì)注入溫度相同，所以在注入井與生產(chǎn)井之間必定存在溫度梯度，在均勻介質(zhì)滲流場中，兩井之間必存在流速最快的最優(yōu)通路。在最優(yōu)通路間，溫度梯度介于注入溫度和熱突破溫度之間，在最優(yōu)通路外其溫度梯度上限與熱儲初始溫度有關(guān)。因此，隨著Tm的升高，最優(yōu)通路外的熱儲溫度升高；隨著Tin的升高，整個熱儲溫度梯度下限升高；隨著SSa的增加，熱儲與圍巖接觸面積增加，熱補償?shù)臒崃吭龆啵紩?dǎo)致Tml的升高。此外，其他條件相同時，V的增加使得熱突破時低溫區(qū)域增大，Tml逐漸減小。

由圖5（f）、（g）、（h） 可知：受圍巖熱補償?shù)挠绊懀?5年后 ，Qt隨著EGS運行時間逐漸增加；隨著Tm、V和SSa的增大，Qt也逐漸增加；對于Tin對Qt的影響，由圖5（f）分析可以發(fā)現(xiàn)，除熱儲初始溫度處于485.66 K時，Qt隨熱媒介質(zhì)注入溫度逐漸增加，其余工況均無明顯規(guī)律。分別取第27年和第29年Qt的值，利用極差分析法將Tm與Tin對Qt的影響進行比較，如表5和6所示，發(fā)現(xiàn)熱儲初始溫度Tm對Qt的影響要遠遠大于Tin對Qt的影響，可近似將Tin對Qt的影響忽略。因此，可以確定，Qt主要受到Tm、Tin、V和SSa的影響。

表5 第27年Tm與Tin極差分析Tab. 5 Range analysis of Tm and Tin in the 27th year

表6 第29年Tm與Tin極差分析Tab. 6 Range analysis of Tm and Tin in the 29th year

由圖5（g）、（h）可知：隨著V的增大，α表現(xiàn)出先減小后增大的趨勢，其最低點大致出現(xiàn)在邊長為400 m的正六面體熱儲附近；隨著SSa的增大，α表現(xiàn)出先增大后減小的趨勢，其最高點出現(xiàn)在正六面體熱儲附近。

經(jīng)上述分析，將表1中的影響因素予以修正，如表7所示，同時將式（11）修正為式（12）：

表7 修正后的Tml、Qt、、α影響因素匯總Tab. 7 Summary of revised influencing factors of Tml, Qt,,α

表7 修正后的Tml、Qt、、α影響因素匯總Tab. 7 Summary of revised influencing factors of Tml, Qt,,α

參數(shù) 影響因素Tml Tm、Tin、V、SSa Qt Tm、V、SSa、t Tout Tm、Tin、V、SSa α V、SSa

為得到式（12）的具體表現(xiàn)形式，需要對上述計算結(jié)果繼續(xù)分析。

3.1 熱儲終溫Tml 變化規(guī)律

為探究表7中4個因素與Tml的定量關(guān)系，以Tin為353.16 K時為基準，確定出Tml353.16隨Tm的變化；然后，基于此分別討論Tin、SSa與V對Tml的增量的變化，如圖6所示；最后，將上述4項進行線性疊加，如式（13）所示：

如圖6（a）所示，極易確定出Tml353.16隨Tm變化公式。如圖6（b）所示，不同Tm的ΔTmlTin隨著Tin呈一致的線性負相關(guān)關(guān)系，因此，ΔTmlTin不受Tm的影響。取6種不同Tm的ΔTmlTin的平均值，得到ΔTmlTin隨Tin的變化公式，如圖6（c）所示。由于比表面積在不同體積時難以一致表示，因此，預(yù)測比表面積引起的熱儲終溫增量ΔTmlS時，應(yīng)以在同體積下與比表面積變化規(guī)律一致的量—各表面積與相應(yīng)體積下正六面體表面積的比值S/SR表示，以消除熱儲體積對ΔTmlS的影響。如圖6（d）所示，得到了ΔTmlS隨S/SR的變化規(guī)律。考慮ΔTmlV隨不同體積熱儲與邊長為500 m的正六面體熱儲體積比值V/V500的變化規(guī)律，用以量化體積對Tml的影響，如圖6（e）所示。所以得到式（13）中各分項的具體表達式為：

圖6 熱儲終溫Tml計算圖Fig.6 Calculation of final storage temperature (Tml)

3.2 大地熱補償熱量 Qt變化規(guī)律

首先，得到前25年的補償熱量隨熱儲初始溫度Tm的變化量Qt25，基于此分別討論Qt隨t和SSa的增量變化ΔQtt和ΔQtS，將上述3個量進行線性疊加。如圖5（g）所示，當熱儲體積為0時，Qt也一定為0。令不同體積得到的補償熱量與邊長為500 m的正六面體熱儲體積得到的補償熱量的比值為ηV，將ηV作為體積系數(shù)與Qt25、ΔQtt和ΔQtS疊加得到的補償熱量相乘，以考慮體積因素對Qt的影響，如式（15）所示：

Qt25與Tm的關(guān)系、ΔQtS與S/SR的關(guān)系及ηV與V/V500的關(guān)系如圖7（a）、（b）和（d）所示。由圖7（e）可知：ΔQtt隨著時間的延長不斷增加，近似呈線性變化；但不同Tm對應(yīng)的ΔQtt的增量不同，且隨著Tm的上升，ΔQtt隨著時間變化的線性曲線斜率不斷增加，截距不斷減小。故在考慮ΔQtt隨著時間的變化時，還要考慮Tm的影響。分別計算出斜率c和截距d隨Tm的變化，如圖7（c）所示。最終得到式（15）中各分項的具體表達式為：

圖7 大地熱補償熱量Qt計算Fig.7 Calculation of geothermal compensation heat (Qt)

3.3 平均產(chǎn)出溫度 Tout變化規(guī)律

式中，Tout353.16、ΔToutS和ΔToutV極易通過數(shù)值計算結(jié)果通過擬合曲線得到，如圖8（a）、（c）和（d）。如圖8（e）所示，隨著Tin的減小，ΔToutTin逐漸增加，可近似看做線性曲線，且不同Tm的ΔToutTin在同一注入溫度的增量不同，隨著注入溫度的上升而不斷增加。因此，在計算ΔToutTin時，需要考慮Tm的影響。與ΔQtt計算方法類似，如圖8（b）所示，分別計算出ΔToutTin隨Tin變化線性曲線的斜率a和截距b隨Tm的變化，即可得到ΔToutTin的計算值，則式（17）中各分項的具體表達式為：

圖8 平均產(chǎn)出溫度 計算Fig.8 Calculation of average output temperature

3.4 系數(shù) α的變化規(guī)律

由表7可知，系數(shù)α只與V和SSa有關(guān)。為得到三者間的定量關(guān)系，將α分為αV和ΔαS。首先，得到αV隨V的變化，如圖9（a）所示；然后，基于αV計算由SSa引起的增量ΔαS，如圖9（b）所示；最后，將αV與ΔαS進行線性疊加得到系數(shù)α的計算式，如式（19）所示：

如圖9所示，最終得到式（19）中各分項的具體表達式為：

圖9 系數(shù) α的計算圖Fig.9 Calculation of coefficientsα

式（21）中各分式可由式（14）、（16）、（18）、（20）表示。將式（21）代入式（5）中，結(jié)合勘測到的熱儲滲透率即可得到Tin、P、t三者關(guān)系，指導(dǎo)EGS建設(shè)。

3.5 計算方程的適用性驗證

目前，大多數(shù)研究主要利用工質(zhì)循環(huán)流量表示熱媒介質(zhì)注入的量及熱媒介質(zhì)產(chǎn)出量，利用注采壓差P與已有研究比較來驗證計算公式適用性的方法難以直接實現(xiàn)。根據(jù)Dupuit公式可知，通過熱儲尺寸和滲透率可將注采壓差P換算為熱媒介質(zhì)在生產(chǎn)井的平均流速，由生產(chǎn)井的尺寸即可得到熱媒介質(zhì)產(chǎn)出量，即工質(zhì)循環(huán)流量。為驗證計算方程的可靠性，將其應(yīng)用至文獻[15]的模型中。將文獻[15]中的熱儲尺寸、熱儲初始溫度、熱媒介質(zhì)注入溫度和不同EGS運行時間代入式（5）中，得到不同EGS運行時間對應(yīng)的注采壓差P；然后，僅使用 COMSOL Multiphysics中達西定律模塊計算得到工質(zhì)循環(huán)流量，所得結(jié)果如圖10所示。研究發(fā)現(xiàn)，利用式（5）計算結(jié)果與文獻[15]所得曲線結(jié)果一致，說明式（5）具有良好的適用性。

圖10 計算公式適用性驗證Fig.10 Applicability verification of calculation formula

4 結(jié) 論

利用理論分析與數(shù)值模擬相結(jié)合，描述EGS地下取熱過程，主要結(jié)論如下：

1）將吸放熱公式與Dupuit公式相結(jié)合，得到了雙井EGS產(chǎn)能與壽命控制方程，揭示了EGS產(chǎn)能隨熱儲滲透率、EGS壽命、熱媒介質(zhì)注入溫度、注采壓差和熱儲形狀系數(shù)的變化關(guān)系。

3）將控制方程和預(yù)測公式與現(xiàn)有研究進行對比，控制方程和預(yù)測公式具有良好的適用性，可有效指導(dǎo)EGS建設(shè)，加強EGS運行期間的實時調(diào)控。