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

（中國科學技術大學熱科學和能源工程系，安徽合肥 230026）

文章編號：1001?246X（2015）05?0586?09

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

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

（中國科學技術大學熱科學和能源工程系，安徽合肥 230026）

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

有限分析算法；多相流動；多孔介質；非均勻介質；多尺度模擬

0 引言

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

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

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

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

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

1 二維不可壓兩相流動中的有限分析算法

二維不可壓油水兩相流動的控制方程可表述為

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

圖1 方程（2）的離散示意圖，陰影部分表示點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

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

c

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

wc定，根據(jù)奇點區(qū)域性質不同，分為兩種情況：①角域壓力呈冪律分布（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陰影部分所示，絕對滲透率k在四個網(wǎng)格內(nèi)有不同的值ki，j，ki＋1，j，ki，j＋1和ki＋1，j＋1，和單相流動一樣，壓力梯度在網(wǎng)格公共點P0附近是發(fā)散的［14－15］.該奇點鄰域（圖2中的小圓）的解析解是有限分析法的核心.對于方程（2），由于奇點附近的飽和度分布無法確定，所以找不到嚴格的奇點鄰域解析解.如果考慮到流動的非定常性，情況將更加復雜.因此，首先要簡化方程（2），以得到近似的奇點鄰域解析解.

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

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

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

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

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

上述公式的詳細推導參見文獻［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é)點附近都是有限值．因此，－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）．

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

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

其中Qw表示穿過控制體邊界的水相流量，與－krwk/μwΔPw相對應（見圖3），由于Qpw與－（krw/μw＋kro/μo）k ΔPw對應，所以Qw和Qpw之間存在一個簡單的關系：

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

4）更新得到tn＋1時刻離散飽和度場的第k＋1次迭代值．重復步驟2）和3），直到得到收斂的．

2 數(shù)值算例

2.1 滲透率棋盤式分布算例

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

設定棋盤式結構中的無量綱絕對滲透率分別為k1＝0.1，1，10和k2＝100，對應k2/k1＝10，100，1 000的情況；進出口的無量綱水相壓力和無量綱飽和度設為：Pin＝3.0，sw in＝1.0和Pout＝1.0，swout＝0.3；相對滲透率率曲線取為：krw＝/（2－）和kro＝2（1－）/（2－）；無量綱粘度設為μw＝1和μo＝2．這種情況下，分流量函數(shù)f（sw）＝（krw/μw）/（krw/μw＋kro/μo） ＝．此時f″（sw）≥0，所以流動屬于活塞式驅替［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下，計算出的飽和度場分布.可以看出，k2/k1＝10情況下，非均勻性不是很強，有限分析法和IMPES方法都收斂的較快，有限分析法的鋒面移動速度略快于IMPES；隨著非均勻性增強，k2/k1＝100和1 000時，隨著網(wǎng)格細分參數(shù)n的增加，有限分析法收斂的速度比傳統(tǒng)IMPES方法要快得多，有限分析法的計算結果在n＝8和n＝16時已基本沒有區(qū)別；而對于傳統(tǒng)IMPES方法，n＝8和n＝16的計算結果都非常不準確，與有限分析法相比，傳統(tǒng)IMEPS方法嚴重低估了鋒面的移動速度.為了能達到和有限分析法相近的準確度，IMPES方法需要將計算網(wǎng)格劃分得非常密，從而需要很大的計算成本，以至不可行.換言之，有限分析方法則可以在較粗的計算網(wǎng)格條件下獲得很高的準確度，而且不依賴介質的非均勻強度.

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

圖8 不同k2/k1情況下FAM和傳統(tǒng)IMPES方法計算的無量綱見水時間隨著網(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 滲透率對數(shù)正態(tài)分布（Log?normal）算例

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

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

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

3 結論和討論

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

圖10 Log?noram l情況下FAM和傳統(tǒng)IMPES方法計算的無量綱見水時間隨著網(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)格加密過程中，有更快的收斂速度，且不依賴于介質的非均勻強度.相反的，如果用傳統(tǒng)的數(shù)值算法去求解，為得到相對準確的結果，往往需要對原始網(wǎng)格進行大幅度加密.

另一個問題是關于多相流粗化.從本文的數(shù)值模擬結果來看，鋒面幾乎處處存在，飽和度場出現(xiàn)波動式分布（見圖5和8）.此時，如果將9×9的網(wǎng)格合并成一個粗網(wǎng)格，并定義一個等效飽和度值來替代該區(qū)域內(nèi)的波動飽和度場，這種做法可能并不合適.換言之，如果用等效飽和度值來計算粗化后網(wǎ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

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

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

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