羅海山,王曉宏,施安峰
(中國(guó)科學(xué)技術(shù)大學(xué)熱科學(xué)和能源工程系,安徽合肥230026)
裂縫-孔隙油藏蒸汽輔助重力泄油過(guò)程的自適應(yīng)網(wǎng)格數(shù)值計(jì)算
羅海山,王曉宏,施安峰
(中國(guó)科學(xué)技術(shù)大學(xué)熱科學(xué)和能源工程系,安徽合肥230026)
為提高裂縫-孔隙雙重介質(zhì)油藏蒸汽輔助重力泄油(SAGD)過(guò)程數(shù)值模擬的計(jì)算速度,根據(jù)溫度、各相飽和度隨空間變化的快慢在不同區(qū)域采用不同尺度網(wǎng)格,提出雙孔雙滲模型下的自適應(yīng)網(wǎng)格算法及裂縫-基質(zhì)方程的解耦方法。對(duì)給出的算例分別采用全精細(xì)網(wǎng)格及自適應(yīng)網(wǎng)格算法進(jìn)行計(jì)算。結(jié)果表明:自適應(yīng)網(wǎng)格法處理這類帶峰面非線性問(wèn)題有很好的計(jì)算精度和計(jì)算速度;裂縫-孔隙雙重介質(zhì)油藏SAGD開(kāi)采可能較早出現(xiàn)蒸汽竄頂這一不利于蒸汽能量利用率的現(xiàn)象。
蒸汽輔助重力驅(qū)油;稠油;滲流;多相流;自適應(yīng)網(wǎng)格法;裂縫;雙重介質(zhì)
對(duì)于稠油油藏蒸汽輔助重力泄油(SAGD)熱采過(guò)程,油藏中存在著一個(gè)不斷膨脹的蒸汽腔,其界面附近相變非常劇烈,溫度和各相飽和度梯度很大[1],在變化劇烈的峰面附近,通常需要非常精細(xì)的網(wǎng)格,如果對(duì)全場(chǎng)都采用均勻精細(xì)網(wǎng)格,將需要很長(zhǎng)的計(jì)算時(shí)間。為了提高計(jì)算速度,僅在需要的區(qū)域采用精細(xì)網(wǎng)格,而在其他區(qū)域采用較粗網(wǎng)格的自適應(yīng)網(wǎng)格法得到了廣泛的應(yīng)用[2-4],并逐步推廣到多孔介質(zhì)多相流及油藏?cái)?shù)值模擬領(lǐng)域[5-9]。筆者對(duì)裂縫-孔隙雙重介質(zhì)稠油油藏SAGD過(guò)程的自適應(yīng)網(wǎng)格計(jì)算方法進(jìn)行研究,提出雙孔雙滲模型下的自適應(yīng)網(wǎng)格算法及裂縫-基質(zhì)方程解耦方法,以提高數(shù)值模擬的計(jì)算效率。
裂縫-孔隙雙重介質(zhì)由高滲透率低體積比的裂縫和低滲透率高體積比的基質(zhì)巖石構(gòu)成。在大尺度下,為使問(wèn)題得到簡(jiǎn)化,通用的Warren-Root模型將雙重介質(zhì)描述成由正交的裂縫網(wǎng)絡(luò)構(gòu)成并將基質(zhì)分割開(kāi)來(lái)[10]。基質(zhì)之間的流動(dòng)通常很少,但基質(zhì)與裂縫之間在壓差作用下的流體交換(亦稱竄流)通常是裂縫-孔隙雙重介質(zhì)油藏中石油被采出的主要原因。根據(jù)基質(zhì)之間是否忽略流動(dòng)和熱傳導(dǎo),W-R模型又可分為雙孔單滲模型[10-11]和雙孔雙滲模型[12]。由于雙孔單滲模型忽略了基質(zhì)巖石之間的熱傳導(dǎo),考慮到巖石的熱傳導(dǎo)是SAGD開(kāi)采的一個(gè)重要因素,因此本文中使用雙孔雙滲模型。
裂縫和基質(zhì)的水組分質(zhì)量守衡方程為
裂縫和基質(zhì)的油組分質(zhì)量守衡方程為
裂縫和基質(zhì)的能量守衡方程:
其中
式中,f和m分別表示裂縫和基質(zhì);下角標(biāo)w、g和o分別表示水、蒸汽和稠油;φ(β)為孔隙度;ρ、H、U分別為密度、焓和內(nèi)能。
α相流體流動(dòng)速度由達(dá)西定律確定為
式中,k和Krα分別為絕對(duì)和相對(duì)滲透率;對(duì)于稠油,黏度μα隨溫度變化非??臁?/p>
qα和E分別表示裂縫與基質(zhì)之間的竄流量和熱交換[10-14],其表達(dá)式為
在熱平衡方程(3)中,^Uα和^Hα取上游方向的值,即如果
相類似,方程(5)中“upstream”表示上游方向的值。σ是基質(zhì)巖石的形狀因子,可以根據(jù)基質(zhì)巖石的特征長(zhǎng)度L*估算為
根據(jù)Kazemi公式[15],常數(shù)γ可簡(jiǎn)單取為12。井口源項(xiàng)Qα由經(jīng)典的Peaceman公式[16]求得。
裂縫和基質(zhì)中熱傳導(dǎo)系數(shù)為
油、水和氣等各相壓力之間的關(guān)系滿足
其中
式中,pcow為油水毛管力;pcog為油氣毛管力;S(β)l為液相總飽和度。
當(dāng)水和蒸汽共存時(shí),多孔介質(zhì)內(nèi)飽和蒸汽壓力-溫度關(guān)系[17]滿足
其中,氣體常數(shù)R=8.31 J/(K·mol),水-汽毛管力
自適應(yīng)網(wǎng)格法包括網(wǎng)格結(jié)構(gòu)、網(wǎng)格細(xì)化和粗化、網(wǎng)格細(xì)分準(zhǔn)則以及自適應(yīng)網(wǎng)格有限差分方程的求解步驟。
自適應(yīng)網(wǎng)格方法(AMR)是將整個(gè)流場(chǎng)的網(wǎng)格結(jié)構(gòu)劃分為若干層,層次由最粗層(k=1)到最細(xì)層(k=kmax)。相鄰層的空間步長(zhǎng)相差2倍或其他倍數(shù)(本文采用2倍)。對(duì)于均勻裂縫-孔隙雙重介質(zhì)的自適應(yīng)網(wǎng)格法的流場(chǎng)計(jì)算,為使計(jì)算簡(jiǎn)便,將裂縫網(wǎng)格結(jié)構(gòu)與基質(zhì)網(wǎng)格結(jié)構(gòu)綁定在一起,讓它們的網(wǎng)格結(jié)構(gòu)完全重合。
為動(dòng)態(tài)刻畫峰面位置,須根據(jù)油藏物理量的變化情況判斷是否需要重新劃分網(wǎng)格結(jié)構(gòu),這就涉及到了網(wǎng)格的細(xì)化與網(wǎng)格細(xì)分準(zhǔn)則。網(wǎng)格的細(xì)化是采用線性或高階曲線插值以確定粗網(wǎng)格所覆蓋的最精細(xì)網(wǎng)格上的物理量。本文討論的SAGD過(guò)程問(wèn)題,油藏中存在溫度和各相飽和度峰面,則網(wǎng)格細(xì)分準(zhǔn)則為對(duì)溫度和各相飽和度設(shè)置一個(gè)閾值(εT和εSα),如果任一粗網(wǎng)格內(nèi)溫度和各相飽和度的變化大于閾值,網(wǎng)格將予以細(xì)分。由于本文研究的介質(zhì)含有裂縫和基質(zhì)兩個(gè)系統(tǒng),因此設(shè)置的閾值同時(shí)包括裂縫與基質(zhì)的溫度和各相飽和度。
新的網(wǎng)格結(jié)構(gòu)確定后,需要確定新出現(xiàn)的粗網(wǎng)格上的物理量。對(duì)于新被激活的粗網(wǎng)格,由其覆蓋的所有最精細(xì)網(wǎng)格上物理值求平均得到粗網(wǎng)格上的值(針對(duì)裂縫-孔隙油藏)。在本文中,孔隙度、壓力和溫度取算術(shù)平均,飽和度取加權(quán)算術(shù)平均(權(quán)值是孔隙度),滲透率和熱傳導(dǎo)系數(shù)取調(diào)和平均。
自適應(yīng)網(wǎng)格算法中,不同層次的時(shí)間步長(zhǎng)是不同的。讓不同層之間的Courant數(shù)νk=Δtk/Δxk相等,即時(shí)間步長(zhǎng)與空間步長(zhǎng)成正比[2]。這樣,越粗的網(wǎng)格其時(shí)間步長(zhǎng)越大,最粗層的時(shí)間步長(zhǎng)是最細(xì)層時(shí)間步長(zhǎng)的2n-1(n是總層數(shù))倍。完整的一輪計(jì)算步驟從最細(xì)層到最粗層逐層推進(jìn),見(jiàn)圖1。計(jì)算完最粗層后該輪計(jì)算結(jié)束。此后,須判斷網(wǎng)格是否需要重新劃分,在確定好網(wǎng)格結(jié)構(gòu)后,進(jìn)入下一輪計(jì)算,如此循環(huán)。
圖1 自適應(yīng)網(wǎng)格算法各層計(jì)算步驟Fig.1 Adaptive mesh refinement(AMR)calculation procedure for advancing solution on different level grids
在自適應(yīng)網(wǎng)格算法中,只有被激活的網(wǎng)格參與計(jì)算,被激活網(wǎng)格的父網(wǎng)格或子網(wǎng)格不直接參與計(jì)算。計(jì)算是分層求解的,當(dāng)計(jì)算某一層次時(shí),假設(shè)當(dāng)前的層次下被激活網(wǎng)格共有n個(gè),使用有限體積法對(duì)方程組(1)~(3)進(jìn)行全隱式離散,由于含有裂縫和基質(zhì)兩個(gè)系統(tǒng),因此可分別得到兩個(gè)3n階的非線性離散方程組,其數(shù)學(xué)形式簡(jiǎn)化為
式中,F(xiàn)1對(duì)應(yīng)裂縫方程組;F2對(duì)應(yīng)基質(zhì)方程組;Xf和Xm分別為裂縫方程組和基質(zhì)方程組的待求主變量。
對(duì)非線性方程組(13)和(14)使用牛頓迭代法,分別得到線性方程組
在雙孔單滲模型下,由于忽略了相鄰基質(zhì)之間的相互影響,根據(jù)方程組(16),任一被激活網(wǎng)格(用i標(biāo)號(hào))對(duì)應(yīng)的基質(zhì)方程組的主變量可用其對(duì)應(yīng)的裂縫網(wǎng)格方程組的主變量表示為
將式(17)帶入式(15)即可消去所有基質(zhì)變量,從而只須求解一個(gè)3n階的線性方程組。這樣裂縫-基質(zhì)方程組的解耦方法減少了一半的變量,大大減少了矩陣運(yùn)算計(jì)算量[11]。然而對(duì)于雙孔雙滲模型,相鄰基質(zhì)網(wǎng)格之間存在著關(guān)聯(lián),為將雙孔單滲模型下成立的方程組解耦方法推廣到雙孔雙滲模型中,在牛頓迭代過(guò)程中讓基質(zhì)方程組(17)中的對(duì)流項(xiàng)與熱傳導(dǎo)項(xiàng)在迭代計(jì)算過(guò)程中都暫時(shí)“顯式化”,即使用上一迭代步的值,從而實(shí)現(xiàn)裂縫-基質(zhì)方程組的解耦。
考慮一個(gè)長(zhǎng)300 m、寬64 m、高32 m的裂縫-孔隙雙重介質(zhì)稠油油藏,長(zhǎng)、寬、高方向分別用Z、X、Y表示。在油藏中間距離底部5 m的位置沿Z方向打一口300 m長(zhǎng)的水平采油井,在采油井正上方5 m處平行打一口等長(zhǎng)的注汽井。油藏兩邊取周期邊界條件,上、下取無(wú)流邊界條件(但不絕熱)。油藏首先進(jìn)行蒸汽吞吐,當(dāng)注氣井與生產(chǎn)井連通后開(kāi)始進(jìn)入SAGD過(guò)程。假設(shè)蒸汽在注氣井內(nèi)流動(dòng)沿程沒(méi)有損失,且油藏在Z方向上所有橫截面沒(méi)有區(qū)別,則該三維油藏模型可以簡(jiǎn)化為X、Y方向的二維油藏模型。油藏的參數(shù)如下:注蒸汽速度Qg(in)=100 t/d;注入蒸汽溫度Tin=523 K;注入蒸汽干度95%;油藏初始溫度T|t=0=320 K;油藏頂部初始?jí)毫Ξa(chǎn)油井井底壓力pwell=2 MPa;裂縫孔隙度φ0=0.01;基質(zhì)孔隙度φ0=0.1;裂縫絕對(duì)滲透率10-13m2;基質(zhì)絕對(duì)滲透率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],各相相對(duì)滲透率可近似等于各相飽和度[18]。本文中,裂縫內(nèi)毛管力忽略不計(jì)。對(duì)于裂縫-孔隙雙重介質(zhì)流體的流動(dòng),基質(zhì)內(nèi)的毛管力卻是一個(gè)很重要的角色[14]。這是由于基質(zhì)毛管力通常比較大,會(huì)造成基質(zhì)各相流體壓力的較大不同。只要基質(zhì)巖石是親油的,基質(zhì)就會(huì)在毛管力作用下吸入裂縫中的水并向裂縫排出油,這是裂縫-孔隙雙重介質(zhì)油藏中重要的驅(qū)油機(jī)制之一。使用經(jīng)典的Brooks-Corey模型[19]計(jì)算基質(zhì)內(nèi)各相流體之間的毛管力和相對(duì)滲透率曲線,束縛水飽和度、束縛氣飽和度、殘余油飽和度分別取為0.2、0、0.2;啟動(dòng)油水毛管力、啟動(dòng)油氣毛管力分別取20、50 kPa。由于是多相流,油相相對(duì)滲透率由Stone II模型[20]得到。
油藏模型方程離散所用的最精細(xì)網(wǎng)格邊長(zhǎng)為1 m,最精細(xì)網(wǎng)格對(duì)應(yīng)的時(shí)間步長(zhǎng)最高取0.5 d。自適應(yīng)網(wǎng)格總層次取4,網(wǎng)格結(jié)構(gòu)的劃分準(zhǔn)則使用裂縫和基質(zhì)的溫度和各相飽和度,閾值均取εT=5 K和εSα=0.05。采用預(yù)處理共軛梯度法求解代數(shù)方程組。圖2為全精細(xì)網(wǎng)格與自適應(yīng)網(wǎng)格法同一時(shí)刻的計(jì)算結(jié)果對(duì)比。
從圖3可以看出,產(chǎn)油量在穩(wěn)定上升,累積油汽比約為0.17。圖2和圖3的計(jì)算結(jié)果對(duì)比顯示,采用全精細(xì)網(wǎng)格求解與自適應(yīng)網(wǎng)格法求解的結(jié)果符合得很好。自適應(yīng)網(wǎng)格結(jié)構(gòu)隨峰面位置而變。圖4中給出第360 d時(shí)的自適應(yīng)網(wǎng)格結(jié)構(gòu)。
由圖4可見(jiàn),在遠(yuǎn)離峰面的位置都是粗網(wǎng)格,而在已經(jīng)被峰面掃過(guò)的區(qū)域其細(xì)網(wǎng)格可以逐漸恢復(fù)成粗網(wǎng)格。
在同一臺(tái)微機(jī)(CPU主頻奔四1.6 GHz,內(nèi)存2G)上,當(dāng)程序計(jì)算到油藏SAGD開(kāi)采的第360 d時(shí),全精細(xì)網(wǎng)格的運(yùn)算耗費(fèi)時(shí)間為25 min,而自適應(yīng)網(wǎng)格法僅用了3 min,自適應(yīng)網(wǎng)格法的計(jì)算速度提高了7倍多。本文中只使用了64×32個(gè)網(wǎng)格數(shù)量即達(dá)到了8倍以上的計(jì)算速度,如果采用更多的網(wǎng)格數(shù)量或使用三維模型計(jì)算,則計(jì)算速度的提高將更加顯著。
圖2 第360 d時(shí)兩種方法計(jì)算的油藏裂縫溫度、壓力、蒸汽飽和度和基質(zhì)油飽和度等值線對(duì)比Fig.2 Comparison of numerical results between AMR solution and fine-grid solution at 360th day.
另外,從圖3還可以看出,在第500 d附近存在一個(gè)拐點(diǎn)。在這之前累積油汽比在緩慢上升,但之后開(kāi)始緩慢下降。累積油汽比的上升說(shuō)明蒸汽能量的利用率可能出現(xiàn)了變化。對(duì)結(jié)果進(jìn)行分析,發(fā)現(xiàn)第500 d前后正好是蒸汽腔開(kāi)始觸及油藏頂部的時(shí)刻,蒸汽腔的膨脹由rising方式變成spreading方式[21],之后蒸汽腔與油藏頂部的接觸面積逐漸擴(kuò)大。圖5顯示了計(jì)算到第1 000 d時(shí)裂縫蒸汽飽和度的等值線。
圖3 累積產(chǎn)油量和累積油汽比對(duì)比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è)向的膨脹程度不高,油藏許多位置還沒(méi)被開(kāi)采。這些現(xiàn)象對(duì)采油的經(jīng)濟(jì)效益將有負(fù)面影響。
裂縫-孔隙雙重介質(zhì)稠油油藏的SAGD開(kāi)采是帶峰面的強(qiáng)非線性問(wèn)題。提出了雙孔雙滲模型解耦的方法。通過(guò)對(duì)算例的數(shù)值求解,將全精細(xì)網(wǎng)格與自適應(yīng)網(wǎng)格法的計(jì)算結(jié)果進(jìn)行了對(duì)比,顯示出自適應(yīng)網(wǎng)格法很好的計(jì)算速度和計(jì)算精度。在裂縫-孔隙雙重介質(zhì)油藏實(shí)施SAGD開(kāi)采可能較早出現(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]韓大匡,陳欽雷,閆存章.油藏?cái)?shù)值模擬基礎(chǔ)[M].北京:石油工業(yè)出版社,1993:254-259.
[8]譚未一,劉福平,李瑞忠,等.滲流方程的滲透率自適應(yīng)權(quán)重網(wǎng)格粗化算法[J].石油大學(xué)學(xué)報(bào):自然科學(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
中國(guó)科學(xué)院知識(shí)創(chuàng)新工程重大項(xiàng)目(KJCX1-YW-21);國(guó)家“863”計(jì)劃項(xiàng)目(2006AA06Z228);國(guó)家自然科學(xué)基金項(xiàng)目(50674081)
羅海山(1982-),男(漢族),福建龍巖人,博士研究生,主要研究方向?yàn)槎嗫捉橘|(zhì)中流動(dòng)理論和數(shù)值計(jì)算以及蒸汽熱采稠油數(shù)值計(jì)算的自適應(yīng)網(wǎng)格法。
1673-5005(2011)01-0140-06