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

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

0 引言

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

裂縫網(wǎng)絡(luò)的表征模型在油藏?cái)?shù)值模擬中就已經(jīng)取得了很多進(jìn)展，近年來相關(guān)模型被越來越多地用于干熱巖儲層的描述中，其中應(yīng)用最廣泛的裂縫模型為離散裂縫模型[7-8]。Sun 等[9]和Shi 等[10]分別使用COMSOL Multiphysics 基于離散裂縫模型分別建立了干熱巖熱能開采過程的二維、三維熱水力三場耦合模型。離散裂縫模型在劃分網(wǎng)格時(shí)，裂縫要與基巖網(wǎng)格相互對齊，在模擬具有復(fù)雜裂縫型儲藏時(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)格剖分過程。Yan 等[12]基于該模型開發(fā)了裂縫巖體的流固耦合模型，計(jì)算結(jié)果與商業(yè)軟件相差不大且計(jì)算速度較快。Karvounis 和Jenny[13]改進(jìn)了嵌入式離散裂縫模型，提出采用自適應(yīng)分層裂縫模型來模擬增強(qiáng)型地?zé)嵯到y(tǒng)的運(yùn)行過程。

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

相比于多尺度有限單元法，多尺度有限容積法不僅可以獲得局部守恒的通量且不會引入過多的求解自由度[19]。目前，多尺度有限容積法已經(jīng)從簡單的多孔介質(zhì)單相流動問題被推廣到了多相流[20]、非牛頓流[21]、流固耦合問題[22]等領(lǐng)域。但是，研究人員發(fā)現(xiàn)對于不均勻性非常強(qiáng)的問題或網(wǎng)格長寬比較大時(shí)，傳統(tǒng)的多尺度有限容積法的求解精度明顯降低，為此學(xué)者提出了迭代多尺度有限容積法[20,23-25]，其通過不斷的迭代求解來改進(jìn)局部問題的邊界條件，從而提高多尺度方法的求解精度。Hajibeygi 等[20]提出了類似于幾何多重網(wǎng)格的迭代多尺度有限容積法，并將其應(yīng)用于不可壓縮多孔介質(zhì)單相和多相流動問題，其中細(xì)尺度光順采用線松弛技術(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)多尺度有限容積求解方法不同，該方法在求解裂縫型流動換熱問題時(shí)主要有兩點(diǎn)優(yōu)勢：①在求解溫度基函數(shù)時(shí)無需單獨(dú)計(jì)算各個(gè)粗網(wǎng)格上的局部瞬態(tài)對流換熱問題，而是通過細(xì)尺度上的離散矩陣及預(yù)先定義的延拓算子（基函數(shù)），通過雅克比迭代即可求得局部溫度基函數(shù)；②通過兩步迭代求解策略，可逐步消除細(xì)網(wǎng)格上的高頻誤差和粗網(wǎng)格上的低頻誤差，對于復(fù)雜裂縫流動換熱問題，精度更高。HT-MsFVM 方法主要實(shí)施過程如下：①基于嵌入式離散裂縫模型，分別采用兩套控制方程描述基巖和裂縫介質(zhì)的流動和換熱過程，其中基巖和裂縫介質(zhì)之間通過質(zhì)量和能量交換項(xiàng)進(jìn)行耦合；②采用傳統(tǒng)有限容積法將兩套方程在細(xì)尺度上分別離散為稀疏矩陣形式；③通過定義多尺度網(wǎng)格系統(tǒng)及限定和延拓算子，將細(xì)尺度上的離散系統(tǒng)分別在細(xì)網(wǎng)格上光順和粗網(wǎng)格上迭代求解，達(dá)到最終所需要的求解精度。

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

1.1 嵌入式離散裂縫模型

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

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

1.2 流動方程

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

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

1.3 換熱方程

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

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

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

或矢量形式：

式中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ì)算得到的速度場v 只能保證在粗網(wǎng)格上的質(zhì)量守恒。這是由于在計(jì)算局部基函數(shù)時(shí)忽略了粗網(wǎng)格內(nèi)邊界上的通量傳遞[19]。為了得到細(xì)網(wǎng)格上守恒的速度場v'，還需要計(jì)算新的壓力場p'：

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

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

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

3 數(shù)值算例

3.1 交叉裂縫算例

本文基于嵌入式離散裂縫模型建立的熱開采模型的有效性已經(jīng)經(jīng)過實(shí)踐驗(yàn)證[30-31]。為了驗(yàn)證多尺度有限容積法的精度，設(shè)計(jì)了圖4-a 所示的“十”字交叉裂縫算例。模型邊長為100 m，水平裂縫和垂直裂縫的長度均為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ù)相對偏差：

式中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ú)求和后的云圖?？梢钥闯?，無論是壓力還是溫度基函數(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 年后得到的壓力場，圖6-c為采用傳統(tǒng)多尺度有限容積法計(jì)算開采10 年后得到的壓力場。從圖中可以看出圖6-a 和圖6-b 幾乎呈現(xiàn)出了一致的壓力分布，根據(jù)式（21）算得的相對偏差為2.41×10－3。圖6-a與圖6-c之間的相對偏差為0.758。圖7 為采用兩種多尺度方法計(jì)算得到的溫度解和細(xì)尺度溫度解間的比較，從圖中可以看出無論是基巖區(qū)域還是裂縫區(qū)域，圖7-a 和圖7-b 中的兩組溫度結(jié)果均非常的接近，相對偏差為2.30×10－3，而采用傳統(tǒng)多尺度有限體積法計(jì)算得到的溫度解（圖7-c）則與細(xì)尺度溫度解（圖7-a）差距較大，相對偏差為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。為了表征儲層的非均質(zhì)性，筆者將滲透率、孔隙度以及熱導(dǎo)率設(shè)定為高斯隨機(jī)分布（圖9）。

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

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

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

多尺度有限容積法所需要的整體計(jì)算時(shí)間主要消耗在粗網(wǎng)格劃分、基函數(shù)獲取和粗細(xì)網(wǎng)格間的迭代求解上。針對特定的干熱巖儲層，進(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）針對采用傳統(tǒng)有限容積法模擬干熱巖儲層采熱過程效率低的問題，基于嵌入式離散裂隙模型，提出了一體化的流動與換熱多尺度有限容積法求解技術(shù)（HT-MsFVM）。相比于傳統(tǒng)的粗化降階方法，該方法采用延拓算子和限定算子建立粗細(xì)尺度間的交互關(guān)系，在提高數(shù)值模擬計(jì)算效率的同時(shí)，保證了計(jì)算精度。

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

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

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

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

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

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

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