二維非均勻多孔介質(zhì)中不可壓兩相驅(qū)替的有限分析算法

鄭曉磊，劉志峰，王曉宏，施安峰

（中國(guó)科學(xué)技術(shù)大學(xué)熱科學(xué)和能源工程系，安徽合肥 230026）

文章編號(hào)：1001?246X（2015）05?0586?09

二維非均勻多孔介質(zhì)中不可壓兩相驅(qū)替的有限分析算法

鄭曉磊，劉志峰，王曉宏，施安峰

（中國(guó)科學(xué)技術(shù)大學(xué)熱科學(xué)和能源工程系，安徽合肥 230026）

為提高油藏?cái)?shù)值模擬算法的計(jì)算效率，在求解單向穩(wěn)態(tài)滲流的有限分析算法基礎(chǔ)上，構(gòu)建二維非均勻多孔介質(zhì)中不可壓兩相滲流的有限分析算法.算法中，網(wǎng)格界面上的平均滲透率不是簡(jiǎn)單地取為相鄰網(wǎng)格滲透率的調(diào)和平均值，而是通過奇點(diǎn)鄰域解析解積分求得.相比于傳統(tǒng)的數(shù)值算法，有限分析算法隨著網(wǎng)格的加密，能夠很快地收斂（僅需將原始網(wǎng)格細(xì)分至2×2或3×3），并且其計(jì)算精度和收斂性不依賴于介質(zhì)的非均勻強(qiáng)度，從而計(jì)算效率得到提高.

有限分析算法；多相流動(dòng)；多孔介質(zhì)；非均勻介質(zhì)；多尺度模擬

0 引言

多孔介質(zhì)中的多相流動(dòng)在很多科學(xué)和工程領(lǐng)域廣泛存在，比如石油開采工業(yè)中常見的水驅(qū)油和蒸汽驅(qū)油等.其它的例子還包括干燥過程、多相流化床反應(yīng)器，熱管中的二元混合物流動(dòng)、以及地?zé)醿?chǔ)能中的鹽水混合物流動(dòng)等等.多相流問題的強(qiáng)非線性使得理論分析異常困難，從而數(shù)值模擬成為研究相關(guān)問題的必要手段.以石油工業(yè)為例，在實(shí)踐中廣泛使用的油藏?cái)?shù)值模擬技術(shù)主要用于預(yù)測(cè)儲(chǔ)層性能，比較各種采收方案，以及測(cè)試不同的運(yùn)營(yíng)策略等等，其目的是實(shí)現(xiàn)企業(yè)利潤(rùn)的最大化.出于這個(gè)原因，已經(jīng)有大量發(fā)展較為成熟的數(shù)值算法在實(shí)際中運(yùn)用［1］，如聯(lián)立求解法（SS）和隱式壓力－顯示飽和度解法（IMPES）［2］.SS方法［3－4］是把油相和水相壓力作為直接的待求變量求解，而飽和度則作為毛管力的函數(shù)，根據(jù)毛管力得到飽和度.IMPES方法［5］的最基本思想是假定毛管力梯度在每個(gè)時(shí)間步里是不變的，利用描述多相流動(dòng)的方程組得到其中某一相的壓力方程，隱式求解該壓力方程從而進(jìn)一步去更新每個(gè)時(shí)間步內(nèi)的飽和度值.如果不假定毛管力梯度在每個(gè)時(shí)間步里是定值，飽和度將隱式求解，就可以得到一種全隱式的方法—順序求解法（SEQ）［2］.

在以上提到的諸方法中，網(wǎng)格節(jié)點(diǎn)間的絕對(duì)滲透率通常都定義為相鄰網(wǎng)格絕對(duì)滲透率的調(diào)和平均值，相對(duì)滲透率則由上游節(jié)點(diǎn)的飽和度確定［2］.眾所周知，在非均勻性很強(qiáng)的情況下，調(diào)和平均算法會(huì)嚴(yán)重低估網(wǎng)格間的界面流量.Cordazzo給出了一個(gè)4×4的單相流動(dòng)算例，如采用調(diào)和平均算法，需要將原始網(wǎng)格細(xì)分至700×700的密網(wǎng)格，其計(jì)算結(jié)果才能與真值相近［6－7］.Romeu和Noetinger建立了一種理論分析方法來研究等效滲透率計(jì)算的準(zhǔn)確性.他們證實(shí)了傳統(tǒng)格式的計(jì)算結(jié)果會(huì)存在很大偏差，而且隨著網(wǎng)格的加密，其向真值的收斂速度很慢.簡(jiǎn)言之，對(duì)于非均勻性很強(qiáng)的介質(zhì)，傳統(tǒng)格式如果沒有進(jìn)行足夠地加密，計(jì)算結(jié)果會(huì)導(dǎo)致很大的誤差［8］.

實(shí)驗(yàn)研究方面，Dave研究了滲透率或潤(rùn)濕性存在差異的2×2方型多孔區(qū)域中的多相流動(dòng)問題［9］.實(shí)驗(yàn)結(jié)果表明，即使?jié)B透率很小的差異也會(huì)導(dǎo)致流線很大的彎曲.流體從四個(gè)象限的中心角點(diǎn)處快速通過形成“竄流”.“竄流”現(xiàn)象導(dǎo)致流線彎曲，從而保持流動(dòng)在高滲透率區(qū)域之間進(jìn)行，并且在角點(diǎn)區(qū)域形成一個(gè)很大的壓力梯度.Dave認(rèn)為，如何用少量的網(wǎng)格模擬出竄流現(xiàn)象，對(duì)于傳統(tǒng)數(shù)值算法是一個(gè)很大的挑戰(zhàn).

在國(guó)內(nèi)，楊權(quán)一等提出了在滲透率間斷附近的自適應(yīng)加密算法［10］；段志田等提出了改進(jìn)的兩相滲流有限元數(shù)值算法［11］；林剛等利用薩曼斯基技巧改造牛頓迭代方法［12］；梁棟根據(jù)混溶和不混溶情況，分別提出了相應(yīng)的粘性分離算法和迎風(fēng)算法［13］；這些方法在一定程度上有效地提高了兩相滲流的計(jì)算效率和精度.

已經(jīng)知道，對(duì)于二維單相流動(dòng)，在不同滲透率區(qū)域交界的角點(diǎn)附近，壓力呈冪律分布，而且其梯度是發(fā)散的.這個(gè)特殊的規(guī)律是導(dǎo)致傳統(tǒng)格式計(jì)算效果差的原因.為了提高計(jì)算效率，在角點(diǎn)冪律解析解的基礎(chǔ)上，提出了針對(duì)單向流動(dòng)的有限分析算法［14－15］.該算法在2×2或3×3的網(wǎng)格細(xì)分條件下，就可以得到相當(dāng)準(zhǔn)確的計(jì)算結(jié)果，并且其準(zhǔn)確性和收斂速度不依賴于介質(zhì)的非均勻強(qiáng)度.本文針對(duì)二維不可壓縮兩相流動(dòng)，進(jìn)一步構(gòu)造相應(yīng)的有限分析算法.

1 二維不可壓兩相流動(dòng)中的有限分析算法

二維不可壓油水兩相流動(dòng)的控制方程可表述為

其中?表示孔隙度；sα表示相飽和度，下標(biāo)α表示水相或者油相；Vα表示各相的速度；k是多孔介質(zhì)的絕對(duì)滲透率，本文作為標(biāo)量處理；krα表示α相的相對(duì)滲透率，它是飽和度sα的函數(shù)；μα是α相的粘度系數(shù)；Pα是α相的壓力；毛管力Pc是兩相壓力之間的差值，也是飽和度的函數(shù)．將（1）式的前兩個(gè)方程相加，并利用sw＋so＝1可以得到

圖1 方程（2）的離散示意圖，陰影部分表示點(diǎn)P0的影響區(qū)域Fig.1 Sketch map of discretization of Eq.（2）on a square cell in 2D case，where the hatched area represents influenced area of grid node P0

本文重點(diǎn)闡述用有限分析方法［16］來數(shù)值求解壓力方程（2）.對(duì)控制容積（見圖1）應(yīng)用流量守恒，可以得到其中Q（pqi＋1/2，j＋1/2）表示穿過控制體相應(yīng)邊界的總流量（包括水相和油相），下標(biāo)pq表示如果以（i＋1/2，j＋1/2）為坐標(biāo)原點(diǎn)，流量是從第p象限流向第q象限．根據(jù)方程（2），（3）中的每個(gè)Q可以分成兩個(gè)部分：Qpw和Qp，分別對(duì)應(yīng)方程（2）中的－（krw/μw＋kro/μo）kΔPw和－krok/μoΔPc，因此方程（3）又可以改寫成

c

有限分析法的關(guān)鍵就是如何計(jì)算界面流量（Qp和Qp）．各界面流量由相應(yīng)的奇點(diǎn)（i＋1/2，j＋1/2）區(qū)域決

wc定，根據(jù)奇點(diǎn)區(qū)域性質(zhì)不同，分為兩種情況：①角域壓力呈冪律分布（ki，jki＋1，j＋1≠ki＋1，jki，j＋1）；②角域壓力呈線性分布（ki，jki＋1，j＋1＝ ki＋1，jki，j＋1）．

1.1 角域壓力呈冪律分布（ki，jki＋1，j＋1≠ki＋1，jki，j＋1）

如圖1陰影部分所示，絕對(duì)滲透率k在四個(gè)網(wǎng)格內(nèi)有不同的值ki，j，ki＋1，j，ki，j＋1和ki＋1，j＋1，和單相流動(dòng)一樣，壓力梯度在網(wǎng)格公共點(diǎn)P0附近是發(fā)散的［14－15］.該奇點(diǎn)鄰域（圖2中的小圓）的解析解是有限分析法的核心.對(duì)于方程（2），由于奇點(diǎn)附近的飽和度分布無法確定，所以找不到嚴(yán)格的奇點(diǎn)鄰域解析解.如果考慮到流動(dòng)的非定常性，情況將更加復(fù)雜.因此，首先要簡(jiǎn)化方程（2），以得到近似的奇點(diǎn)鄰域解析解.

記奇點(diǎn)處的水相飽和度值為s?，大多數(shù)情況下，水相飽和度在奇點(diǎn)附近是連續(xù)的；僅當(dāng)飽和度鋒面恰通過該點(diǎn)時(shí)才不連續(xù).事實(shí)上，飽和度鋒面會(huì)在迅速通過該奇點(diǎn)，對(duì)tn到tn＋1這個(gè)時(shí)段所造成的影響可以忽略．因此，總可以假定飽和度在奇點(diǎn)附近連續(xù)，其梯度是個(gè)有限值，故與發(fā)散的壓力梯度ΔPw相比，方程（2）中的kkro/μoΔPc可以忽略，故方程（2）可以簡(jiǎn)化為

在實(shí)際計(jì)算中，為保證物理上的正確性，離散的（krw/μw＋kro/μo）值總是采用“迎風(fēng)格式”，即定義在上游，故在角域，方程（5）中（krw/μw＋kro/μo）項(xiàng)可視為常數(shù).也就是說，可以把類拉普拉斯方程（6）在奇點(diǎn)鄰域的解析解作為方程（2）一個(gè)合理近似.

根據(jù)上文分析，在ki，jki＋1，j＋1≠ki＋1，jki，j＋1的情況下，方程（4）中的界面流量（Qpc）21，（Qpc）34，（Qpc）41和（Qpc）32可以忽略；界面流量Qpw和離散的水相壓力之間的關(guān)系，可以根據(jù)方程（6）的奇點(diǎn)鄰域解析解來建立.把ki＋1，j＋1，ki，j＋1，ki，j，ki＋1，j，（Pw）i＋1，j＋1，（Pw）i，j＋1，（Pw）i，j，（Pw）i＋1，j，si＋1，j＋1，si，j＋1，si，j，si＋1，j（見圖2）相應(yīng)地簡(jiǎn)記為k1，k2，k3，k4，P1，P2，P3，P4，s1，s2，s3，s4（見圖3）．界面流量Qpw和離散的水相壓力之間的關(guān)系可以寫成［11］

圖2 方網(wǎng)格下計(jì)算節(jié)點(diǎn)示意圖Fig.2 Plot of grid node in a rectangular grid system

圖3 受奇點(diǎn)P0控制的界面流量Qpw和Qpw示意圖Fige.3 Plot of interface fluxes Qpwand Qpcdominated by a grid node

上述公式的詳細(xì)推導(dǎo)參見文獻(xiàn)［10］或［11］.

1.2 角域壓力呈線性分布（ki，jki＋1，j＋1＝ki＋1，jki，j＋1）

ki，jki＋1，j＋1＝ki＋1，jki，j＋1的情況主要發(fā)生在網(wǎng)格加密過程中，這種情況下，壓力梯度ΔPw和毛管力梯度ΔPc在網(wǎng)格節(jié)點(diǎn)附近都是有限值．因此，－kkro/μoΔPc這部分不能夠忽略，方程（2）在此情況下的離散格式將退化到傳統(tǒng)形式［11］．方程（4）中的界面流量Qpw和Qpc可以寫成

其中（Pc）1＝Pc（s1），（Pc）2＝Pc（s2），（Pc）3＝Pc（s3），（Pc）4＝Pc（s4）．

至此，可以簡(jiǎn)述利用有限分析法求解二維不可壓縮兩相流動(dòng)方程（1）的步驟如下.這里，有限分析算法也可以看作是改進(jìn)的IMPES方法.

1）將tn時(shí)刻飽和度場(chǎng)的離散值作為tn＋1時(shí)刻的迭代初值，即：，其中表示tn時(shí)刻飽和度場(chǎng)的離散值，表示tn＋1時(shí)刻的第k次迭代值.

其中Qw表示穿過控制體邊界的水相流量，與－krwk/μwΔPw相對(duì)應(yīng)（見圖3），由于Qpw與－（krw/μw＋kro/μo）k ΔPw對(duì)應(yīng)，所以Qw和Qpw之間存在一個(gè)簡(jiǎn)單的關(guān)系：

事實(shí)上，Qpw在步驟2）中已經(jīng)求出，所以Qw可以直接得到．求解方程（33）可以得到tn＋1時(shí)刻第k＋1次的迭代值.

4）更新得到tn＋1時(shí)刻離散飽和度場(chǎng)的第k＋1次迭代值．重復(fù)步驟2）和3），直到得到收斂的．

2 數(shù)值算例

2.1 滲透率棋盤式分布算例

如圖4所示的二維棋盤式結(jié)構(gòu)在理論和數(shù)值上都引起研究者們的廣泛興趣，經(jīng)常被用來檢驗(yàn)一個(gè)數(shù)值算法的準(zhǔn)確性.本文考慮9×9的方網(wǎng)格系統(tǒng)，分別用有限分析法和IMPES方法計(jì)算不同的k2/k1情況下活塞型驅(qū)替算例.為了得到更準(zhǔn)確的結(jié)果，每一個(gè)原始網(wǎng)格被細(xì)分至n×n的計(jì)算網(wǎng)格.分別給出了k2/k1＝10，100，1 000情況下有限分析法和傳統(tǒng)IMPES方法的計(jì)算結(jié)果，以作比較.

設(shè)定棋盤式結(jié)構(gòu)中的無量綱絕對(duì)滲透率分別為k1＝0.1，1，10和k2＝100，對(duì)應(yīng)k2/k1＝10，100，1 000的情況；進(jìn)出口的無量綱水相壓力和無量綱飽和度設(shè)為：Pin＝3.0，sw in＝1.0和Pout＝1.0，swout＝0.3；相對(duì)滲透率率曲線取為：krw＝/（2－）和kro＝2（1－）/（2－）；無量綱粘度設(shè)為μw＝1和μo＝2．這種情況下，分流量函數(shù)f（sw）＝（krw/μw）/（krw/μw＋kro/μo） ＝．此時(shí)f″（sw）≥0，所以流動(dòng)屬于活塞式驅(qū)替［17］.

圖4 二維棋盤式分布示意圖Fig.4 A sketchmap of 2D checkerboard problem

圖5～7給出了k2/k1＝10、k2/k1＝100及k2/k1＝1 000情況下有限分析法和傳統(tǒng)IMPES方法在不同的加密參數(shù)n下，計(jì)算出的飽和度場(chǎng)分布.可以看出，k2/k1＝10情況下，非均勻性不是很強(qiáng)，有限分析法和IMPES方法都收斂的較快，有限分析法的鋒面移動(dòng)速度略快于IMPES；隨著非均勻性增強(qiáng)，k2/k1＝100和1 000時(shí)，隨著網(wǎng)格細(xì)分參數(shù)n的增加，有限分析法收斂的速度比傳統(tǒng)IMPES方法要快得多，有限分析法的計(jì)算結(jié)果在n＝8和n＝16時(shí)已基本沒有區(qū)別；而對(duì)于傳統(tǒng)IMPES方法，n＝8和n＝16的計(jì)算結(jié)果都非常不準(zhǔn)確，與有限分析法相比，傳統(tǒng)IMEPS方法嚴(yán)重低估了鋒面的移動(dòng)速度.為了能達(dá)到和有限分析法相近的準(zhǔn)確度，IMPES方法需要將計(jì)算網(wǎng)格劃分得非常密，從而需要很大的計(jì)算成本，以至不可行.換言之，有限分析方法則可以在較粗的計(jì)算網(wǎng)格條件下獲得很高的準(zhǔn)確度，而且不依賴介質(zhì)的非均勻強(qiáng)度.

圖5 k2/k1＝10情況下FAM和傳統(tǒng)IMPES方法在不同的網(wǎng)格加密參數(shù)n下計(jì)算的水相飽和度場(chǎng)Fig.5 Saturation fields of a 2D checkerboard calculated with FAM and traditional IMPES under different grid refinement parameter n as k2/k1＝10

圖6 k2/k1＝100情況下FAM和傳統(tǒng)IMPES方法在不同的網(wǎng)格加密參數(shù)n下計(jì)算的水相飽和度Fig.6 Saturation fields of a 2D checkerboard calculated with FAM and traditional IMPES under different grid refinement parameter n as k2/k1＝100

圖7 k2/k1＝1 000情況下FAM和傳統(tǒng)IMPES方法在不同的網(wǎng)格加密參數(shù)n下計(jì)算的水相飽和度Fig.7 Saturation fields of a 2D checkerboard calculated with FAM and traditional IMPES under different grid refinement parameter n as k2/k1＝1 000

圖8給出了k2/k1＝10，100，1 000的情況下有限分析法和傳統(tǒng)IMPES方法計(jì)算得到的出口邊界見水時(shí)間.這里定義出口邊界見水時(shí)間為出口邊界任意位置出現(xiàn)水相飽和度大于0.35的情況.如圖8所示，隨著介質(zhì)非均勻性的增強(qiáng)，在不同的網(wǎng)格加密條件下，有限分析法計(jì)算的見水時(shí)間幾乎是不變的，而對(duì)傳統(tǒng)IMPES方法，隨著網(wǎng)格的加密，計(jì)算的見水時(shí)間會(huì)逐漸減短，而且隨著介質(zhì)非均勻性增強(qiáng)，收斂越來越困難，甚至無法收斂，這也就意味著傳統(tǒng)IMPES方法在粗網(wǎng)格情況下會(huì)嚴(yán)重低估見水時(shí)間.

圖8 不同k2/k1情況下FAM和傳統(tǒng)IMPES方法計(jì)算的無量綱見水時(shí)間隨著網(wǎng)格加密參數(shù)n的變化Fig.8 Water breakthrough time of a 2D checkerboard calculated with FAM and traditional IMPESunder different grid refinement parameter n as k2/k1＝10，100，1 000

2.2 滲透率對(duì)數(shù)正態(tài)分布（Log?normal）算例

地層中的滲透率常滿足對(duì)數(shù)正態(tài)分布，它的概率密度函數(shù)為

圖9 Log?normal分布情況下FAM和傳統(tǒng)IMPES方法在不同的網(wǎng)格加密參數(shù)n下計(jì)算的水相飽和度場(chǎng)Fig.9 Saturation fields calculated with FAM and traditional IMPESmethod under different grid refinement parameter n for the case of lognormally distributed permeability

從圖中可以看出，滲透率對(duì)數(shù)正態(tài)分布情況下，傳統(tǒng)IMPES方法依然嚴(yán)重低估了見水時(shí)間，且隨著網(wǎng)格加密，計(jì)算值收斂得很慢；而有限分析法即使在粗網(wǎng)格情況下也能計(jì)算出相對(duì)準(zhǔn)確的飽和度場(chǎng)，且見水時(shí)間在粗網(wǎng)格條件下也收斂.這個(gè)算例再一次證實(shí)了有限分析算法的高效性.

3 結(jié)論和討論

針對(duì)二維非均勻多孔介質(zhì)中的不可壓兩相流動(dòng)，提出一種有限分析算法.當(dāng)四個(gè)相鄰網(wǎng)格的絕對(duì)滲透率滿足ki，jki＋1，j＋1≠ki＋1，jki，j＋1時(shí)，水相和油相壓力梯度在網(wǎng)格奇點(diǎn)處是發(fā)散的，而毛管力梯度是有限值.因此毛管力梯度在奇點(diǎn)鄰域可以忽略，類拉普拉斯方程（6）在奇點(diǎn)鄰域的冪律解析解可以作為兩相流動(dòng)的一個(gè)近似解，并可以基于此構(gòu)造有限分析格式.而對(duì)于ki，jki＋1，j＋1＝ki＋1，jki，j＋1的情況，相壓力梯度和毛管力梯度都是有限值，這時(shí)的數(shù)值格式自動(dòng)退化為傳統(tǒng)形式.

圖10 Log?noram l情況下FAM和傳統(tǒng)IMPES方法計(jì)算的無量綱見水時(shí)間隨著網(wǎng)格加密參數(shù)n的變化Fig.10 Water breakthrough time calculated with FAM and traditional IMPESmethod under different grid refinement parameter n for the case of lognormally distributed permeability

和傳統(tǒng)數(shù)值算法相比，有限分析法在網(wǎng)格加密過程中，有更快的收斂速度，且不依賴于介質(zhì)的非均勻強(qiáng)度.相反的，如果用傳統(tǒng)的數(shù)值算法去求解，為得到相對(duì)準(zhǔn)確的結(jié)果，往往需要對(duì)原始網(wǎng)格進(jìn)行大幅度加密.

另一個(gè)問題是關(guān)于多相流粗化.從本文的數(shù)值模擬結(jié)果來看，鋒面幾乎處處存在，飽和度場(chǎng)出現(xiàn)波動(dòng)式分布（見圖5和8）.此時(shí)，如果將9×9的網(wǎng)格合并成一個(gè)粗網(wǎng)格，并定義一個(gè)等效飽和度值來替代該區(qū)域內(nèi)的波動(dòng)飽和度場(chǎng)，這種做法可能并不合適.換言之，如果用等效飽和度值來計(jì)算粗化后網(wǎng)格間的各相流量，可能會(huì)帶來相當(dāng)大的誤差.

［1］ Allen M B.Numericalmodeling ofmultiphase flow in porousmedia［J］.Adv Water Resource，1995，8：162－187.

［2］ Aziz K，Setttari A.Petroleum reservoir simulation［M］.Applied Science，1979.

［3］ Breitenbach E A，Thurnau D H，van Poolen H K.Solution of the immiscible fluid flow simulation equation［J］.Society of Petroleum Engineers Journal，1969：155－169.

［4］ Chen C J，Chen H C.Finite analytic numericalmethod for unsteady two?dimensional Navier?Stokes equations［J］.JComput Phys，1984，53：209－226.

［5］ Coats K H，Nielsen R L，Terhune M H，Weber A G.Simulation of three?dimensional two?phase flow in oil and gas reservoirs ［J］.Society of Petroleum Engineers Journal，1967：377－388.

［6］ Cordazzo J，Maliska C R，Romeu R K.A Considerations about the internodal permeability evaluation in reservoir simulation ［C］.The 2nd Brazilian Congress on R＆D in Petroleum and Gas，Rio de Janeiro，June，15－18，2003.

［7］ Cordazzo J，Hurtado F SV，Maliska C R，Da Silva A F C.Numerical techniques for solving partial computationalmethods in engineering［C］.2003.

［8］ Dawe R A，Grattoni C A.Experimental displacement patterns in a 2×2 quadrant block with permeability and wettability heterogeneities—problems for numericalmodeling［J］.Trans Porous Media，2008，71：5－22.

［9］ Douglas J J，Peaceman D W，Jr Rachford H H.A method for calculating multidimensional immiscible displacement［J］. Metallurgical and Petroleum Engineers，1959，216：297－306.

［10］ Yang Q Y，Liu F P，Yang C C，Liu L F，Zhuo L.Adaptive nonuniform grid upscalingmethod of three?dimensional transient heterogeneous fluids in porousmedia［J］.Chinese JComput Phys，2003，20：193－198.

［11］ Duan Z T，Zhang D M.Finite element algorithms for simulating water flow in variably saturated porousmedia［J］.Chinese J Comput Phys，1992，9：15－22.

［12］ Lin G，Shi JM，Lin C B，Lv T，Lin A M.Gas reservoirs simulation and numerical comparison of results［J］.JComput Phys，1991，8：286－278.

［13］ Liang D.New numericalmethods and their theoretical analysis for the two phase displacement problems［J］.JComput Phys，1992，9：286－525.

［14］ Liu Z F，Wang X H.Finite analytic numericalmethod for two?dimensional fluid flow in heterogeneous porous media［J］.J Comput Phys，2013，235：286－301.

［15］ Qu Z X，Liu Z F，Wang X H，Zhao P.Finite analytic numericalmethod for solving two?dimensional quasi?Laplace equation ［J］.Numerical Methods for Partial Differential Equations，2014，DOI 10.1002/num.21863.

［16］ Romeu R K，Noetinger B.Calculation of internodal transmissibilities in finite differencemodels of flow in heterogeneous porous media［J］.Water Resour Res，1995，31：943－959.

［17］ Wu B，Liu Z F，Wang X H.A study on the numerical algorithm for the non?piston?like displacement in oil?water two?phase flows［J］.Journal On Numerical Methods and Computer Applications，2012，33：274－282.

Finite Analytic Numerical M ethod for Two?Phase Incom pressible Flow in 2D Heterogeneous Porous M edia

ZHENG Xiaolei，LIU Zhifeng，WANG Xiaohong，SHIAnfeng
（Department of Thermal Science and Energy Engineering，University ofScience and Technology of China，Hefei，Anhui 230026，China）

A finite analytic numerical scheme is constructed for two?dimensional two?phase flow in heterogeneous porous media. Compared with traditional numerical methods，F(xiàn)AM makes convergence faster as refinement parameter increases，and accuracy is independentwith heterogeneity.In contrast，as using traditional numerical schemes to simulate flow through a strong heterogeneous porousmedium，refinement ratio for grid cell needs to increase dramatically to get an accurate result.Compared with the proposed scheme，traditional numerical scheme underestimates greatly breakthrough time under coarse grids.However，to get a saturation distribution with high resolution even employing the proposed scheme，relative fine grids are still needed.This is different from FAM for solving a single phase flow，where coarse grids provide rather accurate results.

finite analytic method；multi?phase flow；fluid flows in porousmedia；heterogeneousmedia；multi?scale simulation

O362

A

2014－12－04；

2015－01－22

國(guó)家自然科學(xué)基金（11172295，11202205），中國(guó)科學(xué)技術(shù)大學(xué)青年創(chuàng)新基金（WK2090130017））及CNPC?CAS科技合作資助項(xiàng)目

鄭曉磊（1987－），男，安徽六安，博士生，從事油藏?cái)?shù)值模擬研究，E?mail：lzf123＠ustc.edu.cn

Received date： 2014－12－04；Revised date： 2015－01－22