伍銳東， 張春光， 楊志興， 石美雪， 姜瑞忠

( 1. 中海石油(中國(guó))有限公司 上海分公司，上海 200335； 2. 中國(guó)石油大學(xué)(華東) 石油工程學(xué)院，山東 青島 266580 )

0 引言

低滲氣藏常采用酸化等措施改善近井區(qū)物性條件，加劇儲(chǔ)層的平面非均質(zhì)性，呈現(xiàn)復(fù)合氣藏的滲流特征。研究低滲復(fù)合氣藏的滲流特征需考慮低滲區(qū)氣體的非達(dá)西滲流、儲(chǔ)層的應(yīng)力敏感性和各向異性等因素[1-4]，采用傳統(tǒng)滲流模型描述此類氣藏存在局限性。

低滲氣藏非達(dá)西滲流曲線由非線性段和擬線性段組成，非線性段對(duì)開發(fā)的影響不可忽略[5]。馮文光[6]分析流固表面作用、孔喉結(jié)構(gòu)等因素對(duì)天然氣非達(dá)西滲流的影響。李允等[7]基于實(shí)驗(yàn)研究啟動(dòng)壓力梯度對(duì)氣藏開發(fā)的影響，并建立擬啟動(dòng)壓力梯度滲流模型?？紤]啟動(dòng)壓力梯度和應(yīng)力敏感性，馮青[8]建立氣水兩相數(shù)值試井模型，并對(duì)相關(guān)參數(shù)進(jìn)行敏感性分析。劉廣峰等[9]應(yīng)用冪律方程及擬線性方程的分段模型，描述流體的非達(dá)西滲流。歐陽(yáng)偉平等[10]基于三參數(shù)連續(xù)非線性滲流方程建立致密氣藏?cái)?shù)值試井模型，表征氣體的非線性滲流特征。低滲非達(dá)西滲流模型可歸結(jié)為3類：擬啟動(dòng)壓力梯度模型，忽略滲流的非線性段，放大滲流阻力；分段模型，涉及臨界點(diǎn)判斷，應(yīng)用復(fù)雜；連續(xù)非線性滲流模型，多基于實(shí)驗(yàn)擬合，表征形式多樣。為研究復(fù)合氣藏的滲流特征，許峰等[11]與孟凡坤等[12]分別建立雙重和三重介質(zhì)復(fù)合氣藏水平井不穩(wěn)定滲流模型，謝飛[13]與陳軍等[14]進(jìn)一步研究復(fù)合油氣藏水平井、垂直裂縫井的滲流模型。孫高飛[15]與黃雨等[16]分別考慮儲(chǔ)層的應(yīng)力敏感性和啟動(dòng)壓力梯度進(jìn)行低滲復(fù)合氣藏試井分析。

目前，常規(guī)復(fù)合氣藏的研究較為成熟，涉及低滲復(fù)合氣藏斜井非達(dá)西滲流的研究較為鮮見，常采用簡(jiǎn)化的擬啟動(dòng)壓力梯度模型進(jìn)行表征，未綜合考慮儲(chǔ)層的應(yīng)力敏感性和各向異性，不能充分體現(xiàn)低滲復(fù)合氣藏的滲流特征?；诘蜐B氣藏非達(dá)西滲流和應(yīng)力敏感性微觀機(jī)理，筆者考慮滲透率的各向異性，完善滲流運(yùn)動(dòng)方程，建立低滲復(fù)合氣藏斜井非達(dá)西滲流數(shù)學(xué)模型，采用有限元方法進(jìn)行Matlab編程求解，繪制井底壓力動(dòng)態(tài)曲線，對(duì)比不同模型滲流規(guī)律并對(duì)相關(guān)參數(shù)進(jìn)行敏感性分析，應(yīng)用模型進(jìn)行實(shí)際斜井動(dòng)態(tài)壓力擬合并解釋相關(guān)儲(chǔ)層參數(shù)，驗(yàn)證模型實(shí)用性。

1 試井模型建立與求解

1.1 物理模型

在低滲氣藏開發(fā)過程中，受工藝措施或氣藏特點(diǎn)影響，易形成物性差異較大的復(fù)合氣藏。

建立低滲復(fù)合氣藏斜井物理模型(見圖1)。假設(shè)條件：(1)氣藏頂?shù)追忾]不滲透，氣藏半徑為re，地層等厚且厚度為h；(2)儲(chǔ)層劃分為內(nèi)、外區(qū)，內(nèi)區(qū)半徑為r1，外區(qū)半徑為r2，地層初始擬壓力為mi；(3)斜井貫穿整個(gè)氣藏，井斜角為θ，井長(zhǎng)度為L(zhǎng)，井中心的縱向深度為zw，以恒定地面產(chǎn)量qsc生產(chǎn)；(4)內(nèi)、外區(qū)交界面不存在附加壓力降，等溫滲流，流體為單相微可壓縮，忽略重力和毛管力；(5)考慮井筒儲(chǔ)集系數(shù)、表皮因數(shù)及滲透率的各向異性，氣藏水平方向初始滲透率為Knhi，縱向初始滲透率為Knvi(n=1,2，1代表內(nèi)區(qū)，2代表外區(qū))；(6)內(nèi)區(qū)流體遵循達(dá)西滲流，外區(qū)流體考慮非達(dá)西滲流及儲(chǔ)層應(yīng)力敏感性。

圖1 低滲復(fù)合氣藏斜井物理模型Fig.1 Physical model of inclined well in the low-permeability composite gas reservoir

1.2 運(yùn)動(dòng)方程

低滲氣藏滲透率低、束縛水飽和度高，氣相吸附作用增強(qiáng)，孔道壁出現(xiàn)邊界層，氣源壓力較低時(shí)氣體出現(xiàn)明顯的非達(dá)西滲流。低滲氣藏非線性滲流模式見圖2，其中δ為表邊界層厚度，r0為毛細(xì)管半徑，ν為通過巖心的流速，ζ*為最小啟動(dòng)壓力梯度，ζ為擬啟動(dòng)壓力梯度，p為壓力。低滲儲(chǔ)層流體滲流速度與壓力梯度的關(guān)系曲線見圖2(b)，由非線性段ζ*α與擬線性段αβ組成，α為臨界點(diǎn)，η為臨界點(diǎn)對(duì)應(yīng)的壓力梯度，β為滲流曲線的末端點(diǎn)。

圖2 低滲氣藏非線性滲流模式Fig.2 The nonlinear flow model of low-permeability gas reservoir

從微觀角度解釋非達(dá)西滲流的產(chǎn)生原因，考慮剪切應(yīng)力及邊界層厚度的影響，修正Hagen-Poisseuille方程[17]：

(1)

式中：K為巖心滲透率；μ為流體黏度；τ0為流體屈服應(yīng)力；δ/r0表征邊界層的影響，8τ0/3r0表征流體屈服應(yīng)力的影響。

通過微觀實(shí)驗(yàn)擬合，式(1)變形簡(jiǎn)化，得到宏觀上的連續(xù)非線性滲流運(yùn)動(dòng)方程：

(2)

式中：c1、c2為通過實(shí)驗(yàn)擬合得到的參數(shù)。

為增強(qiáng)連續(xù)非線性滲流運(yùn)動(dòng)方程調(diào)整的靈活性，令a=-c2/c1，b=1/c1，則式(2)簡(jiǎn)化為

(3)

式中：a為非線性因數(shù)，反映流體屈服應(yīng)力及邊界層對(duì)滲流的影響；b為啟動(dòng)壓力梯度的倒數(shù)。

當(dāng)0

低滲儲(chǔ)層具有較強(qiáng)的應(yīng)力敏感性，隨氣藏開采的進(jìn)行，巖石的有效應(yīng)力不斷增加，滲透率不斷減小。為描述滲透率隨壓力的變化關(guān)系，采用PEDROSA O A定義的滲透率模量[18]：

(4)

式(4)積分可得

K=Kie-γ(pi-p),

(5)

式中：Ki為儲(chǔ)層初始滲透率；γ為滲透率模量；pi為原始地層壓力。

綜上，考慮低滲氣藏應(yīng)力敏感性和非達(dá)西滲流的運(yùn)動(dòng)方程為

(6)

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

1.3.1 模型建立

基于物理模型假設(shè)，聯(lián)立非達(dá)西運(yùn)動(dòng)方程、質(zhì)量守恒方程及狀態(tài)方程，并引入擬壓力函數(shù)可得低滲復(fù)合氣藏斜井非達(dá)西滲流無(wú)因次微分方程：

(7)

(8)

式(7-8)中：m1D、m2D分別為量綱一的內(nèi)區(qū)和外區(qū)壓力；xD、yD、zD為量綱一的方向坐標(biāo)；tD為量綱一的時(shí)間；rD、r1D分別為量綱一的滲流半徑和內(nèi)區(qū)半徑；γmD為按擬壓力定義的量綱一的滲透率模量；η1,2為內(nèi)外區(qū)導(dǎo)壓系數(shù)比；λmlD為量綱一的非線性滲流中間變量，l=x,y,z，分別代表x、y、z方向變量。

初始條件為

m1D|tD=0=m2D|tD=0。

(9)

定產(chǎn)內(nèi)邊界條件為

(10)

式中：zwD為量綱一的斜井中心縱向深度；εD為量綱一的微變量。

頂?shù)追忾]外邊界條件為

(11)

水平封閉外邊界條件為

(12)

式中：reD為量綱一的最大邊界半徑。

水平定壓外邊界條件為

m2D|rD=reD=0。

(13)

界面連接條件為

m1D(r1D,tD)=m2D(r1D,tD),

(14)

(15)

式中：M1,2為內(nèi)外區(qū)流度比。

1.3.2 模型求解

非線性滲流數(shù)學(xué)模型解析求解較復(fù)雜，采用有限元法進(jìn)行數(shù)值求解[19]。應(yīng)用Galerkin法得到不存在源匯項(xiàng)時(shí)外區(qū)的有限元方程為

(16)

式中：Ni為形函數(shù)的分量。

通過格林公式分部積分，得到內(nèi)部單元和具有封閉條件外邊界單元的有限元方程：

(17)

將有限元方程變換成矩陣形式：

(18)

對(duì)有限元方程的矩陣形式進(jìn)一步簡(jiǎn)化得

(19)

(20)

(21)

式中：Ke為非線性系數(shù)矩陣；Fe為單元載荷向量。

式(19)為該模型外區(qū)的有限元單元平衡方程，采用相同方法可得到內(nèi)區(qū)的有限元單元平衡方程。

無(wú)限導(dǎo)流模型更符合斜井的實(shí)測(cè)壓力分布，但求解困難，通常選取均勻線源模型等價(jià)壓力點(diǎn)表征斜井壓力分布[20-21]。采用CINCO L H等[20]提出的方法描述井底壓力，井底等效壓力點(diǎn)取值為

xD=0.3LDsinθw,yD=1,zD=-0.2LDcosθw,

(22)

應(yīng)用Delta函數(shù)可得到與井位置有關(guān)的線源/匯公式，把斜井所在節(jié)點(diǎn)考慮成單元內(nèi)源匯項(xiàng)，對(duì)井節(jié)點(diǎn)所在單元積分并無(wú)因次化，得到單元源匯項(xiàng)的有限元方程為

(23)

式中：d為線源劃分的節(jié)點(diǎn)數(shù)。

通過Laplace變換并根據(jù)杜哈美原理，得到考慮井筒儲(chǔ)集效應(yīng)與表皮因數(shù)影響的井底壓力解為

(24)

式中：S為表皮因數(shù)；s為L(zhǎng)aplace變量；CD為井筒儲(chǔ)集系數(shù)。

采用Stehfest數(shù)值反演，得到考慮井筒儲(chǔ)集效應(yīng)與表皮因數(shù)的斜井井底壓力解為

(25)

(26)

2 滲流規(guī)律分析

2.1 流態(tài)劃分

利用Matlab編程進(jìn)行模型求解，繪制模型擬壓力動(dòng)態(tài)曲線，并與達(dá)西滲流模型和擬啟動(dòng)壓力梯度模型進(jìn)行對(duì)比[22-23](見圖3)。模型基本參數(shù)：CD=100，S=1，hD=200，θ=60°，r1D=15，η1,2=0.2，γmD=0。非線性滲流模型考慮非線性特征和啟動(dòng)壓力梯度，取a=0.5，bmlD=100；擬啟動(dòng)壓力梯度模型僅考慮啟動(dòng)壓力梯度，取a=0，bmlD=100。兩個(gè)模型的量綱一的擬啟動(dòng)壓力梯度為0.01。

根據(jù)擬壓力動(dòng)態(tài)曲線特征劃分7個(gè)流動(dòng)階段：①井筒儲(chǔ)存，擬壓力及其導(dǎo)數(shù)曲線呈斜率為1的重合直線；②表皮效應(yīng)過渡，擬壓力導(dǎo)數(shù)曲線呈駝峰狀；③井斜角控制，隨井斜角的增大，逐漸呈水平井的滲流特征，井斜角大于60° 后逐漸出現(xiàn)早期垂向徑向流和中期線性流，導(dǎo)數(shù)曲線分別呈水平直線和斜率為0.5的直線；④內(nèi)區(qū)擬徑向流，擬壓力導(dǎo)數(shù)曲線呈值為0.5的水平直線；⑤內(nèi)外區(qū)擬徑向流過渡，擬壓力導(dǎo)數(shù)曲線斜率與內(nèi)外區(qū)導(dǎo)壓系數(shù)比有關(guān)；⑥外區(qū)擬徑向流，不考慮非達(dá)西滲流與應(yīng)力敏感性等因素，擬壓力導(dǎo)數(shù)曲線呈水平直線；⑦邊界影響，受封閉邊界影響，擬壓力及其導(dǎo)數(shù)曲線劇烈上翹(見圖3)。

圖3 不同模型擬壓力動(dòng)態(tài)曲線Fig.3 Dynamic pseudo-pressure curves of different models

僅考慮外區(qū)的低滲特性，擬啟動(dòng)壓力梯度模型和非線性滲流模型的擬壓力及其導(dǎo)數(shù)曲線在外區(qū)擬徑向流階段發(fā)生上翹，非線性滲流模型的上翹幅度明顯低于擬啟動(dòng)壓力梯度模型的。

2.2 敏感性分析

分別對(duì)非線性因數(shù)a、井斜角θ、內(nèi)區(qū)半徑r1D及內(nèi)外區(qū)導(dǎo)壓系數(shù)比η1,2進(jìn)行敏感性分析[24-26]。

2.2.1 非線性因數(shù)

非線性因數(shù)對(duì)擬壓力動(dòng)態(tài)曲線的影響見圖4。由圖4可知，受非線性滲流和擬啟動(dòng)壓力梯度的影響，擬壓力動(dòng)態(tài)曲線在外區(qū)擬徑向流階段發(fā)生上翹。隨非線性因數(shù)的增大，非線性段曲線的斜率減小，最小啟動(dòng)壓力梯度逐漸趨于0，非線性滲流帶來的附加滲流阻力減小，曲線上表現(xiàn)為外區(qū)擬徑向流階段上翹時(shí)間延遲，上翹幅度減小。非線性因數(shù)對(duì)低滲儲(chǔ)層生產(chǎn)動(dòng)態(tài)的影響明顯，忽略非線性滲流特征的達(dá)西滲流模型或擬啟動(dòng)壓力梯度模型，難以客觀描述儲(chǔ)層壓力的動(dòng)態(tài)特征。

圖4 非線性因數(shù)對(duì)擬壓力動(dòng)態(tài)曲線的影響Fig.4 Effect of the nonlinear parameter on dynamic pseudo-pressure curve

2.2.2 井斜角

井斜角對(duì)擬壓力動(dòng)態(tài)曲線的影響見圖5。由圖5可知，在井斜角控制階段，壓力波尚未傳遠(yuǎn)，井斜角對(duì)擬壓力動(dòng)態(tài)曲線的影響較為明顯。井斜角較小時(shí)，擬壓力動(dòng)態(tài)曲線呈近似直井的特征；隨井斜角增大，擬壓力動(dòng)態(tài)曲線逐漸呈近似水平井的特征。當(dāng)井斜角大于60°后逐漸出現(xiàn)早期垂向徑向流和中期線性流，其擬壓力導(dǎo)數(shù)曲線分別呈一條水平線和斜率為0.5的直線。隨井斜角的增大，井長(zhǎng)度逐漸增加，早期徑向流持續(xù)時(shí)間變長(zhǎng)，定產(chǎn)生產(chǎn)時(shí)產(chǎn)量均勻分布，井底壓力降低，擬壓力及其導(dǎo)數(shù)曲線下移。不考慮儲(chǔ)層應(yīng)力敏感性和非達(dá)西滲流時(shí)，擬壓力導(dǎo)數(shù)曲線在內(nèi)區(qū)和外區(qū)的擬徑向流階段匯聚成一條直線。

圖5 井斜角對(duì)擬壓力動(dòng)態(tài)曲線的影響Fig.5 Effect of the well angle on dynamic pseudo-pressure curve

2.2.3 內(nèi)區(qū)半徑

內(nèi)區(qū)半徑對(duì)擬壓力動(dòng)態(tài)曲線的影響見圖6。由圖6可知，內(nèi)區(qū)半徑影響內(nèi)區(qū)擬徑向流的持續(xù)時(shí)間和內(nèi)外區(qū)擬徑向流過渡階段的開始時(shí)間：內(nèi)區(qū)半徑越小，內(nèi)區(qū)擬徑向流的持續(xù)時(shí)間越短，過渡階段出現(xiàn)的時(shí)間越早，當(dāng)內(nèi)區(qū)半徑小到一定程度時(shí)，壓力波迅速傳播到內(nèi)邊界，曲線上不能體現(xiàn)內(nèi)區(qū)擬徑向流的特征。2.2.4 內(nèi)外區(qū)導(dǎo)壓系數(shù)比

圖6 內(nèi)區(qū)半徑對(duì)擬壓力動(dòng)態(tài)曲線的影響Fig.6 Effect of the inner radius on dynamic pseudo-pressure curve

內(nèi)外區(qū)導(dǎo)壓系數(shù)比對(duì)擬壓力動(dòng)態(tài)曲線的影響見圖7。由圖7可知，導(dǎo)壓系數(shù)比主要影響內(nèi)外區(qū)擬徑向流過渡階段曲線的斜率：內(nèi)外區(qū)導(dǎo)壓系數(shù)比小于1.0時(shí)，即外區(qū)滲透性好于內(nèi)區(qū)的，內(nèi)外區(qū)擬徑向流過渡階段的擬壓力導(dǎo)數(shù)曲線呈下凹的特征；內(nèi)外區(qū)導(dǎo)壓系數(shù)比大于1.0時(shí)，即內(nèi)區(qū)滲透性好于外區(qū)的，過渡階段擬壓力導(dǎo)數(shù)曲線呈上凸的特征。隨內(nèi)外區(qū)導(dǎo)壓系數(shù)比的增大，外區(qū)的滲流能力相對(duì)逐漸變差，滲流阻力相對(duì)增大，內(nèi)外區(qū)擬徑向流過渡階段的擬壓力導(dǎo)數(shù)曲線斜率逐漸增大，外區(qū)擬徑向流階段的擬壓力及其導(dǎo)數(shù)曲線位置也越高。

圖7 內(nèi)外區(qū)導(dǎo)壓系數(shù)比對(duì)擬壓力動(dòng)態(tài)曲線的影響Fig.7 Effect of the pressure conduction coefficient ratio on dynamic pseudo-pressure curve

3 應(yīng)用實(shí)例對(duì)比

東海XH氣田A2井的測(cè)壓數(shù)據(jù)表現(xiàn)出明顯的低滲復(fù)合氣藏斜井的滲流特征，分別應(yīng)用常規(guī)試井模型和非達(dá)西試井模型對(duì)A2井進(jìn)行解釋(見圖8)。由圖8可知，常規(guī)試井模型僅考慮氣藏的徑向復(fù)合特征，不能體現(xiàn)外區(qū)擬徑向流階段的低滲特征；非達(dá)西試井模型考慮儲(chǔ)層的低滲特性，精細(xì)刻畫各流動(dòng)階段的滲流特征，模型與實(shí)測(cè)壓力數(shù)據(jù)擬合效果更好，增加非線性因數(shù)、啟動(dòng)壓力梯度等對(duì)擬合結(jié)果的約束，有效解決儲(chǔ)層參數(shù)多解性強(qiáng)的問題。A2井的基礎(chǔ)參數(shù)及試井解釋結(jié)果見表1。

圖8 A2井壓力動(dòng)態(tài)擬合曲線Fig.8 Dynamic pressure fitting curves of well A2

表1 A2井基礎(chǔ)參數(shù)及試井解釋結(jié)果

由A2井的壓力動(dòng)態(tài)曲線可知，該井的井斜角較小，無(wú)早期徑向流階段，類似直井的滲流特征。滲流后期壓力傳播到斷層邊界，壓力及其導(dǎo)數(shù)曲線發(fā)生劇烈上翹。考慮低滲氣藏非線性滲流和儲(chǔ)層應(yīng)力敏感性，改善該井外區(qū)擬徑向流的動(dòng)態(tài)壓力擬合效果。非達(dá)西試井模型解釋外區(qū)的滲透率為5.9×10-3μm2，表皮因數(shù)為8.3；常規(guī)試井模型解釋的外區(qū)滲透率為3.2×10-3μm2，表皮因數(shù)為7.9。

基于產(chǎn)能試井實(shí)測(cè)數(shù)據(jù)，采用二項(xiàng)式產(chǎn)能方程計(jì)算該井的無(wú)阻流量為47.85×104m3/d；基于常規(guī)試井模型和非達(dá)西試井模型解釋的參數(shù)，計(jì)算該井的無(wú)阻流量分別為32.49×104m3/d和45.44×104m3/d，非達(dá)西試井模型解釋的參數(shù)更貼近儲(chǔ)層實(shí)際(見圖9)，其中Qg為測(cè)試流量。忽略低滲氣藏滲流特性的常規(guī)試井模型對(duì)滲透率的解釋與氣藏實(shí)際的有效滲透率偏差較大；非達(dá)西試井模型的建立與應(yīng)用可有效改進(jìn)低滲復(fù)合氣藏參數(shù)解釋的精度，對(duì)低滲氣藏的壓力與產(chǎn)能評(píng)價(jià)具有重要意義。

圖9 模型解釋結(jié)果對(duì)比Fig.9 Comparison of different model interpretation results

4 結(jié)論

(1)基于低滲氣藏非達(dá)西滲流和應(yīng)力敏感性的微觀機(jī)理，完善滲流運(yùn)動(dòng)方程，建立低滲復(fù)合氣藏斜井非達(dá)西滲流數(shù)學(xué)模型；利用等價(jià)壓力點(diǎn)處理內(nèi)邊界條件，采用有限元方法求解井底壓力。

(2)繪制并對(duì)比不同模型的井底壓力動(dòng)態(tài)曲線，根據(jù)曲線特征劃分7個(gè)流動(dòng)階段，即井筒儲(chǔ)存、表皮效應(yīng)過渡、井斜角控制、內(nèi)區(qū)擬徑向流、內(nèi)外區(qū)擬徑向流過渡、外區(qū)擬徑向流及邊界影響階段。

(3)非線性因數(shù)增大，擬壓力動(dòng)態(tài)曲線在外區(qū)擬徑向流階段上翹幅度增大；井斜角大于60°后，逐漸出現(xiàn)早期垂向徑向流段，且隨井斜角增大持續(xù)時(shí)間延長(zhǎng)；內(nèi)區(qū)半徑增大，內(nèi)區(qū)擬徑向流階段持續(xù)時(shí)間延長(zhǎng)；內(nèi)外區(qū)導(dǎo)壓系數(shù)比大于1.0時(shí)，過渡流段斜率為正數(shù)且隨導(dǎo)壓系數(shù)比的增大而增大。

(4)低滲復(fù)合氣藏非達(dá)西滲流模型相較于傳統(tǒng)試井模型壓力動(dòng)態(tài)曲線擬合更好，增加非線性因數(shù)等對(duì)解釋結(jié)果的約束，降低參數(shù)的多解性；兩種模型解釋的滲透率相差明顯，非達(dá)西滲流模型可以合理評(píng)估低滲氣藏的產(chǎn)能。