干熱巖流動換熱多尺度有限容積法

宇 波 李庭宇 韓東旭 孫東亮 楊福勝 魏進家

1. 北京石油化工學(xué)院機械工程學(xué)院 2. 西安交通大學(xué)化學(xué)工程與技術(shù)學(xué)院

0 引言

干熱巖一般指埋藏于地下3 ～10 km，內(nèi)部不存在流體（或僅有少量流體）、溫度高于150 ℃且具有開發(fā)價值的熱巖體，是一種儲量巨大的清潔能源[1-2]。干熱巖熱能開采過程與致密油氣田的開采有相似之處，均需要通過壓裂的方式增強儲層的導(dǎo)流性能。但是油氣開采主要為“物質(zhì)”的采出，而干熱巖開采為“熱能”的采出，二者有本質(zhì)的區(qū)別。數(shù)值模擬是研究干熱巖熱能開發(fā)過程的重要手段之一，采用該方法研究干熱巖儲層中的熱能提取規(guī)律，對于干熱巖開發(fā)方案的制訂、運行優(yōu)化以及采熱能力評價具有重要意義[3-6]。對干熱巖儲層進行數(shù)值模擬主要有兩個難點：①裂縫模型的表征和流動—換熱模型的構(gòu)建；②高效的數(shù)值求解方法。

裂縫網(wǎng)絡(luò)的表征模型在油藏數(shù)值模擬中就已經(jīng)取得了很多進展，近年來相關(guān)模型被越來越多地用于干熱巖儲層的描述中，其中應(yīng)用最廣泛的裂縫模型為離散裂縫模型[7-8]。Sun 等[9]和Shi 等[10]分別使用COMSOL Multiphysics 基于離散裂縫模型分別建立了干熱巖熱能開采過程的二維、三維熱水力三場耦合模型。離散裂縫模型在劃分網(wǎng)格時，裂縫要與基巖網(wǎng)格相互對齊，在模擬具有復(fù)雜裂縫型儲藏時必須使用大量的非結(jié)構(gòu)化網(wǎng)格，嚴重影響數(shù)值計算效率。為了克服離散裂縫模型的缺點，Lee 等[11]發(fā)展了嵌入式離散裂縫模型，該模型將裂縫型計算區(qū)域分成兩種介質(zhì)——基巖和裂縫，前者可以采用結(jié)構(gòu)化網(wǎng)格劃分，而后者可以直接嵌入基巖結(jié)構(gòu)化網(wǎng)格系統(tǒng)中，避免了傳統(tǒng)離散裂縫模型復(fù)雜的非結(jié)構(gòu)化網(wǎng)格剖分過程。Yan 等[12]基于該模型開發(fā)了裂縫巖體的流固耦合模型，計算結(jié)果與商業(yè)軟件相差不大且計算速度較快。Karvounis 和Jenny[13]改進了嵌入式離散裂縫模型，提出采用自適應(yīng)分層裂縫模型來模擬增強型地?zé)嵯到y(tǒng)的運行過程。

實際儲層往往復(fù)雜多變，天然裂縫和人工裂縫的長度可以從幾毫米到數(shù)十米不等。熱儲層內(nèi)的物性參數(shù)，如滲透率、孔隙度、導(dǎo)熱系數(shù)等均呈現(xiàn)出不均勻性。因此，采用傳統(tǒng)的有限單元法[9-10]、有限容積法[13]或者二者的混合求解法[12]進行求解時，需要劃分較細的網(wǎng)格來實現(xiàn)不均勻性的描述，從而導(dǎo)致計算耗時非常大。為此，研究人員發(fā)展了一些以多尺度求解方法為代表的模型降階方法[14]，如多尺度有限單元法[15-16]、多尺度有限容積法[17-18]等。多尺度方法的核心思想是通過局部細尺度（細網(wǎng)格）的求解獲得當(dāng)?shù)鼗瘮?shù)，然后基于獲得的基函數(shù)在粗尺度（粗網(wǎng)格）上對控制方程離散并求解得到粗尺度的解，最后采用局部基函數(shù)重構(gòu)得到細尺度的解。由于不需要直接求解細尺度的離散方程，多尺度方法可以顯著提高模擬效率，成為近年來的研究熱點。

相比于多尺度有限單元法，多尺度有限容積法不僅可以獲得局部守恒的通量且不會引入過多的求解自由度[19]。目前，多尺度有限容積法已經(jīng)從簡單的多孔介質(zhì)單相流動問題被推廣到了多相流[20]、非牛頓流[21]、流固耦合問題[22]等領(lǐng)域。但是，研究人員發(fā)現(xiàn)對于不均勻性非常強的問題或網(wǎng)格長寬比較大時，傳統(tǒng)的多尺度有限容積法的求解精度明顯降低，為此學(xué)者提出了迭代多尺度有限容積法[20,23-25]，其通過不斷的迭代求解來改進局部問題的邊界條件，從而提高多尺度方法的求解精度。Hajibeygi 等[20]提出了類似于幾何多重網(wǎng)格的迭代多尺度有限容積法，并將其應(yīng)用于不可壓縮多孔介質(zhì)單相和多相流動問題，其中細尺度光順采用線松弛技術(shù)；Lunati 等[24]基于多尺度算子提出了一種健壯性較好的多尺度有限容積法，使局部問題的邊界條件更為精確。M?yner 等[25]提出了模擬多孔介質(zhì)流動的多尺度投影限定光滑基方法。該方法預(yù)先定義基巖和裂縫粗網(wǎng)格內(nèi)多尺度基函數(shù)的支撐區(qū)域，并通過迭代光順法使得多尺度基函數(shù)的求解更加便捷。隨后Shah 等[26]將該方法推廣到了求解裂縫型多孔介質(zhì)的流動問題。

受投影限定光滑基模型的啟發(fā)，筆者提出了模擬干熱巖熱儲流動和換熱過程一體化多尺度有限容積求解方法（Thermal-Hydraulic Coupling Multi-scale Finite Volume Method, HT-MsFVM）。與傳統(tǒng)多尺度有限容積求解方法不同，該方法在求解裂縫型流動換熱問題時主要有兩點優(yōu)勢：①在求解溫度基函數(shù)時無需單獨計算各個粗網(wǎng)格上的局部瞬態(tài)對流換熱問題，而是通過細尺度上的離散矩陣及預(yù)先定義的延拓算子（基函數(shù)），通過雅克比迭代即可求得局部溫度基函數(shù)；②通過兩步迭代求解策略，可逐步消除細網(wǎng)格上的高頻誤差和粗網(wǎng)格上的低頻誤差，對于復(fù)雜裂縫流動換熱問題，精度更高。HT-MsFVM 方法主要實施過程如下：①基于嵌入式離散裂縫模型，分別采用兩套控制方程描述基巖和裂縫介質(zhì)的流動和換熱過程，其中基巖和裂縫介質(zhì)之間通過質(zhì)量和能量交換項進行耦合；②采用傳統(tǒng)有限容積法將兩套方程在細尺度上分別離散為稀疏矩陣形式；③通過定義多尺度網(wǎng)格系統(tǒng)及限定和延拓算子，將細尺度上的離散系統(tǒng)分別在細網(wǎng)格上光順和粗網(wǎng)格上迭代求解，達到最終所需要的求解精度。

1 基于嵌入式離散裂縫模型的儲層流動換熱描述

1.1 嵌入式離散裂縫模型

在干熱巖開發(fā)過程中，儲層中的裂縫是流體的主要流動通道和換熱場所，因此裂縫的刻畫至關(guān)重要。嵌入式離散裂縫模型是目前較為流行的裂縫表征方式。該模型有兩個關(guān)鍵點：①采用降維的方法描述裂縫網(wǎng)絡(luò)，在二維問題中，裂縫被看作是一維線段（圖1），而在三維問題里裂縫則是一個二維平面；②將基巖和裂縫看成兩個獨立的介質(zhì)，根據(jù)局部守恒定律構(gòu)建裂縫、基巖和裂縫與裂縫之間的質(zhì)量和能量交換。嵌入式離散裂縫模型最大的優(yōu)點是基巖和裂縫可以通過結(jié)構(gòu)化網(wǎng)格單獨劃分，避免了局部網(wǎng)格的加密（圖1）。

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

1.2 流動方程

如圖1 所示，嵌入式離散裂縫模型將基巖和裂縫看成是兩套獨立的介質(zhì)?；鶐r與裂縫、裂縫與裂縫之間通過竄流函數(shù)進行耦合?；鶐r系統(tǒng)質(zhì)量守恒方程為：

式中TI 表示相交裂縫間的交換系數(shù)；bi表示裂縫段的開度；Li表示裂縫段的長度。為了提高模擬效率，式（5）的推導(dǎo)過程中忽略了相交裂縫處的網(wǎng)格尺寸。

1.3 換熱方程

基于嵌入式離散裂縫模型的思想，仍然采用兩套能量方程來分別描述基巖系統(tǒng)和裂縫系統(tǒng)，同時需要考慮二者之間的能量交換，并假設(shè)流體和巖石之間處于局部熱平衡狀態(tài)。對于基巖系統(tǒng)有

對于第i 條裂縫介質(zhì)，采用以下能量方程進行描述：

采用有限容積法可將上述流動和換熱方程分別離散為如下矩陣形式：

或矢量形式：

式中A 表示離散的系數(shù)矩陣；B 表示方程右端項；x表示待求變量，x=p 或x=T。

式（10）或（11）可以采用傳統(tǒng)的迭代法進行求解。但當(dāng)計算區(qū)域較大，離散裂縫較多時，數(shù)值模擬較耗時，因此筆者采用多尺度有限容積法進行求解。

2 多尺度有限容積法

2.1 多尺度網(wǎng)格系統(tǒng)

圖2 多尺度網(wǎng)格系統(tǒng)圖

2.2 多尺度求解方法

根據(jù)圖2 所示的網(wǎng)格系統(tǒng)，需要定義粗細網(wǎng)格之間的映射關(guān)系：

2）粗網(wǎng)格修正：消除粗網(wǎng)格上的低頻誤差，同時保證粗網(wǎng)格上的質(zhì)量和能量守恒。

值得注意的是，利用多尺度壓力解p 計算得到的速度場v 只能保證在粗網(wǎng)格上的質(zhì)量守恒。這是由于在計算局部基函數(shù)時忽略了粗網(wǎng)格內(nèi)邊界上的通量傳遞[19]。為了得到細網(wǎng)格上守恒的速度場v'，還需要計算新的壓力場p'：

根據(jù)新的壓力場p'可得到守恒的速度場v'：

利用守恒的速度場v'即可通過方程（17）和（18）計算多尺度溫度解。壓力和溫度的多尺度循環(huán)求解過程見圖3：首先，根據(jù)方程（17）和（18）計算多尺度壓力解；其次，根據(jù)方程（19）和（20）重構(gòu)溫度方程所需要的速度場；最后，根據(jù)方程（17）和（18）計算多尺度溫度解。重復(fù)上述步驟直至達到所規(guī)定的多尺度循環(huán)次數(shù)（IC）。

圖3 壓力和溫度多尺度循環(huán)求解過程偽代碼圖

3 數(shù)值算例

3.1 交叉裂縫算例

本文基于嵌入式離散裂縫模型建立的熱開采模型的有效性已經(jīng)經(jīng)過實踐驗證[30-31]。為了驗證多尺度有限容積法的精度，設(shè)計了圖4-a 所示的“十”字交叉裂縫算例。模型邊長為100 m，水平裂縫和垂直裂縫的長度均為40 m。模型四周為非滲透絕熱邊界條件。左下角藍點為注入井的位置，注入流量為500 m3/d，注入工質(zhì)的溫度為20 ℃。圖4 右上角紅點為開采井的位置，開采壓力為1 MPa，其他物性參數(shù)如表1 所示。圖4-b 為多尺度計算網(wǎng)格。細網(wǎng)格總數(shù)為10 000，裂縫細網(wǎng)格總數(shù)為82?；鶐r粗網(wǎng)格個數(shù)為25，裂縫粗網(wǎng)格個數(shù)為4。

圖4 交叉裂縫算例圖

表1 模型物性參數(shù)表

為了定量比較多尺度解和細尺度解之間的差距，定義了如下二范數(shù)相對偏差：

式中xfs表示細尺度壓力解或溫度解；xms表示多尺度壓力解或溫度解。

圖5 為壓力和溫度基函數(shù)的空間分布情況。圖5-a中①和②及圖5-b 中①和②分別為第13 個基巖粗網(wǎng)格和第3 個裂縫粗網(wǎng)格內(nèi)的壓力和溫度基函數(shù)空間分布。從基巖的基函數(shù)可以看出，相交裂縫間的導(dǎo)流作用比導(dǎo)熱作用更明顯。裂縫支撐域內(nèi)的壓力和溫度基函數(shù)的空間分布差別不大。圖5-a 中③和④及圖5-b 中③和④分別為基巖和裂縫粗網(wǎng)格壓力和溫度基函數(shù)各自單獨求和后的云圖?？梢钥闯?，無論是壓力還是溫度基函數(shù)，單獨將基巖或裂縫基函數(shù)求和并不滿足單位分解特性，這是由于根據(jù)式（15）迭代計算得到的基函數(shù)包含了基巖和裂縫的耦合影響，但將圖5-a 中③和④及圖5-b 中③和④分別求和后即可使壓力或溫度基函數(shù)滿足單位分解。

圖5 壓力和溫度基函數(shù)空間分布圖

圖6-a 為采用傳統(tǒng)有限容積法在細尺度上計算開采10 年后得到的壓力解。圖6-b 為采用本文提出的多尺度有限容積法計算開采10 年后得到的壓力場，圖6-c為采用傳統(tǒng)多尺度有限容積法計算開采10 年后得到的壓力場。從圖中可以看出圖6-a 和圖6-b 幾乎呈現(xiàn)出了一致的壓力分布，根據(jù)式（21）算得的相對偏差為2.41×10－3。圖6-a與圖6-c之間的相對偏差為0.758。圖7 為采用兩種多尺度方法計算得到的溫度解和細尺度溫度解間的比較，從圖中可以看出無論是基巖區(qū)域還是裂縫區(qū)域，圖7-a 和圖7-b 中的兩組溫度結(jié)果均非常的接近，相對偏差為2.30×10－3，而采用傳統(tǒng)多尺度有限體積法計算得到的溫度解（圖7-c）則與細尺度溫度解（圖7-a）差距較大，相對偏差為0.186。

3.2 復(fù)雜裂縫算例

為了進一步驗證多尺度有限容積法的計算性能，設(shè)計了圖8-a 所示的復(fù)雜裂縫算例。該算例中，除裂縫條數(shù)為101，其余計算參數(shù)和前述交叉裂縫算例相同。如圖8-b 所示，基巖和裂縫的細尺度網(wǎng)格數(shù)分別為10 000 和2 226?；鶐r和裂縫的粗網(wǎng)格數(shù)分別為25 和59。為了表征儲層的非均質(zhì)性，筆者將滲透率、孔隙度以及熱導(dǎo)率設(shè)定為高斯隨機分布（圖9）。

圖10 為開采20 a 時的細尺度壓力解和分別采用兩種多尺度方法求得的壓力解間的比較。首先從圖10-a 中可以看出，由于裂縫的滲透率較大，導(dǎo)流作用明顯，致使裂縫附近的壓力空間分布與遠場基巖的壓力分布幾乎保持一致，同時可以看出：由本文提出的多尺度方法計算得到的壓力解（圖10-b）與細尺度壓力解（圖10-a）吻合較好，而由傳統(tǒng)多尺度方法得到的壓力解（圖10-c）誤差較大。圖11 為開采20 a 時的采用不同求解方法得到的溫度解對比。從圖11-a 可以看出，由于儲層的非均質(zhì)性和裂縫的存在，使得溫度空間分布呈現(xiàn)出各向異性。從圖11-b 可以看出，本文方法得到的多尺度溫度解能夠非常細致地捕捉到溫度變化的各個細節(jié)，而傳統(tǒng)多尺度計算方法則失效（圖11-c）。

不同開采時間下HT-MsFVM 多尺度解和細尺度解間的二范數(shù)相對偏差如表2 所示。從表2 中可以看出，壓力和溫度的相對偏差隨著開采時間的增加略有增多，這主要是前期數(shù)值計算累積誤差的影響，但總體精度完全可以滿足實際需要。表3 為在開采10 a 時，不同多尺度循環(huán)次數(shù)（即圖3 中Ic 值不同）下多尺度解和細尺度解間的相對偏差。從表中可以看出，隨著循環(huán)次數(shù)的增加多尺度求解精度增加。這主要是循環(huán)次數(shù)較少時，細尺度上的高頻誤差殘留較多，從而導(dǎo)致累積誤差也會增多。

圖6 交叉裂縫算例壓力解比較圖

多尺度有限容積法所需要的整體計算時間主要消耗在粗網(wǎng)格劃分、基函數(shù)獲取和粗細網(wǎng)格間的迭代求解上。針對特定的干熱巖儲層，進行大量工況的數(shù)值模擬時，網(wǎng)格劃分只需進行一次。因此，為了最大程度的利用多尺度方法的計算性能，可以采用離線和在線相結(jié)合的方式。多尺度網(wǎng)格系統(tǒng)的劃分可以采用離線的方式一次計算完成，而基函數(shù)計算和多尺度迭代求解則采用在線的方式根據(jù)不同的工況的條件多次求解。圖12 為兩個算例分別采用傳統(tǒng)有限容積法和本文提出的HT-MsFVM 方法的在線耗時比較?？梢钥闯?，多尺度方法可以明顯提高數(shù)值計算效率，加速比可達4 ～5 倍。

圖7 交叉裂縫算例溫度解比較圖

圖8 復(fù)雜裂縫算例圖

圖9 物性參數(shù)隨機分布圖

圖10 復(fù)雜裂縫算例壓力解比較圖

4 結(jié)論

1）針對采用傳統(tǒng)有限容積法模擬干熱巖儲層采熱過程效率低的問題，基于嵌入式離散裂隙模型，提出了一體化的流動與換熱多尺度有限容積法求解技術(shù)（HT-MsFVM）。相比于傳統(tǒng)的粗化降階方法，該方法采用延拓算子和限定算子建立粗細尺度間的交互關(guān)系，在提高數(shù)值模擬計算效率的同時，保證了計算精度。

圖11 復(fù)雜裂縫算例溫度解比較圖

2）多尺度有限容積法的壓力場和溫度場求解精度會隨著模擬時間的增長而略有降低，主要是由于誤差累積導(dǎo)致的，但總體得到的結(jié)果完全可以滿足工程需求。增加多尺度循環(huán)次數(shù)可以提高計算精度，但循環(huán)次數(shù)過多會降低多尺度的加速效果。對于本文的算例而言，多尺度循環(huán)次數(shù)介于50 ～100 時，相對誤差的量級介于10－3～10－4，數(shù)值計算效率可提高4 ～5 倍。

表2 不同開采時間下HT-MsFVM 多尺度解和細尺度解間的相對偏差表

表3 不同多尺度循環(huán)次數(shù)下HT-MsFVM 多尺度解和細尺度解間的相對偏差表

圖12 多尺度求解方法和細尺度計算方法時間比較圖

3）在干熱巖的熱開采模擬過程中，需要數(shù)值計算大量的算例，從而實現(xiàn)熱開采參數(shù)的優(yōu)化和不確定性評估。對于通過改變邊界條件或物性參數(shù)等進行大量數(shù)值實驗的問題，采用本文提出的多尺度計算方法可有效提高計算效率。

4）為了方便比較，本文暫未考慮物性參數(shù)隨溫度變化帶來的物理場空間分布的影響，因此得到的細尺度壓力和溫度解會與實際情況有一定的偏差。