裂縫－孔隙油藏蒸汽輔助重力泄油過程的自適應(yīng)網(wǎng)格數(shù)值計算

羅海山，王曉宏，施安峰

(中國科學(xué)技術(shù)大學(xué)熱科學(xué)和能源工程系，安徽合肥230026)

裂縫－孔隙油藏蒸汽輔助重力泄油過程的自適應(yīng)網(wǎng)格數(shù)值計算

羅海山，王曉宏，施安峰

(中國科學(xué)技術(shù)大學(xué)熱科學(xué)和能源工程系，安徽合肥230026)

為提高裂縫－孔隙雙重介質(zhì)油藏蒸汽輔助重力泄油(SAGD)過程數(shù)值模擬的計算速度，根據(jù)溫度、各相飽和度隨空間變化的快慢在不同區(qū)域采用不同尺度網(wǎng)格，提出雙孔雙滲模型下的自適應(yīng)網(wǎng)格算法及裂縫－基質(zhì)方程的解耦方法。對給出的算例分別采用全精細網(wǎng)格及自適應(yīng)網(wǎng)格算法進行計算。結(jié)果表明:自適應(yīng)網(wǎng)格法處理這類帶峰面非線性問題有很好的計算精度和計算速度;裂縫－孔隙雙重介質(zhì)油藏SAGD開采可能較早出現(xiàn)蒸汽竄頂這一不利于蒸汽能量利用率的現(xiàn)象。

蒸汽輔助重力驅(qū)油;稠油;滲流;多相流;自適應(yīng)網(wǎng)格法;裂縫;雙重介質(zhì)

對于稠油油藏蒸汽輔助重力泄油(SAGD)熱采過程，油藏中存在著一個不斷膨脹的蒸汽腔，其界面附近相變非常劇烈，溫度和各相飽和度梯度很大［1］，在變化劇烈的峰面附近，通常需要非常精細的網(wǎng)格，如果對全場都采用均勻精細網(wǎng)格，將需要很長的計算時間。為了提高計算速度，僅在需要的區(qū)域采用精細網(wǎng)格，而在其他區(qū)域采用較粗網(wǎng)格的自適應(yīng)網(wǎng)格法得到了廣泛的應(yīng)用［2-4］，并逐步推廣到多孔介質(zhì)多相流及油藏數(shù)值模擬領(lǐng)域［5-9］。筆者對裂縫-孔隙雙重介質(zhì)稠油油藏SAGD過程的自適應(yīng)網(wǎng)格計算方法進行研究，提出雙孔雙滲模型下的自適應(yīng)網(wǎng)格算法及裂縫-基質(zhì)方程解耦方法，以提高數(shù)值模擬的計算效率。

1 模型與控制方程

裂縫－孔隙雙重介質(zhì)由高滲透率低體積比的裂縫和低滲透率高體積比的基質(zhì)巖石構(gòu)成。在大尺度下，為使問題得到簡化，通用的Warren-Root模型將雙重介質(zhì)描述成由正交的裂縫網(wǎng)絡(luò)構(gòu)成并將基質(zhì)分割開來［10］?；|(zhì)之間的流動通常很少，但基質(zhì)與裂縫之間在壓差作用下的流體交換(亦稱竄流)通常是裂縫－孔隙雙重介質(zhì)油藏中石油被采出的主要原因。根據(jù)基質(zhì)之間是否忽略流動和熱傳導(dǎo)，W-R模型又可分為雙孔單滲模型［10-11］和雙孔雙滲模型［12］。由于雙孔單滲模型忽略了基質(zhì)巖石之間的熱傳導(dǎo)，考慮到巖石的熱傳導(dǎo)是SAGD開采的一個重要因素，因此本文中使用雙孔雙滲模型。

裂縫和基質(zhì)的水組分質(zhì)量守衡方程為

裂縫和基質(zhì)的油組分質(zhì)量守衡方程為

裂縫和基質(zhì)的能量守衡方程:

其中

式中，f和m分別表示裂縫和基質(zhì);下角標w、g和o分別表示水、蒸汽和稠油;φ(β)為孔隙度;ρ、H、U分別為密度、焓和內(nèi)能。

α相流體流動速度由達西定律確定為

式中，k和Krα分別為絕對和相對滲透率;對于稠油，黏度μα隨溫度變化非?？?。

qα和E分別表示裂縫與基質(zhì)之間的竄流量和熱交換［10－14］，其表達式為

在熱平衡方程(3)中，^Uα和^Hα取上游方向的值，即如果

相類似，方程(5)中“upstream”表示上游方向的值。σ是基質(zhì)巖石的形狀因子，可以根據(jù)基質(zhì)巖石的特征長度L*估算為

根據(jù)Kazemi公式［15］，常數(shù)γ可簡單取為12。井口源項Qα由經(jīng)典的Peaceman公式［16］求得。

裂縫和基質(zhì)中熱傳導(dǎo)系數(shù)為

油、水和氣等各相壓力之間的關(guān)系滿足

其中

式中，pcow為油水毛管力;pcog為油氣毛管力;S(β)l為液相總飽和度。

當水和蒸汽共存時，多孔介質(zhì)內(nèi)飽和蒸汽壓力－溫度關(guān)系［17］滿足

其中，氣體常數(shù)R=8.31 J/(K·mol)，水－汽毛管力

2 自適應(yīng)網(wǎng)格法

自適應(yīng)網(wǎng)格法包括網(wǎng)格結(jié)構(gòu)、網(wǎng)格細化和粗化、網(wǎng)格細分準則以及自適應(yīng)網(wǎng)格有限差分方程的求解步驟。

2.1 網(wǎng)格結(jié)構(gòu)

自適應(yīng)網(wǎng)格方法(AMR)是將整個流場的網(wǎng)格結(jié)構(gòu)劃分為若干層，層次由最粗層(k=1)到最細層(k=kmax)。相鄰層的空間步長相差2倍或其他倍數(shù)(本文采用2倍)。對于均勻裂縫－孔隙雙重介質(zhì)的自適應(yīng)網(wǎng)格法的流場計算，為使計算簡便，將裂縫網(wǎng)格結(jié)構(gòu)與基質(zhì)網(wǎng)格結(jié)構(gòu)綁定在一起，讓它們的網(wǎng)格結(jié)構(gòu)完全重合。

2.2 網(wǎng)格的細化與網(wǎng)格細分準則

為動態(tài)刻畫峰面位置，須根據(jù)油藏物理量的變化情況判斷是否需要重新劃分網(wǎng)格結(jié)構(gòu)，這就涉及到了網(wǎng)格的細化與網(wǎng)格細分準則。網(wǎng)格的細化是采用線性或高階曲線插值以確定粗網(wǎng)格所覆蓋的最精細網(wǎng)格上的物理量。本文討論的SAGD過程問題，油藏中存在溫度和各相飽和度峰面，則網(wǎng)格細分準則為對溫度和各相飽和度設(shè)置一個閾值(εT和εSα)，如果任一粗網(wǎng)格內(nèi)溫度和各相飽和度的變化大于閾值，網(wǎng)格將予以細分。由于本文研究的介質(zhì)含有裂縫和基質(zhì)兩個系統(tǒng)，因此設(shè)置的閾值同時包括裂縫與基質(zhì)的溫度和各相飽和度。

2.3 網(wǎng)格的粗化

新的網(wǎng)格結(jié)構(gòu)確定后，需要確定新出現(xiàn)的粗網(wǎng)格上的物理量。對于新被激活的粗網(wǎng)格，由其覆蓋的所有最精細網(wǎng)格上物理值求平均得到粗網(wǎng)格上的值(針對裂縫－孔隙油藏)。在本文中，孔隙度、壓力和溫度取算術(shù)平均，飽和度取加權(quán)算術(shù)平均(權(quán)值是孔隙度)，滲透率和熱傳導(dǎo)系數(shù)取調(diào)和平均。

2.4 自適應(yīng)網(wǎng)格有限差分方程的求解步驟

自適應(yīng)網(wǎng)格算法中，不同層次的時間步長是不同的。讓不同層之間的Courant數(shù)νk=Δtk/Δxk相等，即時間步長與空間步長成正比［2］。這樣，越粗的網(wǎng)格其時間步長越大，最粗層的時間步長是最細層時間步長的2n－1(n是總層數(shù))倍。完整的一輪計算步驟從最細層到最粗層逐層推進，見圖1。計算完最粗層后該輪計算結(jié)束。此后，須判斷網(wǎng)格是否需要重新劃分，在確定好網(wǎng)格結(jié)構(gòu)后，進入下一輪計算，如此循環(huán)。

圖1 自適應(yīng)網(wǎng)格算法各層計算步驟Fig.1 Adaptive mesh refinement(AMR)calculation procedure for advancing solution on different level grids

3 方程的解耦

3.1 方程離散

在自適應(yīng)網(wǎng)格算法中，只有被激活的網(wǎng)格參與計算，被激活網(wǎng)格的父網(wǎng)格或子網(wǎng)格不直接參與計算。計算是分層求解的，當計算某一層次時，假設(shè)當前的層次下被激活網(wǎng)格共有n個，使用有限體積法對方程組(1)～(3)進行全隱式離散，由于含有裂縫和基質(zhì)兩個系統(tǒng)，因此可分別得到兩個3n階的非線性離散方程組，其數(shù)學(xué)形式簡化為

式中，F(xiàn)1對應(yīng)裂縫方程組;F2對應(yīng)基質(zhì)方程組;Xf和Xm分別為裂縫方程組和基質(zhì)方程組的待求主變量。

3.2 雙孔雙滲模型下裂縫－基質(zhì)方程的解耦

對非線性方程組(13)和(14)使用牛頓迭代法，分別得到線性方程組

在雙孔單滲模型下，由于忽略了相鄰基質(zhì)之間的相互影響，根據(jù)方程組(16)，任一被激活網(wǎng)格(用i標號)對應(yīng)的基質(zhì)方程組的主變量可用其對應(yīng)的裂縫網(wǎng)格方程組的主變量表示為

將式(17)帶入式(15)即可消去所有基質(zhì)變量，從而只須求解一個3n階的線性方程組。這樣裂縫－基質(zhì)方程組的解耦方法減少了一半的變量，大大減少了矩陣運算計算量［11］。然而對于雙孔雙滲模型，相鄰基質(zhì)網(wǎng)格之間存在著關(guān)聯(lián)，為將雙孔單滲模型下成立的方程組解耦方法推廣到雙孔雙滲模型中，在牛頓迭代過程中讓基質(zhì)方程組(17)中的對流項與熱傳導(dǎo)項在迭代計算過程中都暫時“顯式化”，即使用上一迭代步的值，從而實現(xiàn)裂縫－基質(zhì)方程組的解耦。

4 算例

考慮一個長300 m、寬64 m、高32 m的裂縫－孔隙雙重介質(zhì)稠油油藏，長、寬、高方向分別用Z、X、Y表示。在油藏中間距離底部5 m的位置沿Z方向打一口300 m長的水平采油井，在采油井正上方5 m處平行打一口等長的注汽井。油藏兩邊取周期邊界條件，上、下取無流邊界條件(但不絕熱)。油藏首先進行蒸汽吞吐，當注氣井與生產(chǎn)井連通后開始進入SAGD過程。假設(shè)蒸汽在注氣井內(nèi)流動沿程沒有損失，且油藏在Z方向上所有橫截面沒有區(qū)別，則該三維油藏模型可以簡化為X、Y方向的二維油藏模型。油藏的參數(shù)如下:注蒸汽速度Qg(in)=100 t/d;注入蒸汽溫度Tin=523 K;注入蒸汽干度95%;油藏初始溫度T|t=0=320 K;油藏頂部初始壓力產(chǎn)油井井底壓力pwell=2 MPa;裂縫孔隙度φ0=0.01;基質(zhì)孔隙度φ0=0.1;裂縫絕對滲透率10－13m2;基質(zhì)絕對滲透率K(m)=1.97×10－15m2;裂縫初始各相飽和度|t=0=0;基質(zhì)初始各相飽和度裂縫壓縮系數(shù)10－8Pa－1;基質(zhì)壓縮系數(shù)稠油比熱容cpo=2.0 kJ/(kg·℃);基質(zhì)巖石比熱容ρRcpR=2.39 MJ/(m3·℃)。

通常，裂縫的毛管力很小［14］，各相相對滲透率可近似等于各相飽和度［18］。本文中，裂縫內(nèi)毛管力忽略不計。對于裂縫－孔隙雙重介質(zhì)流體的流動，基質(zhì)內(nèi)的毛管力卻是一個很重要的角色［14］。這是由于基質(zhì)毛管力通常比較大，會造成基質(zhì)各相流體壓力的較大不同。只要基質(zhì)巖石是親油的，基質(zhì)就會在毛管力作用下吸入裂縫中的水并向裂縫排出油，這是裂縫－孔隙雙重介質(zhì)油藏中重要的驅(qū)油機制之一。使用經(jīng)典的Brooks-Corey模型［19］計算基質(zhì)內(nèi)各相流體之間的毛管力和相對滲透率曲線，束縛水飽和度、束縛氣飽和度、殘余油飽和度分別取為0.2、0、0.2;啟動油水毛管力、啟動油氣毛管力分別取20、50 kPa。由于是多相流，油相相對滲透率由Stone II模型［20］得到。

油藏模型方程離散所用的最精細網(wǎng)格邊長為1 m，最精細網(wǎng)格對應(yīng)的時間步長最高取0.5 d。自適應(yīng)網(wǎng)格總層次取4，網(wǎng)格結(jié)構(gòu)的劃分準則使用裂縫和基質(zhì)的溫度和各相飽和度，閾值均取εT=5 K和εSα=0.05。采用預(yù)處理共軛梯度法求解代數(shù)方程組。圖2為全精細網(wǎng)格與自適應(yīng)網(wǎng)格法同一時刻的計算結(jié)果對比。

從圖3可以看出，產(chǎn)油量在穩(wěn)定上升，累積油汽比約為0.17。圖2和圖3的計算結(jié)果對比顯示，采用全精細網(wǎng)格求解與自適應(yīng)網(wǎng)格法求解的結(jié)果符合得很好。自適應(yīng)網(wǎng)格結(jié)構(gòu)隨峰面位置而變。圖4中給出第360 d時的自適應(yīng)網(wǎng)格結(jié)構(gòu)。

由圖4可見，在遠離峰面的位置都是粗網(wǎng)格，而在已經(jīng)被峰面掃過的區(qū)域其細網(wǎng)格可以逐漸恢復(fù)成粗網(wǎng)格。

在同一臺微機(CPU主頻奔四1.6 GHz，內(nèi)存2G)上，當程序計算到油藏SAGD開采的第360 d時，全精細網(wǎng)格的運算耗費時間為25 min，而自適應(yīng)網(wǎng)格法僅用了3 min，自適應(yīng)網(wǎng)格法的計算速度提高了7倍多。本文中只使用了64×32個網(wǎng)格數(shù)量即達到了8倍以上的計算速度，如果采用更多的網(wǎng)格數(shù)量或使用三維模型計算，則計算速度的提高將更加顯著。

圖2 第360 d時兩種方法計算的油藏裂縫溫度、壓力、蒸汽飽和度和基質(zhì)油飽和度等值線對比Fig.2 Comparison of numerical results between AMR solution and fine-grid solution at 360th day.

另外，從圖3還可以看出，在第500 d附近存在一個拐點。在這之前累積油汽比在緩慢上升，但之后開始緩慢下降。累積油汽比的上升說明蒸汽能量的利用率可能出現(xiàn)了變化。對結(jié)果進行分析，發(fā)現(xiàn)第500 d前后正好是蒸汽腔開始觸及油藏頂部的時刻，蒸汽腔的膨脹由rising方式變成spreading方式［21］，之后蒸汽腔與油藏頂部的接觸面積逐漸擴大。圖5顯示了計算到第1 000 d時裂縫蒸汽飽和度的等值線。

圖3 累積產(chǎn)油量和累積油汽比對比Fig.3 Comparison of cumulative oil production and cumulative oil-steam ratio

圖4 第360 d的自適應(yīng)網(wǎng)格結(jié)構(gòu)Fig.4 AMR grid structure at 360th day

圖5 第1000 d裂縫蒸汽飽和度等值線Fig.5 Fracture steam saturation contour line at 1000th day

從圖5中可觀察到，蒸汽腔已經(jīng)竄頂，熱量從頂部外泄到外界，形成熱損失，而油藏側(cè)向的膨脹程度不高，油藏許多位置還沒被開采。這些現(xiàn)象對采油的經(jīng)濟效益將有負面影響。

5 結(jié)束語

裂縫－孔隙雙重介質(zhì)稠油油藏的SAGD開采是帶峰面的強非線性問題。提出了雙孔雙滲模型解耦的方法。通過對算例的數(shù)值求解，將全精細網(wǎng)格與自適應(yīng)網(wǎng)格法的計算結(jié)果進行了對比，顯示出自適應(yīng)網(wǎng)格法很好的計算速度和計算精度。在裂縫－孔隙雙重介質(zhì)油藏實施SAGD開采可能較早出現(xiàn)蒸汽竄頂這一不利于注入蒸汽的能量利用率現(xiàn)象。

［1］BUTLER R M，STEPHENS D J.The gravity drainage of steam heated to parallel horizontal wells［J］.JCPT，1981，90-96.

［2］BERGER M J，OLIGER J.Adaptive mesh refinement for hyperbolic partial differential equations［J］.J Comput Phys，1984，53:484-512.

［3］BERGER M J，COLELLA P.Local adaptive mesh refinement for shock hydrodynamics［J］.J Comput Phys，1989，82:64-84.

［4］VENUTURUMILLI R，CHEN L D.Numerical simulation using adaptive mesh refinement for laminar jet diffusion flames［J］.Numer Heat Transfer B，2004，46(2):101-120.

［5］TRANGENSTEIN J A.Multi-scale iterative techniques and adaptive mesh refinement for flow in porous media［J］.Adv in Water Resources，2002，25:1175-1213.

［6］NILSSON J，GERRITSEN M，YOUNIS R.Towards an adaptive，high-resolution simulator for steam injection processes［R］.SPE 93881，2005.

［7］韓大匡，陳欽雷，閆存章.油藏數(shù)值模擬基礎(chǔ)［M］.北京:石油工業(yè)出版社，1993:254-259.

［8］譚未一，劉福平，李瑞忠，等.滲流方程的滲透率自適應(yīng)權(quán)重網(wǎng)格粗化算法［J］.石油大學(xué)學(xué)報:自然科學(xué)版，2004，28(4):52-55.TAN Wei-yi，LIU Fu-ping，LI Rui-zhong，et al.Adaptive grid roughening algorithm for average permeability of fluid equations［J］.Journal of the University of Petroleum，China(Edition of Natural Science)，2004，28(4):52-55.

［9］WANG X H，QUINTARD M，DARCHE G.Adaptive mesh refinement for one dimensional three-phase flow with phase change in porous media［J］.Numer Heat Transfer B，2006，50:231-268.

［10］WARREN J E，ROOT P J.The behavior of naturally fractured reservoirs［J］.SPE J，1963，3:245-255.

［11］THOMAS L K，DIXON T N，PIERSON R G.Fractured reservoir simulation［J］.SPE J，1983，23:42-54.

［12］HILL A C，THOMAS G W.A new approach for simulating complex fractured reservoirs［C］.SPE Middle East Oil Tech Conference and Exhibition，Bahrain，11-14 March，c1985.

［13］PRUESS K，NARASIMHAN T N.A practical method for modelling fluid and heat flow in fractured porous media［J］.SPE J，1985，25(1):14-26.

［14］OBALLA V，COOMBE D A，BUCHANAN W L.Fac-tors affecting the thermal response of naturally fractured reservoirs［J］.The Journal of Canadian Petroleum Technology，1993，32(8):31-42.

［15］KAZEMI K，MERRILL L S JR，PORTERFIELD K P，et al.Simulation of water-oil flow in naturally fractured reservoirs［J］.SPE J，1976，16(6):317-326.

［16］PEACEMAN D W.Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability［R］.SPE 10528，1982.

［17］UDELL K S.Heat transfer in porous media heated from above with evaporation，condensation，and capillary effects［J］.Journal of Heat Transfer，1983，105:485-492.

［18］HONARPOUR，M M，KOEDERITZ F，HERBERT A.Relative permeability of petroleum reservoirs［M］.Boca Raton，F(xiàn)L:CRC Press Inc，1986:41.

［19］BROOKS R H，COREY A T.Properties of porous media affecting fluid flow［J］.Journal of Irrigation and Drainage Division，Proceedings of the American Society of Civil Engineers，1966，92(IR2):61-88.

［20］STONE H L.Estimation of three-phase relative permeability and residual oil data［J］.Journal of Petroleum Technology，1973，12(4):52-61.

［21］RENARD P，MARSILY G D.Calculating the equivalent permeability:a review［J］.Advances in Water Resources，1987，20:253-278.

［22］CHOW L，BUTLER R M.Numerical simulation of the steam-assisted gravity drainage process(SAGD)［J］.JCPT，1996，35:55-62.

(編輯 沈玉英)

Adaptive mesh refinement technique for numerical simulation of SAGD process in fracture-pore reservoirs

LUO Hai-shan，WANG Xiao-hong，SHI An-feng
(Department of Thermal Science and Energy Engineering，University of Science and Technology of China，Hefei 230026，China)

To improve the numerical simulation calculation speed for steam assisted gravity drainage(SAGD)process in fracture-pore dual-porosity reservoirs，the adaptive mesh refinement(AMR)algorithm as well as the fracture-matrix equation decomposition approach under a dual-porosity(DK)model was proposed by using different grid sizes in different regions based upon the spatial variations of the temperature and phase-saturations.The numerical results show that the proposed AMR technique dealing with such nonlinear problems with the front is fast and can give good precision compared with the fine-grid solution.It is also shown that the steam chamber could reach the top of the reservoir，which is harmful to the thermal efficiency of steam energy for the SAGD process in fracture-pore dual-porosity reservoirs.

steam-assisted gravity drainage;heavy oil;seepage;multiphase flow;adaptive mesh refinement;fracture;dual-porosity media

TE 345

A

10.3969/j.issn.1673-5005.2011.01.028

2010－06－22

中國科學(xué)院知識創(chuàng)新工程重大項目(KJCX1－YW－21);國家“863”計劃項目(2006AA06Z228);國家自然科學(xué)基金項目(50674081)

羅海山(1982－)，男(漢族)，福建龍巖人，博士研究生，主要研究方向為多孔介質(zhì)中流動理論和數(shù)值計算以及蒸汽熱采稠油數(shù)值計算的自適應(yīng)網(wǎng)格法。

1673-5005(2011)01-0140-06