流線數(shù)值模擬中的無(wú)效網(wǎng)格處理方法

王 宇 李治平

1.陜西延長(zhǎng)石油（集團(tuán)）有限責(zé)任公司研究院，陜西 西安 710075；2.中國(guó)地質(zhì)大學(xué)（北京），北京 100083

0 前言

流線模擬技術(shù)能將三維模擬模型轉(zhuǎn)化為一系列的一維流線模型，實(shí)踐表明流線方法是一種快速有效的數(shù)值模擬計(jì)算方法，流線模擬技術(shù)能更好地認(rèn)識(shí)地下流體的分布、運(yùn)移，改善油田開(kāi)發(fā)效果和提高油田最終采收率［1］。在實(shí)際油田開(kāi)采過(guò)程中，油田內(nèi)往往存在無(wú)油非連接區(qū)域，有些局部區(qū)域?yàn)榫植繑鄬佑绊懙乃畡?dòng)力滯留區(qū)或?yàn)槟鄮r，此時(shí)在油藏?cái)?shù)值模擬過(guò)程中對(duì)應(yīng)為不可流動(dòng)區(qū)域，網(wǎng)格化后即為無(wú)效網(wǎng)格。無(wú)效網(wǎng)格的存在會(huì)影響流線數(shù)值模擬的壓力場(chǎng)計(jì)算以及流線追蹤過(guò)程，導(dǎo)致流線數(shù)值模擬失敗。

另外，油藏形狀各式各樣，進(jìn)行網(wǎng)格劃分后往往是非規(guī)則形狀，形成統(tǒng)一的差分過(guò)程，對(duì)角矩陣過(guò)程復(fù)雜，而將油藏劃分為規(guī)則的矩形，對(duì)其內(nèi)部流體不可流動(dòng)區(qū)域按無(wú)效網(wǎng)格處理，形成流程化計(jì)算，計(jì)算過(guò)程更為高效。因此，無(wú)效網(wǎng)格研究對(duì)于流線數(shù)值模擬計(jì)算具有重要意義。

1 壓力數(shù)學(xué)模型與無(wú)效網(wǎng)格的處理

1.1 壓力方程

在二維兩相的前提下，考慮流體與巖石的壓縮性，忽略重力及毛管壓力，考慮維數(shù)因子和產(chǎn)量項(xiàng)，根據(jù)達(dá)西定律、質(zhì)量守恒定律以及流體狀態(tài)方程可得到二維兩相流數(shù)學(xué)模型［2-3］：

式中：K 為油藏的絕對(duì)滲透率，D；μo、μw分別為油、水相黏度，Pa·s；kro、krw分別為油、水相的相對(duì)滲透率，無(wú)因次；p 為油藏壓力，MPa；qo、qw分別為油、水相在單位時(shí)間、單位體積巖石注入（或采出）的質(zhì)量流量，kg/（m3·s）；ρo、ρw分別為油、水相的密度，kg/m3；So、Sw為油、水相飽和度，%；φ為孔隙度，%；t 為時(shí)間，s。

由式（1）～（2）可得到壓力控制公式［4-5］：

對(duì)式（3）五點(diǎn)隱式差分，把引起非線性的系數(shù)作顯式處理后得：

如果把式（4）寫(xiě)成矩陣的形式，即AU=g，A稱為該線性方程組的系數(shù)矩陣，矩陣結(jié)構(gòu)為五對(duì)角形式。如果在某一個(gè)網(wǎng)格處有井存在，則把它作為點(diǎn)源或點(diǎn)匯來(lái)處理，在對(duì)網(wǎng)格建立的差分方程中增加一個(gè)產(chǎn)量項(xiàng)。

1.2 差分過(guò)程無(wú)效網(wǎng)格的處理方法

差分過(guò)程無(wú)效網(wǎng)格的處理方法有三種。第一種方法是直接計(jì)算法，對(duì)研究區(qū)域網(wǎng)格化處理后，只對(duì)有效網(wǎng)格節(jié)點(diǎn)差分，而與無(wú)效網(wǎng)格相連的有效網(wǎng)格得到差分方程后，求解五對(duì)角矩陣方程，此方法會(huì)因無(wú)效網(wǎng)格數(shù)及所處模擬區(qū)域位置的變化而需重新進(jìn)行差分和計(jì)算五對(duì)角系數(shù)矩陣，這種方法無(wú)法形成統(tǒng)一、有規(guī)律的計(jì)算方法，尤其不適應(yīng)程序化、軟件化數(shù)值模擬。

第二種方法是鏡像反映法［6］，將無(wú)效網(wǎng)格作為周邊與其連接的可流動(dòng)網(wǎng)格的鏡像，使其壓力等于周邊網(wǎng)格的壓力，此時(shí)鏡像網(wǎng)格在計(jì)算過(guò)程中無(wú)流體發(fā)生滲流。塊中心網(wǎng)格見(jiàn)圖1，從圖1 可看出，對(duì)無(wú)效網(wǎng)格i、j 點(diǎn)差值后，采用鏡像反映法處理時(shí)，可使此時(shí)在i+1/2，i-1/2，j+1/2，j-1/2 網(wǎng)格邊界處，無(wú)流量通過(guò)，（i，j）網(wǎng)格內(nèi)無(wú)流量變化，這時(shí)（i，j）網(wǎng)格即為無(wú)效網(wǎng)格。

圖1 塊中心網(wǎng)格

鏡像反映法盡管不需隨著無(wú)效網(wǎng)格的改變而重新計(jì)算方程的系數(shù)矩陣，能形成一個(gè)統(tǒng)一的算法，但需要每一個(gè)時(shí)間步長(zhǎng)內(nèi)設(shè)置無(wú)效網(wǎng)格的壓力，使之與周邊網(wǎng)格壓力相等，計(jì)算過(guò)程繁雜，模擬過(guò)程的不斷賦值，大大降低了計(jì)算效率。

第三種方法是網(wǎng)格設(shè)置法，可直接設(shè)置無(wú)效網(wǎng)格的φ及Ki，j為0，即無(wú)效網(wǎng)格，此時(shí)網(wǎng)格就成為無(wú)效網(wǎng)格。因綜合壓縮系數(shù)為式（3）右邊的綜合壓縮系數(shù)意義已不適用，應(yīng)以展開(kāi)形式帶入公式，以避免存在無(wú)意義情況。另式（4）中的λTi+1/2等系數(shù)，均按“上游權(quán)”的取值原則。由于計(jì)算過(guò)程中網(wǎng)格的滲透率取“上游權(quán)”原則，對(duì)于無(wú)效網(wǎng)格則不需要此判斷過(guò)程，只要確定無(wú)效網(wǎng)格邊界滲透率值為0，見(jiàn)圖1。（i，j）網(wǎng)格為無(wú)效網(wǎng)格點(diǎn)，其邊界滲透率kri+1/2、kri-1/2、krj+1/2和krj-1/2均為0。

由于無(wú)效網(wǎng)格的存在，對(duì)應(yīng)網(wǎng)格點(diǎn)的滲透率為0，差分后的五對(duì)角矩陣行列式為0，無(wú)法求解方程組，因此需要將矩陣中為0 的行、列除去。壓力方程求解完畢后，得到對(duì)應(yīng)網(wǎng)格在tn+1時(shí)刻的壓力值，即壓力場(chǎng)。此方法可直接將不規(guī)則形狀油藏劃為矩形油藏進(jìn)行計(jì)算，并形成了統(tǒng)一的系數(shù)矩陣計(jì)算方法，其系數(shù)矩陣計(jì)算規(guī)律簡(jiǎn)單，計(jì)算速度快。

2 流線數(shù)學(xué)模型與無(wú)效網(wǎng)格的處理

流線數(shù)值模擬方法作為新興的數(shù)值模擬技術(shù)方法，與傳統(tǒng)的數(shù)值模擬方法相比有很多優(yōu)點(diǎn)，能為生產(chǎn)實(shí)踐提供指導(dǎo)和技術(shù)支持。當(dāng)模擬區(qū)塊外部、內(nèi)部存在不可流動(dòng)區(qū)域時(shí)，無(wú)效網(wǎng)格處理方法能夠有效解決數(shù)值模擬過(guò)程中的壓力求解、流線追蹤以及流線飽和度計(jì)算等因存在無(wú)效網(wǎng)格而出現(xiàn)的問(wèn)題。

2.1 流線追蹤與無(wú)效網(wǎng)格的處理

在存在無(wú)效網(wǎng)格的流線數(shù)值模擬中，流線應(yīng)避開(kāi)無(wú)效網(wǎng)格區(qū)域。Pollock 提出的追蹤流線軌跡方法是在網(wǎng)格系統(tǒng)中壓力場(chǎng)已知的情況下，應(yīng)用達(dá)西定律建立流體真實(shí)流動(dòng)速度場(chǎng)，然后在此基礎(chǔ)上追蹤流線。

考慮二維情況，根據(jù)達(dá)西定律計(jì)算每個(gè)網(wǎng)格界面上的速度分量［7-8］：

圖2 中采用Pollock 方法進(jìn)行流線追蹤。x 方向的速度定義為：

圖2 流線穿過(guò)網(wǎng)格過(guò)程

對(duì)于無(wú)效網(wǎng)格，應(yīng)沒(méi)有流體流入或流出，網(wǎng)格邊界封閉，無(wú)效網(wǎng)格邊界不存在流體速度，因此無(wú)效網(wǎng)格邊界速度應(yīng)取為0。

為了避免速度為0 而出現(xiàn)無(wú)意義情況，可將無(wú)效網(wǎng)格邊界速度在一定精度下取值無(wú)窮小，這樣可模擬流線在遇到無(wú)效網(wǎng)格時(shí)的情況，此時(shí)流線流動(dòng)到無(wú)效網(wǎng)格邊界的時(shí)間無(wú)窮大，即不可到達(dá)。完整的流線追蹤過(guò)程可得到流線流動(dòng)的完整路徑［9］。

2.2 流線飽和度模型與無(wú)效網(wǎng)格處理

要采用流線方法求解生產(chǎn)階段的滲流數(shù)學(xué)模型，首先要將基于網(wǎng)格的二維或三維模型轉(zhuǎn)化為沿流線的一維模型。

定義沿著流線的傳播時(shí)間（Timeof Flight，TOF）［10］為：，處理后最終得到：即流線飽和度方程。

無(wú)效網(wǎng)格內(nèi)不存在流體流動(dòng)、變化，因此流線不經(jīng)過(guò)無(wú)效網(wǎng)格，飽和度方程僅跟流線所經(jīng)過(guò)的網(wǎng)格有關(guān)，因此無(wú)效網(wǎng)格不影響其他網(wǎng)格飽和度的計(jì)算。

3 應(yīng)用分析

選取長(zhǎng)慶定邊姬塬油田羅1 井區(qū)長(zhǎng)8 油藏內(nèi)的4口采油井、1 口注水井作為應(yīng)用對(duì)象，根據(jù)所研究區(qū)域儲(chǔ)層主要參數(shù)及儲(chǔ)層流體特性資料作為流線數(shù)值模擬的基本參數(shù)。

網(wǎng)格系統(tǒng)采用塊中心網(wǎng)格系統(tǒng)，縱向不劃分網(wǎng)格，平面上不考慮儲(chǔ)層的非均值性，將平面劃分為19×19×1 網(wǎng)格系統(tǒng)。無(wú)效網(wǎng)格區(qū)域分布在3 個(gè)區(qū)域。模擬區(qū)域中央為注水井，4 口采油井分布在儲(chǔ)層四邊，無(wú)效網(wǎng)格分布在下面三個(gè)部分，分別為左下邊的單個(gè)無(wú)效網(wǎng)格、2 個(gè)連接無(wú)效網(wǎng)格以及右下邊6 個(gè)相連的無(wú)效網(wǎng)格。

根據(jù)前面的數(shù)學(xué)模型及第三種無(wú)效網(wǎng)格處理方法，編寫(xiě)了流線數(shù)值模擬計(jì)算程序，模擬3 月后，得到了模擬區(qū)域1 注4 采井網(wǎng)形式的流線分布特征，見(jiàn)圖3。

圖3 存在無(wú)效網(wǎng)格流線分布

從圖3 中可以看出流線從中間1 口注水井出發(fā)，流向周邊的4 口生產(chǎn)井，流線碰到無(wú)效網(wǎng)格區(qū)域會(huì)明顯出現(xiàn)變化，而三部分不可流區(qū)域無(wú)流線通過(guò)，流線繞過(guò)不可流動(dòng)區(qū)域，流向了注水井，流線明顯改變了流線原有形狀，證明了無(wú)效網(wǎng)格處理的有效性及合理性。

由于不可流動(dòng)區(qū)域的存在，靠近不可流動(dòng)區(qū)域，流線分布較密，而流線分布密集地區(qū)反映了注水流量大，驅(qū)油效果好。

4 結(jié)論

a）差分過(guò)程無(wú)效網(wǎng)格處理方法中，網(wǎng)格設(shè)置法計(jì)算最為簡(jiǎn)便，方程系數(shù)矩陣算法最容易得出，直接計(jì)算法無(wú)法形成統(tǒng)一、有規(guī)律的計(jì)算方法，尤其不適應(yīng)程序化、軟件化數(shù)值模擬；鏡像反映法計(jì)算過(guò)程繁雜，模擬過(guò)程的不斷賦值，大大降低了計(jì)算的效率。網(wǎng)格設(shè)置法彌補(bǔ)了直接計(jì)算法和鏡像反映法的不足，既形成了統(tǒng)一的系數(shù)矩陣計(jì)算方法，又有相對(duì)較快的計(jì)算速度。

b）網(wǎng)格設(shè)置法可將不規(guī)則形狀油藏按規(guī)則的矩形油藏進(jìn)行數(shù)值模擬計(jì)算，并形成統(tǒng)一的系數(shù)矩陣計(jì)算方法，簡(jiǎn)化了模擬計(jì)算過(guò)程。

c）Pollock 的追蹤流線軌跡方法在進(jìn)行追蹤流線軌跡過(guò)程中，因無(wú)效網(wǎng)格的存在，在處理網(wǎng)格邊界時(shí)，會(huì)影響流線變化，改變飽和度的計(jì)算過(guò)程。

d）應(yīng)用無(wú)效格處理方法可解決含啟動(dòng)壓力梯度的情況下油藏動(dòng)邊界問(wèn)題，因油水兩相方程動(dòng)邊界很難給出解析形式，所以需要應(yīng)用數(shù)值方法求取確定，通過(guò)小時(shí)間步長(zhǎng)計(jì)算確定動(dòng)邊界范圍，再將動(dòng)邊界外部確定為無(wú)效網(wǎng)格，這樣可解決因啟動(dòng)壓力梯度的存在，動(dòng)邊界外部流體滲流速度為負(fù)的情況。

［1］Batycky R P A.Three-Dimensional Two-Phase Scale Streamline Simulator［D］.Department of Petroleum Engineering，School of Earth Science，Stanford University，1997.

［2］劉慧卿.油藏?cái)?shù)值模擬方法專題［M］.東營(yíng)：石油大學(xué)出版社，2001.Liu Huiqing.Topic of Numerical Simulation of Oil Reservoir［M］.Dongying：China University of Petroleum Press，2001.

［3］郭鳴黎，程?hào)|風(fēng)，李大勇.文25 東復(fù)雜斷塊油藏剩余油分布研究［J］.江漢石油學(xué)院學(xué)報(bào)，2003，25（3）：82-83.Guo Mingli，Cheng Dongfeng，Li Dayong.Remaining Oil Distribation of Complex Fault Block Reserovir in the East of Block Wen 25［J］.Journal of Jianghan Petroleum Institute，2003，25（3）：82-83.

［4］尹 虎，王新海，劉 洪，等.考慮啟動(dòng)壓力梯度的頁(yè)巖氣藏?cái)?shù)值模擬［J］.天然氣與石油，2012，30（4）：43-45.Yi Hu，Wang Xinhai，Liu Hong，et al.Numerical Simulation of Shale Gas Reservoir Considering Startup Pressure Gradient［J］.Natural Gas and Oil，2012，30（4）：43-45.

［5］王 平，姜瑞忠，王公昌.高含水期水驅(qū)狀況影響因素?cái)?shù)值模擬研究［J］.天然氣與石油，2012，30（4）：36-38.Wang Ping，Jiang Ruizhong，Wang Gongchang.Numerical Simulation of Factors Affecting Water-Flooding in High Water-Cut Stage［J］.Natural Gas and Oil，2012，30（4）：36-38.

［6］韓大匡，陳欽雷，閆存章.油藏?cái)?shù)值模擬基礎(chǔ)［M］.北京：石油工業(yè)出版社，1993.87-93.Han Dakuang，Chen Qinlei，Yan Chunzhang.Fundamentals of Numerical Reservoir Simulation［M］.Beijing：Petroleum Industry Press，1993.87-93.

［7］Portella R C M，Hewett T A.Fast 3-D Reservoir Simulation and Applications Using Streamlines［C］.SPE 39061，1998.

［8］姚 軍，吳明錄，戴衛(wèi)華，等.流線數(shù)值試井解釋模型［J］.石油學(xué)報(bào)，2006，27（3）：96-99.Yao Jun，Wu Minglu，Dai Weihua，et al.Streamline Numerical Well Test Interpretation Model［J］.Acta Petrolei Sinica，2006，27（3）：96-99.

［9］Batycky R P.Martin J B，Marco R T.A 3D Field Scale Streamline Simulator With Gravity and Changing Well Conditions［C］.SPE 36726，1996.

［10］吳明錄，姚 軍.不規(guī)則污染井的流線數(shù)值試井解釋模型及其應(yīng)用［J］，新疆石油地質(zhì)，2009，60（3）：373-375.Wu Minglu，Yao Jun.Streamline Numerical Well -Testing Interpretation Model for Irregularly Contaminated Wells with Application［J］，Xinjiang Petroleum Geology，2009，60（3）：373-375.