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

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

0 引言

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

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

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

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

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

1.1 基質(zhì)非線性滲流

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

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

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

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

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

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

1.2 裂縫分形特征

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

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

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

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

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

式中:Cf為裂縫孔隙壓縮系數(shù)，MPa-1；pi為原始地層壓力，MPa；pf為裂縫系統(tǒng)壓力，MPa。

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

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

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

2.1 物理模型

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

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

2.2 數(shù)學模型

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

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

2.3 模型求解

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

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

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

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

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

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

其中Ke，F(xiàn)e的表達形式為

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

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

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

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

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

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

3 滲流規(guī)律

圖2 為碳酸鹽巖油藏水平井的達西、非線性以及分形非線性模型的壓力動態(tài)曲線。各基本參數(shù)取值如下: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。假設(shè)非線性參數(shù)cD=c1D=c2D，分形指數(shù)β=2+θ-Df，則各模型其他參數(shù)的取值如下:達西模型cD=0，θ=0，Df=2；非線性流模型cD=0.3，θ=0，Df=2；分形非線性模型cD=0.5，θ=0.3，Df=1.9。在考慮儲層應力敏感性的基礎(chǔ)上，可以根據(jù)實際情況適當調(diào)整各模型的αKD值。

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

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

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

4 敏感性分析

首先對非線性參數(shù)和分形指數(shù)進行敏感性分析，各基本參數(shù)取值如下: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 分別為非線性參數(shù)cD和分形指數(shù)β對壓力動態(tài)曲線的影響。

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

圖3 非線性參數(shù)對壓力動態(tài)曲線的影響Fig.3 Effect of nonlinear parameters on dynamic pressure curves

圖4 分形指數(shù)對壓力動態(tài)曲線的影響Fig.4 Effect of fractal index on dynamic pressure curves

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

圖5 與圖6 分別為水平井長度LD和縱向位置zwD對壓力動態(tài)曲線的影響，各基本參數(shù)取值如下: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。研究表明，水平井越長，壓力及其導數(shù)曲線越低，早期徑向流越明顯，從而后續(xù)流動階段出現(xiàn)的時間延遲。當水平井長度小到一定程度時，早期徑向流易被井筒儲集效應所掩蓋?？紤]儲層分形和應力敏感性特征，不同長度水平井的壓力動態(tài)曲線并不在后期匯聚在一起，而是呈平行分布。在其他參數(shù)不變的情況下，水平井越靠近油層中部，即zwD的值越接近0.5，早期徑向流階段出現(xiàn)得越早，同時早期線性流持續(xù)時間越長，但水平井縱向上位置的變化并不影響其他流動階段。

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

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

5 應用實例分析

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

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

表2 S63 井基礎(chǔ)參數(shù)取值Table 2 Values of basic parameters of well S63

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

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

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

6 結(jié)論

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

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

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

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