煤層氣藏多分支水平井非穩(wěn)態(tài)產(chǎn)能模型

姜瑞忠，劉秀偉，王 星，郜益華，高 岳，黃巖松

（1.中國石油大學(xué)（華東）石油工程學(xué)院，山東青島 266580；2.中海油研究總院有限責(zé)任公司，北京 100028；3.延長油田股份有限公司志丹采油廠，陜西延安 717500）

根據(jù)國際能源機(jī)構(gòu)（IEA）的統(tǒng)計(jì)資料顯示，世界上煤層氣儲量約為260×1013m3，其中，中國的煤層氣儲量超過10×1013m3［1-2］。煤層氣作為一種優(yōu)質(zhì)的清潔能源，是國家大力提倡開發(fā)的新能源之一，其開發(fā)利用對緩解中國的能源緊張形勢具有重大的戰(zhàn)略意義。由于煤層氣藏儲層的低孔、低滲透特性，常規(guī)直井開采效果不好，而多分支水平井能夠增加井筒與煤層氣藏儲層的接觸面積，增加各個(gè)分支的波及面積和泄氣面積，從而大大提高單井產(chǎn)量和產(chǎn)氣速度［3-5］，又由于鉆井技術(shù)的飛速發(fā)展和成本的降低，多分支水平井技術(shù)在煤層氣藏中的應(yīng)用較為廣泛。

諸多中外學(xué)者依據(jù)不同井型對煤層氣藏的生產(chǎn)動態(tài)開展研究。張先敏等基于分形理論，考慮了低滲透煤層氣藏啟動壓力梯度的影響，建立了分形煤層氣藏直井滲流數(shù)學(xué)模型，分析了分形維數(shù)等對生產(chǎn)動態(tài)的影響［6］。宋洪慶等基于滲流力學(xué)基本理論，推導(dǎo)了耦合煤層氣藏解吸效應(yīng)和低速非達(dá)西滲流特征的壓力和產(chǎn)能特征方程，求解給出了徑向流解析解［7］。羅毅等采用裂縫波及區(qū)和未波及區(qū)等復(fù)合兩區(qū)模型，推導(dǎo)了煤層氣藏壓裂水平井的產(chǎn)能公式，分析了水平長度、儲層厚度以及縫網(wǎng)波及寬度等對產(chǎn)能的影響［8］。ZHAO 等建立了煤層氣藏?zé)o限導(dǎo)流能力壓力直井的滲流數(shù)學(xué)模型，同樣將整個(gè)煤層氣藏分為復(fù)合兩區(qū)，包括受壓裂井影響形成的壓裂改造的內(nèi)區(qū)和未壓裂改造的外區(qū)，求解繪制了產(chǎn)量遞減動態(tài)曲線圖，同時(shí)對比了擬穩(wěn)態(tài)和非穩(wěn)態(tài)的影響［9］。ZHANG 等建立了復(fù)合兩區(qū)煤層氣藏水平井滲流數(shù)學(xué)模型，通過求解井底壓力和產(chǎn)量，分別繪制了井底壓力動態(tài)曲線和產(chǎn)量動態(tài)曲線圖，進(jìn)行了相關(guān)的敏感性參數(shù)分析［10］。

雖然在煤層氣藏開采中多分支水平井技術(shù)應(yīng)用廣泛，但是相關(guān)煤層氣藏多分支水平井非穩(wěn)態(tài)產(chǎn)能模型的研究較少［11］。同時(shí)，大量研究發(fā)現(xiàn)，煤層氣藏儲層滲透率往往存在著較強(qiáng)的應(yīng)力敏感特性［12-16］，在進(jìn)行煤層氣藏生產(chǎn)動態(tài)研究中，應(yīng)考慮到該特性。為此，筆者建立了煤層氣藏多分支水平井非穩(wěn)態(tài)產(chǎn)能模型，耦合了割理系統(tǒng)的滲透率應(yīng)力敏感特性和基質(zhì)系統(tǒng)的非穩(wěn)態(tài)擴(kuò)散，求解繪制出產(chǎn)量遞減動態(tài)曲線，并分析了相關(guān)因素的影響規(guī)律。

1 模型建立

1.1 物理模型

假設(shè)在無限大頂?shù)走吔绶忾]的地層中心，存在著以定產(chǎn)量生產(chǎn)的一口多分支水平井（圖1），該多分支水平井由1 個(gè)水平井筒和N個(gè)分支井筒組成，每個(gè)分支井筒的長度為La，各個(gè)分支井筒之間的距離為d，各個(gè)分支井筒與水平井筒之間的夾角為χa。考慮煤層氣藏的吸附解吸和擴(kuò)散特性，分別用LANGMUIR 等溫吸附理論［17］和FICK 第二擴(kuò)散定律進(jìn)行表征，氣體流動途徑如圖2所示。同時(shí)，忽略重力、毛細(xì)管力、黏滯力以及溫度變化，假設(shè)分支井筒和水平井筒導(dǎo)流能力無限大。

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

圖2 氣體流動途徑示意Fig.2 Schematic of gas flow path

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

為了方便模型推導(dǎo)和求解，在此引入一系列無因次變量（表1），從而建立點(diǎn)源數(shù)學(xué)模型進(jìn)行求解。

1.2.1 割理系統(tǒng)

結(jié)合質(zhì)量守恒方程、狀態(tài)方程、運(yùn)動方程以及PEDROSA［18］提出的表征滲透應(yīng)力敏感特性的滲透率模量表達(dá)式，推導(dǎo)得到割理系統(tǒng)的流動方程為：

結(jié)合擬壓力和按擬壓力定義的滲透率模量，按照表1 中的無因次變量進(jìn)行無因次化之后，得到割理系統(tǒng)的無因次化流動方程為：

表1 無因次變量Table1 Dimensionless variables

1.2.2 基質(zhì)系統(tǒng)

采用FICK 第二定律來描述基質(zhì)系統(tǒng)的煤層氣非穩(wěn)態(tài)擴(kuò)散現(xiàn)象，假設(shè)基質(zhì)塊為球形，球形半徑為R，徑向距離為rm（圖3）［19-20］。

圖3 割理和基質(zhì)系統(tǒng)示意Fig.3 Schematic of cleat and matrix systems

結(jié)合表1 中的無因次變量，基質(zhì)系統(tǒng)中的非穩(wěn)態(tài)擴(kuò)散方程可以表示為：

初始條件為：

外邊界條件為：

球形對稱內(nèi)邊界條件為：

2 數(shù)學(xué)模型求解

對于數(shù)學(xué)模型的求解，整體思路為：先求解出基質(zhì)系統(tǒng)的氣體體積濃度表達(dá)式，然后代入割理系統(tǒng)的流動方程中，通過一系列數(shù)學(xué)方法進(jìn)行求解。

2.1 基質(zhì)系統(tǒng)求解

將基質(zhì)系統(tǒng)擴(kuò)散方程進(jìn)行Laplace變換后，可以得到：

為了求解方便，在此使用變量替換法，令：

則（12）式變?yōu)椋?/p>

（14）式為二階常系數(shù)微分方程組，可以求解得到其通解為：

2.2 割理系統(tǒng)求解

由于割理系統(tǒng)流動方程的非線性較強(qiáng)，在此采用Pedrosa變換和攝動變換消除非線性，其表達(dá)式分別為：

通過Pedrosa 變換和攝動變換后，進(jìn)行Laplace變換，可以得到：

將（17）式代入（22）式，可得：

對（23）式進(jìn)行有限余弦變換，則：

（25）式的通解為：

由外邊界條件可得A=0，則（27）式可寫為：

結(jié)合外邊界條件可以求得：

進(jìn)而點(diǎn)源函數(shù)解為：

結(jié)合有限余弦反演得到Laplace 空間下的點(diǎn)源函數(shù)解為：

為了求解得到井底壓力，在此利用數(shù)值離散和疊加原理的方法，將每個(gè)分支井筒劃分為M個(gè)微元段，同時(shí)微元段線流量密度為qa，b。這樣第a個(gè)分支的第b個(gè)微元段在任意一點(diǎn)形成的壓力均可寫為：

這樣，所有微元段在第α個(gè)分支第β個(gè)微元段所形成的壓力為井底壓力（因?yàn)榧僭O(shè)分支井筒和水平井筒導(dǎo)流能力無限大）：

另外，因各微元段流量總和為水平井筒總流量，有：

（33）式和（34）式可以構(gòu)成一個(gè)矩陣方程組，求解可得井底壓力，進(jìn)而結(jié)合杜哈美原理可得耦合表皮系數(shù)和無因次井筒儲集系數(shù)的井底壓力解為：

然后，通過Stehfest數(shù)值反演［21］和（35）式可以得到真實(shí)空間下的井底壓力解為：

依據(jù)VAN EVERDINGEN 等的研究成果［22］，可以得到Laplace空間下產(chǎn)量的表達(dá)式為：

3 結(jié)果與討論

3.1 模型驗(yàn)證

將本文模型中的單相氣換為單相油，將模型中的f（u）替換為文獻(xiàn)［23］中的f（u），并與其中的模型進(jìn)行對比（圖4）?；緟?shù)取值為：S=0.01，CD=1×10-6，LaD=1，dD=10，N=4，ω=0.1，λ=1×10-4，χa=60°，γpD=0。從圖4 中可以看出，本文模型計(jì)算結(jié)果與文獻(xiàn)［23］模型結(jié)果吻合較好，從而驗(yàn)證了本文模型的可靠性。

圖4 本文模型與文獻(xiàn)［23］中模型的對比Fig.4 Comparison between model in this paper and model in Reference［23］

3.2 典型曲線分析

煤層氣藏多分支水平井?dāng)M壓力曲線和產(chǎn)能曲線（圖5）的基本參數(shù)取值為：S=0.2，CD=5×10-6，LaD=1，dD=8，N=4，ω=0.1，ωad=5，λ=1×10-9，χa=60°，γpD=0。根據(jù)曲線特征，將其劃分為9個(gè)流動階段：①井筒儲集階段。擬壓力和擬壓力導(dǎo)數(shù)曲線重合，產(chǎn)能和產(chǎn)能導(dǎo)數(shù)曲線重合。②表皮效應(yīng)階段。擬壓力導(dǎo)數(shù)曲線呈駝峰狀，產(chǎn)能和產(chǎn)能導(dǎo)數(shù)曲線繼續(xù)下降。③分支井筒第一徑向流階段。擬壓力導(dǎo)數(shù)呈水平線，產(chǎn)能導(dǎo)數(shù)曲線呈下降趨勢。④分支井筒第一線性流階段，擬壓力導(dǎo)數(shù)曲線斜率為0.5，產(chǎn)能導(dǎo)數(shù)曲線呈上升趨勢。⑤分支井筒的第二徑向流階段。擬壓力導(dǎo)數(shù)呈水平線，產(chǎn)能與產(chǎn)能導(dǎo)數(shù)曲線平行。⑥多分支水平井的整體第二線性流階段。擬壓力導(dǎo)數(shù)曲線斜率為0.5，產(chǎn)能導(dǎo)數(shù)呈水平線。⑦割理系統(tǒng)擬徑向流階段。擬壓力導(dǎo)數(shù)呈水平線，產(chǎn)能與產(chǎn)能導(dǎo)數(shù)曲線平行。⑧基質(zhì)系統(tǒng)向割理系統(tǒng)的竄流階段。擬壓力導(dǎo)數(shù)和產(chǎn)能導(dǎo)數(shù)曲線均呈現(xiàn)出一個(gè)寬的凹子。⑨整體徑向流階段，此時(shí)整個(gè)系統(tǒng)的流動達(dá)到平衡狀態(tài)。

圖5 煤層氣藏多分支井?dāng)M壓力曲線和產(chǎn)能曲線Fig.5 Pseudo-pressure curves and productivity curves of multi-branched wells in CBM reservoir

從考慮非穩(wěn)態(tài)擴(kuò)散和擬穩(wěn)態(tài)擴(kuò)散的產(chǎn)能曲線對比（圖6）中可知：前期，兩者的產(chǎn)能曲線重合；而后期，擬穩(wěn)態(tài)擴(kuò)散呈現(xiàn)出的凹子窄且深，而非穩(wěn)態(tài)擴(kuò)散呈現(xiàn)的凹子寬且淺，采用擬穩(wěn)態(tài)擴(kuò)散模型時(shí)，竄流強(qiáng)度較高，所以就計(jì)算產(chǎn)能而言精確度不高。

圖6 非穩(wěn)態(tài)擴(kuò)散和擬穩(wěn)態(tài)擴(kuò)散的對比Fig.6 Comparison between unsteady diffusion and pseudo-steady diffusion

3.3 敏感性參數(shù)分析

無因次滲透率模量 從圖7 中可以看出，無因次滲透率模量主要影響產(chǎn)能曲線最后部分。其原因?yàn)?，在生產(chǎn)初期井底壓力和地層原始壓力相差較小，割理系統(tǒng)滲透率的應(yīng)力敏感性較差。隨著時(shí)間的推移，應(yīng)力敏感性逐漸增加，割理系統(tǒng)滲透率下降，氣體流動阻力增強(qiáng)，從而對產(chǎn)能造成影響。無因次滲透率模量取值越小，后期產(chǎn)能較高；相反，若取值越大，則產(chǎn)能較低。

圖7 無因次滲透率模量對產(chǎn)能曲線的影響Fig.7 Effect of dimensionless permeability modulus on productivity curve

分支井筒間距 從圖8 中可以看出，分支井筒間距造成的影響在產(chǎn)能曲線后部分，隨著分支井筒間距的減小，產(chǎn)能和產(chǎn)能導(dǎo)數(shù)曲線下降。其原因?yàn)?，在生產(chǎn)前期，各個(gè)分支井筒之間的干擾并未形成；隨著生產(chǎn)的進(jìn)行，分支井筒之間的干擾開始形成，同時(shí)分支井筒間距越小，干擾發(fā)生的時(shí)間越早，干擾程度越大，從而較小的分支井筒間距會造成產(chǎn)能下降的時(shí)間提前且下降的速度變快。

圖8 分支井筒間距對產(chǎn)能曲線的影響Fig.8 Effect of branch spacing on productivity curve

分支井筒個(gè)數(shù) 分支井筒個(gè)數(shù)反映多分支水平井的生產(chǎn)規(guī)模，個(gè)數(shù)越多，流入水平井筒的氣體越多，從而生產(chǎn)能力越強(qiáng)。從圖9中可以看出，隨著分支井個(gè)數(shù)的增加，在產(chǎn)能曲線中間段部分，產(chǎn)能越大，且產(chǎn)能降低速度越小，在中期整體表現(xiàn)為產(chǎn)能和產(chǎn)能導(dǎo)數(shù)曲線的上升，但在后期可以看出增加分支井筒個(gè)數(shù)的增產(chǎn)效果不明顯。

圖9 分支井個(gè)數(shù)對產(chǎn)能曲線的影響Fig.9 Effect of number of branches on productivity curve

竄流系數(shù) 從竄流系數(shù)對產(chǎn)能曲線的影響（圖10）可以看出，竄流系數(shù)越小，說明基質(zhì)系統(tǒng)和割理系統(tǒng)之間的物理性質(zhì)差別越大，基質(zhì)系統(tǒng)解吸出來的氣體越難擴(kuò)散竄流進(jìn)入割理系統(tǒng)中被產(chǎn)出。在產(chǎn)能曲線上表現(xiàn)為竄流系數(shù)越小，竄流開始的時(shí)間越晚，凹子右移且下移。

圖10 竄流系數(shù)對產(chǎn)能曲線的影響Fig.10 Effect of crossflow coefficient on productivity curve

吸附系數(shù) 從圖11 中可以看出，吸附系數(shù)越大，竄流階段的凹子就越深且越寬。這是因?yàn)槲较禂?shù)代表了基質(zhì)系統(tǒng)以吸附形式存在的氣體體積。因此，吸附系數(shù)越大，解吸出來的氣體擴(kuò)散竄流至割理系統(tǒng)的量越多，從而產(chǎn)能曲線的位置越高，竄流階段的凹子越明顯。

圖11 吸附系數(shù)對產(chǎn)能曲線的影響Fig.11 Effect of adsorption coefficient on productivity curve

彈性儲容比 從彈性儲容比對產(chǎn)能曲線的影響（圖12）可以看出，彈性儲容比對產(chǎn)能曲線中期部分影響較大。彈性儲容比越大，表明割理系統(tǒng)的儲集能力越強(qiáng)，氣體量越大，在生產(chǎn)壓差不變的情況下，產(chǎn)出的氣體越多。具體表現(xiàn)為竄流階段的凹子越淺越窄，同時(shí)凹子出現(xiàn)時(shí)間推遲，產(chǎn)能曲線中部位置抬升越高。

圖12 彈性儲容比對產(chǎn)能曲線的影響Fig.12 Effect of storage ratio on productivity curve

4 結(jié)論

綜合考慮割理系統(tǒng)的應(yīng)力敏感特性及基質(zhì)系統(tǒng)的非穩(wěn)態(tài)擴(kuò)散建立了煤層氣藏多分支水平井非穩(wěn)態(tài)產(chǎn)能數(shù)學(xué)模型，利用一系列數(shù)學(xué)方法求解得到產(chǎn)能表達(dá)式。通過模型對比驗(yàn)證了所提出模型的可靠性，并且繪制了典型的擬壓力曲線以及產(chǎn)能曲線。

分別考慮擬穩(wěn)態(tài)擴(kuò)散和非穩(wěn)態(tài)擴(kuò)散的產(chǎn)能曲線對比，說明了非穩(wěn)態(tài)擴(kuò)散模型就產(chǎn)能計(jì)算而言精確度更高。無因次滲透率模量的增加會降低產(chǎn)能；分支井筒間距和分支井筒個(gè)數(shù)主要影響產(chǎn)能曲線的中間位置，分支井筒間距越小產(chǎn)能曲線位置越低，分支井筒個(gè)數(shù)的增加會引起產(chǎn)能的增加；竄流系數(shù)和吸附系數(shù)主要影響竄流階段；彈性儲容比除了影響竄流階段之外，其取值越大，產(chǎn)能越大。此外，基于本文模型還可以進(jìn)行考慮煤層氣藏非線性滲流的產(chǎn)能研究。

符號解釋

La——分支井筒的長度，m；d——分支井筒間距，m；χa——分支井筒與水平井筒之間夾角，（°）；h——儲層厚度，m；x,y,z——笛卡爾坐標(biāo)，m；m——擬壓力，Pa2/（Pa·s）；下標(biāo)D——無因次；Kfhi——水平方向割理系統(tǒng)初始滲透率，m2；Tsc——標(biāo)準(zhǔn)狀態(tài)下的溫度，K；psc——標(biāo)準(zhǔn)狀態(tài)下的壓力，Pa；Q——多分支水平井產(chǎn)量，m3/s；T——溫度，K；mi——初始條件下的擬壓力，Pa2/（Pa·s）；γp——擬滲透率模量，mPa·s/Pa2；Vm——基質(zhì)系統(tǒng)吸附氣體體積濃度，m3/m3；VEi——標(biāo)準(zhǔn)狀態(tài)下的氣體平衡體積濃度，m3/m3；Cm——基質(zhì)系統(tǒng)下的氣體體積濃度，m3/m3；Cmi——基質(zhì)系統(tǒng)下的初始?xì)怏w體積濃度，m3/m3；t——時(shí)間，s；Λ——無因次中間變量；L——參考長度，m；q——點(diǎn)源產(chǎn)量，m3/s；λ——竄流系數(shù)；D——擴(kuò)散系數(shù)，m2/s；R——基質(zhì)單元半徑，m；ωad——吸附系數(shù)；mL——Langmuir 擬壓力，Pa2/（Pa·s）；mf——割理系統(tǒng)擬壓力，Pa2/（Pa·s）；ω——彈性儲容比；?f——割理系統(tǒng)下的孔隙度，f；Ctf——割理系統(tǒng)壓縮系數(shù)，Pa-1；μ——黏度，Pa·s；r——徑向距離，m；κ——無窮小垂向距離，m；σ——無窮小徑向距離，m；Kfvi——垂直方向割理系統(tǒng)初始滲透率，m2；rm——基質(zhì)單元徑向距離，m；γ——滲透率模量，Pa-1；pi——儲層初始壓力，Pa；pf——割理系統(tǒng)壓力，Pa；Z——?dú)怏w壓縮因子，小數(shù)；Ct——儲層綜合壓縮系數(shù)，Pa-1；zwD——無因次點(diǎn)源縱向坐標(biāo)；——Laplace 空間下無因次基質(zhì)系統(tǒng)氣體體積濃度；u——Laplace 因子；——Laplace 空間下無因次割理系統(tǒng)擬壓力——Laplace 空間下無因次基質(zhì)系統(tǒng)吸附氣體體積濃度；——Laplace空間下無因次點(diǎn)源產(chǎn)量；τfD——Pedrosa 變換后的無因次割理系統(tǒng)擬壓力；——Laplace 空間下無因次割理系統(tǒng)擬壓力——有限余弦變換后的Laplace 空間下無因次割理系統(tǒng)擬壓力；gn——無因次中間變量；n——取0，1，2，3，…；A——中間系數(shù)；I0——第一類零階虛宗量Bessel函數(shù)；B——中間系數(shù)；K0——第二類零階虛宗量Bessel 函數(shù)；I1——第一類一階虛宗量Bessel 函數(shù)；K1——第二類一階虛宗量Bessel 函數(shù)；M——每個(gè)分支井筒離散微元段數(shù)；qa,b——第a個(gè)分支第b個(gè)微元段的流量，m3/s；——Laplace 空間下第a個(gè)分支第b個(gè)微元段的壓力響應(yīng)；ΓDa,b——第a個(gè)分支第b個(gè)微元段的無因次積分長度；——第α個(gè)分支第β個(gè)微元段無因 次中點(diǎn)坐標(biāo)；N——分支井筒個(gè)數(shù)；——Laplace 空間下第a個(gè)分支第b個(gè)微元段的線流量密度；CD——無因次井筒儲集系數(shù)；mwD——真實(shí)空間下無因次井底壓力；——第a個(gè)分支第b個(gè)微元段的無因次長度；S——表皮系數(shù)；——Laplace空間下的無因次井底壓力。