隨機(jī)裂隙網(wǎng)絡(luò)儲(chǔ)層與井筒熱流耦合數(shù)值模擬

單丹丹 閆 鐵 李 瑋 孫士慧 逯廣東 趙 歡

1.東北石油大學(xué)石油工程學(xué)院 2.提高油氣采收率教育部重點(diǎn)實(shí)驗(yàn)室·東北石油大學(xué)

0 引言

干熱巖是一種深層熱巖體，其資源儲(chǔ)量極其豐富[1-2]。增強(qiáng)型地?zé)嵯到y(tǒng)（Enhanced Geothermal System，EGS）是目前開(kāi)采干熱巖資源最有效的技術(shù)手段[3-4]。EGS的設(shè)計(jì)與優(yōu)化主要依據(jù)兩方面重要因素——出力與壽命[5]。對(duì)系統(tǒng)出力與壽命進(jìn)行評(píng)價(jià)，要反復(fù)模擬EGS熱流耦合過(guò)程。針對(duì)這一耦合過(guò)程，建立相應(yīng)的數(shù)學(xué)模型，并進(jìn)行井筒—隨機(jī)裂隙網(wǎng)絡(luò)儲(chǔ)層耦合數(shù)值模擬，可以預(yù)測(cè)EGS運(yùn)行狀態(tài)、利用效率與使用壽命，對(duì)系統(tǒng)進(jìn)行采熱評(píng)價(jià)[6]。對(duì)于熱儲(chǔ)內(nèi)裂隙巖體熱流耦合的數(shù)學(xué)模型可分為等效連續(xù)介質(zhì)模型與離散裂隙網(wǎng)絡(luò)模型兩大類(lèi)[7]，后者不僅與現(xiàn)實(shí)更接近，且能更好地模擬水熱遷移過(guò)程，因而受到廣泛的關(guān)注[8-14]。顯式模擬裂隙的離散裂隙網(wǎng)絡(luò)模型雖接近于真實(shí)情況但計(jì)算量較大[15]。由于裂隙寬度方向相對(duì)于其他兩個(gè)方向，其尺度小到可以忽略，Juanes等[16]提出在描述裂隙的過(guò)程中可以采用低一維的單元來(lái)進(jìn)行。在裂隙網(wǎng)絡(luò)生成上，大部分研究都依蒙特卡羅方法隨機(jī)生成[17]。關(guān)于隨機(jī)裂隙網(wǎng)絡(luò)儲(chǔ)層的熱流耦合研究很少考慮到井筒部分的熱損失，而在對(duì)井筒—熱儲(chǔ)耦合模擬的研究中，大多是將儲(chǔ)層看作是單孔隙等效多孔介質(zhì)模型。Jiang等[18]建立了井筒與熱儲(chǔ)循環(huán)通道熱流耦合瞬態(tài)模型，其中儲(chǔ)層被視為單一孔隙度的等效多孔介質(zhì)。Zeng等[19]提取了DP23-1井的地質(zhì)資料，研究了干熱巖儲(chǔ)層的熱潛能，其儲(chǔ)層也被看作是等效多孔介質(zhì)。曹文炅等[20]在局部非熱平衡假設(shè)的基礎(chǔ)上，用理想的垂直井模型模擬了EGS的采熱能力，熱儲(chǔ)部分同樣以多孔介質(zhì)來(lái)代替，忽略了裂隙的存在。長(zhǎng)深井筒的熱損失及隨機(jī)裂隙網(wǎng)絡(luò)儲(chǔ)層的熱流耦合過(guò)程在循環(huán)系統(tǒng)的整體采熱評(píng)價(jià)中都占據(jù)著主導(dǎo)作用，為此有必要開(kāi)展二者的熱流耦合模擬。由于COMSOL軟件在多物理場(chǎng)耦合與裂隙細(xì)微結(jié)構(gòu)、井筒細(xì)長(zhǎng)結(jié)構(gòu)的多尺度耦合上具有強(qiáng)大優(yōu)勢(shì)，故選用其進(jìn)行二維有限元數(shù)值模擬，可形象展現(xiàn)系統(tǒng)滲流傳熱過(guò)程，在此基礎(chǔ)上分析熱物性參數(shù)與幾何參數(shù)變化對(duì)裂隙面溫度、采出溫度、熱開(kāi)采速率等的影響，得出影響系統(tǒng)產(chǎn)能與壽命的各項(xiàng)因素，以期為相關(guān)領(lǐng)域設(shè)計(jì)與施工提供理論依據(jù)。

1 井筒—隨機(jī)裂隙網(wǎng)絡(luò)儲(chǔ)層熱流耦合數(shù)學(xué)模型

1.1 基本假定

人工熱儲(chǔ)可以簡(jiǎn)化為由基質(zhì)巖塊和裂隙組成的雙重介質(zhì)模型，基質(zhì)巖塊視為孔隙介質(zhì)，相比于裂隙，其滲透率極低，而人工壓裂形成的裂隙網(wǎng)絡(luò)才是熱儲(chǔ)層中主要的滲流通道[21]。由于高壓作用，液態(tài)水不可能發(fā)生氣化，導(dǎo)致系統(tǒng)內(nèi)的流動(dòng)為單相水流。井筒、開(kāi)孔與裂隙網(wǎng)絡(luò)構(gòu)成水循環(huán)的主要通道，在注入井與采出井井口壓力不變的條件下，這一通道內(nèi)的流體流動(dòng)相對(duì)于幾十年的系統(tǒng)運(yùn)行時(shí)間而言，很快會(huì)達(dá)到流速穩(wěn)定，因此在EGS模擬中，流動(dòng)過(guò)程是穩(wěn)態(tài)的而傳熱過(guò)程是瞬態(tài)的，這種近似對(duì)于系統(tǒng)的長(zhǎng)期運(yùn)行而言是合理的。對(duì)于穩(wěn)定的滲流系統(tǒng)，井筒、開(kāi)孔及裂隙網(wǎng)絡(luò)中流體的流速應(yīng)相同，均服從達(dá)西定律[22]。巖體的熱量主要通過(guò)對(duì)流換熱及熱傳導(dǎo)的形式向外傳遞。不考慮巖體及裂隙網(wǎng)絡(luò)的體積變化，認(rèn)為系統(tǒng)整體不存在力學(xué)場(chǎng)，僅對(duì)井筒與熱儲(chǔ)耦合的滲流場(chǎng)和溫度場(chǎng)進(jìn)行模擬。

1.2 數(shù)學(xué)模型

水在EGS中的流動(dòng)可分為4個(gè)性質(zhì)不同的子區(qū)域：①具有井筒壁的細(xì)長(zhǎng)流道—注入井井筒、采出井井筒，其與基巖之間只有熱量傳遞，無(wú)質(zhì)量傳遞；②無(wú)井筒壁的開(kāi)孔位置，視為多孔介質(zhì)，與熱儲(chǔ)裂隙及基巖既有傳熱過(guò)程，也有傳質(zhì)過(guò)程；③服從達(dá)西定律的人工熱儲(chǔ)內(nèi)裂隙網(wǎng)絡(luò)；④極低滲透性基質(zhì)巖塊（包括除井筒、開(kāi)孔、裂隙網(wǎng)絡(luò)外的一切地層巖石）。模型基于前述假定條件，流動(dòng)為單相水流，既不考慮循環(huán)水與巖石的化學(xué)作用及水的物性變化，也不考慮巖石熱應(yīng)變引起的裂隙網(wǎng)絡(luò)孔隙率及滲透率的變化。在以上假設(shè)基礎(chǔ)上，分別得出井筒與基巖之間的熱傳導(dǎo)模型、裂隙網(wǎng)絡(luò)內(nèi)滲流模型、裂隙網(wǎng)絡(luò)中的傳熱模型及開(kāi)孔、基巖滲流與傳熱模型。

井筒與基巖單位長(zhǎng)度的徑向熱傳導(dǎo)為：

式中α表示地層熱擴(kuò)散系數(shù)，m2/s；t表示加熱（或冷卻）時(shí)間，h；r表示井筒半徑，m。

裂隙網(wǎng)絡(luò)內(nèi)滲流平衡方程為：式中df表示裂隙寬度，m；εw表示裂隙的孔隙率，無(wú)量綱；ρw表示水的密度，kg/m3；Kf表示裂隙內(nèi)的滲透率，m2；μ表示水的動(dòng)力黏度，Pa·s；Qf表示基巖與裂隙面的流量交換；n表示裂隙面法向；表示沿裂隙α切向求導(dǎo)；p表示壓強(qiáng)場(chǎng)，為一矢量。

裂隙網(wǎng)絡(luò)溫度場(chǎng)方程為：

其中

式中cw表示水的比熱容，J/(kg·K)；Tf表示裂隙內(nèi)水的溫度，K；λf表示水的熱傳導(dǎo)系數(shù)，W/(m·K)；Wf表示裂隙表面水從基巖吸收的熱量，W/m2；h表示裂隙水與基巖邊界處的對(duì)流換熱系數(shù)，W/(m2·K)。

開(kāi)孔、基巖滲流場(chǎng)方程為：

創(chuàng)客是堅(jiān)守創(chuàng)新、堅(jiān)持實(shí)踐、樂(lè)于分享并且追求美好生活的一群人，是把興趣與愛(ài)好努力變成現(xiàn)實(shí)的人，是社會(huì)迎來(lái)新一輪的“科技社會(huì)化”浪潮，是一場(chǎng)快速由工業(yè)社會(huì)向信息社會(huì)過(guò)渡的運(yùn)動(dòng)。創(chuàng)客空間的普及發(fā)展，使分布式、數(shù)字化、個(gè)性化、定制化的電腦網(wǎng)絡(luò)制造方式取代傳統(tǒng)的工廠加工制作方式。作為未來(lái)人才培養(yǎng)基地的學(xué)校應(yīng)該培養(yǎng)更多的創(chuàng)客，打造“創(chuàng)客校園”。

其中

式中ε表示開(kāi)孔、基巖的孔隙率，無(wú)量綱；K表示開(kāi)孔、基巖的滲透率，m2；Qm表示滲流的源匯項(xiàng)。

開(kāi)孔、基巖溫度場(chǎng)控制方程為：

式中ceff表示開(kāi)孔、基巖的等效比熱容，J/(kg·K)；λeff表示開(kāi)孔、基巖的等效熱傳導(dǎo)系數(shù)，W/(m·K)；Tr表示開(kāi)孔、基巖的溫度，K。

2 數(shù)值模型與邊界條件

2.1 概念模型

EGS循環(huán)系統(tǒng)包括注入井、人工熱儲(chǔ)、采出井以及基質(zhì)巖塊4個(gè)部分，所建立的概念模型見(jiàn)圖1。

2.2 二維裂隙網(wǎng)絡(luò)的生成

圖1 EGS概念模型圖[23]

Dershowitz和Einstein[24]總結(jié)了幾種巖石節(jié)理體系且提出對(duì)其描述所需的參數(shù)。復(fù)雜的三維裂隙網(wǎng)絡(luò)不但會(huì)使運(yùn)算速度緩慢，且其依據(jù)參數(shù)隨機(jī)生成，和二維計(jì)算結(jié)果相差并不很大，因此采用二維隨機(jī)生成的裂隙網(wǎng)絡(luò)系統(tǒng)進(jìn)行模擬，所要用到的參數(shù)有裂隙的密度、跡長(zhǎng)、開(kāi)度及產(chǎn)狀。在裂隙的生成上以Monte-Carlo方法為基礎(chǔ)，擬定兩組方向裂隙，分別與水平方向呈30°和150°，且具明顯的各向異性。裂隙網(wǎng)絡(luò)中各條裂隙的跡長(zhǎng)服從正態(tài)分布，其長(zhǎng)度平均值為60 m，方差為20 m，裂隙的數(shù)目為800個(gè)，裂隙的范圍是 500 m×500 m。

2.3 有限元模型

借助數(shù)學(xué)軟件MATLAB進(jìn)行代碼編寫(xiě)并運(yùn)行得出裂隙網(wǎng)絡(luò)模型，之后將生成的腳本文件導(dǎo)入到AutoCAD中形成圖形文件，去掉伸出裂隙邊界的孤立段，再進(jìn)行井筒及周?chē)鶐r的繪制，完成建模過(guò)程，最終導(dǎo)入到COMSOL軟件中進(jìn)行模擬計(jì)算。依據(jù)上述方法建立的有限元模型如圖2所示。

圖2 有限元模型圖

有限元模型中，假定采用以下的初始條件及邊界條件。其中熱物性參數(shù)數(shù)據(jù)參見(jiàn)表1。模型的滲流場(chǎng)、溫度場(chǎng)及井筒壁（保溫層）邊界條件如下：

表1 熱物性參數(shù)表

1）滲流場(chǎng)邊界條件：注入井井口壓力，16 MPa；采出井井口壓力，10 MPa；熱儲(chǔ)外的花崗巖滲透率非常低，故取不透水邊界，即q=0。

2）溫度場(chǎng)邊界條件：注入井井口溫度為20 ℃，地溫梯度為50 ℃/km，熱流邊界條件為井筒壁（保溫層）與周?chē)貙拥膫鳠徇^(guò)程用多孔介質(zhì)傳熱模塊里的薄層來(lái)設(shè)置，厚度取0.01 m。

3）井筒壁（保溫層）邊界條件：在注、采井井筒左右兩邊添加線段，代表兩井的井筒壁及保溫層邊界，1、2、3、4代表注入井井筒左右兩邊的井筒壁及保溫層邊界，5、6、7、8則代表采出井井筒左右兩邊的井筒壁及保溫層邊界（圖3）。

3 模擬結(jié)果與分析

基于假定的熱物性參數(shù)與幾何參數(shù)，研究這些參數(shù)變化對(duì)采出溫度與熱開(kāi)采速率的影響。擬定的幾組算例如表2所示。

3.1 對(duì)算例1進(jìn)行模擬與分析

圖3 井筒壁（保溫層）示意圖

以算例1為參考條件，在給定的已知參數(shù)條件下對(duì)其進(jìn)行數(shù)值模擬，得出系統(tǒng)運(yùn)行5年、10年、15年、20年的溫度分布云圖（圖4）。可以看出，運(yùn)行初期注入井周?chē)貙佑捎跓醾鲗?dǎo)作用而降溫的范圍較小，即影響半徑較小，隨著運(yùn)行時(shí)間的增長(zhǎng)，被井筒內(nèi)水流改變的地層溫度場(chǎng)范圍逐漸增大。裂隙儲(chǔ)層內(nèi)的溫度場(chǎng)分布也有類(lèi)似的性質(zhì)，即時(shí)間越長(zhǎng)，井筒周?chē)貙蛹盁醿?chǔ)內(nèi)基巖被冷卻的范圍越大。這是由于水流流經(jīng)系統(tǒng)時(shí)帶走了地層中的熱量，然而水循環(huán)沒(méi)有停止，而裂隙儲(chǔ)層的溫度卻越來(lái)越低，導(dǎo)致流入采出井的熱水溫度也越來(lái)越低，從而使開(kāi)采溫度逐年降低。

表2 數(shù)值模擬算例表

圖4 溫度分布云圖

3.2 裂隙面與采出井溫度

在裂隙面上選取1～5號(hào)點(diǎn)（圖5-a），分析其溫度隨時(shí)間的變化規(guī)律（圖6-a）：裂隙面上各點(diǎn)溫度都隨時(shí)間推移而降低，且越靠近注入井溫度越低，越靠近采出井溫度越高。這是由于注入井井筒內(nèi)水溫較低而采出井附近基巖溫度較高，注入井內(nèi)的低溫水就近改變周?chē)严秲?chǔ)層，高溫的基巖就近改變附近裂隙水所致。在采出井井筒內(nèi)選取A～E共5個(gè)點(diǎn)（圖5-b），分析其溫度隨時(shí)間的變化規(guī)律，得到如圖6-b所示的變化曲線。可以得出：溫度總體上都隨時(shí)間的增加先升高而后降低，這是由于基巖溫度較裂隙水高，在其發(fā)生熱交換時(shí)，水流溫度升高而基巖溫度降低，隨著時(shí)間的增加，基巖降低的溫度來(lái)不及補(bǔ)給從而導(dǎo)致水流溫度下降，采出井內(nèi)水溫降低。

圖5 裂隙面及井筒取點(diǎn)位置簡(jiǎn)圖

3.3 井筒壁與保溫層

圖7 為算例1、4、5、6的采出井采出溫度變化曲線。由圖可知，熱突破時(shí)間隨開(kāi)孔長(zhǎng)度（L0）增大而提前，且相應(yīng)的冷尾跡效應(yīng)也更明顯。若以采出井采出溫度相對(duì)最高采出溫度下降20%作為系統(tǒng)運(yùn)行壽命的指標(biāo)，則L0=300 m、400 m、500 m時(shí)的最高采出溫度分別為114.00 ℃（第5年），116.98 ℃（第4年）與117.81 ℃（第3年），系統(tǒng)的壽命分別為12年、10年、8年。三者之中L0=300 m時(shí)系統(tǒng)壽命最長(zhǎng)，但最高采出溫度較L0=400 m時(shí)低了2.98 ℃，且晚了1年時(shí)間。綜合出力與壽命兩項(xiàng)因素，選取L0=400 m作為最佳開(kāi)孔長(zhǎng)度，不僅可以得到較高的開(kāi)采速率，還能保持一定的系統(tǒng)壽命。

圖6 各點(diǎn)溫度隨時(shí)間變化曲線圖

圖7 開(kāi)孔長(zhǎng)度影響采出溫度曲線圖

圖8 井筒壁熱傳導(dǎo)系數(shù)影響溫度曲線圖

圖8給出了算例1、7的采出井井筒內(nèi)溫度變化曲線。橫坐標(biāo)代表標(biāo)號(hào)為D～H各點(diǎn)，這5點(diǎn)坐標(biāo)分別為 H（500.1 m，950 m）、G（500.1 m，1 400 m）、F（500.1 m，1 850 m）、D（500.1 m，2 300 m）、E（500.1 m，2 750 m）。可以看出，開(kāi)采初期（1年），溫度沿井筒向上逐漸降低，且隨熱傳導(dǎo)系數(shù)增大而降低得更快，這是由于熱傳導(dǎo)系數(shù)增大促進(jìn)了與基巖的換熱過(guò)程。開(kāi)采前、中期（3、5、7年），溫度沿井筒向上降低的速度減緩，甚至中期不降反升，其原因在于儲(chǔ)層溫度不斷降低，導(dǎo)致井內(nèi)流體溫度降低，井內(nèi)流體從周?chē)鶐r吸熱使得溫度升高，熱傳導(dǎo)系數(shù)越大越有利于熱量吸收，從而溫度升高更快。進(jìn)入開(kāi)采后期（9年），基巖被井內(nèi)流體冷卻得較徹底，致使二者溫度相差不大，這時(shí)井內(nèi)流體溫度與井筒壁熱傳導(dǎo)系數(shù)關(guān)系不大。

圖9給出了算例1、8的溫度變化曲線。有保溫層的井筒壁與地層之間的換熱受阻，會(huì)使得采出溫度很快升高到最大值，即熱突破時(shí)間提前，但溫度卻很快下降，這與減小井筒壁熱傳導(dǎo)系數(shù)的情況一樣。說(shuō)明保溫層能提高初、前期的采出溫度，進(jìn)而提高這期間的開(kāi)采速率。

圖9 保溫層影響產(chǎn)出溫度曲線圖

3.4 采出溫度及熱開(kāi)采速率

圖10 、11為裂隙滲透率及裂隙寬度影響采出溫度與熱開(kāi)采速率的變化曲線圖。當(dāng)裂隙滲透率（Kf）增大至 1.5×106mD，裂隙寬度（df）增加至 0.30 cm時(shí)，最高采出溫度均超過(guò)120 ℃，且二者對(duì)采出溫度及熱開(kāi)采速率的影響表現(xiàn)出相同的規(guī)律性，即隨參數(shù)值的增大而產(chǎn)生較早熱突破，冷尾跡也更明顯，最高采出溫度和最大熱開(kāi)采速率也都有所提高。原因在于二者的物理效應(yīng)是等同的，裂隙寬度增加會(huì)導(dǎo)致滲透率增大，從而使冷卻時(shí)間提前，熱開(kāi)采速率提高，縮短開(kāi)采壽命。實(shí)際工程中，裂隙寬度會(huì)不斷變大，這是由于循環(huán)水會(huì)對(duì)裂隙產(chǎn)生擠壓變形，使流過(guò)區(qū)域?qū)挾仍黾樱M(jìn)而擴(kuò)大換熱面積，降低了開(kāi)采年限。為保持長(zhǎng)時(shí)間的熱能開(kāi)采，要想改善這種狀況，應(yīng)停止開(kāi)采一段時(shí)間，待儲(chǔ)層溫度恢復(fù)后再次提取地?zé)崮堋?/p>

4 結(jié)論

1）考慮井筒流動(dòng)換熱的井筒—熱儲(chǔ)耦合模擬研究，不僅實(shí)現(xiàn)了對(duì)增強(qiáng)型地?zé)嵯到y(tǒng)的完整性評(píng)價(jià)，還得出一項(xiàng)重要結(jié)論——注、采井的開(kāi)孔長(zhǎng)度對(duì)系統(tǒng)產(chǎn)能與壽命會(huì)產(chǎn)生重要影響，根據(jù)4種開(kāi)孔長(zhǎng)度的算例結(jié)果，L0=400 m時(shí)是最佳開(kāi)孔長(zhǎng)度，此時(shí)系統(tǒng)具有最佳出力與壽命，且在井筒壁上加保溫材料可以有效提高開(kāi)采初、前期采出溫度，以減少熱損失，提高開(kāi)采速率。

圖10 裂隙滲透率影響采出溫度及熱開(kāi)采速率曲線圖

圖11 裂隙寬度影響采出溫度及熱開(kāi)采速率曲線圖

2）開(kāi)采初期，滲流主要在裂隙中，不連通裂隙對(duì)滲流及傳熱影響較小，使得注入井周?chē)霈F(xiàn)明顯低溫區(qū)，隨時(shí)間增加，低溫區(qū)逐漸沿裂隙通道向采出井方向推移，并且范圍越來(lái)越大，到后期，低溫區(qū)蔓延至采出井，不連通裂隙對(duì)巖石的傳熱也已經(jīng)進(jìn)行。溫度的傳播過(guò)程始于裂隙通道，進(jìn)而擴(kuò)展到整個(gè)人工熱儲(chǔ)層區(qū)域，這時(shí)應(yīng)停止開(kāi)采一段時(shí)間，待儲(chǔ)層熱量恢復(fù)后，才能繼續(xù)維持較長(zhǎng)時(shí)間的地?zé)豳Y源提取。

3）裂隙滲透率、裂隙寬度等參數(shù)對(duì)開(kāi)采速率的影響都呈現(xiàn)正相關(guān)性，即隨參數(shù)值增大，達(dá)到最大采出溫度的時(shí)間縮短，提高了熱開(kāi)采速率，降低了開(kāi)采年限。