宇 波 李庭宇 韓東旭 孫東亮 楊福勝 魏進(jìn)家

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

0 引言

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

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

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

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

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

1 基于嵌入式離散裂縫模型的儲(chǔ)層流動(dòng)換熱描述

1.1 嵌入式離散裂縫模型

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

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

1.2 流動(dòng)方程

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

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

1.3 換熱方程

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

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

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

或矢量形式：

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

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

2 多尺度有限容積法

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

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

2.2 多尺度求解方法

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

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

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

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

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

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

3 數(shù)值算例

3.1 交叉裂縫算例

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

圖4 交叉裂縫算例圖

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

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

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

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

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

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

3.2 復(fù)雜裂縫算例

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

表2 不同開采時(shí)間下HT-MsFVM 多尺度解和細(xì)尺度解間的相對(duì)偏差表

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

圖12 多尺度求解方法和細(xì)尺度計(jì)算方法時(shí)間比較圖

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

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