縫洞型巖溶熱儲流動傳熱耦合數值模擬

姚 軍 張 旭 黃朝琴 鞏 亮 楊文東 李 陽

中國石油大學（華東）油氣滲流研究中心

0 引言

碳酸鹽巖巖溶熱儲是一種典型的地熱儲層，具有出水量大且地熱利用后的尾水易于回灌的優勢，是我國最具開發利用潛力的主力地熱儲層[1]。我國碳酸鹽巖的分布總面積占陸地面積的三分之一，裸露面積約為90×104km2，隱伏面積達250×104km2以上[2]。處于較高熱流值背景下的幾大主要區域有：渤海灣盆地、鄂爾多斯盆地、蘇北盆地、江漢盆地、楚雄盆地、蘭坪思茅盆地、昂拉仁錯盆地和羌塘盆地[2]。碳酸鹽巖熱儲資源量占我國水熱型地熱資源總量的70%～80%[3]。由于區域構造運動和巖溶作用，巖溶熱儲內發育大量的裂縫及溶洞，裂縫長度跨越從厘米級到千米級，天然裂縫網絡處于連通或不連通狀態；溶洞尺寸跨越從厘米級到米級，空間分布呈現與裂縫相連或分離特征[4-7]。因此，與孔隙型或裂縫型熱儲相比，縫洞型熱儲的儲集空間類型多樣（孔、縫、洞），具有顯著的多尺度、強非均質特征，流體流動既有多孔介質區域（巖石基質及裂縫系統）的滲流又有溶洞區域的自由流[8]。目前，針對孔隙型和裂縫型熱儲的熱采動態已進行大量數值模擬研究[9-30]，但由于縫洞型熱儲內具有復雜的多尺度、強非均質、多流態特征，縫洞型熱儲熱采過程中涉及的流動、傳熱特征及熱采動態尚不清晰，亟需建立一套適用于縫洞型巖溶熱儲流動傳熱模擬的理論和方法。

當前國內外學者在縫洞型碳酸鹽巖滲流理論方面做了大量的工作，發展了等效連續介質和離散介質數學模型，其中等效連續介質模型又可分為等效單重介質模型和三重介質模型。等效單重介質模型是將整個縫洞型介質視為一個連續體，通過等效參數表示其非均質性[31-32]。因此，該模型具有計算效率高、參數求取簡單等優點。原理上，等效單重介質模型僅適用于縫洞系統發育程度高、連通性好的儲層。另外，無法考慮裂縫網絡連通性和其他裂縫屬性對流動過程的影響。三重介質模型是雙重介質模型的擴展和延伸，將縫洞型介質劃分為三個平行的連續介質系統，即低滲透基巖、高滲透裂縫和溶洞，通過竄流函數有機耦合各系統[33-34]。三重介質模型能夠有效刻畫縫洞型介質中的優先流現象，并充分考慮基巖、裂縫和溶洞系統間的物質交換，比較符合實際。但該模型對縫洞網絡幾何形態過于簡化，且竄流函數也難以確定，不適用于離散分布的大尺度裂縫和溶洞的模型，具有一定的局限性。離散介質模型是將裂縫和溶洞的幾何結構進行顯式表征，能夠準確描述縫洞型介質中的流體流動過程[8,35-36]。目前，國內外僅對縫洞型介質內流體流動進行了研究，對于其內部流動換熱機理及熱采特性的研究，未見相關的研究報道。對此，本研究建立了基于離散縫洞網絡方法的熱流耦合模型，研究了裂縫及溶洞的多尺度、強非均質及多流態特征對縫洞型熱儲內流動和傳熱過程的影響，對標裂縫型熱儲開發動態，系統分析縫洞型巖溶熱儲開發特征，以期為其高效開發提供理論基礎。

1 縫洞型巖溶熱儲熱流耦合模型

1.1 縫洞網絡模型

該研究中，通過直線段和圓分別表征裂縫和溶洞幾何形狀[32,37]，如圖1所示。

圖1 縫洞型巖溶熱儲及離散縫洞網絡示意圖

天然縫洞系統通常呈現顯著的多尺度特征，裂縫長度和溶洞直徑服從冪律分布，其概率密度函數n(lf)、n(lc)表達式[38-42]：

式中l表示裂縫長度或溶洞直徑，m；下標f和c分別表示裂縫和溶洞；α是密度項；a表示冪律尺寸指數。

該模型唯一的特征長度參數是裂縫或溶洞最小尺寸和最大尺寸，即lmin和lmax。大量現場數據表明，裂縫和溶洞的長度指數a分布范圍分別介于1.5～3.5和2.2～2.6[41-42]。

1.1.1 裂縫密度與連通參數

式中γf表示裂縫密度，1/m，為單位模型面積內裂縫長度之和，與裂縫長度密度分布函數有關[43-44]；L表示模型尺寸，m；pf表示連通參數，表征裂縫網絡連通程度[37]。

pf值越大，系統連通程度越高。當大于連通閾值（pfc）時，裂縫網絡是連通的；反之，則不連通。對于二維隨機裂縫網絡，pfc具有尺度無關性，pfc值在較小的范圍內變化，即5～7之間，也就是pfc=6±1。

1.1.2 溶洞孔隙度

溶洞孔隙度（φc）表示單位模型面積內溶洞面積之和[35, 45]：

式中Ac表示模型內溶洞總面積，m2；U表示求和運算符；Nc表示溶洞總數。

1.2 流動傳熱耦合模型

基于離散縫洞網絡方法，建立適用于模擬縫洞型巖溶熱儲內流體流動及熱量傳遞的熱流耦合模型。離散縫洞網絡方法將縫洞型巖溶熱儲視為由巖石基質、裂縫系統和溶洞系統組成的儲層，其中裂縫和溶洞嵌套于巖石基質中，相互分離或連接形成縫洞網絡（圖1）。巖石基質和裂縫系統視為多孔介質區域，溶洞系統視為自由流區域。在流動模擬中，考慮不同流動區域間的流動耦合，即多孔介質區域滲流與溶洞區域自由流耦合，滲流區域滿足Darcy方程，自由流區域滿足Navier-Stokes (N-S)方程，兩區域交界面采用Beavers-Joseph-Saffman(簡稱BJS)條件進行流動耦合[46-48]。在傳熱模擬中，考慮熱傳導和熱對流兩種傳熱方式，熱巖石基質與冷流體在縫洞界面上熱量交換滿足局部非熱平衡效應[9]。

1.2.1 流動方程

基質內滲流方程[9]：

式中ρf表示流體密度，kg/m3；Sm表示基質儲水系數，1/Pa；pp表示多孔介質區域流體壓力，Pa；t表示時間，s；um表示基質內滲流速度，m/s；km表示基質滲透率，mD；ηf表示流體動力黏度，Pa·s；g表示重力加速度，m/s2；z表示沿重力方向單位向量。

裂縫內滲流方程[9]：

式中df表示裂縫開度，m；Sf表示裂縫內儲水系數，1/Pa；uf表示裂縫內滲流速度，m/s表示沿裂縫切向方向的梯度算子；kf表示裂縫滲透率，mD。

溶洞內流動方程[8,35]：

式中uc表示溶洞內流動速度，m/s；σf表示溶洞內流體應力，Pa；I表示單位張量表示流體應變張量；pc表示溶洞內流體壓力，Pa。

多孔介質區域與自由流區域交界面耦合流動滿足BJS邊界條件[46-48]：

式中up表示多孔介質內滲流速度，m/s；K表示多孔介質內滲透率張量，mD；n、τ分別表示多孔介質與溶洞界面單位法向向量和單位切向向量；β表示切向阻力系數。

基于Darcy定律估算縫洞型巖溶熱儲的等效滲透率（keq）：

式中Qout表示出口端的體積流量，m2/s；累加項（∑）表示出口端裂縫的流量計算；積分項（∫）表示出口端巖石基質和溶洞的流量計算；pin、pout分別表示注入端和出口端壓力，Pa。

1.2.2 傳熱方程

基質內傳熱方程[9]：

式中 (ρC)eff表示有效熱焓，J/(m3·K)；λeff表示有效導熱系數，W/(m·K)；下標s和f分別表示固體和流體；ε表示基質孔隙度；C表示比熱容，J/(kg·K)；λ表示導熱系數，W/(m·K)；q表示基質與流體在裂縫或溶洞界面上的換熱量，W/m3。

裂縫內傳熱方程[9]：

式中qf=h(Ts－Tf)為牛頓換熱公式，表示單位面積上基質與流體在裂縫面上的換熱量，W/m2；h表示對流換熱系數，W/(m2·K)。

流體在裂縫與溶洞相交處滿足溫度相等。溶洞內傳熱方程：

式中Tc表示溶洞內流體溫度，℃；qc表示單位面積上基質與流體在溶洞界面上的換熱量，W；A表示溶洞單元與基質單元接觸面積，m2；vc表示與基質單元接觸的溶洞單元體積，m3。

出口端平均溫度（Tout）計算公式為[15]：

整體采熱量（Wtotal）為[18]：

式中Qm表示出口端質量流量，kg/(m·s)；hpro、hinj分別表示開采流體和注入流體比焓，J/kg。

可采熱量（Wr）為：

式中Vp、Vc分別表示多孔介質區域體積、溶洞區域體積，m3；Ti、Tinj分別表示儲層初始溫度、注入流體溫度，℃。

熱儲整體熱采率（γglobal）為：

筆者基于有限元軟件COMSOL Multiphysics進行二次開發實現縫洞型巖溶熱儲熱流耦合模型的數值求解。基于Matlab編程生成滿足規定的概率分布的離散縫洞網絡模型，然后基于COMSOL Mutiphysics with Matlab將縫洞網絡幾何模型導入有限元軟件進行求解。采用無厚度平面單元模擬裂縫，實體單元模擬基質和溶洞。基于達西流和蠕動流模塊求解多孔介質和溶洞區域內流體流動；基于多孔介質傳熱模塊求解多孔介質區域換熱過程，自定義偏微分方程模塊求解裂縫和溶洞內換熱過程；進一步考慮流體與基質巖塊在裂縫和溶洞界面上的質量和能量交換，編寫相應的程序代碼，并嵌套在軟件中實現基于離散縫洞網絡方法的熱流耦合模型求解。與連續介質模型相比，離散介質方法可以通過流量等效等方法，對裂縫進行降維處理，簡化為1維線單元，并將裂縫開度耦合到控制方程中，這樣可以避免在裂縫區域進行大量的精細網格剖分，降低計算量，提高計算效率。

1.2.3 流動模型驗證

該研究中滲流—自由流耦合模型采用經典的Beavers—Joseph實驗模型進行驗證，如圖2-a所示。考慮多孔介質區域滲流及自由流區域流動為穩定層流，注入壓力（pin）設為0.5 Pa，開采壓力（pout）設為0 Pa，模型長度（L）設為0.5 m，高度（h）設為0.1 m；多孔介質區域孔隙度設為0.2，滲透率設為1 000 mD。流體的密度設為1 kg/m3，黏度設為1 Pa·s，在交界面上BJS切向阻力系數（β）設為1.0。圖2-b為注入端速度分布結果，表明數值解與解析解[49]基本一致，驗證了滲流—自由流耦合模型的正確性。

圖2 滲流—自由流耦合模型及注入段速度分布對比圖

1.2.4 熱流耦合模型驗證

通過與等效介質方法模擬縫洞介質熱流耦合過程的結果對比，驗證該研究中離散縫洞方法的準確性。如圖3-a所示，模型大小為1 m×1 m，多孔介質為均質，直徑0.2 m的溶洞及開度為0.15 mm的貫穿裂縫位于模型中心。縫洞介質初始飽和水，初始溫度和壓力分別為200 ℃和0 MPa。模型上下邊界為不滲透絕熱邊界，模型左右側為注入端和采出端，壓力分別為pin= 0.2 MPa和pout= 0 MPa，注入溫度為Tin= 60 ℃。詳細的模型參數如下：基質孔隙度0.05，基質滲透率為1 mD，裂縫滲透率為1 000 mD，流體的密度為1 000 kg/m3，黏度為0.001 Pa·s，流體導熱系數為 0.5 W/(m·K)，流體比熱容為4 200 J/(kg·K－1)，基質的密度為2 600 kg/m3，基質導熱系數為3.0 W/(m·K)，基質比熱容為 1 000 J/(kg·K－1)，基質儲水系數為1×10－9Pa－1，裂縫內儲水系數為1×10－10Pa－1，在交界面上BJS切向阻力系數（β）為1.0，對流換熱系數為3 000 W/(m2·K)。對于連續介質模型，裂縫顯式表征為2維薄單元層（圖3-b）；對于離散介質模型，裂縫表征為1維線單元（圖3-b）。

圖3 熱流耦合模型設置及模型網格剖分圖

圖4給出了縫洞介質在熱采1天后，縫洞介質內流速分布、流體壓力分布及溫度分布情況，結果表明基于連續介質模型和離散介質模型的結果基本一致。另外，圖5給出了沿中間橫軸縫洞介質內流速分布、沿中間縱軸縫洞介質內流速分布及出口平均溫度隨時間變化情況，可以發現裂縫內的流速遠高于溶洞和基質內流速，而溶洞內流速遠高于基質內流速，因此，與裂縫型介質類似，縫洞型介質中裂縫也是流體高速流動的主要通道。離散介質模型結果與連續介質模型結果的一致性，驗證了該研究中縫洞介質熱流耦合數值模型的準確性。

圖4 熱采1天后的流速、壓力、溫度分布圖

圖5 熱流耦合數值模擬的流速、溫度變化圖

2 模型設置與結果分析

2.1 模型設置

該研究基于長為L=100 m的正方形區域（水平面），建立了一系列2維隨機縫洞網絡。為分析縫洞尺寸分布的影響，假設縫洞位置和走向是隨機分布的。首先基于縫洞幾何統計信息，建立裂縫網絡和溶洞網絡，然后對二者進行疊加，生成離散縫洞網絡。裂縫長度（lf）服從冪律分布，上限和下限分別設置為lf,max=50L和lf,min=L/50。研究了3種裂縫長度指數情況：af=1.5、2.5和 3.5；3種裂縫密度情況：γf=0.2 m－1、0.4 m－1和 0.6 m－1；針對每一種裂縫參數組合下的裂縫網絡，計算了相應的裂縫網絡連通參數（pf）。溶洞直徑上限和下限分別設置為lc,max=L/5和lc,min=L/25。因為溶洞尺寸冪律指數變化范圍比較小（2.2～2.6），所以我們只考慮了溶洞尺寸指數ac=2.4的情況。該研究中溶洞孔隙度為φc=0.2。生成的裂縫網絡涵蓋了pf低于、接近和高于連通閾值（pfc）的情況，也就是存在不連通、處于過渡階段和連通的裂縫網絡。當af較大（如af=3.5）時，系統由小裂縫（lf＜L）控制；當af較小（如af=1.5）時，系統受大裂縫（lf與模型尺寸L相當）控制；當af=2.5時，系統受大裂縫和小裂縫同時控制。

縫洞型介質熱流耦合模擬參數已在1.2.4節描述，其儲層埋深為2 500 m，初始儲層為飽和水狀態，初始壓力為25 MPa。基于水平對井（即1口注入水平井和1口開采水平井）進行開采，左側注入井和右側開采井的壓力分別為27.5 MPa和25 MPa，保證流體在儲層中的循環流動。儲層初始溫度為100 ℃，注入邊界流體注入溫度為30 ℃。流動和傳熱邊界條件設置如圖3-a所示。

2.2 結果分析

2.2.1 流動過程

為了研究縫洞型熱儲內流體流動特征，本文對比分析了不同縫洞網絡參數（即af、γf、pf和φc）下，離散裂縫網絡和離散縫洞網絡的流速分布（圖6）以及等效滲透率變化情況（圖7）。為了方便對比和分析流體在具有相同裂縫幾何參數下的裂縫網絡和縫洞網絡內的流動，基于裂縫網絡內平均流速，對相應的裂縫網絡和縫洞網絡內的裂縫流速進行無量綱化（圖6）。

圖6 離散裂縫與離散縫洞網絡無量綱化流速分布圖

針對裂縫網絡內流體流動情況，如圖6、7所示，①對于af=3.5的裂縫網絡，流速相對較高的流動通道并不能貫穿注采兩端，所以流體只能通過低滲基質和不連通的裂縫網絡流動到出口端，系統的等效滲透率低；②對于af=1.5的裂縫網絡，注采兩端出現多條連通的高速流動通道，系統的等效滲透率高；③對于af=2.5的裂縫網絡，當裂縫密度較低時（如γf=0.2 m－1、0.4 m－1），系統內不存在貫穿注采端的高速流動通道，系統等效滲透率較低；而當裂縫密度增大時，如γf=0.6 m－1，系統內出現幾條貫通的高速流動通道，系統等效滲透率較高。

針對存在溶洞、縫洞網絡內流體流動情況時，我們發現溶洞對流體流動起重要作用，主要分為兩種作用：①溶洞的存在可以增多系統內貫穿注采端的高速流通通道或者使系統從不連通狀態變為連通狀態，例如，與af=2.5，γf=0.2 m－1的裂縫網絡相比（圖 6-a），由于溶洞的存在，系統從不連通狀態轉變為連通狀態，也就是相應的縫洞網絡內存在貫穿注采端的高速流動通道（圖6-b），系統的等效滲透率增大了將近5倍（圖7）；②溶洞的存在可以增大系統內與溶洞相連接的裂縫上的流速，相比于對應的離散裂縫網絡，造成系統的等效滲透率出現增大的情況（圖7）。

如圖7-a所示，系統連通參數（pf）對等效滲透率（keq）起控制作用，具體表現為keq隨pf的增加而增大。當pf＜pfc時，keq較低，這是因為系統內不存在貫通的高速流動通道，流體流動主要受低滲的基質控制。當pf＞pfc時，keq較高，這是因為系統內存在大量的貫通的高速流動通道，流體流動阻力低；當pf≈pfc時，keq出現突然增大，這是因為系統在該區域從不連通狀態變為了連通狀態，流動阻力顯著降低。另外，我們發現，溶洞的存在會增大系統的等效滲透率，特別是當系統處于連通閾值（pfc）附近時，系統等效滲透率出現顯著增大（圖7-b）。這是因為在pfc附近，系統由于溶洞的存在從不連通變為連通的概率大。

圖7 離散縫洞與離散裂縫模型的等效滲透率及比值隨連通參數變化圖

2.2.2 傳熱過程

為了研究縫洞型熱儲內傳熱特征，筆者對比分析了不同縫洞網絡參數（即af、γf、pf和φc）下，離散裂縫網絡和離散縫洞網絡的溫度分布（圖8），以及熱采效率評價指標變化情況（圖9）。

圖8 熱采1 800天后離散裂縫網絡及離散縫洞網絡內溫度分布圖

圖9 整體熱采率為30%時對應的熱采時間及其熱采時間比值隨連通參數變化圖

結合圖6-a、8-a分析裂縫網絡內流體流動下的傳熱情況，①對于af=3.5的裂縫網絡，冷卻前沿推進較慢，溫度分布較均勻；這是因為系統內不存在貫穿注采端的高速流動通道，流體流動速度較低，換熱過程由傳熱效率低的熱傳導方式控制。②對于af=1.5的裂縫網絡，冷卻前緣推進較快，冷流體波及面積較大，且溫度分布呈現不均勻狀態；這是因為系統內存在貫通的高速流通通道，換熱過程由傳熱效率高的熱對流方式控制。③對于af=2.5的裂縫網絡，當裂縫密度較低時（如γf=0.2 m－1、0.4 m－1），與af=3.5的裂縫網絡內的溫度分布特征一致；而當裂縫密度增大時（如γf=0.6 m－1），存在較明顯的冷卻前緣錐進現象；這是因為系統內僅存在少數幾條貫通的高速流動通道。

而針對存在溶洞、縫洞網絡內傳熱情況時（圖6-b、8-b），①發現對于af=3.5的裂縫網絡，溶洞的存在對熱儲傳熱過程的影響較小，主要是因為溶洞的存在只能改變局部流體流動速度，而不能在系統內形成貫通注采端的高速流動通道，換熱過程主要還是由熱傳導控制。②對于af=1.5的裂縫網絡，溶洞的存在使冷卻前緣推進加速，冷流體波及面積更大；這是因為溶洞的存在使系統內流體流速加快，熱對流過程增強。③對于af=2.5的裂縫網絡，當裂縫密度較大時（如γf=0.6 m－1），溶洞對系統傳熱的作用與對af=1.5的裂縫網絡的傳熱作用一致；而當裂縫密度較低時（如γf=0.2 m－1、0.4 m－1），從均勻的冷卻前緣變化為明顯的錐進的冷卻前緣，冷流體波及面積增大；這是因為溶洞的存在使系統從不連通狀態變化為連通狀態，系統換熱過程從效率低的熱傳導過程變為效率高的熱對流過程。

為了研究熱儲的換熱效率，定義了熱采效率評價指標：熱儲整體熱采率達到裂縫型熱儲可采熱量的30%，所對應的熱采時間（即tγglobal=30%）根據式（20）計算得到。熱采時間越短，表明在開采相同的熱儲能量情況下，熱采速度和熱采效率更高。如圖9所示，系統連通參數（pf）對熱采時間（即熱采效率）tγglobal起控制作用，具體表現為tγglobal隨pf的增加而降低。當pf＜pfc時，熱采時間長，熱采效率低；這是因為系統的傳熱主要受熱傳導控制，導致換熱效率較低。當pf＞pfc時，熱采時間短，熱采效率較高；這是因為系統內流動速度較大，傳熱過程主要受換熱效率高的熱對流方式控制。另外，我們發現，相比于裂縫網絡，溶洞的存在會增大系統的熱采效率，并且整體上，隨著pf的增大，溶洞對系統熱采效率的提高作用越大。特別地，當系統處于pfc附近時，由于溶洞的存在，系統的熱采效率出現顯著增大。

3 結論

1）針對縫洞型巖溶熱儲的多尺度、強非均質和多流態的特點，提出了基于離散縫洞網絡方法的熱流耦合數值模型，其中多孔介質滲流區采用達西定律描述、溶洞自由流區域采用N—S方程描述，兩區域間采用BJS邊界條件進行耦合，并基于局部非熱平衡理論，描述冷流體與熱基巖在縫洞界面上的熱量交換。通過與解析解和連續介質模型解對比，驗證了離散縫洞熱流耦合模型的準確性。

2）裂縫網絡連通性是控制和評價縫洞型熱儲流動傳熱效果的關鍵參數，而溶洞的存在對熱儲內的流動傳熱效果起重要影響。對于連通性較差的系統，流體流動速度小、等效滲透率低，熱采效果差；對于連通性好的系統，流體流動速度大、等效滲透率高，熱采效果好。溶洞主要通過兩個作用影響熱儲內的流體流動和傳熱：一是溶洞的存在會增加系統內貫穿的高速流動通道數量，甚至使系統從不連通變為連通；二是溶洞的存在會增大系統內局部流動通道速度。因此，對于基質滲透率較低，不存在連通或連通性較差的熱儲，能夠基于裂縫網絡連通性指標，指導縫洞型碳酸鹽巖熱儲進行酸化或者壓裂，形成具有開發價值的人工熱儲。另外，在酸化或壓裂造縫過程中，盡量使改造裂縫與溶洞相交，提高熱采效率。