黃朝琴, 周 旭, 劉禮軍, 黃 濤, 姚 軍, 王曉光, HERVé Jourde

(1.中國石油大學(xué)(華東)油氣滲流研究中心,山東青島 266580; 2.法國國家科學(xué)研究中心水文科學(xué)所,蒙彼利埃 34095)

中國西部碳酸鹽巖油藏油氣資源豐富,其中縫洞型約占2/3,是增儲上產(chǎn)的現(xiàn)實領(lǐng)域[1]。與孔隙型和裂縫型碳酸鹽巖油藏不同,縫洞型油藏從儲集空間到流動規(guī)律都有較大差異[2-5],傳統(tǒng)的Darcy滲流理論已不完全適用[6]。對此,研究者提出了離散縫洞模型[7-8],采用Darcy-N-S(Navier-Stokes)方程來表征滲流-自由流耦合流動。縫洞型碳酸鹽巖油藏埋藏深度超過5 300 m,為深層油氣資源[9]。裂縫和溶蝕孔洞在開發(fā)過程中易發(fā)生變形,具有強應(yīng)力敏感性[10-13]。Murad等[14]在Beavers-Joseph-Saffman(BJS)條件基礎(chǔ)上,考慮巖石變形影響,將BJS條件擴展至流固耦合研究中。Badia等[15]應(yīng)用擴展的BJS條件,基于有限元對Biot-N-S(Navier-Stokes)方程進行了數(shù)值研究,Ambartsumyan等[16]則通過引入Biot-N-S方程拉格朗日乘子來構(gòu)造穩(wěn)定有限元格式。Zhang等[17]針對縫洞介質(zhì)提出了一種流固耦合計算格式,應(yīng)用FLAC3D軟件[18]開展了數(shù)值模擬研究。但其研究忽略了滲流-自由流界面上的切向流動;同時裂縫未采用幾何降維處理,導(dǎo)致計算量大難以實際應(yīng)用。對此,筆者建立離散縫洞流固耦合數(shù)學(xué)模型。多孔介質(zhì)滲流區(qū)域采用Biot方程,并對裂縫進行降維處理,采用Galerkin有限元求解;溶洞中采用N-S方程,采用Taylor-Hood混合元方法求解;兩區(qū)域間通過擴展BJS條件耦合。通過兩個經(jīng)典數(shù)值算例驗證了模型和方法的正確性。在此基礎(chǔ)上,通過實際縫洞型油藏算例對開采過程中的裂縫和溶洞變形及其對產(chǎn)量的影響進行研究。

1 離散縫洞流固耦合數(shù)學(xué)模型建立

縫洞型碳酸鹽油藏一般經(jīng)歷了多期構(gòu)造運動和巖溶作用[1],導(dǎo)致介質(zhì)中發(fā)育著大量裂縫和溶蝕孔洞,如圖1(a)所示。因此縫洞型儲層介質(zhì)為一復(fù)雜的離散縫洞模型,如圖1(b)所示。在該模型中,基質(zhì)巖塊和裂縫組成滲流系統(tǒng),為經(jīng)典的離散裂縫模型,其中的流動滿足Darcy定律;溶洞系統(tǒng)為自由流區(qū)域,流體流動采用N-S方程;兩區(qū)域間采用BJS條件[19-20]耦合。

圖1 縫洞型油藏巖心及離散縫洞模型示意圖Fig.1 Typical core of fractured vuggy carbonate reservoirs and discrete fracture-vug model

1.1 滲流區(qū)域流固耦合模型

為進一步考慮滲流區(qū)域介質(zhì)的變形影響,假設(shè)多孔介質(zhì)滲流區(qū)域滿足彈性小變形假設(shè),并考慮微可壓縮流體的等溫單相流動過程,其流固耦合過程由以下方程予以描述[12]。

(1)平衡方程為

(1)

式中,σp為應(yīng)力,以拉伸為正,下標(biāo)p表示多孔介質(zhì)區(qū)域;ρp、ρs和ρf分別為多孔介質(zhì)、固體骨架和流體的密度;φ為多孔介質(zhì)孔隙度;g為重力加速度。

(2)本構(gòu)方程為

σ′=σp+αppI=C:εs.

(2)

式中,σ′為巖石有效應(yīng)力;α為Biot有效應(yīng)力系數(shù);pp為孔隙流體壓力;I為單位張量;C為固體骨架的彈性模量;εs為固體骨架的應(yīng)變張量,滿足以下幾何協(xié)調(diào)方程。

(3)幾何方程為

(3)

式中,us為固體骨架位移,εv為相應(yīng)的體積應(yīng)變。

(4)Darcy方程為

(4)

式中,v為多孔介質(zhì)滲流速度;k為多孔介質(zhì)滲透率;μ為流體黏度。

(5)連續(xù)性方程為

(5)

其中

式中,q為源匯項;Kf為孔隙流體的體積彈性系數(shù)(流體壓縮系數(shù)的倒數(shù));Ks為固體骨架的體積彈性系數(shù);M為Biot模量,描述多孔介質(zhì)體積不變的情況下,在流體壓力作用下進入多孔介質(zhì)中的流體[21-22]。

在彈性小變形條件下,彈性模量C和Biot有效應(yīng)力系數(shù)α均可視為常數(shù)?？紫抖圈蘸蜐B透率k一般是有效應(yīng)力的函數(shù)[23]。研究中所有的微裂縫(包括節(jié)理)和微溶孔均包含于基質(zhì)巖塊系統(tǒng)中,可視為等效的連續(xù)介質(zhì)模型。對于宏觀裂縫,如圖2所示。

圖2 實際粗糙裂縫簡化示意圖Fig.2 Schematic of a real rough fracture and the simplified conceptual fracture

由于裂縫壁面粗糙且相互接觸,與此同時裂縫通常被其他物質(zhì)充填[24];因此本文中將宏觀裂縫視為狹長的高滲透區(qū)域,仍屬于多孔介質(zhì)滲流區(qū)域。然而,相對于整體縫洞介質(zhì)研究區(qū)域的特征尺度,裂縫可簡化為一個彈性薄層。為避免裂縫區(qū)域的大量精細(xì)網(wǎng)格剖分,本文中在幾何上對裂縫采用降維處理[25],建立相應(yīng)的離散裂縫流固耦合模型。

1.2 溶洞自由流區(qū)域流動方程

考慮溶洞中微可壓縮流體的低雷諾數(shù)等溫層流流動問題,滿足下述方程。

(1)Navier-Stokes方程為

(6)

式中,uf和σf分別為溶洞中的流體速度和應(yīng)力。

(2)本構(gòu)方程為

σf=-pfI+2μεf.

(7)

式中,εf為流體應(yīng)變張量,滿足以下幾何協(xié)調(diào)方程。

(3)幾何方程為

(8)

(4)連續(xù)性方程為

(9)

1.3 滲流-自由流耦合界面條件

如圖1(b)所示,在滲流-自由流耦合交界面Γpc上,需要引入相適應(yīng)的界面條件。對于單一流場問題,一般采用BJS條件來耦合Darcy方程和N-S方程。對于流固耦合問題,需要在界面處進一步考慮多孔介質(zhì)骨架變形的影響,根據(jù)質(zhì)量和動量守恒定律可采用拓展的BJS界面條件[14,26-27],

(10)

式中,n和τ分別為滲流和自由流耦合交界面Γpc的單位法向量(指向滲流區(qū)域)和單位切向量;β為耦合交界面上的BJS速度滑移系數(shù)。

2 有限元數(shù)值求解

采用有限元方法求解上述離散縫洞流固耦合數(shù)學(xué)模型。其中滲流區(qū)域采用標(biāo)準(zhǔn)的Galerkin有限元方法,自由流區(qū)域采用Taylor-Hood混合元方法,裂縫采用降維處理。

2.1 滲流區(qū)域數(shù)值離散格式

對于滲流區(qū)域的Biot多孔彈性力學(xué)方程,數(shù)值計算格式以節(jié)點位移和壓力為未知量,即us-pp格式[28]。采用標(biāo)準(zhǔn)的Galerkin有限元方法進行離散,則多孔介質(zhì)區(qū)域流固耦合方程(1)和(5)離散后的有限元計算格式為

(11)

其中

式中,{us}為節(jié)點固體骨架位移列陣;{pp}為節(jié)點孔隙壓力列陣;{fp}為多孔介質(zhì)區(qū)域體積力列陣;{qp}為多孔介質(zhì)區(qū)域流體體積力和源匯項列陣;變量符號上面的點表示對時間的一階導(dǎo)數(shù);Kp為多孔介質(zhì)固體骨架的標(biāo)準(zhǔn)剛度矩陣;B為應(yīng)變-位移矩陣;D為排水條件下的巖石固體骨架材料常數(shù)矩陣;L為多孔彈性介質(zhì)的流固耦合矩陣;對于三維問題α=α[1,1,1,0,0,0]T,二維問題α=α[1,1,0,0]T,稱為Biot有效應(yīng)力系數(shù)列陣;H為滲透矩陣;N為形函數(shù)向量;Ms為多孔介質(zhì)流體儲集矩陣。

上述方程系統(tǒng)可采用標(biāo)準(zhǔn)的向后差分時間格式進行離散和求解,詳細(xì)的推導(dǎo)和求解過程可參考文獻[28]。

2.2 自由流區(qū)域數(shù)值離散格式

對于自由流區(qū)域,需要聯(lián)立求解方程(6)和(9),由方程(6)可知速度插值函數(shù)最好比壓力插值函數(shù)高一階,以達到較好的計算精度[29]。對此應(yīng)用加權(quán)殘量法推導(dǎo)有限元方程,采用經(jīng)典的Taylor-Hood混合單元(Pn-Pn-1)構(gòu)造插值函數(shù)[28-29],具體表達式為

(12)

式中,I和J分別為速度與壓力單元節(jié)點數(shù);Nu和Np分別為速度和壓力的形函數(shù)。

對于微可壓縮流體的等溫流動過程,可忽略流體密度的變化,因此在自由流區(qū)域可采用不可壓縮流體流動的Navier-Stokes方程,離散后的數(shù)值計算格式為

(13)

其中

式中,{uf}為節(jié)點的流體速度列陣;{pf}為節(jié)點的流體壓力列陣;{ff}為自由流區(qū)域載荷列陣,包括體積力和外邊界給定的應(yīng)力邊界條件;T為自由流區(qū)域給定的應(yīng)力邊界條件。

上述方程系統(tǒng)(13)式可采用標(biāo)準(zhǔn)的向后差分時間格式進行離散和求解,詳細(xì)的推導(dǎo)和求解過程可參考文獻[29]和[30]。

2.3 離散裂縫數(shù)值離散格式

如前所述,實際裂縫可視為一狹長的高滲透多孔介質(zhì)區(qū)域。在數(shù)值處理中,對于滲流場可通過等效流量原理對裂縫進行降維處理,目前已有大量相關(guān)研究[3-5]。對于應(yīng)力場,裂縫單元開度d和單元長度l相比,往往相差多個數(shù)量級。此時,類似于離散裂縫滲流模型中的等效流量原理,可將應(yīng)力場中的裂縫也進行降維處理,如圖3所示。此時,相應(yīng)的變形剛度矩陣為

(14)

其中下標(biāo)th表示薄層裂縫。當(dāng)d/l的比值取值在0.01～0.1時可取得較高的計算精度[32-33]。對于方程(12)和(13)中的其他項可進行類似的處理,對于三維問題,裂縫降維后為一裂縫面,將出現(xiàn)兩個方向的單元長度l。此時,裂縫單元開度d應(yīng)和單元長度l較小的方向取比值進行分析。

(15)

圖3 裂縫薄層單元簡化示意圖Fig.3 Schematic of a simplified fracture thin-layer model

多孔介質(zhì)流固耦合方程(11)和自由流區(qū)域方程(13)通過引入交界面邊界條件方程(10)即可形成縫洞型介質(zhì)的流固耦合數(shù)值計算格式。對兩方程采用相同的時間離散格式,本文中采用向后差分格式,可得到最終的有限元計算格式。

3 數(shù)值算例分析

3.1 數(shù)值驗證算例

3.1.1 流固耦合驗證算例

考慮如圖4所示的一維固結(jié)問題,該物理過程可分解為兩個階段:第一階段,在多孔介質(zhì)上表面瞬時加載均勻載荷,為非排水過程,導(dǎo)致孔隙壓力升高;第二階段為排水階段,壓力逐漸降低,垂向位移持續(xù)增加。本算例模型為1 m×1 m×100 m(x×y×z),采用結(jié)構(gòu)化網(wǎng)格剖分,網(wǎng)格尺寸為0.5 m×0.5 m×0.5 m(x×y×z)。其他計算參數(shù)取自于文獻[12]7.5節(jié)中的Berea砂巖試驗數(shù)據(jù),上覆施加載荷為0.1 MPa。由于排水和非排水過程存在一定的差異,兩者之間的主要差別在于是否考慮流體流出,前者沒有流體流出,因此測量得到的泊松比并不一樣,具體參數(shù)取值為:剪切模量G=6 GPa,泊松比ν=0.2(排水),泊松比ν=0.33(不排水),孔隙度φ=0.19,滲透率k=1.875×10-13m2,彈性模量E=14.4 GPa,拉梅常數(shù)λ=4 GPa,Biot模量M=12.25 GPa,Biot有效應(yīng)力系數(shù)α=0.79,流體可壓縮系數(shù)Cf=3.5×10-10Pa-1,流體黏度μ=1.0 mPa·s。

圖5為第二階段排水過程中本文數(shù)值解與解析解的對比。顯然,兩者具有很好的一致性,驗證了本文中多孔介質(zhì)區(qū)域流固耦合數(shù)值計算的正確性。

圖4 一維固結(jié)物理過程示意圖Fig.4 Schematic of 1D consolidation process

圖5 一維固結(jié)問題數(shù)值解與解析解的比較Fig.5 Comparison of pore pressure and vertical displacement between numerical results and analytical solutions

3.1.2 滲流自由流耦合驗證算例

該驗證模型采用經(jīng)典的Beavers-Joseph試驗?zāi)Ｐ蚚18,32],如圖6(a)所示。考慮自由流-滲流穩(wěn)定耦合層流問題,pin=0.5 Pa,pout=0 Pa,模型長度L=0.5 m,高度h=0.1 m,兩區(qū)域的壓力梯度為1.0 Pa/m,滲流區(qū)域的孔隙度為0.2,滲透率k=1×10-12m2。

為簡化計算,流體的密度取值為1 kg/m3,黏度為1 Pa·s,在交界面上BJS速度滑移系數(shù)β取值為1.0。有限元計算采用如圖6(b)所示的結(jié)構(gòu)化網(wǎng)格剖分,并在自由流-滲流交界面處網(wǎng)格加密。圖7為本文數(shù)值解和解析解的對比,兩者基本吻合,驗證了滲流自由流耦合數(shù)值計算的正確性。詳細(xì)的解析解形式可參考文獻[20]和[21]。

圖6 滲流-自由流耦合流動Fig.6 Schematic of coupling porous-free flow system and its structured meshing

圖7 本文數(shù)值解與解析解對比Fig.7 Comparison between analytical solution and numerical solution

3.2 縫洞介質(zhì)流固耦合算例

考慮如圖8所示的碳酸鹽巖縫洞型介質(zhì)模型,該模型整體處于奧陶系,頂深為5500 m,右端溶洞與兩條正交裂縫相連通,模型的右下端存在底水能量補充,生產(chǎn)井鉆遇左端溶洞并進行定壓生產(chǎn)(45 MPa)。模型的初始流體壓力、上覆壓力和右側(cè)水平側(cè)壓均由圖9地應(yīng)力深度分布曲線計算得到,其中右側(cè)水平側(cè)壓取為最大水平主應(yīng)力。圖9中的數(shù)據(jù)取自現(xiàn)場測試數(shù)據(jù)和文獻中數(shù)據(jù)[35-38]。多孔介質(zhì)區(qū)域的左端和底部均為滑動邊界。

實際生產(chǎn)過程為油水兩相流動,本文中模型僅考慮了單相流動,具體計算參數(shù)為:流體密度ρf=870 kg/m3,固體骨架密度ρs=2 650 kg/m3,流體黏度μf=20 mPa·s,固體骨架彈性模量Em=48 GPa,裂縫彈性模量Ef=24 GPa,泊松比ν=0.27,基質(zhì)巖塊孔隙度φm=0.09,裂縫孔隙度φf=0.5,基質(zhì)巖塊滲透率km=1×10-14m2,裂縫滲透率kf=1×10-10m2,BJS速度滑移系數(shù)β=1.0,Biot有效應(yīng)力系數(shù)α=0.79,流體可壓縮系數(shù)Cf=5×10-10Pa-1,剪切強度τn=80 MPa。其中基質(zhì)巖塊包含了微小裂縫和溶蝕孔洞,其滲透率取值大于實際純基質(zhì)的滲透率。圖10為流固耦合效應(yīng)對產(chǎn)量的影響,計算中采用三角形單元進行網(wǎng)格剖分(圖11(a)),包含19 072個單元,其中多孔介質(zhì)區(qū)域采用二次插值函數(shù)單元,溶洞區(qū)域采用P2-P1型Taylor-Hood單元。

圖8 縫洞型介質(zhì)幾何模型及約束加載圖Fig.8 Geometry and constraints of the fractured vuggy medium

圖9 地應(yīng)力與孔隙壓力隨深度變化曲線Fig.9 Variation of stress and pore pressure with depth

圖10 不同模型的產(chǎn)量對比Fig.10 Comparison of productions between different models

如圖10所示,由于井打在溶洞上,生產(chǎn)井日產(chǎn)量初期很高,但由于能量補充不足迅速下降;對于流固耦合模型,由于充分考慮了固體骨架的變形作用,使油藏的供液能力在初期有所提高,因此其初期日產(chǎn)量和累積產(chǎn)量相比于不考慮流固耦合效應(yīng)的純流動模型均有所提高。但隨著開采的進行,由于孔隙壓力的降低,固體骨架進一步變形導(dǎo)致多孔介質(zhì)區(qū)域的孔隙度及滲透率降低,使井的日產(chǎn)量和累積產(chǎn)量均低于不考慮流固耦合效應(yīng)的純流動模型。計算中,基質(zhì)巖塊的孔隙度和滲透率變化與平均應(yīng)力相關(guān)[39-40],表達式為

(16)

對于裂縫,孔隙度仍采用式(16),但滲透率與裂縫開度變化相關(guān)[40],表達式變?yōu)?/p>

(17)

式中,kf0為裂縫初始滲透率;c1為經(jīng)驗參數(shù)(本算例取值0.071 2);Δd為裂縫開度變化量。

圖11(b)中給出了耦合模型計算達到穩(wěn)定后的壓力云圖和速度矢量分布,計算結(jié)果顯示流體在溶洞中的流動速度較快,可近似處理為等勢體;多孔介質(zhì)區(qū)域中裂縫作為高導(dǎo)流通道對壓力和速度也有顯著影響。圖11(c)中給出了圍壓狀態(tài)下最大主應(yīng)力和最小主應(yīng)力的差值云圖,可看到在該荷載工況下,裂縫尖端具有明顯的應(yīng)力集中,但溶洞的應(yīng)力集中部位并非都出現(xiàn)在其頂部和底部,采用Tresca強度準(zhǔn)則進一步分析其剪切破壞區(qū)域,

(18)

式中,σ1和σ3分別為最大和最小主應(yīng)力;τ為巖石所受的剪切應(yīng)力。

破壞程度F計算公式為

(19)

式中,F為破壞程度,F≥1,表示巖石已發(fā)生剪切破壞或臨界破壞區(qū)域,F<1,表示未發(fā)生破壞;τn為巖石的剪切強度。

圖11 縫洞型介質(zhì)流固耦合算例結(jié)果Fig.11 Results of hydro-mechanical process for fractured vuggy medium example

根據(jù)上述判定準(zhǔn)則可得到巖石剪切破壞區(qū)域分布,如圖11(d)所示。破壞區(qū)域主要集中于裂縫尖端和溶洞尖點區(qū)域。計算中宏觀裂縫的初始開度均為0.05 m,壓力場穩(wěn)定后,垂直裂縫開度平均增加0.026 m;水平裂縫開度平均減少0.036 m,閉合程度較高。該結(jié)果表明,水平裂縫與左端溶洞破壞后相連通,這對于溝通縫洞結(jié)構(gòu)體有利,但裂縫的閉合對其連通程度也有降低作用,同時破壞區(qū)域的儲滲參數(shù)變化也有重要影響,因此破壞區(qū)域?qū)α鲃拥挠绊懶杈C合考慮和研究。

4 結(jié) 論

(1)針對深層縫洞型油藏的強應(yīng)力敏感性,創(chuàng)建了離散縫洞模型的流固耦合數(shù)學(xué)模型,其中多孔介質(zhì)滲流區(qū)域采用Biot方程,溶洞區(qū)域視為自由流區(qū)域采用N-S方程,兩區(qū)域間采用擴展的Beavers-Joseph-Saffman邊界條件進行耦合。采用有限元法對Biot-N-S(Navier-Stokes)流固耦合數(shù)學(xué)模型進行了數(shù)值離散:滲流區(qū)域采用經(jīng)典Galerkin有限元方法;自由流區(qū)域采用Taylor-Hood混合元方法。裂縫采用幾何降維處理以提高計算效率:在流場中,裂縫為高滲流通道;在應(yīng)力場中,裂縫處理為薄層單元。通過一維固結(jié)和Beavers-Joseph實驗?zāi)Ｐ蛢蓚€數(shù)值算例驗證了模型和數(shù)值方法的正確性。

(2)溶洞中的壓力傳播速度較快,相對于滲流區(qū)域可視為一流動等勢體。在降壓生產(chǎn)過程中,裂縫尖端和溶洞附近區(qū)域易發(fā)生較大面積的破壞,而較高的流體壓力對于溶洞和裂縫壁面具有一定的支撐作用;因此在此類油氣藏的開發(fā)中應(yīng)適當(dāng)采取保壓措施,以避免溶洞坍塌和裂縫閉合。