各向異性油藏水平井多角度人工裂縫線性單元計算方法

方思冬，程林松，辛一男，何聰鴿

（1.中國石油大學(xué)（北京）石油工程學(xué)院，北京 102249；2.中國石油勘探開發(fā)研究院，北京 100083）

文章編號：1001?246X（2015）05?0595?08

各向異性油藏水平井多角度人工裂縫線性單元計算方法

方思冬1，程林松1，辛一男1，何聰鴿2

（1.中國石油大學(xué)（北京）石油工程學(xué)院，北京 102249；2.中國石油勘探開發(fā)研究院，北京 100083）

考慮各向異性油藏滲透率張量的表征，利用Green函數(shù)和拉普拉斯變換建立裂縫流動的一維單元，裂縫的流量分布采用節(jié)點(diǎn)線性插值，裂縫內(nèi)的流動處理為線性積分，耦合地層與人工裂縫的流動，建立有限導(dǎo)流裂縫井底壓力的求解方法.結(jié)果表明：多裂縫壓裂水平井存在壓裂裂縫線性流、地層線性流、系統(tǒng)徑向流3種流動形態(tài)，壓裂裂縫條數(shù)越多，相同的生產(chǎn)時間，無因次井底壓降越??；裂縫條數(shù)對流動影響明顯.隨著裂縫條數(shù)的增加，壓降變化減?。涣芽p長度和導(dǎo)流能力有相似的變化.人工裂縫與井筒角度越大，產(chǎn)能越大，當(dāng)裂縫垂直于井筒時，產(chǎn)量最大；地層最大滲透率方向垂直于人工裂縫時產(chǎn)量最大，平行裂縫時產(chǎn)量最小.當(dāng)人工裂縫垂直于井筒，并同時垂直于地層最大滲透率方向時，達(dá)到最大產(chǎn)量值.

各向異性油藏；壓裂水平井；試井；點(diǎn)源函數(shù)；多角度裂縫

0 引言

目前，非常規(guī)油氣藏，如致密油氣藏以及頁巖氣藏，普遍采用壓裂水平井進(jìn)行開發(fā)［1－2］.國外一部分學(xué)者主要利用Gringarten［3－4］提出的源函數(shù)方法求解壓力動態(tài)，如Evans和Guo［5］利用實(shí)空間的點(diǎn)源函數(shù)和Newman乘積原理得到了有多條裂縫的水平井壓力動態(tài).Rbeawi［6］利用Green函數(shù)求解均質(zhì)油藏中具有多條傾斜裂縫壓裂水平井的壓力分布.Valko［7－8］和Zhu［9］建立了體積源函數(shù)，并通過該方法得到了多條裂縫壓裂水平井的半解析解.其他一些學(xué)者主要采用拉普拉斯變換和數(shù)值反演求解不同情況下的壓力動態(tài)，如Ozkan［10－12］建立了拉氏空間內(nèi)考慮不同的井型、邊界類型以及不同油藏類型的試井模型，Raghavan［13－14］利用盒狀地層拉氏空間點(diǎn)源解，采用疊加原理得到了多條裂縫壓裂水平井拉氏空間解.Brown［15］提出利用三線性流模型研究非常規(guī)油氣藏的壓裂水平井動態(tài)，該模型計算簡便，物理過程清晰，但無法準(zhǔn)確表征壓裂水平井各個流動階段.而對于各向異性油藏的研究，劉月田［16］利用坐標(biāo)變換，將各向異性油藏轉(zhuǎn)化為各向同性油藏進(jìn)行求解.何應(yīng)付［17］將各向異性坐標(biāo)變換應(yīng)用到邊界元方法中求解封閉地層各向異性油藏壓力動態(tài).姚軍［18］基于滲透率張量理論與有限元和有限差分方法進(jìn)行了各向異性油藏的數(shù)值模擬研究.廉培慶［19］利用半解析解求取各向異性油藏壓裂水平井產(chǎn)能.以前研究多裂縫壓裂水平井主要是研究與井筒正交裂縫的情況，而與井筒存在不同角度的裂縫研究很少，各向異性油藏的研究主要簡化成最大和最小滲透率兩個方向，而缺乏對實(shí)際油藏各向異性情況的準(zhǔn)確描述.本文吸收有限元線性單元的思想，結(jié)合源函數(shù)方法和滲透率張量對各向異性的表征，建立各向異性油藏壓裂水平井壓力動態(tài)的求解方法，劃分流動階段，分析各向異性和人工裂縫參數(shù)對壓力動態(tài)的影響.

1 裂縫流動模型的建立和求解

1.1 物理模型

水平井穿過多條與水平井筒存在一定夾角的裂縫.整個油藏為水平等厚、均質(zhì)、各向異性油藏，且最大滲透率方向與笛卡爾坐標(biāo)軸方向存在一定夾角；油藏頂?shù)追忾]，只考慮無限大外邊界；水平井僅在裂縫處射孔生產(chǎn)；水平井定產(chǎn)生產(chǎn)，不考慮井筒摩阻引起的附加壓力降，各條裂縫在井筒處壓力相等；水平井井筒考慮為無限導(dǎo)流而人工裂縫為有限導(dǎo)流；地層和流體微可壓縮，流體為單相，且地層中流動符合達(dá)西定理；不考慮重力對流動的影響.

建立平面二維滲流的各向異性油藏滲透率張量模型，坐標(biāo)軸1，2為滲透率張量的主軸，且軸1為最大滲透率方向，軸2為最小滲透率方向，并假定所取的坐標(biāo)系（x，y）與（1，2）存在θ角度［20］：

式中，K1為最大滲透率，μm2；K2為最小滲透率，μm2；θ為坐標(biāo)系與滲透率最大值夾角，rad；Kxy為二階滲透率張量，μm2.

在平面滲透率張量的基礎(chǔ)上，建立考慮滲透率張量的不穩(wěn)定流滲流模型

式中，kxx，kyx，kyy為滲透率張量的分量，μm2；?為孔隙度；Ct為壓縮系數(shù)，1·atm－1；μ為粘度，mPa·s.引入特征變換將（2）式化簡為

式中，η1＝（kyy－/kxx）/（?μCt）．

1.2 壓裂水平井壓力動態(tài)求解模型

1）油藏向人工裂縫流動模型

Gringaden［2］給出了實(shí)空間下壓裂直井的線源

式中，q為壓裂直井產(chǎn)量（恒定），cm3·s－1；B為體積系數(shù)，xw為井半徑，cm；xf為裂縫半長，cm；h為油藏厚度，cm.

根據(jù)已有壓裂直井裂縫上流量分布規(guī)律［4］，將裂縫劃分成若干單元如圖1所示，每個單元內(nèi)流量分布根據(jù)單元兩個端點(diǎn)線性插值構(gòu)成如圖2所示，基于式（5）得出油藏向裂縫單元流動壓降：

圖1 多級壓裂人工裂縫模型Fig.1 Fractured horizontal wellmodel

圖2 裂縫單元模型Fig.2 Fracture elementmodel

式中，ε＝l/lf；qe1，qe2為端點(diǎn)處流量強(qiáng)度，cm3·（s·cm－1）－1；lf為單元體長度，cm.α為裂縫與x軸夾角，rad；Δp（x，y，t）為油藏中任意一點(diǎn)壓降，atm.將式（6）無因次化為

2）裂縫內(nèi)流動模型

目前有限導(dǎo)流裂縫內(nèi)流動簡化為穩(wěn)定線性流動［8－9，19］，根據(jù)達(dá)西定理由裂縫遠(yuǎn)端向射孔點(diǎn)積分，本文采用流量積分法，但流量在每個單元內(nèi)由端點(diǎn)線性插值得到，即裂縫中任意的壓力表示為

將式（10）無因次化并轉(zhuǎn)換到拉普拉斯空間下

式中，CfD＝2πkxf/（kfwf）；lfD＝lf/xf．

以單條裂縫劃分為4個單元為例，根據(jù)裂縫網(wǎng)格的對稱性，建立裂縫流動附加矩陣

3）油藏滲流與裂縫內(nèi)流動的耦合

數(shù)學(xué)模型中，假設(shè)壓裂的裂縫條數(shù)為n，每條裂縫離散單元個數(shù)為m（i），i＝1∶n，每一條壓裂裂縫都假定擁有不同的特性，并且可以在井筒的任意位置.認(rèn)為壓裂裂縫都在相同的井筒壓力下生產(chǎn)，即為水平井無限導(dǎo)流能力模型.水平井總的流量是每一條壓裂裂縫中流出流量的總和.

拉普拉斯空間下，流量和壓力滿足以下表達(dá)式

式中，nt為壓裂水平井全部一維單元數(shù).在無限導(dǎo)流能力的假設(shè)條件下，水平井井筒中各處保持相等的壓力，

將每個單元對應(yīng)式（13），（14），（15）相應(yīng)的公式形式寫成如下矩陣方程

裂縫流動矩陣（12）和油藏流動矩陣（16）耦合求解，將（16）矩陣中壓力部分系數(shù)矩陣表示為：AqD＝pD；式（12）裂縫的流動表示為pwfD－pD＝CqD；以上兩式可以合并為（A＋C）qD＝PwfD.通過求解耦合矩陣方程組，在Laplace空間計算得到水平井井底壓力之后，再利用Stehfest數(shù)值反演方法［20］可以求得實(shí)空間下壓裂水平井有限導(dǎo)流裂縫井底壓力解及每條裂縫的產(chǎn)量.

2 模型對比及壓力動態(tài)分析

2.1 模型對比

將計算結(jié)果與經(jīng)典有限導(dǎo)流裂縫［21］的計算結(jié)果進(jìn)行對比如圖3所示.選取低無因次導(dǎo)流能力（無因次導(dǎo)流系數(shù)為20）有限導(dǎo)流裂縫為例，該條件下裂縫內(nèi)流動阻力大，流動形態(tài)復(fù)雜，與以往計算結(jié)果對比更有意義，對比不同離散網(wǎng)格數(shù)目下計算結(jié)果發(fā)現(xiàn)需要一定數(shù)量的網(wǎng)格才能滿足求解精度，當(dāng)劃分線性單元數(shù)為14時，計算結(jié)果與Blasingame有限導(dǎo)流壓力和壓力導(dǎo)數(shù)曲線幾乎重合.該模型優(yōu)點(diǎn)在于計算精度是可控的，受網(wǎng)格數(shù)影響.

2.2 壓力動態(tài)分析

圖3 有限導(dǎo)流裂縫壓力對比Fig.3 Finite fracture flow validation

計算前提：地層參數(shù)無因次條件，油藏為各向異性，水平滲透率為垂直滲透率2倍，水平井存在3條與井筒呈不同角度的壓裂裂縫，裂縫兩翼不等長.無因次參考長度選取最長裂縫半長，按照式（7）無因次化，各參數(shù)取值依次為：裂縫1半長100m；夾角60°；裂縫2半長120m；夾角90°；裂縫1半長80 m；夾角70°.分別考慮無限導(dǎo)流裂縫、有限導(dǎo)流裂縫（導(dǎo)流能力較大）、有限導(dǎo)流裂縫（導(dǎo)流能力較?。┤鐖D4、圖5所示.

無限導(dǎo)流和高導(dǎo)流能力裂縫，流動形態(tài)分為3種：①地層線性流動階段，是指在流動初期地層中流體線性的流向各條壓裂裂縫.在雙對數(shù)圖上表現(xiàn)為無因次壓力導(dǎo)數(shù)曲線為1/2斜率的直線段.②過渡流階段，該階段沒有明顯特征.③擬徑向流動階段，是指對于整個油藏，如果生產(chǎn)時間很長，且壓力波未傳播到邊界，則流體以擬徑向流的形式向水平井及壓裂裂縫區(qū)域流動.該流動段在雙對數(shù)診斷圖上表現(xiàn)為無因次壓力導(dǎo)數(shù)曲線為0.5值水平直線段.對于低導(dǎo)流能力裂縫，流動形態(tài)分為5種：①地層裂縫雙線性流動階段，是指在流動初期地層中流體線性的流向各條壓裂裂縫，壓裂裂縫內(nèi)部同時存在線性流動，在雙對數(shù)圖上表現(xiàn)為無因次壓力導(dǎo)數(shù)曲線為1/4斜率的直線段.②和④為過渡流階段，而③為地層線性流動階段，⑤是擬徑向流動階段.

圖4 無限導(dǎo)流與高導(dǎo)流裂縫壓力Fig.4 Pressure of infinite and high conductive fracture

圖5 低導(dǎo)流裂縫壓力Fig.5 Pressure of low conductive fracture

3 實(shí)例計算與參數(shù)敏感性分析

3.1 實(shí)例計算

選取大慶長垣外圍致密油區(qū)塊，部分試驗(yàn)區(qū)采用壓裂水平井開采，選取典型致密油藏試驗(yàn)區(qū)域進(jìn)行實(shí)例計算，壓裂水平井微地震監(jiān)測如圖6所示，油層厚度為3.2m，原始地層壓力18MPa，井底流壓9MPa，基質(zhì)滲透率為0.23×10－3μm2，原油密度0.826g·cm－3，原油體積系數(shù)為1.052，原油粘度為1.45mPa·s，有效水平井長度為1 400m，壓裂16段，每段壓2－3簇縫，人工裂縫平均半長為154m，井筒半徑0.05m.計算其無因次壓降曲線與實(shí)際動態(tài)數(shù)據(jù)對比如圖7所示，計算結(jié)果能夠擬合實(shí)際數(shù)據(jù)點(diǎn).

圖6 壓裂水平井微地震監(jiān)測圖Fig.6 Micro?seismicmap of fractured horizontalwell

圖7 壓力擬合Fig.7 Pressure fitting

3.2 儲層各向異性對壓力動態(tài)的影響

儲層各向異性對壓力動態(tài)的影響在以往的研究中只是反映在滲透率主值上［16－17］，無因次壓降分布如圖8所示，對于實(shí)際地層，最大滲透率方向與井筒呈一定角度，而角度的大小對壓力動態(tài)和產(chǎn)能影響較大，如圖9所示，滲透率主值方向與人工裂縫垂直時無因次壓降最小，定井底流壓產(chǎn)量最大，隨著滲透率主值方向與人工裂縫角度的減小，產(chǎn)量減小.因?yàn)?，裂縫半長很大，垂直裂縫面的滲流區(qū)域很大，而平行裂縫方向流入裂縫的流量所占比例較小.

3.3 裂縫條數(shù)對壓力動態(tài)的影響

隨著壓裂工藝技術(shù)的不斷完善和提高，多裂縫壓裂水平井得到了廣泛的應(yīng)用，特別是針對非常規(guī)油藏的開發(fā)與增產(chǎn).如圖10所示：人工裂縫條數(shù)對壓力動態(tài)影響很大，隨著裂縫條數(shù)的增加，同樣的生產(chǎn)時間時，無因次壓降變小.這是由于增加裂縫的條數(shù)可以使流體更多的流入井筒中，從而減少了滲流阻力.但是隨著裂縫條數(shù)的繼續(xù)增加及生產(chǎn)時間的增加，這種變化會越來越不明顯.雖然裂縫條數(shù)增加能提高水平井產(chǎn)能，但隨著裂縫的增加，產(chǎn)能的增加幅度會減小.

圖8 滲透率主值方向與人工裂縫存在夾角Fig.8 Pressure distribution of different angles between permeabilitymain principal value and fractures

圖9 滲透率主值與裂縫的角度壓力Fig.9 Pressure of different angles between permeability main principal value and fractures

3.4 裂縫半長對壓力動態(tài)的影響

裂縫半長是影響壓裂水平井生產(chǎn)動態(tài)的又一個重要因素.裂縫半長對產(chǎn)能有直接的影響，通過對比不同裂縫半長對應(yīng)井底壓力曲線如圖11所示，可以看出，隨著裂縫長度的增加，無因次壓降變小，但這種變小并不是無限制的，隨著裂縫半長的增加，無因次壓降減小的幅度變小.

圖10 不同裂縫條數(shù)井底壓力Fig.10 Pressure distribution of different number of fractures

圖11 不同裂縫半長井底壓力Fig.11 Pressure distribution of different fractures half length

3.5 裂縫導(dǎo)流能力對壓力動態(tài)的影響

裂縫導(dǎo)流能力的大小直接影響到壓裂水平井最終的產(chǎn)能，通過對比不同裂縫導(dǎo)流能力對應(yīng)井底壓力曲線的影響如圖12所示，可以看出，裂縫導(dǎo)流能力越大，無因次壓降越小，但隨著裂縫導(dǎo)流能力的增加，無因次壓降減小的程度變小，特別是當(dāng)無因次導(dǎo)流能力達(dá)到上百數(shù)量級后，產(chǎn)能增加幅度很小.

3.6 裂縫角度對壓力動態(tài)的影響

目前對裂縫與井筒角度對壓力和產(chǎn)能的影響研究較少，通過分析不同裂縫與井筒夾角井底壓力如圖13所示，裂縫與井筒夾角越大，無因次壓降越小，對應(yīng)產(chǎn)能越大，當(dāng)裂縫與井筒垂直時，產(chǎn)能達(dá)到最大值.通過對比不同裂縫與井筒夾角井底壓力圖還可發(fā)現(xiàn)，角度對產(chǎn)能的影響不如前面幾個因素明顯，對產(chǎn)能的影響程度小于裂縫數(shù)目，半長等.

4 結(jié)論

1）結(jié)合各向異性油藏滲透率張量模型以及偏微分方程變換將各向異性油藏轉(zhuǎn)化為各向同性油藏，形成了相應(yīng)滲透率和空間位置表征，為建立滲流方程提供基礎(chǔ).

圖12 不同裂縫導(dǎo)流能力井底壓力Fig.12 Pressure distribution of different fracture conductivity

圖13 不同裂縫與井筒夾角井底壓力Fig.13 Pressure distribution between different angles of fractures and wellbore

2）基于格林函數(shù)，結(jié)合有限元一維單元插值函數(shù)，建立裂縫流動單元的數(shù)學(xué)表征，推導(dǎo)裂縫單元之間的流動關(guān)系，耦合地層向裂縫流動和裂縫內(nèi)流動方程，得到任意角度多裂縫壓裂水平井定產(chǎn)條件下的井底壓力.

3）計算結(jié)果與已有解析結(jié)果對比，壓力曲線以及壓力導(dǎo)數(shù)曲線幾乎相同，驗(yàn)證模型的可靠性與準(zhǔn)確性.

4）分析多裂縫壓裂水平井壓力特征，存在壓裂裂縫線性流、地層線性流、系統(tǒng)徑向流3種流動形態(tài)，通過參數(shù)敏感性分析發(fā)現(xiàn)，壓裂裂縫條數(shù)越多，相同的生產(chǎn)時間，無因次井底壓降越小；裂縫條數(shù)對流動影響明顯，但隨著壓裂條數(shù)的增加，壓降變化幅度減?。涣芽p長度和導(dǎo)流能力也有相似的變化.

5）對比地層最大滲透率方向與人工裂縫與井筒呈不同角度的情況，得出：人工裂縫與井筒角度越大，產(chǎn)能越大，當(dāng)裂縫垂直于井筒時，產(chǎn)量最大，但角度影響程度不及裂縫長度、導(dǎo)流能力，地層最大滲透率方向垂直于人工裂縫時，產(chǎn)量最大，因此，當(dāng)人工裂縫垂直于井筒，并同時垂直于地層最大滲透率方向時，達(dá)到最大產(chǎn)量.

［1］ Yao Jun，Sun Hai，F(xiàn)an Dongmei，et al.Transportmechanisms and numerical simulation of shale gas reservoirs［J］.Journal of China University of Petroleum，2013，37（1）：91－98.

［2］ Huang Yanzhang.Flow mechanism of low permeability reservoirs［M］.Beijing：Petroleum Industry Press，1998：58－79.

［3］ Gringarten A C，Ramey H J.The use of source and green functions in solving unsteady?flow problems in reservoirs［R］.SPE 3818，1973.

［4］ Gringarten A C，Ramey H J.Unsteady?state pressure distributions created by a well with a single horizontal fracture，partial penetration，or restricted entry［R］.SPE 3819，1974.

［5］ Guo G，Evans R D.Pressure?transient behavior and inflow performance of horizontalwells intersecting discrete fractures［R］. SPE 26446，1993.

［6］ Al Rbeawi，Djebbar T.Transient pressure analysis of a horizontal well with multiple inclined hydraulic fractures using type?curvematching［R］.SPE 149902，2012.

［7］ Valko P，Amini S.The method of distributed volumetric sources for calculating the transient and pseudosteady?state productivity of complex well?fracture configurations［R］.SPE 106279，2007.

［8］ Amini S，ValkóP.Using distributed volumetric sources to predict production from multiple?fractured horizontal wells under non?Darcy?flow conditions［R］.SPE 120110，2010.

［9］ Zhu D，Magalhaes F V.Predicting the productivity ofmultiple?fractured horizontal gas wells［R］.SPE 106280，2007.

［10］ Ozkan E，Raghavan R.Some new solutions to solve problems in well testanalysis：Part1?Analytical considerations［R］.SPE，1998.

［11］ Ozkan E，Raghavan R.Supplement to new solutions for well?test?analysis problems：Part 1?Analytical considerations［R］. SPE，1991.

［12］ Ozkan E，Raghavan R.New solutions for well?test?analysis problems：Part 1?Analytical considerations（includes associated papers 28666 and 29213）［R］.SPE 18615，1991.

［13］ Raghavan R S，Chen C，Agarwal B.An analysis of horizontalwells intercepted bymultiple fractures［R］.SPE 27652，1997. ［14］ Brown M，Ozkan E，Raghavan R.Practical solutions for pressure?transient responses of fractured horizontal wells in unconventional shale reservoirs［R］.SPE 125043，2010.

［15］ Ozkan E，Brown M L，Raghavan R.Comparison of fractured?horizontal?well performance in tight sand and shale reservoirs ［R］.SPE 121290，2010.

［16］ Liu Yuetian，Zhang Jichang.Stable permeating flow and productivity analysis for anisotropic reservoirs in horizontal wellnetworks［J］.Petroleum Exploration and Development，2004，31（1）：94－96.

［17］ He Yinfu，Yin Hongjun，Yang Zhengming，et al.Boundary elementmethod for unstable flow in uniform anisotropic reservoir ［J］.Daqing Petroleum Exploration and Development，2006，25（5）：37－40.

［18］ Yao Jun，Li Yaqin，Huang Chaoqin.Finite flement simulation of anisotropic reservoir considering full permeability tensor［J］. Journal ofWuhan Polytechnic University，2009，28（3）：1－6.

［19］ Lian Peiqing，Cheng Linsong，Cao Renyi，et al.A coupling model of low permeability reservoir and fractured horizontal wellbore in nonsteady state［J］.Chinese Journal of Computational Physics，2010，27（2）：203－209.

［20］ Tong Dengke，Chen Yinlei.A note on Stehfestmethond on Laplace numerical inversion［J］.Acta Petrolei Sinica，2001，22 （6）：91－92.

［21］ Blasingame T A，Poe B D.Semi?analytic solutions for a well with a single finite?conductivity vertical fracture［R］.SPE 26424，1993.

Linear Element M ethod for M ulti?angle Fractured Horizontal W ell in Anisotropic Reservoir

FANG Sidong1，CHENG Linsong1，XIN Yinan1，HE Congge2
（1.College of Petroleum Engineering，China University of Petroleum（Beijing），Beijing 102249，China；
2.Research Institute ofPetroleum Exploration and Development，PetroChina，Beijing 100083，China）

With Green function and Laplace transformation，one?dimensional element is established considering effectof permeability tensor of anisotropic reservoir.Inflow of fracture is obtained by linear interpolation of endpoints and flow in fracture is treated with linear integral of flow rate.Coupling flow in formation and fractures，calculating method for bottom?hole pressure is formed semi?analytically.It shows that there are three flow regimes including fracturing linear flow，formation linear flow and system radial flow. Themore the fractures，the less the dimensionless pressure and number of factures has significant impacton flow rate.With increase of fractures increasing rate is dropping in the same time.Fracture length and conductivity have similar characteristics.Flow rate is improving as angle between fracture and wellbore is increasing.Flow rate reachesmaximum as fracture is perpendicular to wellbore，and vice versa.Angle between maximum permeability and fracture has similar impact on production.In summary，production rate reach maximum as fracture is perpendicular to wellbore and direction ofmaximum permeability.

anisotropic reservoir；fractured horizontalwell；well test；source function；multi?angle fractures

TE312

A

2014－09－13；

2014－12－11

國家自然科學(xué)基金（51174215/E0403）資助項(xiàng)目

方思冬（1988－），男，在讀博士，主要從事低滲透油藏工程和數(shù)值模擬研究，E?mail：jbwolfgang＠163.com

Received date： 2014－09－13；Revised date： 2014－12－11