李彥澤 段 彬 商 琳，2，3 孫彥春 龔麗榮 辛春彥 邢建鵬

（1.中國(guó)石油冀東油田公司，河北 唐山 063000；2.中國(guó)石油冀東油田公司博士后科研工作站，河北 唐山 063000；3.中國(guó)石油勘探開(kāi)發(fā)研究院博士后流動(dòng)站，北京 100083；4.中國(guó)石油大學(xué)（北京）安全與海洋工程學(xué)院，北京 102200）

0 引 言

油藏?cái)?shù)值模擬技術(shù)自20世紀(jì)50年代誕生，到21世紀(jì)的今天已發(fā)展為一項(xiàng)較為成熟的常用技術(shù)，作為一種定量預(yù)測(cè)技術(shù)，具有精度高、三維可視化、能夠直觀揭示油氣流在儲(chǔ)層內(nèi)的滲流過(guò)程等特點(diǎn)。在理論研究上，學(xué)者將其應(yīng)用于探索多孔介質(zhì)儲(chǔ)層的復(fù)雜滲流問(wèn)題；在油田開(kāi)發(fā)中，數(shù)值模擬技術(shù)也廣泛應(yīng)用于開(kāi)發(fā)方案、開(kāi)發(fā)調(diào)整方案、油氣藏提高采收率方案等，為油田開(kāi)發(fā)技術(shù)政策的制定提供了重要支撐［1-6］。目前，Eclipse、CMG、VIP 等是目前應(yīng)用比較廣泛的商業(yè)數(shù)值模擬軟件，都以達(dá)西滲流模型為核心，在中高滲透油藏中成熟應(yīng)用。流體在特低滲透、致密砂巖等多孔介質(zhì)儲(chǔ)層中的流動(dòng)，由于邊界效益表現(xiàn)為存在最小啟動(dòng)壓力梯度的低速非線性滲流，達(dá)西線性滲流難以表征此類儲(chǔ)層中流體的滲流過(guò)程［7-9］。因此，目前市面上常用的商業(yè)數(shù)值模擬軟件在特低滲透、致密砂巖油藏?cái)?shù)值模擬中應(yīng)用誤差相對(duì)增大。

自20世紀(jì)90年代以來(lái)，眾多學(xué)者開(kāi)展了低滲透、致密砂巖儲(chǔ)層非線性滲流規(guī)律研究，也探索了相應(yīng)的數(shù)值模擬方法［10-16］，主要認(rèn)識(shí)有3 個(gè)方面：一是普遍認(rèn)識(shí)到流體在特低滲透、致密砂巖等多孔介質(zhì)儲(chǔ)層中的流動(dòng)，存在最小啟動(dòng)壓力梯度，且為低速非達(dá)西滲流，達(dá)西滲流模型在一定程度上夸大了儲(chǔ)層的滲流能力；二是部分學(xué)者提出的擬啟動(dòng)壓力梯度模型，不能表征低速非線性滲流階段，表現(xiàn)出儲(chǔ)層滲流阻力增大，造成采用擬啟動(dòng)壓力梯度模型的數(shù)值模擬軟件的預(yù)測(cè)結(jié)果誤差較大；三是部分學(xué)者提出變滲透率的數(shù)值模擬方法，該方法可以準(zhǔn)確地描述特低滲透、致密儲(chǔ)層多孔介質(zhì)的非達(dá)西滲流，但未能給出連續(xù)光滑的滲透率狀態(tài)方程，而且現(xiàn)有的變滲透率數(shù)值模擬軟件沒(méi)有能夠給出壓力梯度場(chǎng)和滲透率場(chǎng)的變化規(guī)律，因此變滲透率數(shù)值模擬方法未得到很好地應(yīng)用。

為解決非線性滲流模擬問(wèn)題，通過(guò)理論模型推導(dǎo)，建立非線性滲流模型，并求解非線性滲流系數(shù)。鑒于tNavigator 中的啟動(dòng)壓力梯度限制為定值，無(wú)法體現(xiàn)非線性滲流規(guī)律，本文研究考慮對(duì)滲透率進(jìn)行時(shí)變，即通過(guò)對(duì)滲透率進(jìn)行修正來(lái)反映啟動(dòng)壓力梯度的變化，從而反映非線性滲流?，F(xiàn)有軟件無(wú)法實(shí)現(xiàn)壓力梯度的網(wǎng)格表征，通過(guò)前期網(wǎng)格篩選，在定壓條件下進(jìn)行模擬。

1 非線性滲流模型構(gòu)建

原油成分極其復(fù)雜，并不是單純的牛頓流體，其存在屈服應(yīng)力，而孔喉的微小使得原油的非牛頓性質(zhì)更強(qiáng)，不可以忽略。以經(jīng)典的毛管模型為例，流體勻速通過(guò)毛管時(shí)力學(xué)平衡關(guān)系式為

式中：Δp——毛管兩端的驅(qū)動(dòng)壓差，MPa；

r——毛管半徑，μm；

τ0——流體屈服應(yīng)力，MPa；

μ——流體黏度，mPa·s；

v——流體滲流速度，cm/s；

L——毛管長(zhǎng)度，μm。

Δpπr2是驅(qū)動(dòng)力，是內(nèi)摩擦力，τ0體現(xiàn)出流體的非牛頓性對(duì)流體滲流的影響計(jì)算公式為

式中 ?p——驅(qū)動(dòng)壓力梯度，MPa/m。

原油中的極性成分在滲流通道固相表面的吸附，導(dǎo)致邊界層的存在是必然的。從前述實(shí)驗(yàn)可知，低滲透油藏孔喉細(xì)微，比表面大，邊界層的影響是不可以忽略的，考慮邊界層后流量的計(jì)算公式為

式中：q——單根毛細(xì)管流量，cm3；

δ——邊界層厚度，μm。

δ體現(xiàn)出固體與液體分子、液體與液體分子之間相互作用對(duì)流體滲流的影響；油藏巖心相當(dāng)于一組半徑均為r的平行毛細(xì)管埋置于固體之中。如果與流動(dòng)方向相垂直的每單位橫截面面積上有N根這樣的毛細(xì)管，則通過(guò)該多孔介質(zhì)塊的流量計(jì)算公式為［17-18］

式中：Q——截面流量，cm3；

K——滲透率，10-3μm2；

A——截面積，cm2。

式（4）相對(duì)于達(dá)西定律更符合實(shí)際情況，但形式太過(guò)復(fù)雜，多位學(xué)者研究了邊界流體對(duì)低滲透、致密儲(chǔ)層滲流規(guī)律的影響，認(rèn)為邊界層厚度隨驅(qū)動(dòng)壓力梯度的增大而減小。龍騰［15］測(cè)定了去離子水通過(guò)微管的流量與驅(qū)替壓力梯度之間的關(guān)系，并計(jì)算出去離子水在微管中邊界層厚度隨驅(qū)替壓力梯度的變化曲線。閆慶來(lái)等［5］指出隨著驅(qū)替壓力梯度的增大，邊界層流體厚度呈指數(shù)減薄，說(shuō)明邊界層厚度是驅(qū)替壓力梯度的函數(shù)。理論與實(shí)驗(yàn)證實(shí)驅(qū)替壓力梯度越大，儲(chǔ)層的邊界層越薄。由前述邊界層實(shí)驗(yàn)有邊界層比例與壓力梯度呈冪函數(shù)關(guān)系，與黏度呈線性關(guān)系，對(duì)于單根毛細(xì)管，有

式中a，b——擬合系數(shù)，常量。

式（5）非線性化較高，對(duì)于后面模型推廣與實(shí)用較差，因此，通過(guò)進(jìn)一步結(jié)合實(shí)驗(yàn)數(shù)據(jù)，對(duì)于壓力梯度采用反比例函數(shù)表示，即

由邊界層厚度定量表征可知，在壓力梯度過(guò)大時(shí)，邊界層比例減少較小，趨于穩(wěn)定，并不完全符合整體函數(shù)關(guān)系，去除掉壓力梯度過(guò)大后的點(diǎn)所得到的2 種表征方式系數(shù)與誤差如表1所示，轉(zhuǎn)換為反比例函數(shù)后誤差為15.17%，相對(duì)于冪函數(shù)關(guān)系12.77%的誤差，誤差增加不大，且可極大減小模型復(fù)雜程度。因此，在后續(xù)推導(dǎo)過(guò)程中使用簡(jiǎn)化過(guò)后的邊界層厚度定量表征公式。

表1 邊界層公式簡(jiǎn)化前后誤差對(duì)比Table 1 Error comparison of boundary layer formula before and after simplification

忽略高階無(wú)窮小項(xiàng)后，有

最終有

2 非線性滲流數(shù)值模擬方法

外部軟件計(jì)算模型擬合參數(shù)，非線性滲流模型里的a1、a2兩個(gè)參數(shù)，對(duì)于不同類型儲(chǔ)層，值均不同，通過(guò)現(xiàn)有的巖心數(shù)據(jù)，利用matlab 進(jìn)行非線性擬合。通過(guò)對(duì)研究區(qū)塊進(jìn)行全面、系統(tǒng)、綜合研究，應(yīng)用油氣勘探、開(kāi)發(fā)過(guò)程中獲取的地震、測(cè)井、測(cè)試、鉆井以及開(kāi)發(fā)動(dòng)態(tài)等資料，建立更加貼近于客觀沉積演化規(guī)律的精細(xì)化三維地質(zhì)模型。對(duì)油藏儲(chǔ)層進(jìn)行精細(xì)描述和預(yù)測(cè)，為非線性滲流模擬提供可靠的地質(zhì)依據(jù)。對(duì)于采用的非線性滲流模擬方法，滲透率場(chǎng)和壓力場(chǎng)有直接影響，防止場(chǎng)數(shù)據(jù)的極值化，將滲透率場(chǎng)與壓力場(chǎng)初始化。對(duì)于井周網(wǎng)格進(jìn)行平面方向的逐個(gè)篩選，并將其賦值，便于直接調(diào)用。將建立的非線性滲流模型通過(guò)滲透率時(shí)變的方式導(dǎo)入，輸出考慮非線性滲流的模擬結(jié)果，并與其他模型比較。導(dǎo)入歷史數(shù)據(jù)，通過(guò)不同算法對(duì)模擬結(jié)果進(jìn)行智能歷史擬合，反饋數(shù)據(jù)，調(diào)整非線性滲流參數(shù)（圖1）。

圖1 致密砂巖油藏非線性滲流模擬流程Fig.1 Workflow of nonlinear flow simulation in tight sandstone reservoir

3 應(yīng)用實(shí)例

3.1 模擬流程

由于現(xiàn)場(chǎng)壓力梯度測(cè)量較為困難，因此根據(jù)實(shí)驗(yàn)室已開(kāi)展的巖心驅(qū)替實(shí)驗(yàn)數(shù)據(jù)進(jìn)行非線性滲流參數(shù)的擬合，根據(jù)儲(chǔ)層分類結(jié)果，選取冀東油田南堡403X1 斷塊典型井組進(jìn)行計(jì)算，并選取滲透率在1×10-3～3×10-3μm2內(nèi)的巖心數(shù)據(jù)進(jìn)行非線性擬合。南堡403X1 斷塊Ed2儲(chǔ)層為一套辮狀河三角洲前緣砂體沉積，儲(chǔ)層平均滲透率為2.8×10-3μm2，屬典型的特低滲油藏和層狀構(gòu)造巖性油藏，采用壓裂注水開(kāi)發(fā)技術(shù)進(jìn)行開(kāi)發(fā)，經(jīng)過(guò)近10 a 滾動(dòng)開(kāi)發(fā)，已經(jīng)進(jìn)入高含水階段，注入水沿人工裂縫竄流嚴(yán)重，剩余油分布復(fù)雜，急需探索剩余油分布規(guī)律，為開(kāi)發(fā)調(diào)整提供依據(jù)。結(jié)合鉆井資料、物性參數(shù)、生產(chǎn)動(dòng)態(tài)等資料對(duì)M 區(qū)塊進(jìn)行屬性建模，平均滲透率為1.29×10-3μm2，平均孔隙度為11.66%，平均深度為3 324.2 m，含油飽和度為46.2%，含水飽和度為53.8%。

3.1.1 初始化場(chǎng)圖

根據(jù)推導(dǎo)的非線性滲流模型，在采取滲透率時(shí)變的過(guò)程中，時(shí)變的滲透率是每一個(gè)時(shí)間步在上一個(gè)時(shí)間步的基礎(chǔ)上進(jìn)行倍乘，不加以控制會(huì)導(dǎo)致滲透率無(wú)限倍乘出現(xiàn)極值，需對(duì)其進(jìn)行初始化，防止極值情況出現(xiàn)。利用滲透率場(chǎng)初始化、壓力場(chǎng)初始化模型（圖2），使每步調(diào)用初始滲透率參與時(shí)變計(jì)算。

圖2 初始化模型Fig.2 Initialization model

3.1.2 井周網(wǎng)格篩選

tNavigator 無(wú)法做到壓力梯度的表征，因此根據(jù)泄油半徑距離，通過(guò)篩選網(wǎng)格的方式進(jìn)行壓力梯度的模擬，通過(guò)Voronoi 圖指定半徑的方式進(jìn)行網(wǎng)格篩選?？蓪⒕芊秶鷥?nèi)的網(wǎng)格進(jìn)行篩選，便于之后對(duì)壓力梯度進(jìn)行表征。以網(wǎng)格長(zhǎng)度為最小單位步長(zhǎng)進(jìn)行篩選（圖3）。

圖3 井周網(wǎng)格篩選Fig.3 Grid screening around the well

將井附近網(wǎng)格篩選過(guò)后，對(duì)于井周圍不同距離的網(wǎng)格可近似表示壓力梯度。

3.1.3 滲透率時(shí)變

tNavigator 中的啟動(dòng)壓力梯度限制為定值，無(wú)法體現(xiàn)非線性滲流規(guī)律，采用對(duì)滲透率進(jìn)行時(shí)變的方法，即通過(guò)對(duì)滲透率進(jìn)行修正來(lái)反映啟動(dòng)壓力梯度的變化，從而反映非線性滲流。對(duì)于井周不同距離的網(wǎng)格進(jìn)行非線性滲流模型的導(dǎo)入，改變滲透率，在井底流壓大于0 MPa 的情況下，將推導(dǎo)的非線性滲流模型導(dǎo)入，結(jié)合實(shí)際生產(chǎn)情況，對(duì)井周不同范圍內(nèi)的網(wǎng)格進(jìn)行滲透率時(shí)變。

3.2 結(jié)果分析

歷史擬合情況及對(duì)比，對(duì)導(dǎo)入考慮非線性滲流影響的模型作智能歷史擬合，采用拉丁超立方實(shí)驗(yàn)算法對(duì)日產(chǎn)油量、累計(jì)產(chǎn)油量進(jìn)行擬合（圖4、圖5），選取擬合程度較高且參數(shù)均科學(xué)合理的擬合模型，并與現(xiàn)有的擬啟動(dòng)壓力梯度模型進(jìn)行比較。

圖4 不同模型日產(chǎn)油量對(duì)比Fig.4 Comparison of daily oil production of different models

圖5 不同模型累計(jì)產(chǎn)油量對(duì)比Fig.5 Comparison of cumulative oil production of different models

模擬結(jié)果表明，考慮非線性滲流的模型，模擬日產(chǎn)油量與實(shí)際產(chǎn)油量相關(guān)性、趨勢(shì)性較好。累計(jì)產(chǎn)油量由非線性模型擬合結(jié)果為12.085 9×104m3，與實(shí)際產(chǎn)油量11.424 8×104m3相比，誤差為5.47%，擬啟動(dòng)模型產(chǎn)油量為9.826 8×104m3，誤差為13.99%。由此不難發(fā)現(xiàn)，非線性模型擬合日產(chǎn)油量和累計(jì)產(chǎn)油量與實(shí)際生產(chǎn)情況更接近，符合率更高，由此表明該方法適用于特低滲、致密砂巖油藏?cái)?shù)值模擬。

4 結(jié) 論

（1）基于非線性滲流理論及啟動(dòng)壓力實(shí)驗(yàn)、邊界層實(shí)驗(yàn)，推導(dǎo)了不同類型儲(chǔ)層的非線性滲流新模型，綜合考慮了啟動(dòng)壓力梯度、邊界層厚度對(duì)非線性滲流的影響，新模型考慮了流體黏度的影響。

（2）基于非線性滲流模型，開(kāi)展商業(yè)軟件二次開(kāi)發(fā)，通過(guò)篩選井周網(wǎng)格及對(duì)滲透率場(chǎng)與壓力場(chǎng)初始化，結(jié)合實(shí)驗(yàn)數(shù)據(jù)，對(duì)滲透率進(jìn)行時(shí)變反映非線性滲流對(duì)儲(chǔ)層的影響，形成了非線性滲流數(shù)值模擬方法。

（3）利用本文方法對(duì)典型井組歷史擬合，通過(guò)對(duì)比來(lái)看擬啟動(dòng)壓力梯度模擬方法夸大了地層阻力，數(shù)值模擬結(jié)果低于實(shí)際產(chǎn)量，非線性滲流模型更符合實(shí)際生產(chǎn)情況，誤差更小，由此證明該方法適用于特低滲致密砂巖油藏?cái)?shù)值模擬。