李瑞軒，黃云龍，李 源

（西安石油大學(xué)石油工程學(xué)院，陜西西安 710300）

隨著油氣工業(yè)的快速發(fā)展，我國(guó)油氣田開發(fā)已經(jīng)進(jìn)入瓶頸期，從世界經(jīng)濟(jì)環(huán)境來看，面對(duì)日益嚴(yán)峻的資源枯竭、能源危機(jī)等問題，國(guó)家正在加大對(duì)于油氣資源的開發(fā)和利用。我國(guó)油氣資源消費(fèi)量與進(jìn)口量逐年增加，這一現(xiàn)象表現(xiàn)出對(duì)國(guó)外資源依賴程度明顯上升，這已經(jīng)影響到國(guó)家的經(jīng)濟(jì)安全。與世界其他富油地區(qū)相比，中國(guó)的石油資源有限且開采技術(shù)難度大，面臨許多世界級(jí)技術(shù)難題，若想實(shí)現(xiàn)我國(guó)油氣自給自足，則必須加大技術(shù)投入，降低生產(chǎn)成本，增加經(jīng)濟(jì)效益。通過建立油藏滲流模型，可以預(yù)測(cè)未來開采時(shí)段油氣水流動(dòng)情況和儲(chǔ)層變化情況，能夠幫助油藏工程師合理的調(diào)整生產(chǎn)方案，經(jīng)濟(jì)高效的開采油氣資源。

文章根據(jù)我國(guó)油藏?cái)?shù)值模擬技術(shù)的發(fā)展歷程，將國(guó)內(nèi)應(yīng)用于油藏的滲流模型劃分為三類，總結(jié)前人在具體模型發(fā)展上做出的成果與貢獻(xiàn)，旨在分析各自模型發(fā)展的內(nèi)在聯(lián)系與區(qū)別，探究我國(guó)油藏滲流基礎(chǔ)理論與技術(shù)研究進(jìn)展，把握未來油藏滲流模型發(fā)展趨勢(shì)，為我國(guó)油藏的高效開發(fā)提供實(shí)際依據(jù)。

1 油藏滲流模型簡(jiǎn)述

根據(jù)國(guó)內(nèi)學(xué)者研究的領(lǐng)域和得到模型所運(yùn)用的數(shù)學(xué)方法，大致可以分為三類：解析模型、半解析模型和數(shù)值模型。每種模型的代表性成果總結(jié)如下：

1.1 解析模型

解析法通常采用一系列假設(shè)簡(jiǎn)化實(shí)際油藏，求解的模型較為理想，采用數(shù)學(xué)推導(dǎo)、演繹而求解油藏滲流模型的方法，可取得簡(jiǎn)化模型的精確解。

1989 年戴榕菁等[1]采用Weber 變換，得出了有越流的無(wú)限大、多層油藏定壓生產(chǎn)的產(chǎn)量相關(guān)函數(shù)和油藏壓力分布的解析解。1990 年張洪亮等[2]考慮了油井表皮系數(shù)，研究了層間越流現(xiàn)象的n 層恒定井底壓力，考慮了層內(nèi)垂向滲透率模型得出三層模型的解析解。1995 年伊?xí)詵|等[3]提出了均質(zhì)油藏?zé)o限大N 層地層試井模型，得出Laplace 空間解和實(shí)空間函數(shù)的精確解，分別用拉普拉斯空間解和實(shí)空間解對(duì)結(jié)果進(jìn)行驗(yàn)證。2000 年劉月田等[4]討論了橢圓形區(qū)域含中心點(diǎn)源的位勢(shì)流動(dòng)，采用相似流動(dòng)替換法得出解析解。文章采用位勢(shì)疊加原理并考慮到相互平行的等壓直線邊界間的點(diǎn)源到邊界直線的流動(dòng)，方法具有其獨(dú)創(chuàng)性。2001 年向開理等[5]引入分形幾何理論，考慮到井筒續(xù)流、井筒相分布和油井表皮等的影響，建立了具有不穩(wěn)定滲流效應(yīng)的分形油藏有效井徑模型，并采用Laplace 變換得到了該模型的解析解。2008 年吳小慶等[6]通過對(duì)溶洞型油藏的研究，并采用正交變換法建立了相關(guān)試井模型，同時(shí)得到其解析解，他們考慮到方程邊界條件的特點(diǎn)，忽略掉重力作用效果將求解偏微分方程初值轉(zhuǎn)化成求解常微分方程初值，對(duì)其模型進(jìn)行了一系列完善和求解。

實(shí)際油藏研究中需要考慮的影響因素較多，較為復(fù)雜，而解析法通常不能考慮開發(fā)時(shí)間對(duì)油藏各項(xiàng)參數(shù)的影響，且求解方法較為單一；再者未能考慮天然裂縫網(wǎng)絡(luò)對(duì)油藏滲流的影響，且只能解決穩(wěn)態(tài)單相達(dá)西滲流問題，在實(shí)際油藏滲流中計(jì)算誤差較大。由此油田工作者們引出了半解析模型和更為精確的數(shù)值模型，進(jìn)一步精確刻畫油田實(shí)際情況，對(duì)油藏開發(fā)做出更細(xì)致的模擬。

1.2 半解析模型

半解析法通常運(yùn)用點(diǎn)源函數(shù)，將油藏劃分成每一小部分，用函數(shù)描述每一小部分油藏的流動(dòng)狀態(tài)，之后對(duì)每一個(gè)劃分的小部分進(jìn)行疊加，所求出其方程的半解析解。但是半解析法不能很精準(zhǔn)的描述油藏整體滲流。

1990 年虞紹永等[7]根據(jù)在雙滲地層中的段塞流的壓力特征分析，提出了相應(yīng)的滲流模型，對(duì)相應(yīng)的模型應(yīng)用了Laplace 變換和數(shù)值反演的方法，并考慮幾種極限狀態(tài)，求出了相對(duì)應(yīng)的半解析解。2004 年韓國(guó)慶等[8]建立了分支井的半解析產(chǎn)能預(yù)測(cè)模型。其中利用格林函數(shù)處理油藏壓力模型響應(yīng)，求解出三個(gè)瞬時(shí)點(diǎn)源解，之后沿油井軌跡進(jìn)行積分，求出壓力分布。2006 年郝斐等[9]針對(duì)低滲透油藏滲流規(guī)律，引入玻爾茲曼變換，考慮邊界條件，求出了地層壓力分布的函數(shù)。2008 年安永生等[10]根據(jù)油藏的各向異性，利用坐標(biāo)變換的方法將油藏的參數(shù)轉(zhuǎn)化成為同性空間進(jìn)行求解，再與已有的半解析模型進(jìn)行耦合，求出其新的預(yù)測(cè)模型。2009年鄧玉珍[11]將不穩(wěn)定滲流的過程劃分成每個(gè)小時(shí)間點(diǎn)，再把每個(gè)小時(shí)間點(diǎn)的滲流看成穩(wěn)定滲流狀態(tài)。求出每個(gè)小時(shí)間點(diǎn)的滲流流量和壓力，在同一累加就是總的流量和壓力。以及2009 年吳曉東等[12]使用這種方法求出了多分支水平井半解析產(chǎn)能預(yù)測(cè)模型。2015 年楊萬(wàn)有等[13]針對(duì)各向異性蛇曲井的滲流規(guī)律，更新了井筒流動(dòng)壓力梯度計(jì)算方法，與已有的半解析模型進(jìn)行更好的融合。2015 年孫玉豹等[14]根據(jù)紐曼乘積方法確定源函數(shù)，再利用Laplace 變換、數(shù)值反演等方法求解了油藏半解析解。2015 年賈品等[15]引用Laplace 空間源函數(shù)求得油藏流動(dòng)解析解。2016 年任宗孝等[16]基于Valko體積源函數(shù)描述縫網(wǎng)的基礎(chǔ)上，建立了致密油藏體積壓裂水平井解析解模型，但此方法假設(shè)縫網(wǎng)太過理想，于是任宗孝[17]在2018 年更新了模型，基于Ozkan 體積源函數(shù)描述縫網(wǎng)的基礎(chǔ)上，又利用“星-三角形”變換方法，建立了致密油藏體積壓裂水平井解析解模型。2017年王家航等[18]利用內(nèi)外邊界Green 函數(shù)建立了致密油藏直井體積壓裂非穩(wěn)態(tài)壓力分析半解析模型。2017 年王家航等[19]還是利用了之前的方法，建立了非均質(zhì)油藏多段壓裂水平井非穩(wěn)態(tài)壓力分析半解析模型。2018年侯紹繼等[20]建立了徑向低滲非Darcy 滲流模型，并利用Laplace 變換結(jié)合無(wú)窮級(jí)數(shù)法對(duì)模型進(jìn)行了求解。

由此可見，國(guó)內(nèi)大多數(shù)學(xué)者都是采用Laplace 變換對(duì)油藏滲流方程進(jìn)行變換并求解。也有一些學(xué)者利用源函數(shù)對(duì)油藏滲流方程進(jìn)行求解。但使用的源函數(shù)不能精準(zhǔn)描述邊界復(fù)雜油氣層儲(chǔ)層。但是半解析法不能用于油藏?cái)?shù)值模擬。

1.3 數(shù)值模型

對(duì)于油藏?cái)?shù)值模擬技術(shù)，通過文獻(xiàn)調(diào)研發(fā)現(xiàn)，國(guó)內(nèi)外基本都采用兩類模擬油藏的滲流過程的方法，一種是基于多重連續(xù)介質(zhì)滲流理論的方法，另外一種是基于離散裂縫模型的方法。第一種方法將真實(shí)裂縫處理成連續(xù)介質(zhì)，將油藏的流動(dòng)空間采用多重連續(xù)介質(zhì)描述；第二種方法將裂縫處理為離散介質(zhì)，形成離散裂縫模型，此方法的好處在于建立的裂縫能夠最大程度還原真實(shí)狀態(tài)。而且連續(xù)介質(zhì)模型還具有處理微小裂縫的優(yōu)點(diǎn)。

隨著我國(guó)油氣田開采技術(shù)的快速發(fā)展和對(duì)油藏的深入研究，越來越需要更加貼合實(shí)際情況的油藏滲流模型。采用的數(shù)值求解方法是將偏微分方程離散后進(jìn)行計(jì)算，通過數(shù)學(xué)運(yùn)算，進(jìn)而計(jì)算出油藏區(qū)域上每個(gè)具體坐標(biāo)點(diǎn)上的值。雖然數(shù)值模擬對(duì)方程進(jìn)行近似處理的時(shí)候會(huì)存在一定的誤差，得到的是近似解，但這種方法卻比解析模型和半解析模型更能真實(shí)地反映出復(fù)雜條件下地層的滲流規(guī)律。早在1953 年，Peacema 等率先將數(shù)值模擬方法應(yīng)用于氣相流動(dòng)的問題。在此之后，這種模擬方法便大量用于研究油藏滲流規(guī)律。時(shí)至今日，各國(guó)學(xué)者對(duì)于解析模型和數(shù)值模型及其求解方法都進(jìn)行了大量的創(chuàng)新和改進(jìn)，使得油藏?cái)?shù)值模擬領(lǐng)域愈漸先進(jìn)。

2003 年姜漢橋等[21]得出了油藏流體滲流的不完全可逆變形介質(zhì)模型,考慮到Richtmyer 線性化和隱式處理狀態(tài)方程，得出了其滲流模型的數(shù)值解。2004 年何書梅等[22]通過研究強(qiáng)非均質(zhì)油藏板9 斷塊的各項(xiàng)開發(fā)特征，得出了該斷塊油藏的滲流模型，采用邊界控制不規(guī)則網(wǎng)格系統(tǒng)優(yōu)化網(wǎng)格系統(tǒng)和隱式差分方法對(duì)該斷塊進(jìn)行了數(shù)值模擬研究。2004 年王正茂等[23]通過考慮骨架砂剝離情況、體積應(yīng)變現(xiàn)象以及可動(dòng)砂在孔隙中的沉積對(duì)儲(chǔ)層各項(xiàng)物性參數(shù)的影響，導(dǎo)出了出砂儲(chǔ)層物性參數(shù)的動(dòng)態(tài)模型，隨后采取有限差分法對(duì)模型進(jìn)行求解，并采用有限方法求解了巖土平衡方程，將二者聯(lián)合后求得了模型的數(shù)值解。2004 年劉偉[24]綜合考慮了堵劑驅(qū)油過程中的物理特性，確定了凍膠類堵劑與地層以及儲(chǔ)層內(nèi)流體之間關(guān)系的局部加密網(wǎng)格多輪次調(diào)剖滲流模型，采用行索引存儲(chǔ)方案和預(yù)處理雙共軛梯度算法進(jìn)行了模型求解。2006 年高慧梅、姜漢橋等[25]通過研究滲透率與壓力的指數(shù)增長(zhǎng)規(guī)律，利用邊界元法得出油藏內(nèi)相關(guān)點(diǎn)的壓力，并采用正則攝動(dòng)法處理方程，得出雙重介質(zhì)油藏的任意形狀應(yīng)力敏感性不穩(wěn)定滲流模型。2012 年張?jiān)省⒀α羀26]對(duì)其致密油藏?cái)?shù)值模擬模型進(jìn)行有限元離散，通過離散連續(xù)介質(zhì)和采取三角形網(wǎng)格劃分油藏區(qū)域，進(jìn)而提高了其模型精度。2012年姚軍、劉文超[27]對(duì)非達(dá)西滲流低滲透油藏動(dòng)邊界數(shù)學(xué)模型進(jìn)行坐標(biāo)變換，將其轉(zhuǎn)變?yōu)槎ㄟ吔鐢?shù)學(xué)模型，隨后采取有限差分法和全隱式求解法，得到了定流量?jī)?nèi)邊界與定井底壓力內(nèi)邊界的模型數(shù)值解。2013 年陳忠[28]通過全隱式差分對(duì)低滲透油藏油水兩相非線性滲流模型方程進(jìn)行離散，構(gòu)造了非線性滲流數(shù)值模型，提出擬牛頓法討論數(shù)值模型，并驗(yàn)證了算法的收斂性，此方法對(duì)于低滲透油藏中小規(guī)模數(shù)值模擬收斂速度較快。通過研究低滲透多孔介質(zhì)滲流邊界模型采取空間坐標(biāo)變換法處理模型邊界，并選擇隱式有限差分得出模型的數(shù)值解。2015 年李威威[29]也應(yīng)用隱式有限差分方法對(duì)模型進(jìn)行求解，但應(yīng)用本征正交分解法構(gòu)造函數(shù)對(duì)微分方程進(jìn)行降低方程維數(shù)求解。2018 年劉勇等[30]針對(duì)裂縫性油藏求解，針對(duì)非結(jié)構(gòu)四邊形網(wǎng)格的計(jì)算效率過低的現(xiàn)象，提出了三種提高計(jì)算效率的計(jì)算方法。2019 年黃迎松[31]對(duì)壓裂井組非線性滲流模型進(jìn)行有限差分求解，其中利用Taylor 公式展開對(duì)模型進(jìn)行線性化處理，之后進(jìn)行數(shù)值求解。

2 求解非線性偏微分方程的常見方法與研究

2.1 求解非線性偏微分方程的常見方法

近年來我國(guó)油藏?cái)?shù)值模擬技術(shù)得到了廣泛的研究和應(yīng)用，各種成熟的商業(yè)軟件不斷推出，為我國(guó)油氣田開發(fā)做出了顯著貢獻(xiàn)。而油藏?cái)?shù)值模擬的重點(diǎn)問題在于處理非線性方程的數(shù)值解，現(xiàn)階段常用于求解非線性偏微分方程的數(shù)值近似方法包括有限差分法、有限元法和有限體積法。

有限差分方法（FDM）是采用最早的數(shù)值模擬方法，至今也仍被廣泛應(yīng)用。該方法將求解區(qū)域劃分為性質(zhì)單一的網(wǎng)格，每個(gè)網(wǎng)格用單一節(jié)點(diǎn)來表示，利用單一節(jié)點(diǎn)表示求解區(qū)域。有限差分法利用Taylor 級(jí)數(shù)展開等方法，每一個(gè)控制方程的導(dǎo)數(shù)用每一個(gè)節(jié)點(diǎn)的差商值來代替，對(duì)方程進(jìn)行離散化處理，從而建立代數(shù)方程組。該方法直接將微分問題轉(zhuǎn)變?yōu)榇鷶?shù)問題的近似數(shù)值解法，表達(dá)形式簡(jiǎn)單，是發(fā)展較早且比較成熟的數(shù)值計(jì)算方法。

有限元方法的基礎(chǔ)是變分原理和加權(quán)余量法，其求解思想是把求解區(qū)域劃分為有限個(gè)單元，但每個(gè)單元形狀不一，在每個(gè)單元內(nèi)，選擇合適的節(jié)點(diǎn)作為求解函數(shù)的插值點(diǎn)，將微分方程改寫成插值函數(shù)的線性表達(dá)式，借助于變分原理或加權(quán)余量法，將微分方程離散求解。

有限體積法（FVM）又稱為控制體積法。利用計(jì)算流體力學(xué)中的思想，將計(jì)算區(qū)域劃分為每個(gè)單一的控制體，每一個(gè)網(wǎng)格周圍都是一個(gè)單一的控制體，對(duì)每一個(gè)控制體進(jìn)行積分，得出一組離散方程。其中網(wǎng)格點(diǎn)上的因變量的數(shù)值是未知數(shù)。為了求出控制體積的積分，必須假定數(shù)值在網(wǎng)格點(diǎn)之間的變化規(guī)律。從而對(duì)方程進(jìn)行求解。

通常主要用有限差分法求解微分方程。通過將有限差分網(wǎng)格的疊加，然后利用網(wǎng)格系統(tǒng)對(duì)方程進(jìn)行離散，通常采用結(jié)構(gòu)化網(wǎng)格，如塊中心網(wǎng)格、點(diǎn)中心網(wǎng)格和徑向網(wǎng)格等。其方法簡(jiǎn)單、概念明確，但存在網(wǎng)格取向性嚴(yán)重、處理復(fù)雜邊界較為困難、求解精度低等缺點(diǎn)。

有限元方法在處理偏微分方程中計(jì)算精度高，且易于處理復(fù)雜邊界，將連續(xù)的模擬油藏空間離散為有限多個(gè)單元的組合體，由于單元本身形狀和相互間連接方式的變化，可以對(duì)幾何形狀復(fù)雜的儲(chǔ)層進(jìn)行模擬。并利用片分思想將全域上的未知場(chǎng)函數(shù)由每個(gè)單元內(nèi)近似插值函數(shù)來表示，最后通過整體疊加思想進(jìn)行組裝求解。

為有效利用網(wǎng)格，充分考慮邊界的復(fù)雜性，可見采用有限元方法建立流動(dòng)方程的數(shù)值計(jì)算方式存在明顯優(yōu)勢(shì)，其優(yōu)點(diǎn)可概括為：（1）便于對(duì)含有復(fù)雜邊界的油藏區(qū)域進(jìn)行整體剖分；（2）避免了直接構(gòu)建全域近似函數(shù)；（3）處理偏微分方程的計(jì)算精度較高。

2.2 有限差分法在油藏研究方面的應(yīng)用成果

饒盛文等[32]考慮到低滲透油藏中存在啟動(dòng)壓力梯度，故其滲流規(guī)律不再符合經(jīng)典的達(dá)西定律，在結(jié)合了啟動(dòng)壓力梯度和毛管力對(duì)滲流的影響因素，建立了徑向水驅(qū)油兩相滲流微分方程，并利用有限差分法進(jìn)行數(shù)值求解。此外饒盛文[33]通過油水兩相滲流理論和相關(guān)實(shí)驗(yàn)成果，推導(dǎo)建立了致密油藏油水兩相滲流的一般微分方程，并介紹了低滲透油藏油、水兩相滲流公式及求解方法，即采用隱式聯(lián)立求解的方法進(jìn)行求解。對(duì)低滲油藏油水兩相模型的建立、計(jì)算方法等方面作了一定深度創(chuàng)新。

2.3 有限元法在油藏研究方面的應(yīng)用成果

有限元法采用若干個(gè)互相連接不重合的小單元對(duì)求解區(qū)域進(jìn)行空間離散，后在各個(gè)單元內(nèi)選取結(jié)點(diǎn)，進(jìn)而進(jìn)行插值并獲得近似函數(shù)，通常由未知場(chǎng)變量函數(shù)在各個(gè)單元節(jié)點(diǎn)上的數(shù)值或插值表達(dá)單元中的近似函數(shù)。這種方法將一個(gè)自由度為無(wú)限的連續(xù)區(qū)域轉(zhuǎn)化為有限多個(gè)單元集合體，從而通過合并求解。

張?jiān)实萚26]研究了致密油藏非達(dá)西滲流有限元數(shù)值模擬，考慮到有限單元法對(duì)復(fù)雜邊界的適用性較好、計(jì)算精度較高等特點(diǎn)，研究了油藏非達(dá)西滲流數(shù)值模擬方法，并編程了對(duì)應(yīng)的數(shù)模軟件。

呂杭[34]把有限元方法應(yīng)用于三維油、水兩相裂縫性油藏，并編寫相應(yīng)的模擬軟件。主要考慮了單滲透率雙孔隙度裂縫模型，將油、水井處理成單元邊界，并采用等參元坐標(biāo)變換對(duì)其進(jìn)行求解，拓展了有限元法在裂縫油藏?cái)?shù)值模擬中的應(yīng)用。

程林松等[35]采用六面體單元對(duì)三維直井模型中裂縫性油藏進(jìn)行了空間離散，應(yīng)用隱式壓力隱式飽和度法和有限元法求解得出了數(shù)值結(jié)果。

蔣廷學(xué)等[36]考慮啟動(dòng)壓力梯度、油-水兩相流動(dòng)、變裂縫導(dǎo)流能力等因素，對(duì)人工裂縫基于等效滲流阻力法應(yīng)用有限元法進(jìn)行了離散剖分。

2.4 有限體積法在油藏研究方面的應(yīng)用成果

有限差分方法滿足局部守恒原理，但是對(duì)于復(fù)雜邊界儲(chǔ)層問題靈活性較差，難以適用；有限元方法由于基于非結(jié)構(gòu)網(wǎng)格，在空間離散上的靈活性增強(qiáng)，但僅符合全局網(wǎng)格內(nèi)的物質(zhì)守恒原則，不能實(shí)現(xiàn)局部區(qū)域內(nèi)的流體物質(zhì)守恒，極容易產(chǎn)生數(shù)值解的震蕩現(xiàn)象。而有限體積法結(jié)合了有限差分和有限元法的優(yōu)點(diǎn)，既可確保局部網(wǎng)格內(nèi)的流體流動(dòng)守恒，同時(shí)也能滿足復(fù)雜邊界儲(chǔ)層的問題處理。

呂心瑞[37]在Delaunay 三角形基礎(chǔ)上基于離散裂縫模型構(gòu)建了有限體積控制單元，推導(dǎo)油-水兩相流動(dòng)控制方程的有限體積數(shù)值格式，對(duì)裂縫性油藏油-水兩相流動(dòng)進(jìn)行了數(shù)值模擬研究。

張芮菡[38]通過研究有天然大型裂縫的非均質(zhì)性強(qiáng)、尺度變化多的油藏，結(jié)合了有限元法和有限體積法，得出了離散型裂縫的網(wǎng)格模型的解，隨后考慮到物質(zhì)平衡、運(yùn)動(dòng)方程、狀態(tài)方程以及連續(xù)性方程，進(jìn)而得出了不混溶的微相可壓縮油水兩相滲流方程。

3 結(jié)論

（1）根據(jù)三種不同類型的滲流模型，列舉了近年來我國(guó)學(xué)者所發(fā)表的代表性成果。它們大多數(shù)都致力于模擬最真實(shí)原油滲流方式，有些影響因素考慮的過于理想化并受限于解決條件，其都與真實(shí)油藏裂縫滲流方式有一定的差距。所以后來的學(xué)者應(yīng)當(dāng)不斷優(yōu)化現(xiàn)有模型，才能精確的模擬真實(shí)油藏滲流狀態(tài)。

（2）介紹了常見求解非線性方程的三種方法，包括有限差分法、有限體積法和有限元法；分析了三種方法的相應(yīng)原理，敘述了有限元法在油藏?cái)?shù)值模擬方面的優(yōu)勢(shì)：①便于對(duì)含有復(fù)雜邊界的油藏區(qū)域進(jìn)行整體剖分；②避免了直接構(gòu)建全域近似函數(shù)；③處理偏微分方程的計(jì)算精度較高；并列舉了三種方法在油藏?cái)?shù)值模擬方面的部分應(yīng)用成果。