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

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

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

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

1 嵌入裂縫多相流動模型的建立

1.1 嵌入裂縫多相流動模型

1.1.1 控制方程

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

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

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

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

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

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

1.1.2 方程離散

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 計算實例

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

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

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

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

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

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

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

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

圖9 不含有天然裂縫的多級水力壓裂水平井在初始壓力40 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 不含有天然裂縫的多級水力壓裂水平井在悶井10 d后的地層壓力分布Fig.10 Formation pressure distribution of a multistage hydraulically fractured horizontal well without natural fracturesafter 10 d of shut-in

不含有天然裂縫和含有天然裂縫的多級水力壓裂水平井的井底壓力相對對數(shù)曲線如圖11何圖12所示（Ⅰ為井儲段，Ⅱ為表皮段，Ⅲ為線性流段，Ⅳ為過渡段或線性流段，Ⅴ為線性流段，Ⅵ為邊界流段）。

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

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

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

3 結(jié)論

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

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

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