劉 靜, 鄭黎明, 李?lèi)傡o, 夏軍勇, 蒲春生

(1.中國(guó)石油大學(xué)(華東) 石油工程學(xué)院， 山東 青島 266580； 2.非常規(guī)油氣開(kāi)發(fā)教育部重點(diǎn)實(shí)驗(yàn)室(中國(guó)石油大學(xué)(華東))， 山東 青島 266580; 3.燕山大學(xué) 車(chē)輛與能源學(xué)院， 河北 秦皇島 066004)

0 引言

低頻波動(dòng)采油技術(shù)，是一種利用物理場(chǎng)激勵(lì)油層達(dá)到提高原油采收率目的的增產(chǎn)技術(shù)[1,2].揭示低頻彈性波作用下儲(chǔ)層滲流動(dòng)力學(xué)變化機(jī)制是該技術(shù)應(yīng)用和推廣的基礎(chǔ).然而，關(guān)于低頻波動(dòng)采油機(jī)理的系列研究主要集中于實(shí)驗(yàn)研究和半定量簡(jiǎn)化動(dòng)力學(xué)模型分析[3-6].實(shí)驗(yàn)?zāi)M主要定性揭示了儲(chǔ)層物性受振動(dòng)參數(shù)的影響敏感性，無(wú)法定量解釋波動(dòng)對(duì)儲(chǔ)層滲流過(guò)程的影響規(guī)律.理論模型研究方面，部分研究人員通過(guò)建立半定量簡(jiǎn)化波動(dòng)滲流模型(將波動(dòng)力直接疊加到滲流運(yùn)動(dòng)方程)或者簡(jiǎn)化為單一毛管，研究波動(dòng)對(duì)粘附層、液滴脫落等的影響[7-10].關(guān)于初始時(shí)刻滲流場(chǎng)與波場(chǎng)的耦合作用問(wèn)題，研究人員也開(kāi)展了部分理論分析工作，但主要針對(duì)一維、二維單相流體的作用過(guò)程，對(duì)于油水兩相滲流場(chǎng)與波場(chǎng)耦合模型研究目前較少[11,12].

基于此，借鑒經(jīng)典的孔隙介質(zhì)彈性波傳播理論模型，并引入目前低頻波動(dòng)采油技術(shù)已有實(shí)驗(yàn)機(jī)理認(rèn)識(shí)，聯(lián)立波動(dòng)條件下的流體運(yùn)動(dòng)方程、固體運(yùn)動(dòng)方程、連續(xù)性方程和狀態(tài)方程，利用COMSOL數(shù)值模擬軟件，建立物理場(chǎng)作用下油水兩相滲流的物理模型和數(shù)學(xué)模型，開(kāi)展波動(dòng)耦合滲流規(guī)律的研究，揭示低頻波動(dòng)對(duì)采油過(guò)程的影響及振動(dòng)參數(shù)影響的敏感性.

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

低頻波動(dòng)采油對(duì)實(shí)際地層滲流的影響較為復(fù)雜，為研究方便，本文模擬中假設(shè)儲(chǔ)層巖石滿(mǎn)足線(xiàn)彈性變化關(guān)系，假設(shè)低頻采油激振力較小，不引起巖石本身彈塑性質(zhì)的變化，仍采用線(xiàn)性本構(gòu)關(guān)系.

1.1 運(yùn)動(dòng)方程和連續(xù)方程

油水滲流過(guò)程中，不僅液相在多孔介質(zhì)中發(fā)生運(yùn)動(dòng)，固相巖石受到載荷作用也會(huì)發(fā)生形變，從而產(chǎn)生運(yùn)動(dòng).初始時(shí)刻流體流動(dòng)不考慮慣性作用和流固耦合作用，結(jié)合Darcy定律，基于均勻、各向同性、彈性多孔介質(zhì)的經(jīng)典Biot彈性波傳播模型，考慮流固耦合的固相、液相運(yùn)動(dòng)方程為[13-15]：

(1)

固相、液相的連續(xù)性方程為：

(2)

式(1)、(2)中：下標(biāo)w表示水相，o表示油相；τ為作用在單元體上的總應(yīng)力，MPa；u為固體位移，m；U為水相對(duì)于固體的相對(duì)位移，U=Ui+Uv，Ui為水相初始宏觀(guān)位移，Uv為水相波誘導(dǎo)位移，m；W為油相對(duì)于固體的相對(duì)位移，W=Wi+Wv，Wi為油相初始宏觀(guān)位移，Wv為油相波誘導(dǎo)位移，m；pw為水相壓力，MPa，po為油相壓力，MPa；ρw為水相密度，kg/m3；ρo為油相密度，kg/m3；ρ為多孔介質(zhì)平均密度，kg/m3，ρd為固體顆粒密度，kg/m3；ηw為水相粘度，mPa·s；ηo為油相粘度，mPa·s；kw為水相滲透率，μm2；ko為油相滲透率，μm2；sw為水相飽和度；so為油相飽和度；φ為孔隙度.

1.2 狀態(tài)方程

彈性波作用下流體密度隨壓力變化而變化，在滲流力學(xué)中常用流體壓縮系數(shù)和壓力變化量的函數(shù)表示，此文中用體積模量來(lái)表示液體壓縮方程.

(3)

式(3)中：Kw為水相體積模量，MPa；Ko為油相的體積模量，MPa.

考慮固體壓縮性，用孔隙壓力引起固體密度改變和有效應(yīng)力引起的固體密度變化之和來(lái)表示固體壓縮性[13].

(4)

式(4)中：φ為孔隙度，εv為固體骨架體應(yīng)變，εv=·u；α為Biot系數(shù)；1/Qc為Biot引入的假設(shè)下一截面的代表流固間耦合關(guān)系的系數(shù)，為正值，其中K為巖石體積模量，MPa；Kdry為干燥固體骨架體積模量，MPa.

其中K的取值可利用斑塊飽和模型計(jì)算得到，巖石骨架體積模量和流體模的關(guān)系可表示為：

(5)

式(5)中：Ki(i=w，o)為第i相流體的體積模量，MPa；Kmi為單獨(dú)飽和某一相流體時(shí)的巖石體積模量，MPa.

巖石骨架的變形以及流體壓力的變化會(huì)引起孔隙度發(fā)生改變，根據(jù)Biot理論傳播模型，孔隙度、形變和孔隙壓力的關(guān)系式可表示為[14,15]：

(6)

式(6)中：P為油水兩相滲流系統(tǒng)平均孔隙壓力，用體積平均法求得p=swpw+sopo，MPa.

波動(dòng)下粘度變化方程[16]：

ηo=ηo(ω)≈ηiF(ω)

(7)

孔隙變化方程：

(8)

式(8)中：k為絕對(duì)滲透率，μm2；k0為初始絕對(duì)滲透率，μm2；φ0為初始孔隙度.

相滲模型采用Mualem理論模型，毛管力曲線(xiàn)采用以下方程：

(9)

式(9)中：m、n、l為van Genuchten毛管力曲線(xiàn)隱式方程的擬合系數(shù).

1.3 定解條件

1.3.1 初始?jí)毫?/p>

在該一維儲(chǔ)層中，彈性波作用之前，油水處于滲流平衡狀態(tài).在單相流動(dòng)區(qū)，可根據(jù)液柱重力梯度算出初始?jí)毫Ψ植紶顩r.各相的壓力遵循單相流動(dòng)時(shí)的計(jì)算方法，同一點(diǎn)不同流相間的壓力差為毛管力.初始時(shí)刻壓力可表示為：

pf(x,0)=pe

(10)

式(10)中：pe為油藏壓力，MPa.

1.3.2 初始飽和度

對(duì)于多相滲流問(wèn)題，需要指定油藏的初始飽和度分布，在單相流動(dòng)區(qū)，飽和度為定值；在油水兩相區(qū)，各相飽和度通過(guò)聯(lián)立運(yùn)動(dòng)方程、連續(xù)性方程、狀態(tài)方程等的波動(dòng)耦合滲流模型求得.初始時(shí)刻飽和度可表示為：

sw(x,0)=swc，so(x,0)=1-swc

(11)

式(11)中：swc為束縛水飽和度.

1.3.3 初始位移和速度

假設(shè)研究的初始時(shí)刻，固相位移為0，流體宏觀(guān)滲流位移也為0，但在彈性波作用之后的滲流過(guò)程中同時(shí)存在宏觀(guān)滲流和波誘導(dǎo)滲流，初始時(shí)刻速度不等于0.

流固初始滲流位移為：

u(x,0)=0,U(x,0)=W(x,0)=0

(12)

初始滲流速度為：

(13a)

(13b)

1.4 邊界條件

彈性傳播方向與研究?jī)?chǔ)層滲流方向一致時(shí)，油藏邊界包括頂部、底部、注入井邊界和采出井邊界四部分，其中，低頻波動(dòng)邊界可與注入或采出邊界條件重疊.

油藏邊界處水相壓力等于注入端壓力，油相壓力等于采出端壓力：

pw(0,t)=pin,po(L,t)=pout

(14)

震源邊界采用最簡(jiǎn)單的位移邊界表示，本文中用諧函數(shù)形式表示如下：

u(0,t)=u0(x,t)=u0exp(iωt)

(15)

式(15)中：u0為初始振幅.

無(wú)限遠(yuǎn)或強(qiáng)制位移邊界(數(shù)值模擬過(guò)程中引入吸收邊界條件)處：

u(∞,t)=0

(16)

2 模型求解

2.1 物理模型建立

利用COMSOL進(jìn)行有限元模擬，基于控制方程分別構(gòu)建固體運(yùn)動(dòng)方程、流體運(yùn)動(dòng)方程、壓力-位移方程(水相和油相壓力方程)和飽和度-位移方程的系數(shù)型偏微分方程形式，一維滲流模擬時(shí)，忽略z方向的控制方程分量，并建立彈性波一維驅(qū)替物理模型，彈性波作用下一維驅(qū)替模型中，保證彈性波傳播的主方向與原滲流場(chǎng)方向一致，且沿x方向，如圖1所示.

在研究區(qū)域縱向外圍各拓展L×d/4，右側(cè)外圍拓展d×d/4，作為PML吸收層，以盡可能地降低彈性波在人工邊界處虛假反射，提高數(shù)值模擬的結(jié)果準(zhǔn)確性，且PML吸收層上、下邊界為不滲透邊界.滲流場(chǎng)設(shè)定為定壓開(kāi)采，左側(cè)為注入端，右側(cè)為采出端，初始時(shí)刻Pf=2.0 MPa/m.在計(jì)算區(qū)域左側(cè)加載u=u0exp(iωt)的波場(chǎng).

圖1 計(jì)算區(qū)域與邊界條件

2.2 波動(dòng)耦合滲流與Darcy滲流對(duì)比分析

將波動(dòng)耦合滲流模型模擬結(jié)果與Darcy滲流模型模擬結(jié)果對(duì)比，分析波動(dòng)對(duì)流體滲流過(guò)程的影響，基礎(chǔ)參數(shù)如表1所示.圖2為不加載彈性波的Darcy滲流和加載彈性波耦合滲流兩種情況下含水飽和度和含油飽和度的分布圖.

表1 基礎(chǔ)參數(shù)表

模擬發(fā)現(xiàn)當(dāng)同樣開(kāi)采500 s時(shí)，相較與Darcy滲流水驅(qū)前緣(飽和度在x=4 m處發(fā)生明顯變化)，波動(dòng)耦合滲流中水驅(qū)前緣(含油飽和度在x=8 m處發(fā)生明顯變化)運(yùn)動(dòng)更快，離注入端更遠(yuǎn).說(shuō)明在一維定壓開(kāi)采過(guò)程中，彈性波作用于儲(chǔ)層，給儲(chǔ)層一個(gè)反復(fù)的振動(dòng)作用力，改變了儲(chǔ)層物性，促進(jìn)了油水滲流，低頻波動(dòng)強(qiáng)化采油技術(shù)有較好的增滲增油效果，從理論上說(shuō)明波動(dòng)可以提高水驅(qū)波及范圍和水驅(qū)效率的機(jī)理.

(a)Darcy滲流含油飽和度變化

(b)波動(dòng)耦合滲流含油飽和度變化圖2 作用500 s時(shí)波動(dòng)耦合滲流和 Darcy滲流油水分布變化圖

2.3 波動(dòng)參數(shù)對(duì)滲流的影響

前期研究表明振動(dòng)頻率為固有頻率時(shí)激勵(lì)效果最佳，本模型模擬基于參數(shù)表，令ω=15 Hz，η0=100 mPa·s，k0=50×10-3μm2，φ0=20%，模擬不同振幅u0=100～350μm對(duì)含油分布、油水滲流速度和孔隙壓力的影響，為振動(dòng)采油參數(shù)優(yōu)化提供依據(jù).

2.3.1 振幅對(duì)含油飽和度影響規(guī)律分析

模擬結(jié)果如圖3所示，相同模擬時(shí)間內(nèi)(500 s)，當(dāng)振幅u0在100～350μm之間變化時(shí)，隨著振幅的增大，水驅(qū)前緣向左偏移，即振幅越小時(shí)同一位置處油相優(yōu)先被水相驅(qū)替.當(dāng)u0<300μm時(shí)，水驅(qū)前緣在空間上隨振幅的變化不明顯，一方面可能是因?yàn)轵?qū)替時(shí)間較短，另一方面可能是振幅足夠強(qiáng)時(shí)飽和度變化敏感性變?nèi)酰丛谝欢▋?chǔ)層物性特征條件下可選擇合適的振幅強(qiáng)度作為工作振幅.

(a)u0=100 μm

(b)u0=150 μm

(c)u0=200 μm

(d)u0=250 μm

(e)u0=300 μm

(f)u0=350 μm圖3 500 s時(shí)不同振幅下含油 飽和度的空間分布

對(duì)不同振幅條件下軸線(xiàn)A-A′含水飽和度變化進(jìn)行繪圖，如圖4所示.由于計(jì)算時(shí)間較小，含水飽和度變化差距較??；彈性波場(chǎng)激勵(lì)下飽和度沿軸線(xiàn)A-A′發(fā)生波動(dòng)性變化，呈先急劇減小、后波動(dòng)性增加與減小的趨勢(shì)；根據(jù)急劇減小段(x<10 m)的變化速率可知，振幅愈大，彈性波場(chǎng)激勵(lì)作用愈強(qiáng)，對(duì)含水飽和度的影響愈大，這一定程度上亦反映了實(shí)驗(yàn)揭示的波動(dòng)降低束縛水飽和度的現(xiàn)象.

圖4 500 s時(shí)不同振幅下A-A′線(xiàn)含水 飽和度變化

2.3.2 振幅對(duì)流速影響規(guī)律分析

模擬結(jié)果如圖5所示，在圖5(a)中水相流速隨振幅增大波動(dòng)變化，但變化范圍不大(0～0.5×10-5m/s之間)，當(dāng)振幅u0=150μm時(shí)，水相流速最小.在圖5(b)中油相速度隨振幅的增加先增大后減小，當(dāng)振幅u0=250μm時(shí)油相速度與水相速度最接近.油水相速度越接近，流速越小，水驅(qū)突進(jìn)程度越低，波及面積越大，驅(qū)油效果越明顯.油藏內(nèi)部相對(duì)注入端，振動(dòng)對(duì)油藏內(nèi)部油水相流速的促進(jìn)明顯降低，說(shuō)明隨傳播距離的增加，波動(dòng)引起的滲流速度變化逐漸降低，反映了波衰減情況下對(duì)滲流的影響減弱.綜合水相和油相的速度變化規(guī)律，最佳振幅為150～250μm之間.

(a)水相速度隨振幅的變化規(guī)律

(b)油相速度隨振幅的變化規(guī)律圖5 500 s時(shí)振幅對(duì)A-A′線(xiàn)處液體速度 的影響規(guī)律

隨著距離的增加，從注入水端面到油水交界處(變化的曲面)到純油水區(qū)，油水飽和度逐漸發(fā)生變化，毛管力、油水相對(duì)滲透率亦發(fā)生變化，水相流速出現(xiàn)先增加后減小的情況；論文分析過(guò)程中考慮了波場(chǎng)對(duì)滲流場(chǎng)的激勵(lì)作用，尤其是低滲孔隙介質(zhì)中激勵(lì)作用愈強(qiáng)，流固耦合交錯(cuò)下流速變化幅度愈大，因此波動(dòng)耦合滲流場(chǎng)流速變化比達(dá)西流速變化更為劇烈.另外由于彈性波傳播過(guò)程是一個(gè)周期性的變化過(guò)程，隨距離和時(shí)間均會(huì)出現(xiàn)相位變化(振幅出現(xiàn)正負(fù)之分)，波誘導(dǎo)流速既可實(shí)現(xiàn)正向激勵(lì)，又可表現(xiàn)為反向抑制，當(dāng)然在考慮毛管力作用下實(shí)際孔隙中油水流度以正向激勵(lì)為主，因此在沿距離上流速出現(xiàn)了波動(dòng)性的變化趨勢(shì)，表現(xiàn)為流速速度振動(dòng)劇烈；且可能在一定位置或時(shí)間時(shí)反向波誘導(dǎo)流速大于原滲流場(chǎng)流速情形，使得圖5(b)中流速出現(xiàn)負(fù)值，當(dāng)孔隙介質(zhì)滲透率較高、波誘導(dǎo)流速較小，原滲流場(chǎng)流速較高時(shí)，油水兩相流速即可全部表現(xiàn)為正值.

2.3.3 振幅對(duì)壓力影響規(guī)律分析

圖6為彈性波作用下，500 s時(shí)A-A′上的壓力隨振幅的變化規(guī)律.在距振動(dòng)端較近的位置(x<12 m)，壓力隨振幅的增大先減小后增大，在u0=250μm時(shí)達(dá)到最大，隨后減小，反映了在距波源較近的位置處u0=250μm為最佳振幅.在油藏中部(12 m

綜合以上各物性特征的變化規(guī)律，建議低頻波動(dòng)強(qiáng)化采油最佳振幅為150～300μm，礦場(chǎng)應(yīng)用時(shí)根據(jù)設(shè)備和具體井筒、儲(chǔ)層情況選擇.

圖6 500 s時(shí)A-A′線(xiàn)處孔隙 壓力隨振幅的變化規(guī)律

3 結(jié)論

(1)基于Biot理論和等效流體理論，考慮油水兩相初始Darcy滲流、低頻波動(dòng)下孔隙-應(yīng)變關(guān)系、狀態(tài)方程、連續(xù)性方程和輔助方程，引入實(shí)際開(kāi)發(fā)油藏邊界條件和初始條件，建立低頻波作用下一維油水兩相流固耦合滲流模型.

(2)通過(guò)對(duì)模型求解和敏感性參數(shù)分析，認(rèn)為與常規(guī)油水兩相Darcy滲流規(guī)律相比，波動(dòng)耦合滲流場(chǎng)中儲(chǔ)層受到一個(gè)反復(fù)的振動(dòng)作用力，改變了儲(chǔ)層物性，促進(jìn)了油水滲流，說(shuō)明低頻波動(dòng)強(qiáng)化采油技術(shù)有較好的增滲增油效果.在地層固有頻率15 Hz下，基于振幅對(duì)含油飽和度、水相流速和油相流速、壓力等的影響規(guī)律，優(yōu)選出低頻波動(dòng)強(qiáng)化采油最佳振幅為150～300μm.