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

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

1. 北京石油化工學院機械工程學院 2. 西安交通大學化學工程與技術學院

0 引言

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

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

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

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

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

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

1.1 嵌入式離散裂縫模型

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

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

1.2 流動方程

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

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

1.3 換熱方程

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

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

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

或矢量形式：

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

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

2 多尺度有限容積法

2.1 多尺度網格系統

圖2 多尺度網格系統圖

2.2 多尺度求解方法

根據圖2 所示的網格系統，需要定義粗細網格之間的映射關系：

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

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

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

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

圖3 壓力和溫度多尺度循環求解過程偽代碼圖

3 數值算例

3.1 交叉裂縫算例

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

圖4 交叉裂縫算例圖

表1 模型物性參數表

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

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

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

圖5 壓力和溫度基函數空間分布圖

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

3.2 復雜裂縫算例

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

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

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

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

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

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

圖8 復雜裂縫算例圖

圖9 物性參數隨機分布圖

圖10 復雜裂縫算例壓力解比較圖

4 結論

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

圖11 復雜裂縫算例溫度解比較圖

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

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

表3 不同多尺度循環次數下HT-MsFVM 多尺度解和細尺度解間的相對偏差表

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

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

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