縫洞型碳酸鹽巖油藏水平井分形非線性滲流

姜瑞忠，張春光，郜益華，耿艷宏，余 輝，李昊遠

（1.中國石油大學（華東）石油工程學院，山東青島 266580；2.中海油研究總院有限責任公司，北京 100028；3.西安石油大學地球科學與工程學院，西安 710000）

0 引言

縫洞型油藏的孔洞、基質及裂縫系統相互影響，為了描述其儲層復雜的孔隙結構，出現了三重介質模型的概念［1-2］。描述縫洞型碳酸鹽巖油藏的滲流特征，須充分考慮裂縫系統的分形特征、致密基質塊的非線性滲流特征，以及儲層的應力敏感性［3-5］。

國內外學者針對三重介質模型進行了深入的研究。姚軍等［6］創新變井筒儲存三重介質油藏試井理論，對介質間竄流、彈性儲容比及變井儲參數對壓力響應的影響進行了研究。Camacho-Velazquez等［7］對所建立的三重孔隙模型考慮了溶洞系統與井筒連通和不連通2 種情況下的解，并對樣板曲線進行了詳細的研究。李成勇等［8］建立了三重介質油藏的水平井各向異性滲流數學模型，并對水平井的壓力動態進行了系統的研究。張冬麗等［9］建立了三重介質數值試井模型，很好地反映出了縫洞型儲層的非均質性特征，該數值試井模型可進一步考慮縫洞型油藏中具體的裂縫和溶洞，進而推廣至兩相流試井。Gomez 等［10］利用全局最優化解法對縫洞型儲層試井曲線解釋方法進行研究，提高了計算效率和試井解釋精度。Li 等［11］對不同邊界條件下縫洞型儲層的動態特征進行了研究，形成了一套系統的刻畫縫洞型儲層動態特征的方法。不僅如此，分形理論在儲層滲流領域的應用也取得了較大的進展。Chang 等［12］提出了分形裂縫網絡嵌入基質的表征公式，并通過對擴散方程系數的適當修正來描述分形系統中的單相流動。Acuna 等［13］將分形理論應用于裂縫性儲層試井中，并基于試井解釋明確了分形參數的物理意義。Razminia 等［14］利用分形幾何和分數階微積分的概念分析了兩區復合儲層的數學模型，并推導了徑向復合系統中井底瞬態壓力解析解。張本健等［15］針對裂縫性氣藏試井滲流問題，使用分形網絡模擬氣藏裂縫系統并將分形網絡嵌入到氣藏中，計算出試井模型的流動動態特征典型曲線，結果表明分形網絡能夠很好地模擬裂縫性氣藏的滲流特征。

雖然目前對三重介質模型的研究較多，但仍存在無法精確表征儲層非均質性特征的問題，而分形理論的應用使得非均質儲層，尤其是裂縫的表征更靈活且更貼近實際，但該理論在縫洞型儲層中的應用研究較少且未成體系。此外，對縫洞型儲層致密基質塊非線性滲流特征的研究也不可忽視。研究低滲儲層最常用的滲流模型主要有擬啟動壓力梯度模型、分段模型及連續模型，這3 種模型均由不同數學函數擬合實驗數據所得，缺乏物理背景，表征不夠靈活。

針對縫洞型碳酸鹽巖儲層滲流模型目前存在的問題，綜合考慮裂縫分形特征和基質系統非線性滲流機理，完善非線性滲流理論，建立三重介質油藏水平井分形非線性滲流模型；基于有限元原理求解水平井動態壓力，并對比各模型的壓力動態曲線，總結滲流規律。在分析非線性參數、分形指數及水平井長度等相關參數的敏感性之后，將所建立模型進行實際井動態壓力擬合并解釋相關滲流參數，以期驗證模型的準確性及實用價值。

1 縫洞型碳酸鹽巖油藏滲流機理

1.1 基質非線性滲流

將低滲儲層微觀滲流機理、邊界層流動理論和毛細管滲流模型相結合，簡化并改進非線性滲流模型。考慮邊界層厚度與流體屈服應力值，在Hagen-Poisseuille 定律的基礎上進行修正［16］可得

式中:v為通過巖心的流速，cm/s；K為巖心滲透率，mD；μ為流體黏度，mPa·s；r0為毛管半徑，mm；δ為邊界層厚度，mm；τ0為流體屈服應力，MPa為壓力梯度，MPa/m。

對于相同毛管，δ/r0與壓力梯度成反比［17］，故假設；對于同一種流體，屈服應力值保持不變，假設8τ0/3r0=a2。將式（1）變形并簡化可得

式中:a1，a2均為實驗擬合得到的常數值。

式中:c1和c2是通過對實驗數據進行非線性擬合得到的參數，c1反映了流體的屈服應力以及邊界層對滲流的影響，c2反映了邊界層對滲流的影響。

當c1=0 時，式（3）簡化為達西模型；當c2=0 時，式（3）簡化為擬啟動壓力梯度模型。相當于啟動壓力梯度項。當v＞0 時最小啟動壓力梯度為啟動壓力梯度的存在與否以及取值大小均取決于c1，c2的值，即取決于礦場實際［18］。

1.2 裂縫分形特征

利用分形維數及異常擴散系數來表征分形裂縫中的滲透率和孔隙度，能夠精確地反映地下儲層復雜的裂縫形態。設分形體內流體聚集在某節點處，且節點體積處處相等，則滲流節點的數量可結合節點密度求出，即

式中:Df為分形維數；α為節點的數量；ds為譜維數；rw為井筒半徑，m；Kfw為生產井處的滲透率，mD；φfw為生產井處的孔隙度；Cft為裂縫綜合壓縮系數，MPa-1。

考慮維數為d的歐氏空間內嵌入分形體，可得到用分形維數表征的孔隙度表達式為

式中:B為描述對稱性的幾何常量；V為儲層總體積，m3；Vf為儲層孔隙體積，m3；Vs為儲集流體的節點體積，m3；r為距生產井的距離，m。

裂縫孔隙度主要受孔隙空間聚集方式（用分形維數表征）的影響，且與有效上覆壓力呈指數關系，所以考慮應力敏感性的裂縫分形孔隙度表達式［19-20］為

式中:Cf為裂縫孔隙壓縮系數，MPa-1；pi為原始地層壓力，MPa；pf為裂縫系統壓力，MPa。

裂縫滲透率不僅取決于孔隙空間的聚集方式，還與裂縫之間的連通性（用異常擴散系數表征）密不可分。儲層壓力變化會引起裂縫張合狀態的變化，從而影響儲層滲透率，在對應力敏感地層不穩定試井的研究中發現，滲透率與基巖上覆壓力之間也成指數關系。在分形孔隙度定義的基礎上引入異常擴散系數便可得到考慮應力敏感性的分形滲透率［19-20］，即

式中:αK為滲透率模數，MPa-1；θ為異常擴散系數。

2 縫洞型碳酸鹽巖分形油藏水平井非線性滲流模型

2.1 物理模型

圖1 物理模型示意圖Fig.1 Schematic diagram of the physical model

建立圓形三重介質碳酸鹽巖分形油藏水平井物理模型［21］（圖1），并設該模型符合以下假定條件:①模型由基質、裂縫和溶洞系統組成，基質塊向裂縫和溶洞系統發生擬穩態竄流，溶洞系統向裂縫發生擬穩態竄流；②外邊界為封閉或定壓條件，儲層厚度為h，原始壓力為pi，水平井半長為L，水平井縱向深度為zw，油藏半徑為re；③考慮儲層的表皮系數、井筒儲集系數以及裂縫滲透率的各向異性，裂縫系統水平滲透率為Kfh，垂直滲透率為Kfv，基質系統滲透率為Km，溶洞系統滲透率為Kv；④水平井位于油藏平面中心，產量為Q）；⑤忽略重力與毛管壓力的影響；⑥考慮裂縫系統的應力敏感性及裂縫的分形特征，基質向裂縫系統的竄流為非線性滲流。

2.2 數學模型

在笛卡爾坐標系中，將流體流動的質量守恒方程和狀態方程與非線性運動方程［式（3）］相結合，可得到縫洞型碳酸鹽巖油藏水平井分形非線性滲流模型［22］，即

式中:rD為無因次滲流半徑；LD為無因次水平井半長；hD為無因次儲層厚度；tD為無因次時間；pfD為無因次裂縫壓力；pmD為無因次基質壓力；pvD為無因次溶洞壓力；δLVDf為裂縫的無因次中間變量；δLVDv為溶洞的無因次中間變量；λmf為基質向裂縫竄流系數；λvf為溶洞向裂縫竄流系數；λmv為基質向溶洞竄流系數；ωf為裂縫彈性儲能比；ωv為溶洞彈性儲能比；下標D 代表無因次；下標m，f，v 分別代表基質、裂縫和溶洞系統。

2.3 模型求解

對于縫洞型碳酸鹽巖油藏水平井分形非線性滲流模型，首先采用隱式求解裂縫系統，應用Galerkin法，得到裂縫系統的有限元方程為

式中:Ni為形函數的分量，i=1，2，…，n。

以格林公式為基礎，通過分部積分可得到內部單元和封閉條件外邊界單元的有限元方程為

由此可得，單元有限方程的矩陣形式為

假設E為單位矩陣，其他各矩陣的表達形式為

對矩陣方程進一步簡化可得

其中Ke，Fe的表達形式為

式中:Ke為非線性系數矩陣；Fe為單元載荷向量；上標n表示時刻變量。

式（18）為三重介質分形非線性模型裂縫系統的有限元單元平衡方程，同理可采用顯式求解方法得到基質和溶洞系統的有限元單元平衡方程，且計算結果表明隱式求解裂縫系統、顯式求解基質和溶洞系統的方法符合計算精度。

考慮到無限導流能力模型比均勻流量模型求解困難，通常選取0.7L處為水平井在這2 種模型下的等價壓力點，從而可以借助均勻流量模型評價水平井井底壓力。通過把水平井考慮成單元源匯項并積分，應用Delta 函數將水平井劃分成n個節點，并進行無因次化得到單元源匯項有限元方程。

所求壓力解經過Laplace 變換后，通過Duhamel原理引入井筒儲集與表皮系數:

式中:S為表皮因數；s為Laplace 變量；CD為井筒儲集系數。

采用Stehfest 數值反演，得到考慮井筒儲集效應與表皮因數的水平井無因次井底壓力為

3 滲流規律

圖2 為碳酸鹽巖油藏水平井的達西、非線性以及分形非線性模型的壓力動態曲線。各基本參數取值如下:CD=100，S=0，ωf=0.01，ωv=0.1，λmf=0.001，λmv=0.000 5，λvf=0.1，LD=5，hD=200，αKD=0.02。假設非線性參數cD=c1D=c2D，分形指數β=2+θ-Df，則各模型其他參數的取值如下:達西模型cD=0，θ=0，Df=2；非線性流模型cD=0.3，θ=0，Df=2；分形非線性模型cD=0.5，θ=0.3，Df=1.9。在考慮儲層應力敏感性的基礎上，可以根據實際情況適當調整各模型的αKD值。

不同模型在不同流動階段曲線的形態各異，但壓力動態曲線仍然呈9 個特征流動階段（圖2）。以達西模型為例，具體流動階段劃分及不同流動階段壓力動態曲線特征［23-24］如表1 所列。通過分析和對比非線性、分形非線性模型的壓力動態曲線與達西模型的壓力動態曲線的異同點（圖2）可發現，非線性參數主要影響基質-裂縫竄流強度，而分形指數使得溶洞-裂縫竄流后的流動階段的壓力動態曲線發生上翹。

圖2 不同模型壓力動態曲線對比1.井筒儲存階段；2.表皮效應過渡階段；3.早期徑向流階段；4.早期線性流階段；5.中期徑向流階段；6.溶洞-裂縫竄流階段；7.溶洞-裂縫擬徑向流階段；8.基質-裂縫竄流階段；9.系統擬徑向流階段Fig.2 Comparison of dynamic pressure curves of different models

表1 不同流動階段壓力動態曲線特征Table 1 Characteristics of dynamic pressure curves in different flow stages

4 敏感性分析

首先對非線性參數和分形指數進行敏感性分析，各基本參數取值如下:CD=100，S=0，ωf=0.1，ωv=0.5，λm=0.001，λv=0.1，LD=5，hD=200，αKD=0.02。圖3 與圖4 分別為非線性參數cD和分形指數β對壓力動態曲線的影響。

非線性參數主要對基質向裂縫的竄流段有顯著影響，基本不影響溶洞向裂縫系統的竄流段。隨著非線性參數的增大，基質向裂縫的竄流強度減弱，持續時間變短。在壓力導數曲線上表現為第2個下凹段延遲出現，下凹程度減弱，但竄流段結束的時間不受影響。

圖3 非線性參數對壓力動態曲線的影響Fig.3 Effect of nonlinear parameters on dynamic pressure curves

圖4 分形指數對壓力動態曲線的影響Fig.4 Effect of fractal index on dynamic pressure curves

分形指數對壓力動態曲線的影響主要體現在溶洞向裂縫竄流后的流動階段，而對早期徑向流和早期線性流階段只產生微小的影響。從溶洞-裂縫擬徑向流階段開始，壓力及其導數曲線隨著時間的增加而逐漸上翹，且上翹幅度隨著分形指數的增大而加劇。分形指數并不影響各個竄流段出現的時間以及竄流強度，所以壓力導數曲線上下凹段的深度和持續時間均保持不變。

圖5 與圖6 分別為水平井長度LD和縱向位置zwD對壓力動態曲線的影響，各基本參數取值如下:CD=100，cD=0，β=0.1，S=0，ωf=0.01，ωv=0.1，λm=0.001，λv=0.1，hD=200，αKD=0.02。研究表明，水平井越長，壓力及其導數曲線越低，早期徑向流越明顯，從而后續流動階段出現的時間延遲。當水平井長度小到一定程度時，早期徑向流易被井筒儲集效應所掩蓋。考慮儲層分形和應力敏感性特征，不同長度水平井的壓力動態曲線并不在后期匯聚在一起，而是呈平行分布。在其他參數不變的情況下，水平井越靠近油層中部，即zwD的值越接近0.5，早期徑向流階段出現得越早，同時早期線性流持續時間越長，但水平井縱向上位置的變化并不影響其他流動階段。

圖5 水平井長度對壓力動態曲線的影響Fig.5 Effect of horizontal well length on dynamic pressure curves

圖6 水平井縱向位置對壓力動態曲線的影響Fig.6 Effect of horizontal well position on dynamic pressure curves

5 應用實例分析

為了驗證本文模型的正確性，進一步認識具有分形特征的縫洞型碳酸鹽巖油藏的滲流特征，對四川CN 油田S63 井的滲流參數進行了解釋。為了便于試井解釋，假定裂縫垂向滲透率與水平滲透率相等，即β=1。S63井的基礎參數取值如表2 所列。

各個模型的異同已在滲流規律及敏感性分析部分進行了詳細對比，但為方便對比、分析本文分形非線性模型的合理性，同時應用達西模型和本文模型進行試井擬合（圖7）。2 個模型的擬合結果均顯示壓力導數與壓力曲線分開。實際壓力導數曲線和理論曲線在前期擬合程度相對較低，這是由于實際儲層的井筒儲集系數發生變化所導致的，在實際試井過程中普遍存在類似的問題。在滲流的中后期，雖然實測數據不完整，壓力導數曲線仍明顯出現反映孔洞縫三重介質儲層滲流特征的2 個下凹段，且相對于達西模型，本文模型的擬合結果更加精確，理論計算值和實測值擬合得更好。選取本文的分形非線性模型進行S63 井的試井解釋，表3為該模型所解釋的滲流參數。

表2 S63 井基礎參數取值Table 2 Values of basic parameters of well S63

實際生產曲線的基質向裂縫的竄流段明顯比達西模型所呈現的下凹程度小，且試井解釋下的非線性參數為0.15，基質滲透率為1.52 mD，這都體現了基質確實存在非達西滲流。擬合結果表明，采用本文模型進行滲流參數解釋更加精確，這說明儲層考慮分形特征和非線性滲流后更符合油井的實際情況，可應用本文的分形非線性模型進行實際試井解釋及滲流規律分析。

圖7 實際算例壓力動態擬合曲線Fig.7 Dynamic pressure fitting curves of the actual example

表3 S63 井滲流參數解釋Table 3 Explanation of seepage parameters of well S63

6 結論

（1）引入分形理論描述裂縫系統的分形特征，同時應用非線性滲流新模型描述致密基質塊的滲流特征，在此基礎上建立縫洞型油藏水平井分形非線性滲流模型，并基于有限元方法求解水平井動態壓力。

（2）通過對比達西模型與分形非線性模型的異同點，總結縫洞型碳酸鹽巖油藏水平井的滲流規律，分析各個流動階段的特點，將滲流過程劃分為9 個流動階段，即井筒儲存、表皮效應過渡、早期徑向流、早期線性流、中期徑向流、溶洞向裂縫竄流、溶洞-裂縫擬徑向流、基質向裂縫竄流以及系統擬徑向流。

（3）非線性參數主要影響基質向裂縫系統的竄流段，隨著非線性參數的增大，基質向裂縫系統竄流強度減小，持續時間變短。分形指數不影響竄流出現的時間及強度，它對壓力動態曲線的影響主要體現在溶洞向裂縫竄流后的階段，壓力及導數曲線隨時間增加而逐漸上翹，且上翹幅度隨分形指數的增大而加劇。

（4）利用本文的分形非線性模型進行實際井動態壓力擬合并解釋相關滲流參數，考慮了儲層分形特征和非線性滲流后，模型更貼近礦場實際，可用其進行實際試井解釋及滲流規律分析。