基于Darcy-Stokes耦合模型的縫洞型介質(zhì)等效滲透率分析

李亞軍,姚 軍,黃朝琴,劉永輝

(中國石油大學(xué)石油工程學(xué)院,山東青島266555)

基于Darcy-Stokes耦合模型的縫洞型介質(zhì)等效滲透率分析

李亞軍,姚 軍,黃朝琴,劉永輝

(中國石油大學(xué)石油工程學(xué)院,山東青島266555)

根據(jù)縫洞型介質(zhì)的特點和流體在不同尺度空間的流動特征,建立Darcy-Stokes耦合數(shù)學(xué)模型,推導(dǎo)出考慮基巖滲透率的廣義立方定律和單裂縫多孔介質(zhì)的等效滲透率張量表達式,以此為基礎(chǔ)提出裂縫預(yù)處理技術(shù),將溶洞-裂縫-基巖耦合流動問題轉(zhuǎn)化為溶洞-基巖耦合流動問題,分析不同縫洞結(jié)構(gòu)介質(zhì)體的等效滲透率和滲透特性。結(jié)果表明:縫洞結(jié)構(gòu)的存在不同程度地影響介質(zhì)體的等效滲透率;含圓形溶洞介質(zhì)的等效滲透率與基巖滲透率成線性關(guān)系;裂縫預(yù)處理技術(shù)可以提高計算效率。

縫洞型介質(zhì);Darcy-Stokes耦合模型;連續(xù)介質(zhì);等效滲透率;流動模擬

縫洞型介質(zhì)是基巖孔隙、裂縫和溶洞交疊嵌套分布的多尺度復(fù)雜固體介質(zhì),在油氣田開發(fā)、水利科學(xué)和生物科學(xué)中,許多問題的研究對象都具有縫洞型介質(zhì)的特征[1-3]?？p洞型介質(zhì)具有多尺度性和非均質(zhì)性,其滲透率與多孔介質(zhì)相比存在較大差異?？p洞型介質(zhì)內(nèi)的流動狀態(tài)和流動規(guī)律異常復(fù)雜,縫洞內(nèi)的自由流動和基巖中的滲流共存[4],Saffman[5]在Beavers和Joseph[6]實驗研究的基礎(chǔ)上,從理論上驗證了其實驗結(jié)論,并通過理論修正得到了Beavers-Joseph-Saffman滑移速度邊界條件(簡稱BJS條件),得到廣泛應(yīng)用,在此基礎(chǔ)上形成了Darcy-Stokes耦合模型[7]。耦合模型的數(shù)值解法主要有單域法和兩域法,單域法是指在整個區(qū)域上取統(tǒng)一的速度有限元空間進行分析,如采用標準有限元法[8]或混合有限元法[9],兩域法是指在兩個區(qū)域分別取不同的速度有限元空間,如采用Galerkin法和混合有限元耦合的數(shù)值解法[10]。單域法容易進行理論和計算分析,但其有限元空間構(gòu)造復(fù)雜,不利于實際計算;兩域法在理論和數(shù)值分析時較困難。筆者根據(jù)縫洞型介質(zhì)的特點和流體在不同尺度空間的流動特征,建立Darcy-Stokes耦合數(shù)學(xué)模型,基于非協(xié)調(diào)Crouzeix-Raviart型三角形有限元格式建立模型的數(shù)值解法,分析不同縫洞結(jié)構(gòu)對介質(zhì)體等效滲透率和滲透特性的影響。

1 數(shù)學(xué)模型

如圖1(a)所示,將整個縫洞型介質(zhì)研究區(qū)域Ω分解為多孔介質(zhì)基巖區(qū)域Ωd和自由流動區(qū)域Ωs,其中Ωs包含裂縫區(qū)域Ωf和溶洞區(qū)域Ωv,即Ω=Ωd∪Ωs,Ωs=Ωf∪Ωv,各區(qū)域的邊界表示為?Ωb(b=d,s,f,v);正方形網(wǎng)格塊的邊界稱為外邊界Γj(j=1,2,3,4),nj為相應(yīng)外邊界Γj的單位外法向量。

根據(jù)等效原則和等效滲透率計算原理[11],縫洞型非均質(zhì)離散介質(zhì)體的等效滲透率(k)可由滿足等效原則的均質(zhì)各向異性連續(xù)介質(zhì)體(圖1(b))的滲透率表示。

圖1 縫洞型介質(zhì)等效示意圖Fig.1 Equivalent schematic diagram of fractured-vuggy media

1.1 控制方程

多孔介質(zhì)區(qū)域Ωd和自由流動區(qū)域Ωs的控制微分方程分別表示如下:

其中

式中,μ為流體動力黏度,Pa·s;km為多孔介質(zhì)滲透率,m2;vd為多孔介質(zhì)區(qū)域滲流速度,m/s;pd為多孔介質(zhì)表征單元體中的平均壓力,Pa;f為單位體積力,N/m3;q為源匯項;vs為自由流動區(qū)域流體速度,m/s;ps為自由流動區(qū)域流體壓力,Pa;Dvs為變形速率張量。

1.2 邊界條件

Darcy-Stokes耦合模型的關(guān)鍵在于建立多孔介質(zhì)區(qū)域與自由流動區(qū)域交界面Γ=?Ωd∩?Ωs耦合條件,采用BJS條件并結(jié)合法向流速連續(xù)條件和法向應(yīng)力連續(xù)條件,交界面Γ處的邊界條件表示為

其中,ns為?Ωs的單位外法向量;τ為?Ωs的單位切向量;α為滑移系數(shù),取決于交界面處的多孔介質(zhì)的幾何特征和結(jié)構(gòu)特征,為無量綱常數(shù)。

外邊界條件對于準確求解等效滲透率張量非常重要,采用周期邊界條件可以保證得到符合物理意義的滲透率張量[12]。整個縫洞型介質(zhì)研究區(qū)域Ω的周期邊界條件表示為

1.3 裂縫預(yù)處理方法

對于實際縫洞型介質(zhì)的流動,由于裂縫數(shù)量較多,且裂縫開度尺寸與研究區(qū)域的數(shù)量級相差較大,為得到較精確的解,裂縫區(qū)域剖分的網(wǎng)格須足夠密。針對水平單裂縫及其周圍多孔介質(zhì)(圖2),考慮基巖滲透率,得到廣義立方定律和單裂縫多孔介質(zhì)的等效滲透率,利用等效滲透率張量表征的連續(xù)介質(zhì)體代替原裂縫介質(zhì)進行計算。

圖2 單裂縫多孔介質(zhì)示意圖Fig.2 Schematic diagram of single fractured porousmedia

假定在系統(tǒng)橫向和縱向分別施加單位壓力梯度,則單位黏度流體通過系統(tǒng)后的速度的解析解分別為

其中,式(7)表示x方向裂縫內(nèi)的流體速度分布;式(8)表示x方向基巖內(nèi)的流體速度分布;式(9)表示y方向基巖內(nèi)的流體速度分布,即單裂縫多孔介質(zhì)y方向的平均速度。

根據(jù)平均速度計算公式[11],得到單裂縫多孔介質(zhì)在x方向的平均速度為

此即為考慮基巖滲透率的單裂縫多孔介質(zhì)的廣義立方定律。根據(jù)達西公式的形式,得單裂縫多孔介質(zhì)系統(tǒng)在x方向和y方向的等效滲透率kfx和kfy為

忽略式(11)右端第二項,則變成經(jīng)典連續(xù)條件(僅考慮基巖和溶洞交界面處的法向速度而忽略切向速度)[5]時的解析解。對裂縫方位角θ≠0的情況,則根據(jù)滲透率張量轉(zhuǎn)換公式[13]求得等效滲透率張量為

計算每條裂縫及其周圍多孔介質(zhì)區(qū)域的等效滲透率張量,利用等效滲透率張量表征的連續(xù)介質(zhì)體代替原裂縫介質(zhì),將溶洞-裂縫-基巖耦合流動問題轉(zhuǎn)化為溶洞-基巖耦合流動問題。

令Qd和Qs分別代表多孔介質(zhì)區(qū)域和自由流動區(qū)域的流動控制微分方程(1)和(2),則在原始縫洞型介質(zhì)研究區(qū)域上的積分表達式為

對裂縫進行簡化后,其積分表達式變?yōu)?/p>

基于非協(xié)調(diào)元中自由度最小的Crouzeix-Raviart型三角形元,采用穩(wěn)定化有限元法對Darcy-Stokes耦合問題進行數(shù)值求解,此方法與壓力常數(shù)單元結(jié)合自動滿足inf-sup條件和質(zhì)量守恒條件,Feng等[14]嚴格證明了該方法的穩(wěn)定性和收斂性,在L2范數(shù)下,速度的計算誤差階數(shù)為2,壓力的為1。

2 計算實例

2.1 含水平單裂縫介質(zhì)

如圖2,多孔介質(zhì)中含水平貫穿裂縫,lx=ly=8 cm,km=10×10-3μm2,α=1。利用所給出的數(shù)值解法和式(11)、(12)分別計算不同裂縫開度時裂縫介質(zhì)的等效滲透率張量,并與考慮經(jīng)典連續(xù)條件所計算的解析解進行對比,驗證數(shù)值算法的正確性。計算結(jié)果見表1。

表1 水平單裂縫介質(zhì)的等效滲透率主值Table 1 Principlal value of equivalent permeability of porous media containing single horizontal fracture

比較發(fā)現(xiàn):解析解與數(shù)值解基本一致,證明了數(shù)值解法的正確性;隨著裂縫開度的增大,kxx急劇增大,而kyy變化較為緩慢;經(jīng)典連續(xù)條件和BJS條件下kxx的解析結(jié)果顯示,是否考慮基巖和溶洞交界面處的切向速度對最終等效滲透率的計算結(jié)果影響不大,而在縫洞型介質(zhì)流動模擬的計算中,忽略切向速度將使問題容易處理;在確定垂直于裂縫延伸方向的等效滲透率系數(shù)時,相當于將裂縫視為具有無限滲透能力,其內(nèi)部流體為等勢體;裂縫開度越小,裂縫在整個流動區(qū)域所占比重越小,即純流體流動區(qū)域所占比例較小,則利用Darcy定理描述整個區(qū)域的流體流動誤差越小,數(shù)值計算結(jié)果越接近解析結(jié)果;當整個研究區(qū)域皆為自由流動區(qū)時,相當于利用立方定律描述其滲透性。

2.2 含圓形溶洞介質(zhì)

如圖3,邊長為1.5 cm的正方形區(qū)域中心為一半徑為rvug的圓形溶洞,周圍為滲透率為km的多孔介質(zhì)。由于對稱性,溶洞介質(zhì)的等效滲透率的主值相等,即kxx=kyy。不同尺寸溶洞下基巖滲透率與等效滲透率的關(guān)系見圖4,不同基巖滲透率時溶洞尺寸與等效滲透率的關(guān)系見圖5。

由圖4可看出,對相同結(jié)構(gòu)的溶洞介質(zhì),其等效滲透率隨基巖滲透率的增大而增大,并成線性關(guān)系,且溶洞半徑越大,直線的斜率越大。由圖5可看出:在基巖滲透率保持不變的情況下,當溶洞半徑相對較小時,其等效滲透率變化緩慢,而當溶洞尺寸接近外邊界時,等效滲透率快速增加;基巖滲透率較小時,溶洞的滲透性受到較大限制,與基巖滲透率較大的情形相比,等效滲透率變化較為緩慢。

圖5 溶洞尺寸與等效滲透率的關(guān)系Fig.5 Relation curve between equivalent permeability andrvug

2.3 復(fù)雜縫洞型介質(zhì)

如圖6,邊長l=12 cm的正方形縫洞型介質(zhì),基巖滲透率為100×10-3μm2,內(nèi)含8條裂縫fi(i=0,1,…,7)和3個溶洞,所有裂縫的開度h=1 000 μm,裂縫f0的傾角β=30°,方形溶洞的尺寸為2.5 cm×8 cm,圓形溶洞的半徑為1.5 cm。

圖6 復(fù)雜縫洞型介質(zhì)Fig.6 Complicated fractured-vuggy media

分別計算如下4種介質(zhì)系統(tǒng)的等效滲透率張量:(Ⅰ)介質(zhì)體內(nèi)僅含溶洞;(Ⅱ)含溶洞和裂縫f0;(Ⅲ)含溶洞和裂縫fi(i=0,1,2,3);(Ⅳ)含溶洞和所有裂縫fi(i=0,1,…,7)。對裂縫進行簡化處理時,將每條裂縫區(qū)域擴展至寬度為0.8 cm的含裂縫多孔介質(zhì)區(qū)域,將真實縫洞型介質(zhì)和裂縫簡化處理后的等效滲透率計算結(jié)果進行對比分析,結(jié)果見表2。

由表2可知:裂縫簡化處理后的計算結(jié)果和真實縫洞系統(tǒng)的等效滲透率基本一致,相對誤差均低于2%;經(jīng)過裂縫預(yù)處理,求解的網(wǎng)格單元數(shù)減少一半,節(jié)省了計算量,提高了計算效率;當縫洞在多孔介質(zhì)內(nèi)部離散分布時,縫洞型介質(zhì)的等效滲透率比基巖滲透率略有增加,與基巖滲透率的數(shù)量級相同;當縫洞貫穿整個研究區(qū)域并相互連通時,縫洞型介質(zhì)的等效滲透率提高了4個數(shù)量級,介質(zhì)體的滲透性大大改善。

表2 不同介質(zhì)系統(tǒng)的等效滲透率張量k及網(wǎng)格單元數(shù)nE比較Table 2 Comparison of equivalent permeability tensork and grid-element numbernEof different system s

3 結(jié) 論

(1)含圓形溶洞介質(zhì)的等效滲透率與基巖滲透率成線性關(guān)系。

(2)縫洞結(jié)構(gòu)的存在會在不同程度上影響介質(zhì)體的等效滲透率:嵌套于多孔介質(zhì)內(nèi)部且邊界遠離研究區(qū)域外邊界的縫洞不會顯著影響整體的滲透性,其數(shù)量級與基巖滲透率相同,且其內(nèi)部流體可視為等勢體;邊界靠近研究區(qū)域外邊界的縫洞以及貫穿整個研究區(qū)域并相互連通的縫洞,將會顯著提高整體的滲透性,其等效滲透率提高4個數(shù)量級。

(3)裂縫預(yù)處理技術(shù)將溶洞-裂縫-基巖耦合流動問題轉(zhuǎn)化為溶洞-基巖(自由流動和多孔介質(zhì)滲流)耦合流動問題,能保證計算精度,并且提高了計算效率。

[1] 姚軍,王子勝.縫洞型碳酸鹽巖油藏試井解釋理論與方法[M].東營:中國石油大學(xué)出版社,2007.

[2] 張志才,陳喜,石朋,等.喀斯特流域分布式水文模型及植被生態(tài)水文效應(yīng)[J].水科學(xué)進展,2009,20(6):806-811.ZHANG Zhi-cai,CHEN Xi,SH IPeng,et al.Distributed hydrological model and eco-hydrological effect of vegetation in Karst watershed[J].Advances inWater Science,2009,20(6):806-811.

[3] 宋付權(quán),許友生,吳鋒民.我國生物滲流的研究現(xiàn)狀和展望[J].浙江師范大學(xué)學(xué)報:自然科學(xué)版,2007,30(4):372-376.SONG Fu-quan,XU You-sheng,WU Feng-min.Survey and prospect of researches in biological porous flow in China[J].Journal of Zhejiang Normal University(Natural Sciences),2007,30(4):372-376.

[4] HUANG Z,YAO J,L I Y,et al.Per meability analysis of fractured vuggy porous media based on homogenization theory[J].Science China Technological Sciences,2010,53(3):839-847.

[5] SAFFMAN P G.On the boundary condition at the surface of a porous medium[J].Studies in Applied Mathematics,1971,L(2):93-101.

[6] BEAVERS G S,JOSEPH D D.Boundary conditions at a naturally permeable wall[J].Journal of Fluid Mechanics,1967,30:197-207.

[7] SAL INGER A G,AR IS R,DERBY J J.Finite element formulations for large-scale,coupled flows in adjacent porous and open fluid domains[J].Numerical Methods in Fluids,1994,18(2):1185-1209.

[8] GARTL INGD K,H ICKOX C E,G IVLER R C.Simulation of coupled viscous and porous flow problems[J].Computational Fluid Dynamics,1996,7(1):23-48.

[9] ARBOGAST T,BRUNSON D S.A computational method for approximating a Darcy-Stokes system governing a vuggy porous medium[J].Computational Geosciences,2007,11(3):207-218.

[10] KANSCHAT G,R IV IEREB.A strongly conservative finite element method for the coupling of Stokes and Darcy flow[J].Journal of Computational Physics,2010,229(17):5933-5943.

[11] 李亞軍,姚軍,黃朝琴,等.裂縫性油藏等效滲透率張量計算及表征單元體積研究[J].水動力學(xué)研究與進展:A輯,2010,25(1):1-7.L I Ya-jun,YAO Jun,HUANG Zhao-qin.Calculation of equivalent permeability tensor and study on representative element volume for modeling fractured reservoirs[J].Chinese Journal of Hydro-Dynamics(series A),2010,25(1):1-5.

[12] DURLOFSKY L J.Numerical calculation of equivalent grid block per meability tensors for heterogeneous porous media[J].Water Resources Research,1991,27(5):699-708.

[13] 姚軍,李亞軍,黃朝琴,等.裂縫性油藏等效滲透率張量的邊界元求解方法[J].油氣地質(zhì)與采收率,2009,16(6):80-83.YAO Jun,L I Ya-jun,HUANG Zhao-qin,et al.Calculation of equivalent permeability tensors of fractured reservoirs using boundary element method[J].Petroleum Geology and Recovery Efficiency,2009,16(6):80-83.

[14] FENGM,Q IR,ZHU R,et al.Stabilized Crouzeix-Raviart element for the coupled Stokes and Darcy problem[J].Applied Mathematics and Mechanics,2010,31(3):393-404.

(編輯 劉為清)

Estimating equivalent permeability of fractured-vuggy media based on coupled Darcy-Stokes model

L I Ya-jun,YAO Jun,HUANG Zhao-qin,LIU Yong-hui

(College of Petroleum Engineering in China University of Petroleum,Qingdao266555,China)

A coupled Darcy-Stokes model was developed to estimate equivalent permeability(EP)of fractured-vuggy media(FVM)according to the structure property of FVM and the fluid flow characteristics in different-scale voids.The generalized cubic law considering bedrock permeability and equivalent per meability tensor of single-fractured porous media were derived.The pretreatment technique of fractured porous media(FPM)was developed.Thus the vug-fracture-matrix coupled flow problem was converted into the vug-matrix problem.The equivalent permeability and permeability characteristics of different FVMs were analyzed.The results show that the EP of FVM was affected by the existence of fractures and vugs in various degrees.The EP of porous media with single circle vug is linear with the matrix permeability.The pretreatment technique of FPM can enhance the computation efficiency.

fractured-vuggy media;coupled Darcy-Stokes model;continuous media;equivalent permeability;flow simulation

TE 319;TE 344

A

10.3969/j.issn.1673-5005.2011.02.016

2010-11-08

國家科技重大專項(2008ZX05014-005-03);高等學(xué)校博士點基金項目(20090133110006);中國石油大學(xué)(華東)研究生創(chuàng)新基金項目(Z10-03)

李亞軍(1984-),男(漢族),山東費縣人,博士研究生,研究方向為油氣田開發(fā)理論與系統(tǒng)工程。

1673-5005(2011)02-0091-05