劉 靜, 鄭黎明, 李悅靜, 夏軍勇, 蒲春生

(1.中國石油大學(華東) 石油工程學院， 山東 青島 266580； 2.非常規(guī)油氣開發(fā)教育部重點實驗室(中國石油大學(華東))， 山東 青島 266580; 3.燕山大學 車輛與能源學院， 河北 秦皇島 066004)

0 引言

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

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

1 數(shù)學模型

低頻波動采油對實際地層滲流的影響較為復雜，為研究方便，本文模擬中假設儲層巖石滿足線彈性變化關系，假設低頻采油激振力較小，不引起巖石本身彈塑性質的變化，仍采用線性本構關系.

1.1 運動方程和連續(xù)方程

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

(1)

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

(2)

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

1.2 狀態(tài)方程

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

(3)

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

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

(4)

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

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

(5)

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

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

(6)

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

波動下粘度變化方程[16]：

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

(7)

孔隙變化方程：

(8)

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

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

(9)

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

1.3 定解條件

1.3.1 初始壓力

在該一維儲層中，彈性波作用之前，油水處于滲流平衡狀態(tài).在單相流動區(qū)，可根據液柱重力梯度算出初始壓力分布狀況.各相的壓力遵循單相流動時的計算方法，同一點不同流相間的壓力差為毛管力.初始時刻壓力可表示為：

pf(x,0)=pe

(10)

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

1.3.2 初始飽和度

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

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

(11)

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

1.3.3 初始位移和速度

假設研究的初始時刻，固相位移為0，流體宏觀滲流位移也為0，但在彈性波作用之后的滲流過程中同時存在宏觀滲流和波誘導滲流，初始時刻速度不等于0.

流固初始滲流位移為：

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

(12)

初始滲流速度為：

(13a)

(13b)

1.4 邊界條件

彈性傳播方向與研究儲層滲流方向一致時，油藏邊界包括頂部、底部、注入井邊界和采出井邊界四部分，其中，低頻波動邊界可與注入或采出邊界條件重疊.

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

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

(14)

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

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

(15)

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

無限遠或強制位移邊界(數(shù)值模擬過程中引入吸收邊界條件)處：

u(∞,t)=0

(16)

2 模型求解

2.1 物理模型建立

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

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

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

2.2 波動耦合滲流與Darcy滲流對比分析

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

表1 基礎參數(shù)表

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

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

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

2.3 波動參數(shù)對滲流的影響

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

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

模擬結果如圖3所示，相同模擬時間內(500 s)，當振幅u0在100～350μm之間變化時，隨著振幅的增大，水驅前緣向左偏移，即振幅越小時同一位置處油相優(yōu)先被水相驅替.當u0<300μm時，水驅前緣在空間上隨振幅的變化不明顯，一方面可能是因為驅替時間較短，另一方面可能是振幅足夠強時飽和度變化敏感性變弱，即在一定儲層物性特征條件下可選擇合適的振幅強度作為工作振幅.

(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時不同振幅下含油 飽和度的空間分布

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

圖4 500 s時不同振幅下A-A′線含水 飽和度變化

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

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

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

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

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

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

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

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

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

3 結論

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

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