用于多層油藏層間干擾預測的滲流阻力計算方法

馬子麟，師永民，李宜強，管 錯，師 翔，李文宏，張忠義，左 靜

（1.北京大學地球與空間科學學院，北京 100871；2.中國石油大學（北京）石油工程學院，北京 102249；3.中海油研究總院有限責任公司，北京 100028；4.中國石油長慶油田分公司勘探開發研究院，陜西西安 710018）

多層砂巖油藏在中國廣泛分布，儲量和產量均約占全國總量的50%［1］。由于層間非均質性的影響，多層油藏注水開發過程中普遍存在層間干擾問題，導致低滲透層不能充分動用，十分不利于采收率的提高［2-5］。因此，層間干擾的表征和分析一直是研究的熱點和難點。業界主要采用物理模擬與數值模擬研究層間干擾特征。早期研究常采用采收率表征層間干擾對地層生產狀況的最終影響［6-10］；近年來，許多學者在研究中發現，隨著開發過程的進行，層間干擾程度是動態改變的［11-14］，但目前仍需通過模擬方法或分層測試表征、預測層間干擾，尚無較為便捷的方法。

耿正玲等研究發現，層間干擾是由水驅油過程中層間滲流阻力差異引起的［15-16］。通過計算層間滲流阻力差異的動態變化，可以預測不同時刻各層吸水量和產液量，分析開采過程中的層間干擾程度，為生產指標的劈分、開發措施的調整提供依據［17-20］。針對油水兩相滲流阻力計算，學者們給出了不同形式的求取公式與解法［21-24］，但存在計算復雜或精度不高的問題，不利于推廣應用。為此，基于Buckley-Leverett理論，推導了一維水驅油過程滲流阻力計算公式，并通過對油水總流度倒數與含水率導數變化關系的擬合，建立了兼顧精度與效率的兩相滲流阻力求取方法。利用單層模型數值模擬，驗證該計算方法的準確性，并將該算法應用于多層一維流動層間滲流阻力差異的求取，以期為多層油藏層間干擾的預測以及開發方案調整提供方法和依據。

1 滲流阻力計算方法

1.1 單層滲流阻力計算方法

考慮一維定截面積地層中水驅油過程，模型假設為：①非活塞式驅替，忽略重力和毛管壓力的影響，注采平衡。②油、水相流動均服從達西線性滲流定律。③不考慮巖石和流體壓縮性。

對于兩相區內任一含水飽和度為Sw的截面，通過該截面的瞬時流量可由兩相滲流達西定律得到：

該截面上滲流阻力為：

兩相區內滲流阻力需由（2）式對兩相區長度積分得到：

利用Buckley-Leverett 方程（簡稱B-L 方程），可以將位置轉化為含水飽和度。B-L方程給出任一含水飽和度平面的前進距離［25-26］為：

設兩相區前緣含水飽和度為f′(Sw)=f′(Swf)，則前緣相對于注入端前進距離為：

（4）式與（5）式相除并整理可得：

對（6）式微分可得：

采出端見水前，兩相區兩端含水飽和度分別為：在x=xf處，Sw=Swf；在x=xm處，Sw=1-Sor。將（7）式代入（3）式并變換積分上下限后，可得兩相區內滲流阻力為：

（8）式積分項中的分母部分實際為某一含水飽和度下油水總流度，是某一含水飽和度下地層流體總流動能力的表征，記為：

將（9）式代入（8）式，同時將積分項上下限互換以消去負號，可得：

設從注入端到采出端距離為f′(Sw)=f′(Swf)，在水驅前緣到達采出端前，即xf＜L時，巖心中存在純油相區與油水兩相區，此時總滲流阻力為兩區域滲流阻力之和，對其進行整理，可得：

當采出端見水后，為繼續分析注采井間油水流動情況，可看作水驅前緣穿過采出端繼續向前推進，即將注采井距視為水驅前緣前進距離的一部分（xf≥L）［25-26］。此時注采井間純油相區消失，兩相區滲流阻力即為總滲流阻力；同時前緣含水飽和度平面也消失，因此（11）式中的前緣含水飽和度Swf要相應變為采出端含水飽和度Swe，即：

求取總滲流阻力關鍵在于求出表達式（見水前為（11）式，見水后為（12）式）中的積分項，為便于求解，可將積分項轉化為以含水率導數為自變量的形式：

由于油水總流度倒數1/λt(Sw)與含水率導數f′(Sw)之間為隱函數關系，此類積分通常需要采用數值積分［25］或圖解法［26-27］得到。其中數值積分方法往往需要專業軟件，同時傳統的梯形法在數據點之間采用線性插值，可能存在一定誤差；圖解法則需要人工讀圖，工作量較大且不易保證結果準確。為此，筆者考慮利用多項式擬合得到兩者之間的顯式函數關系，即將二者的關系表示為：

（14）式積分得到：

將（15）式分別代入（11）式與（12）式，同時采用礦場單位，滲流阻力最終計算公式可表示為分段函數，即：

通過多項式擬合得到各擬合系數am值之后，可以由（16）式較為便捷地計算出任意給定含水率導數值對應的滲流阻力值；又由于含水率導數與累積注入量（以注入孔隙體積倍數表示）為倒數關系，最終可得到滲流阻力隨累積注入量的變化曲線。

1.2 層間阻力差異計算方法

基于建立的單層滲流阻力計算方法，可計算多層油藏水驅過程中各層滲流阻力變化。假設：①各層間無竄流。②各層流體性質及相對滲透率曲線相同。對（16）式分析可知，任一時刻單層滲流阻力由該層累積注水量確定，而各層的注水量又需由各層滲流阻力的相對大小確定，因此需采用迭代法計算，計算流程如下：①給定儲層物性參數滲透率K、孔隙度?、油水黏度μo與μw、注采井距L、截面積A及注水量步長qi。②由（5）式計算前緣位置。③比較xf與L，若xf＜L，由（16）式計算見水前各層滲流阻力Rtik，并根據Rtik值劈分注水量qi+1。④重復步驟②—③至xf≥L。⑤見水后，由f′(Swe)=L?A/Qt計算此時采出端含水率導數，由（16）式計算見水后各層滲流阻力Rtik，并根據Rtik值劈分注水量qi+1。⑥重復步驟⑤，直至達到設定的注水量上限，計算結束。

2 應用實例

2.1 單層滲流阻力計算與驗證

采用某中高滲透巖心的相對滲透率數據進行計算，相對滲透率曲線形態如圖1 所示。對應巖心滲透率為200 mD，其余參數設為：注采井距為100 m，地層截面積為1 m2，孔隙度為20%，地層原油黏度為10 mPa·s，注入水黏度為0.5 mPa·s。

圖1 油水相對滲透率曲線Fig.1 Oil-water relative permeability curves

首先，根據相對滲透率數據繪制含水率曲線，由切線法［25-27］求出水驅前緣含水飽和度。然后，取其右側數據點，計算出不同含水飽和度對應的總流度倒數與含水率導數。以含水率導數為橫坐標，以總流度倒數為縱坐標，使用多項式對兩者在平面直角坐標系內確定的散點進行擬合，得到總流度倒數與含水率導數的關系曲線。所選數據采用三次多項式即可達到較好的擬合效果（圖2）。

圖2 總流度倒數與含水率導數關系曲線Fig.2 Relationship between reciprocal of total mobility and derivative of water cut

給定不同含水率導數，即可求出對應的總流度倒數，進而計算出積分項，得到滲流阻力；同時，根據B-L 方程，可由含水率導數求出對應的注入量，得到滲流阻力隨注入量的變化規律。對比本文方法計算結果與數值積分和數值模擬結果（圖3）發現：三者計算得到的滲流阻力變化趨勢基本一致；但本文方法計算結果與數值模擬結果在滲流阻力大小和變化趨勢上均偏差更小，尤其在高-特高含水期。與數值積分方法相比，本文方法簡化了計算過程，并且具有更高的精度。本文方法計算結果較數值模擬結果下降略快，一方面可能是由于本文方法計算時忽略了巖石和流體的壓縮性，另一方面可能是由于水驅前緣讀取及曲線擬合過程中的誤差所致；基于B-L 方程計算得到的水驅指標變化快于數值模擬結果，這在其他研究中亦有發現［27］。同時在注水初期，由于數值模擬中能反映出井間壓力平衡的過程，使得折算出的滲流阻力在注水初期表現為短暫上升再下降的趨勢。但總體來說，本文方法與數值模擬結果較為吻合，相對誤差低于12%，說明本文方法能夠用于滲流阻力變化的預測與監控。

圖3 不同方法計算所得注水過程中滲流阻力變化對比Fig.3 Seepage resistance changes calculated by different methods during water flooding

分析滲流阻力的變化（圖3）發現，滲流阻力整體呈現下降最終逐漸減緩的趨勢，隨著注入量增加，依次分為線性、漸緩和近平緩3個下降階段。線性下降階段對應無水采油期，因為由（16）式可知，當水驅前緣未到達采出端時，f′(Sw)=f′(Swf)，滲流阻力實際只是前緣位置xf的函數，而前緣位置又與注入量成正比。當采出端見水后，隨著含水率上升，滲流阻力下降速率開始減小，進入漸緩下降階段，對應中-高含水期。隨著含水率進一步上升，進入近平緩下降階段，滲流阻力下降速率進一步降低，對應特高含水期。滲流阻力的下降速率與含水率的上升速率具有一定對應性，這是由于兩者都受到B-L 方程所描述的含水飽和度平面移動速率控制。

2.2 多層非均質油藏滲流阻力計算

以雙層介質滲流過程為例，根據計算流程計算水驅過程中各層滲流阻力的變化。同時，引入滲流阻力比值的概念，更加直觀地反映層間阻力差異，分析層間干擾程度的變化特征。

2.2.1 注水過程中滲流阻力變化特征

利用本文方法可以計算得到層間干擾程度的動態變化。設兩層滲透率分別為200 和100 mD，滲透率級差為2，注采井距為100 m，滲流面積為1 m2，孔隙度為20%，地層原油黏度為10 mPa·s，注入水黏度為0.5 mPa·s。分析計算得到的各層滲流阻力與層間滲流阻力比值（圖4）可知，隨著注水過程的進行，各層滲流阻力變化趨勢相近，但變化速率不同，高滲透層阻力低，下降快，低滲透層阻力高，下降慢。滲流阻力比值反映了兩層滲流阻力的差異，表現為一條近拋物線曲線，當注入量為2 PV 時，滲流阻力比值達到最大，由初始值2升高至最大值（約為3.7），層間干擾程度達到最大；當注入量大于2 PV時，滲流阻力比值呈緩慢下降趨勢，而此時高滲透層已進入特高含水期，注入水無效循環嚴重，該層滲流阻力近乎為定值。因此，對多層油藏應及早進行以分層注水為主的調控措施，抑制層間干擾，避免無效注水；而在層間滲流阻力差異達到最大值后，調控措施應以封堵高滲透層為主。

圖4 注水過程中各層滲流阻力及滲流阻力比值變化Fig.4 Changes of seepage resistance and the resistance ratio of each layer during water flooding

2.2.2 滲透率級差對層間阻力差異的影響

滲透率級差是層間干擾的主要影響因素［7-9］。控制低滲透層滲透率為100 mD，依次改變高滲透層滲透率為200，300 和400 mD，形成級差為2，3 和4的3 個模型，計算得到滲流阻力比值變化規律。結果（圖5）表明，滲透率級差越大，滲流阻力比值越大，層間滲流阻力差異越大。滲透率級差為2時，滲流阻力比值最大值約為3.7；滲透率級差為3 時，滲流阻力比值最大值約為8；滲透率級差為4 時，滲流阻力比值最大值超過12。隨著滲透率級差的增加，滲流阻力比值最大值迅速增加，反映層間干擾程度迅速增大，給油藏產油能力帶來不利影響。因此，對于多層合注油藏，及早實施調控措施抑制層間干擾至關重要。

圖5 不同滲透率級差儲層注水過程中滲流阻力比值變化Fig.5 Changes of seepage resistance ratios with different permeability ratios during water flooding in reservoirs

3 結論

基于Buckley-Leverett 理論，推導了兼顧計算效率與精度的一維地層兩相滲流阻力計算方法，并通過數值模擬驗證了該方法的有效性與實用性。結果表明，本文方法與數值模擬的計算結果在各開發階段均具有較好的一致性，相對誤差低于12%，能夠表征水驅油過程中滲流阻力變化規律。

在單層滲流阻力計算的基礎上，形成了基于層間阻力差異計算的多層油藏層間干擾程度快速預測方法，并應用該方法分析了中高滲透油藏層間干擾程度隨注水開發過程及滲透率級差的變化特征，為多層合注油藏開發措施的調整及調控時機的確定提供了依據。

在計算多層油藏層間阻力差異時假定各層間無竄流，隔層發育較為完善的油藏可以滿足這一假設。因此，在應用文中流程進行層間干擾程度預測時，需要加強對隔層發育狀況的認識，以確保預測結果的準確性。

符號解釋

am——多項式擬合系數；

A——地層截面積，m2；

i——注水量步長序號，i=1，2，3，…；

k——多層油藏層編號，k=1，2，3，…；

K——絕對滲透率，D；

Kro——油相相對滲透率，f；

Krw——水相相對滲透率，f；

L——注采井距，m；

m——多項式中擬合系數對應項的指數，m=0，1，2，3，…，n；

n——擬合得到的多項式最高次項指數，正整數；

p——給定時刻任一含水飽和度平面所在位置壓力，MPa；

qi——注水量步長，m3；

Q——兩相區內任一截面上的瞬時流量，m3/ks；

Qt——給定時刻累積注入量，m3；

R——兩相區內滲流阻力，（mPa·s）/（D·m）或MPa/（m3/d）；

Rt——注采井間總滲流阻力，（mPa·s）/（D·m）或MPa/（m3/d）；

Rtik——第i步第k層總滲流阻力，（mPa·s）/（D·m）或MPa/（m3/d）；

Sw——含水飽和度，f；

Swf——兩相區前緣含水飽和度，f；

Sor——殘余油飽和度，f；

Swe——采出端含水飽和度，f；

x——給定時刻任一含水飽和度平面所在位置，m；

xf——給定時刻兩相區前緣（前緣含水飽和度平面）所在位置，m；

xm——給定時刻兩相區末端（最大含水飽和度平面）所在位置，m；

?——地層孔隙度，f；

λt(Sw)——兩相區內油水總流度，f；

μo——地層原油黏度，mPa·s；

μw——注入水黏度，mPa·s。