頁巖油藏裂縫網(wǎng)絡(luò)多相滲流數(shù)值模擬研究

咸玉席，陳超峰，封 猛，郝有志

（1.中國(guó)科學(xué)技術(shù)大學(xué)石油天然氣研究中心，安徽合肥230027；2.中國(guó)石油新疆油田分公司勘探事業(yè)部，新疆克拉瑪依834000；3.中國(guó)石油集團(tuán)西部鉆探工程有限公司試油公司，新疆克拉瑪依834000）

近年來，頁巖油開發(fā)成為熱點(diǎn)[1-2]。頁巖油儲(chǔ)層具有低孔低滲的特點(diǎn)，烴類化合物表現(xiàn)為自生自儲(chǔ)，有機(jī)質(zhì)結(jié)構(gòu)與性質(zhì)顯著影響油氣的流動(dòng)[3]。頁巖致密儲(chǔ)層經(jīng)水平井和壓裂技術(shù)改造，形成尺度較大并與天然裂縫網(wǎng)絡(luò)溝通的裂縫網(wǎng)絡(luò)，可提高其產(chǎn)能。人工形成裂縫的尺度與微納孔隙裂縫相差很大[4]，導(dǎo)致頁巖油在不同尺度裂縫內(nèi)的流動(dòng)機(jī)理不同，存在表征流體單一、裂縫尺度范圍大、劃分網(wǎng)格要求精度高和流體參數(shù)在裂縫界面處不連續(xù)等困難，難以利用簡(jiǎn)單的解析模型表征油在整個(gè)裂縫網(wǎng)絡(luò)中的流動(dòng)，導(dǎo)致無法對(duì)頁巖油藏復(fù)雜裂縫中的多相流體流動(dòng)進(jìn)行數(shù)值模擬表征。目前，廣泛應(yīng)用雙重介質(zhì)模型進(jìn)行試井解釋，但該模型并不適合有大尺度導(dǎo)流裂縫的油藏[5-6]。早期學(xué)者使用單重介質(zhì)模型模擬導(dǎo)流裂縫，但該模型需要極細(xì)密的網(wǎng)格描述裂縫和基質(zhì)，計(jì)算量大、效率低，難以應(yīng)用到實(shí)際的油藏?cái)?shù)值模擬中[7]。后來發(fā)展的離散裂縫模型假定裂縫內(nèi)法向上的物理量保持不變，裂縫節(jié)點(diǎn)與基質(zhì)節(jié)點(diǎn)在網(wǎng)格內(nèi)重合，將裂縫及基質(zhì)方程疊加后計(jì)算，實(shí)現(xiàn)了對(duì)單相流體流動(dòng)的模擬，但不適用于多相流體流動(dòng)。頁巖油藏經(jīng)過水力壓裂后，導(dǎo)流裂縫的幾何形狀較復(fù)雜，需要沿裂縫劃分非規(guī)則網(wǎng)格，但這類離散裂縫模型計(jì)算效率低。李憲文等人[8]利用微地震數(shù)據(jù)，通過反演得到次生裂縫的分布特征，建立了壓裂水平井滲流數(shù)學(xué)模型，定量表征了滲流場(chǎng)的變化規(guī)律。李玉梅等人[9]利用離散元方法，分析了水力裂縫長(zhǎng)度、天然裂縫傾角、內(nèi)摩擦角及施工凈壓力對(duì)縫網(wǎng)的影響，但未考慮基質(zhì)與裂縫間的滲流關(guān)系。同樣，采用數(shù)值模擬方法模擬頁巖油氣藏時(shí)，也是利用傳統(tǒng)試井方法進(jìn)行反演，并未考慮多相滲流在頁巖網(wǎng)絡(luò)裂縫的流動(dòng)[10-12]。

黎水泉等人[13-15]對(duì)離散裂縫模型進(jìn)行了改進(jìn)，提出了基于達(dá)西定律的交換流量耦合裂縫及基質(zhì)方程，利用該方程可以準(zhǔn)確地模擬油藏裂縫中的多相流動(dòng)，但劃分網(wǎng)格時(shí)要求沿著裂縫的走向劃分，計(jì)算量大，求解大型矩陣出現(xiàn)奇異，不具有適用性。S.H.Lee等人[16]把裂縫近似看作基質(zhì)網(wǎng)格中的井源，提出了嵌入裂縫模型，利用Peaceman公式表征裂縫與基質(zhì)網(wǎng)格間的流體交換，該模型網(wǎng)格劃分效率高，但不能準(zhǔn)確表征多相流體在裂縫和基質(zhì)之間的交換。為此，筆者就多相流體在頁巖油藏復(fù)雜裂縫中的流動(dòng)進(jìn)行了理論模型和數(shù)值求解方面的研究，建立了嵌入裂縫多相流體流動(dòng)的控制模型和網(wǎng)格劃分方法，實(shí)現(xiàn)了數(shù)值模擬的高效計(jì)算，為頁巖油藏的產(chǎn)能評(píng)估提供了新的技術(shù)方法。

1 嵌入裂縫多相流動(dòng)模型的建立

1.1 嵌入裂縫多相流動(dòng)模型

1.1.1 控制方程

頁巖油藏中裂縫開度很小，其中流體的流動(dòng)可以被看作一維流動(dòng)。裂縫和基質(zhì)中α相流體的質(zhì)量守恒方程分別為：

式中：wf為裂縫寬度，m；τf為裂縫走向；qfm,α為單位時(shí)間內(nèi)從單位長(zhǎng)度裂縫流入基質(zhì)的α（α=w,o）相流體的流量，m3/（m·s）；vf,α，vm,α分別為裂縫和基質(zhì)中α相流體的流速，m/s；ρα為α相流體的密度，kg/m3；μα為α 相流體的黏度，mPa·s；φf，φm分別為裂縫和基質(zhì)的孔隙度；Kf，Km分別為裂縫和基體的絕對(duì)滲透率，mD；Kf,rα，Km,rα分別為裂縫和基質(zhì)中α相流體的相對(duì)滲透率；pf,α，pm,α分別為裂縫和基質(zhì)中α相流體的壓力，MPa；Sf,α，Sm,α分別為裂縫和基質(zhì)中α相流體的飽和度；δ為克羅尼克函數(shù)；?為拉布拉斯算子；t為時(shí)間，s；r為裂縫走向變量。

頁巖油藏裂縫的開度很小，其中基質(zhì)中的壓力場(chǎng)在跨越裂縫時(shí)保持連續(xù)，多相流體在裂縫-基質(zhì)中存在不同的流動(dòng)特征，導(dǎo)致界面處相飽和度發(fā)生間斷和基質(zhì)的相飽和度場(chǎng)在裂縫處存在不連續(xù)，如圖1所示。

圖1 油藏中垂直于裂縫方向上的油相飽和度分布曲線Fig.1 Saturation distribution curve of oil phase perpendicular to the fracture direction in oil reservoirs

裂縫與兩側(cè)基質(zhì)的物質(zhì)流動(dòng)交換量并不相等，基質(zhì)與裂縫的交換量可以表示為：

式中為裂縫兩側(cè)邊界的外法線方向，其中表示裂縫兩側(cè)。

1.1.2 方程離散

嵌入基質(zhì)的裂縫被規(guī)則的基質(zhì)網(wǎng)格分割為Nf段一維裂縫網(wǎng)格。頁巖油藏裂縫寬度為wf，取第k段裂縫網(wǎng)格為研究對(duì)象，裂縫內(nèi)流體的流動(dòng)被認(rèn)為是一維的，采用一維離散格式把裂縫中α相質(zhì)量守恒方程全隱式離散為：

式中：Δt為時(shí)間步長(zhǎng)，s；n為t=nΔt時(shí)刻的變量；ΔLk為第k段裂縫網(wǎng)格的長(zhǎng)度，m；φf,k為第k段裂縫網(wǎng)格的孔隙度。

式中：各組分流體的相對(duì)滲透率Kf,rα使用上限值。

為實(shí)現(xiàn)裂縫與基質(zhì)之間的流量交換，在裂縫兩側(cè)建立了位于裂縫網(wǎng)格垂直平分線上的參考點(diǎn)M+和M?，其與裂縫的距離分別為d+和d?（見圖2）。由于壓力場(chǎng)在裂縫附近連續(xù)分布，且其導(dǎo)數(shù)在單側(cè)基質(zhì)中也保持連續(xù)，所以可以通過裂縫及基質(zhì)網(wǎng)格 中的壓力單側(cè)插值確定參考點(diǎn)M+和M?上的壓力，裂縫與兩側(cè)基質(zhì)的流量交換可以表示為：

圖2 裂縫網(wǎng)格兩側(cè)的參考點(diǎn)Fig.2 Reference points on both sides of fracture meshes

假設(shè)基質(zhì)網(wǎng)格中的空間步長(zhǎng)為Δx=Δy=h，以空間坐標(biāo)為(i,j)的網(wǎng)格為研究對(duì)象，使用二維九點(diǎn)離散格式，基質(zhì)中α相質(zhì)量守恒方程可以全隱式離散為：

式（12）中，各組分流體的相對(duì)滲透率Km,rα使用上限值，而系數(shù)T滿足：

由于裂縫穿過基質(zhì)網(wǎng)格，導(dǎo)致流體的物理量在網(wǎng)格中發(fā)生間斷，該網(wǎng)格不參與實(shí)際計(jì)算，而是為其他網(wǎng)格的計(jì)算提供邊界條件。裂縫兩側(cè)基質(zhì)中的物理量分別保持連續(xù)，通過單側(cè)插值的方法利用基質(zhì)網(wǎng)格中的壓力pm,i?1,j和pm,i,j+1及 裂縫壓力pf,k確定網(wǎng)格(i,j)中S1部 分的壓力場(chǎng)，得到網(wǎng)格(i,j)進(jìn)入網(wǎng)格(i?1,j)和(i,j+1)的流量（見圖3，圖中的虛線界面可以看作整個(gè)油藏基質(zhì)的內(nèi)邊界，且滿足第二類邊界條件），為網(wǎng)格(i?1,j)和 (i,j+1)的流量計(jì)算提供邊界條件。同理，該類網(wǎng)格與相鄰常規(guī)網(wǎng)格的流量交換都可以使用上述方法求出。

圖3 內(nèi)邊界替代覆蓋裂縫的基質(zhì)網(wǎng)格Fig.3 Replacement of matrix meshes covering fractures by the inner boundary

求解裂縫控制方程時(shí)，若裂縫位于油藏內(nèi)部，則在兩端取絕流邊界條件；若裂縫延伸至油藏邊界，其兩端的邊界條件由油藏的外邊界條件確定。

1.2 多裂縫交叉網(wǎng)絡(luò)多相流動(dòng)模型

頁巖油藏中壓裂裂縫網(wǎng)絡(luò)中存在多裂縫交叉現(xiàn)象（見圖4），在交叉點(diǎn)處多條裂縫相互連通，但由于交叉點(diǎn)體積僅為e2量級(jí)，遠(yuǎn)小于網(wǎng)格尺度，導(dǎo)致無法模擬多相流體在多裂縫交叉網(wǎng)絡(luò)中的流動(dòng)。設(shè)置交叉點(diǎn)體積僅為e2量級(jí)，可以避免在計(jì)算中出現(xiàn)奇異值。

圖4 交叉于一點(diǎn)的3條裂縫Fig.4 Three fractures intersecting at one point

為避免小體積交叉點(diǎn)而導(dǎo)致求解矩陣出現(xiàn)奇異值，利用裂縫網(wǎng)格間的界面?zhèn)鲗?dǎo)率計(jì)算得到交叉處的交換流量表征裂縫網(wǎng)格交叉點(diǎn)的物理量。單相流動(dòng)中，交叉點(diǎn)O處流體滿足質(zhì)量守恒方程：

式中：viO為流體通過i條 裂縫的流速，m/s；l為匯集到交叉點(diǎn)的裂縫網(wǎng)格數(shù)。

由于交叉點(diǎn)O的長(zhǎng)度極小，即LO?Li，通過裂縫i流入交叉點(diǎn)O的流速可以表示為：

將式（14）代入式（15），得到交叉點(diǎn)O處的壓力為：

根據(jù)式（14）和式（15），可得單相流動(dòng)中裂縫網(wǎng)格之間的流速為：

式中：Tij為裂縫網(wǎng)格i與j之間的界面?zhèn)鲗?dǎo)率。

多相流動(dòng)時(shí)，裂縫網(wǎng)格間的流速為：

2 計(jì)算實(shí)例

參考頁巖油藏開發(fā)中的一口實(shí)例井，采用停泵壓降分析軟件反演某頁巖油藏一口多級(jí)壓裂水平井的裂縫參數(shù)。某頁巖油藏區(qū)域?yàn)? 000 m×2 000 m×100 m，水平井的水平段長(zhǎng)1000 m，水力壓裂生成了20條裂縫；井儲(chǔ)系數(shù)為0.2 m3/MPa。地層裂縫區(qū)域初始?jí)毫?0 MPa,滲透率為4.93 mD，孔隙度為8.0%，地層溫度為80℃；油藏日產(chǎn)液量40 m3，生產(chǎn)時(shí)間為365 d。頁巖油藏相對(duì)滲透率曲線如圖5所示。

圖5 頁巖油藏相對(duì)滲透率曲線Fig.5 Relative permeability curve of the shale oil reservoir

地層中存在天然裂縫，但難以確定具體位置和尺寸。為此，假設(shè)水平井附近有單條裂縫、相互交叉2條裂縫和相交于一點(diǎn)的3條裂縫。水平井多級(jí)水力壓裂產(chǎn)生的裂縫與天然裂縫如圖6所示，其中水力壓裂產(chǎn)生的裂縫半長(zhǎng)最長(zhǎng)為110 m，最短為50 m，平均有效裂縫半長(zhǎng)為91.75 m。圖6設(shè)置了與人工裂縫相交及不相交的單條裂縫、與人工裂縫相交的交叉裂縫、與人工裂縫不相交的交于一點(diǎn)的3條裂縫等3種情況，用于表征天然裂縫可能與人工裂縫構(gòu)成的裂縫網(wǎng)絡(luò)形式。

圖6 多級(jí)壓裂水平井裂縫示意Fig.6 Fractures of a multistage fractured horizontal well

為驗(yàn)證建立的模型，計(jì)算多相流體在多裂縫交叉網(wǎng)絡(luò)中流動(dòng)時(shí)，需要從宏觀物理量上反映有天然裂縫和無天然裂縫的區(qū)別，設(shè)以定產(chǎn)液量40 m3/d生產(chǎn)，計(jì)算其井底壓力數(shù)據(jù)變化，驗(yàn)證人為設(shè)置裂縫對(duì)井底壓力的影響。含有單條裂縫、相互交叉2條裂縫和相交于一點(diǎn)的3條裂縫的多級(jí)壓裂頁巖油水平井的網(wǎng)格劃分如圖7所示。人工裂縫和天然裂縫附近網(wǎng)格較密，不含有裂縫的區(qū)域網(wǎng)格稀疏。可以看出，水力壓裂裂縫區(qū)域網(wǎng)格加密，可以實(shí)現(xiàn)水力壓裂裂縫與天然裂縫的相交；多條裂縫相交于一點(diǎn)的裂縫也可以進(jìn)行網(wǎng)格劃分，并實(shí)現(xiàn)與人工裂縫的交叉。

圖7 含有不同類型裂縫的多級(jí)壓裂水平井離散網(wǎng)格Fig.7 Discrete meshes of a multistage fractured horizontal well with different types of fractures

本算例中給出了定產(chǎn)量和初始?jí)毫Γ搸r油藏地層壓力隨著生產(chǎn)時(shí)間增長(zhǎng)而發(fā)生變化，存在生產(chǎn)時(shí)間較短、壓力難以波及到天然裂縫和無法有效反映存在多條裂縫流體的流動(dòng)引起的壓力變化的問題。因此，假設(shè)圖6中的水平段僅有人工裂縫，給定產(chǎn)液40 m3/d和初始?jí)毫?0 MPa, 計(jì)算出初始時(shí)刻的地層壓力（見圖8），可以看出，多級(jí)壓裂水平井裂縫的壓力為基本上相同；壓力波及范圍區(qū)域如圖9所示，可以看出，無論多級(jí)壓裂水平井中的裂縫長(zhǎng)短，縫內(nèi)壓力在40 MPa下的作用下壓力影響區(qū)域邊界大致相同。悶井10 d以后，縫內(nèi)附近的地層壓力下降，其壓力波及范圍快臨近天然裂縫（見圖10）。選此時(shí)壓力35 MPa為計(jì)算壓力，對(duì)不含有天然裂縫和含有天然裂縫的頁巖油藏裂縫網(wǎng)絡(luò)進(jìn)行數(shù)值模擬，比較其壓力變化。

圖8 不含有天然裂縫的多級(jí)水力壓裂地層壓力分布Fig.8 Formation pressure distribution of a multistage hydraulically fractured horizontal well without natural fractures

圖9 不含有天然裂縫的多級(jí)水力壓裂水平井在初始?jí)毫?0 MPa下壓裂改造后的壓力波及區(qū)域Fig.9 Pressure swept zone of a multistage hydraulically fractured horizontal well without natural fractures with an initial pressure of 40 MPa

圖10 不含有天然裂縫的多級(jí)水力壓裂水平井在悶井10 d后的地層壓力分布Fig.10 Formation pressure distribution of a multistage hydraulically fractured horizontal well without natural fracturesafter 10 d of shut-in

不含有天然裂縫和含有天然裂縫的多級(jí)水力壓裂水平井的井底壓力相對(duì)對(duì)數(shù)曲線如圖11何圖12所示（Ⅰ為井儲(chǔ)段，Ⅱ?yàn)楸砥ざ危鬄榫€性流段，Ⅳ為過渡段或線性流段，Ⅴ為線性流段，Ⅵ為邊界流段）。

圖11 不含有天然裂縫的多級(jí)水力壓裂水平井井底壓力相對(duì)對(duì)數(shù)曲線Fig.11 Relative logarithmic curves of the bottomhole pressure of a multistage hydraulically fractured horizontal well without natural fractures

圖12 含有天然裂縫的多級(jí)水力壓裂水平井井底壓力相對(duì)對(duì)數(shù)曲線Fig.12 Relative logarithmic curves of the bottomhole pressure of a multistage hydraulically fractured horizontal well with natural fractures

從圖11和圖12可以看出，不含有天然裂縫和含有天然裂縫的多級(jí)水力壓裂水平井井底壓力在井儲(chǔ)段存在明顯差異，不含有天然裂縫的井底壓力變化較快，而含有天然裂縫的井底壓力變化持續(xù)時(shí)間較長(zhǎng)，說明致密頁巖油藏中存在裂縫溝通的情況，這與圖6中存在天然裂縫與水力壓裂裂縫溝通相一致。表皮段（Ⅱ）、線性流段（Ⅲ）、過渡段或線性流段（Ⅳ）和線性流段（Ⅴ）的壓力變化趨勢(shì)基本一致，其中過渡段或線性流段（Ⅳ）既呈現(xiàn)出過渡段的特征，又具有線性流的特征，在含有天然裂縫地層中的時(shí)間較短，這與天然裂縫與水力壓裂裂縫的溝通有關(guān)。邊界流段（Ⅵ）的壓力變化趨勢(shì)差異較大，不含有天然裂縫的多級(jí)水力壓裂水平井井底壓力相對(duì)對(duì)數(shù)曲線中的壓力導(dǎo)數(shù)斜率較高，而含有天然裂縫的井底壓力導(dǎo)數(shù)斜率較低，這主要是由縫內(nèi)壓力變化波及到天然裂縫區(qū)域，天然裂縫內(nèi)的多相流體流向水力壓裂裂縫。

3 結(jié)論

1）嵌入裂縫多相流動(dòng)模型和多裂縫交叉網(wǎng)絡(luò)多相流動(dòng)模型，可以對(duì)頁巖油藏復(fù)雜裂縫中的多相流體流動(dòng)進(jìn)行數(shù)值模擬表征，并根據(jù)地層壓力變化規(guī)律得到頁巖油藏隨生產(chǎn)時(shí)間的變化情況。

2）利用裂縫與兩側(cè)基質(zhì)的流量交換，實(shí)現(xiàn)復(fù)雜裂縫嵌入的網(wǎng)格劃分，計(jì)算多相流體的流動(dòng)。但油水在致密頁巖油藏中的滲流速度較低，導(dǎo)致天然裂縫附近的壓力變化與人工裂縫附近的壓力變化相比較小。

3）利用井底壓力變化特征可以表征天然裂縫與水力壓裂裂縫的溝通，根據(jù)井儲(chǔ)段、過渡段和邊界流段的差異判斷溝通水力壓裂裂縫的天然裂縫規(guī)模及天然裂縫與水力壓裂裂縫的距離，為評(píng)價(jià)頁巖油藏儲(chǔ)層提供科學(xué)依據(jù)。