基于Matlab PDETOOL的滲流力學可視化教學與上機實驗

蘇關東， 顧 勛， 趙蘭苓， 孫德永， 韓貝宇， 陳林輝

(1.中國石油大學(北京) a.石油工程學院； b.理學院,北京 102249；2.新疆醫科大學 基礎醫學院,烏魯木齊 830011)

基于Matlab PDETOOL的滲流力學可視化教學與上機實驗

蘇關東1a， 顧 勛1a， 趙蘭苓1b， 孫德永1a， 韓貝宇1a， 陳林輝2

(1.中國石油大學(北京) a.石油工程學院； b.理學院,北京 102249；2.新疆醫科大學 基礎醫學院,烏魯木齊 830011)

為提高滲流力學的課堂效率，針對滲流力學學科自身的特點及其在教學實踐中普遍存在的問題，提出了將MatlabPDE工具箱引入課堂的教學方法：介紹了利用Matlab中的PDE工具箱(PDETOOL)對滲流力學中經典數學模型進行求解的方法和步驟，歸納了常見的油氣滲流數學模型和邊界條件，并利用PDETOOL對滲流力學中一源一匯、兩源兩匯、穩定滲流、彈性不穩定滲流等經典的模型進行了數值模擬，實現了各經典滲流模型壓力場、流場的可視化，有利于幫助學生理解流函數、壓力波、導壓系數等抽象概念，加深學生對壓力波傳播，不穩定滲流壓力衰竭等抽象過程的感性認識，提升學生上課的興趣和熱情；同時提出了在滲流力學教學過程中增設上機實驗內容的教學建議，充分發揮可視化教學的優勢，為數值模擬、油藏工程等課程打基礎，使學生對石油工程整個專業的知識體系產生宏觀的認識，提高教學質量。

PDE工具箱； 滲流力學； 偏微分方程； 數值模擬； 可視化； 教學

0 引 言

滲流力學(fluid flow in porous media)是研究流體在多孔介質中運動規律的學科，是流體力學的一個重要分支[1], 主要研究油氣水等在巖石這種多孔介質中的流動規律以及血液在毛細血管及微毛細血管網絡中的流動規律。石油工程專業，滲流力學課程屬于專業主干課和專業基礎課，這門課程主要研究地下油、氣、水及其混合物在地層中的流動規律[1-2]，它在石油工程專業課程體系中處于承上啟下的重要地位，類似于材料力學在土木工程專業知識體系中的地位，它以高等數學、油層物理等課程為基礎，同時又為油藏工程、采油工程、數值模擬等課程提供理論基礎和依據。

滲流力學對石油工程專業知識體系而言，其重要性毋庸置疑，但是它的教學難度也不可忽視：從研究對象來看，由于滲流力學所研究的油氣水賦存于地下，其運動狀態無法觀察，使得學生缺乏對滲流的形象認識；從研究手段來看，滲流力學主要以歐拉方法(場論)來描述地下油氣水的運動，以微分方程為數學工具對地下油氣水的流動進行研究，因為微分方程理論本身就具有一定的學習難度，使得滲流力學的學習難上加難，以至于滲流力學課程具有理論性強、邏輯性嚴密、公式多等特點，在教學實踐中出現了教師教學難、學生學習效果差等教學現狀[3]。

PDETOOL是Matlab中一個基于有限元方法求解偏微分方程的工具箱，由于其具有操作簡便，可視化效果好，而被廣泛用于熱傳導[4-5]、油氣管道集輸[6]、煤層氣擴散[7]、電磁場[8]、應力場[8-11]、地下水滲流[12]等領域的數值模擬。本文將其運用到滲流力學的教學中，不僅可有針對性地從教學方法和教學手段上解決目前滲流力學教學實踐中普遍存在的問題，同時還可以開闊學生在本專業的知識面，培養研究型人才。

在教學形式上， 教師不僅可以在課堂上展示可視化仿真結果，還可以根據課時安排，以上機實驗或者課后大作業的形式進行仿真可視化教學，增加學生對仿真過程的體驗，以期達到最好的教學效果。

1 基于PDETOOL的偏微分方程的解法

偏微分方程數學模型包括：四類偏微分方程模型和兩類邊界條件。這四類偏微分方程分別是橢圓型(Elliptic)、拋物型(Parabolic)、雙曲型(Hyperbolic)，本征值方程(Eigenmodes)，其中雙曲型偏微分方程和拋物型偏微分方程涉及到時間變量，可用于動態的數值模擬；兩類邊界條件包括：狄里克雷邊界條件(Dirichlet)和諾曼邊界條件(Neumann)，偏微分方程模型及其邊界條件的基本形式歸納如表1所示[6,13-15]

表1 偏微分方程基本類型及邊界條件

解偏微分方程的方法主要有3類[4]，包括：有限差分法、有限元法和譜方法。其中，使用有限元法時可以直接利用MATLAB中的PDE工具箱(PDETOOL)進行求解，操作簡單易學，可視化效果好，運用在教學中有利于幫助學生理解各種抽象的滲流場，激發學生們的學習興趣，提高教學質量。

2 油氣滲流的數學模型

流體在多孔介質中的滲流按照流體的壓力場是否隨時間變化可以分為兩大類：穩定滲流和非穩定滲流，分別對應不同的微分方程數學模型，它們的基本形式和在Matlab PDETOOL中的參數設置及模型的適用條件歸納如表2所示[16-17]。表2中p,K,μ,Ct,t分別為壓力、滲透率、黏度、綜合彈性系數和時間。對于具體的問題還需要知道邊界條件和初始條件才可以得到偏微分方程的解。

油藏的邊界條件分為兩類：外邊界條件和內邊界條件，外邊界條件指的是處于油藏外部，一般由地質作用形成的邊界所滿足的水力學條件，由油藏的成藏過程決定，主要分為供給邊界和封閉邊界；而內邊界條件指的是在油藏內部、一般由人工鉆井而形成的邊界所滿足的水動力學條件，由生產制度決定，可分為定壓邊界條件和定產邊界條件。外邊界的數學表達形式以及其對應的地質條件和內邊界的數學表達形式及其對應的生產制度如表3所示。表3中pe,pw,A,Q分別為供給壓力，井底流壓，滲流截面積和產量。

表2 滲流力學的數學模型

表3 滲流力學數學模型的邊界條件

根據滲流力學的基本模型，結合Matlab PDE工具箱的功能特性，并借鑒文獻[6-12]所用的相關方法，設計流程框圖，得到利用PDE工具箱解滲流數學模型的流程圖，如圖1所示。

圖1 PED工具箱解滲流數學模型流程圖

3 算 例

算例部分通過3個典型的案例，進一步介紹使用PDETOOL求解二位滲流場的方法，同時為教師課堂教學實踐提供教案參考。這3個經典案例分別是：一源一匯、兩匯、水平井開發油藏滲流模型，無解析解的滲流模型，五點井網法開發油藏動態模型。

3.1 經典穩定滲流場數值模擬

算例1 如圖3所示，一均質地層供給邊界為半徑1 000 m的圓(可近似看作無限大地層)，在距離圓心400 m處，以圓的直徑為對稱軸，打兩口井，分別求出兩口井都是生產井(兩匯)時，以及一注一采(一源一匯)時的壓力分布。用水平井段為1 000 m的水平井生產壓力分布又如何？

解 用一源一匯、兩匯、水平井3種方式生產油藏的壓力分布圖分別如圖2～4所示。

由圖2(b)可見：在注水井和生產井處分別產生一個高壓峰和低壓谷，對應圖2中流線由注水井指向生產井；圖2(a)中，在近井地帶，等壓線較密，壓力梯度大，與圖2(b)對應，可見近井地帶的壓力漏斗，根據達西定律可推斷兩井中心的連線流量最大，為主流線所在的位置，對應圖2(b)中明顯可見：兩井連線的中軸線上的壓力分布與外邊界的供給壓力一致，因而可等效為供給邊界；在油藏的外圍等壓線較疏，壓力梯度較小，對應圖2(b)中的壓力坡降平緩，甚至幾乎沒有明顯變化。

(a)

(b)

由圖3(b)可知，與以一源一匯的方式相比，當以兩匯的方式生產時，同樣可見壓力漏斗，不同的是出現了兩個低壓谷，同時在兩井連線的中軸線上，出現了一個極值點(高壓中心)，此時在該點附近的流體由于壓力相對較低，勢能相對較小而無法跨越該點，從一口生產井流向另一口生產井，因而使該中軸線成為分流線，等效于斷層，對應圖3(a)中的中軸線附近的流線只沿著中軸線流動而無法穿越中軸線，而處于該點的流體也因此無法被采出而形成死油點。

同時，通過圖2、圖3的對比可以看出：在一源一匯、一注一采情況下，壓力場是以兩井連線的中垂線為中心對稱軸，中心對稱分布的；而在兩匯的情況下，壓力場是以兩井連線的中垂線為對稱軸，軸對稱分布的。這也是在滲流力學中，用鏡像原理來求解滲流數學模型的事實依據。

由圖4與圖3對比可知，采用水平井生產時，井壁與油藏的接觸面積增大，藍色部分的低壓中心區域面積明顯擴大，即泄油面積擴大，從而增大了單井產量；同時，低壓區的擴大也可以解釋當水平井作為注入井時，相對于直井具有較低的注入壓力的原因。

(a)

(b)

(a)

(b)

3.2 無解析解的穩定滲流場數值模擬

算例2 如圖5(a)所示，在一個水平無限大的均質等厚地層中， 直線斷層邊界和直線供給邊界夾角約為60°，對稱于兩個邊界位置上有一口生產井，生產井距離直線斷層的距離約為200 m，泄油區供給半徑re為5 km，井半徑rw為0.1 m，供給邊界上壓力pe為20.0 MPa，井底流壓pw為4.0 MPa，地層滲透率K=0.8 μm2,地層厚度h=10 m，液體黏度μ=4 mPa·s，求地層壓力場及流場分布。

(a) (b)

圖5 油藏地質模型及其鏡像原理解法示意圖

解法1(鏡像原理)

根據鏡像原理：以封閉邊界(斷層)為對稱軸作的井為同號井，以供給邊界為對稱軸作的井為異號井(正號井為生產井，負號井為注水井)，按順時針方向依次做鏡像井，結果如圖5(b)所示，1號井(生產井)按順時針鏡像得到的2號井為注水井，而按照順時針鏡像得到的2號井應該為生產井，自相矛盾，說明鏡像原理——用于求解壓力分布解析解的方法在這種情況下不適用，即該微分方程在此邊界條件下沒有解析解。

通過本例，教師即可以讓學生認識到解析方法在實際生產應用中的不足，同時還可以順勢簡單地介紹一下油藏數值模擬這門課程的基本任務和基本方法以及常用的油藏數值模擬軟件例如：Eclipse，CMG等，既提高了學生的學習興趣，又能使課堂的內容更貼近生產實際，還可以使學生在學習油藏數值模擬這門課程時更加容易接受；在教學安排上，教師既可以將本例在課堂上用解法2演示一遍，也可以提供基本的思路將本題布置成作業，供學生課后探究。

解法2(用Matlab的PDETOOL進行求解)

按照題中給出的已知條件設置PDE參數，并將得到的數值解可視化如圖6所示。由圖6(a)可見，供給邊界與流線垂直，與等壓線和等勢線相切(平行)，說明供給邊界為流體提供地下滲流的能量，“推”著流體向井口運動；而封閉邊界則與流線相切(平行)，與等勢線等壓線垂直，說明流體只能夠貼著封閉邊界運動而不能穿過封閉邊界；同時，封閉邊界垂直截斷等勢線也說明了壓力波無法穿透封閉邊界，這也是流體無法穿過封閉邊界的動力學原因。

(a)

(b)

近井地帶的等壓線比較緊密，遠井地帶的等壓線比較稀疏，說明在生產井附近，壓力下降較快，生產井對近井地帶的壓力場影響較大，對應圖6(b)中可見壓力漏斗。

3.3 不穩定滲流場數值模擬

算例3(五點井網法) 如圖7(a)所示，有一足夠大水平均質油藏(可視為無限的地層)，以五點井網法對其進行開發，注水井壓力pe=10.0 MPa，生產井井底流壓為pw=7.5 MPa，油井半徑rw=0.1 m，地層厚度h=5 m，滲透率K=1 μm2，黏度μ=4 mPa·s，注水開發周期為T，求該過程中地層壓力場和速度場分布？

解 為了對油藏滲流模型進行求解，從井網中選取一個小單元[18]，如圖7(b)所示。

圖7 五點井網及其單元示意圖

由圖8可知，隨著時間的推移，,注水井處的壓力波緩慢向采油井方向傳遞，高壓中心擴散；采油井附近，隨著地下油的采出對油藏泄壓，低壓中心向四周擴散。隨著時間增長，壓力波推進的速度越來越慢，最后幾乎靜止，達到準靜態(擬穩態)。在最終驅替完畢時，在兩井中軸線中心的區域壓力值不變，這就是賦存剩余油的死油區(見圖8(c))，此時應當調整注采關系(注采井網)以提高波及效率，從而提高采收率。

(c) t=T

圖8 五點井網單元開發動態圖

本例很好地展示了油藏開發從非穩態階段到擬穩態階段的過程，同時又和油藏工程數值模擬緊密聯系，有利于學生理解油氣田開發的動態過程。

特別說明：由于PDETOOL本身的不完善性，在教學實踐中大可不必嚴格拘泥于物理模型和數值模型的幾何比例關系，可以為了達到較好的可視化效果而對其進行適當調整。

4 結論與展望

(1) 歸納總結了用Matlab自帶的PDE工具箱對偏微分方程的求解方法和基本流程以及常見的油氣滲流數學模型及邊界條件。

(2) 用PDETOOL對滲流力學中經典的一源一匯，兩匯，水平井等滲流模型以及五點井網法開發油藏的動態過程進行了數值模擬，有利于加深學生對滲流力學理論的理解，同時為教師課堂教學實踐提供教案參考。

(3) 提出了相關的教學建議：除了在課堂上演示滲流場的可視化過程外，還可以在滲流力學教學實踐中，根據具體情況，增設上機實驗的內容，既可以上機實驗課的形式進行， 也可以課后報告的形式進行，豐富課堂形式，使學生投入到可視化教學的體驗中，以充分發揮可視化教學的優勢。

(4) 將PDE工具箱引入滲流力學的教學，不僅有利于學生對抽象的滲流理論產生形象化的認識，提高教學質量；同時有利于學生提前了解油藏工程、數值模擬等課程，對石油工程專業知識體系有全面的認識，提高畢業生的質量。

致謝 感謝中國石油大學(北京)石油工程學院滲流力學本科教研組的支持與幫助。衷心感謝中國石油大學(北京)理學院趙蘭苓老師的指導以及在論文修改上的幫助。特別感謝中國石油大學(北京)石油工程學院龐占喜老師、涂斌老師和李春蘭老師給予的幫助。特別感謝中國石油大學(北京)石油工程學院隋微波老師、朱舟元老師和東曉虎老師提供的部分圖片。衷心感謝中國石油大學(北京)鄭力會教授、張鵬教授給予的鼓勵。

[1] 宋洪慶，姜琳婧，朱維耀，等. 滲流力學課程研究型教學探索[J]. 中國冶金教育, 2014(2): 5-7.

[2] 袁迎中，戚志林，嚴文德，等. “滲流力學”課程創新思維教學方法探索[J]. 中國地質教育, 2014(2): 60-62.

[3] 陳軍斌. “油氣滲流力學”課程教學改革的研究與實踐[J]. 石油教育, 2016(3): 32-35.

[4] 張哲任，徐 政，徐 韜，等. 采用有限體積法的特高壓直流輸電系統接地極穩態溫度場仿真分析[J]. 高電壓技術, 2012，38(2): 328-334.

[5] 葉長燊. MATLAB在穩態與動態導熱過程分析中的應用[J]. 計算機與應用化學,2006, 23(10): 986-990.

[6] 張純靜，申龍涉，杜義朋，等. MATLAB PDE工具箱在稠油輸送管道中的應用[J]. 當代化工, 2012, 41(1): 48-49.

[7] 馬 勝，虞青松. MATLAB在模擬煤層氣擴散中的應用[J]. 山西焦煤科技, 2012(1): 44-45.

[8] 黃作英，闕沛文. 基于MATLAB的電磁場數值分析[J]. 計算機工程與應用, 2004(36): 196-198.

[9] 李 彤，劉阿龍. 基于MATLAB的圓環承受均布壓力分析[J]. 機械工程與自動化,2008(2): 14-16.

[10] 潘東輝，馬崇武. MATLAB/PDE在彈性力學可視化教學中的應用[J]. 力學與實踐, 2014, 36(4): 500-504.

[11] 徐 勇，李順才，佟軍民. 基于MATLAB的平面應力問題研究[J]. 煤礦機械, 2003(10): 60-62.

[12] 牛 宏，梁 杏，尼勝楠，等. MATLAB PDE工具箱在地下水科學實驗教學中的應用[J]. 安全與環境工程, 2013, 20(6): 6-11.

[13] 馮桂蓮. 偏微分方程的MATLAB數值解法及可視化[J]. 計算機技術與發展, 2013, 23(12): 120-123.

[14] 李 明. 偏微分方程的MATLAB解法[J]. 湖南農機, 2010, 37(5): 89-91.

[15] 田 兵. 用MATLAB解偏微分方程[J]. 陰山學刊(自然科學版), 2006, 20(4): 12-13.

[16] 劉慧卿.油氣滲流力學基礎[M]. 北京:中國石油大學出版社, 2013:30-47.

[17] 程林松.滲流力學[M]. 北京:石油工業出版社, 2011:1-23.

[18] 董平川，韓德金，牛彥良，等. 油藏多相滲流的面向對象有限元程序設計[J]. 巖土力學, 2009, 30(4): 1115-1121.

Visual Teaching and Computer-based Experiment of Fluid Flow in Porous Media with PDETOOL of Matlab

SU Guandong1a， GU Xun1a， ZHAO Lanling1b， SUN Deyong1a， HAN Beiyu1a， CHEN Linhui2

(1a. College of Petroleum Engineering, 1b. College of Science, China University of Petroleum (Beijing), Beijing 102249 China; 2. School of Basic Medical Sciences，Xinjiang Medical University, Urumqi 830011 China)

In order to improve the classroom efficiency and the teaching quality of the course of fluid flow in porous media, the teaching method of introducing Matlab PDETOOL into the classroom is put forward, and the method and steps to solve the classical mathematical model in fluid flow in porous media with the PDE toolbox in Matlab are introduced, according to the characteristics of the subject of fluid flow in porous media and the problems in teaching practice. The PDETOOL is used to simulate the classical models of fluid flow in porous media such as source-sink, source-sink, steady-flow, and elastic-unsteady flow in porous media. It helps students understand the abstract concepts of stream function, pressure wave and conductivity coefficient, enhances students’ perceptual knowledge of the process of the propagation of pressure wave and unsteady flow in porous media. At the same time, it is suggested to add computer-based experiment to the teaching process of fluid flow in porous media. On conclusion, this study will lay a foundation for the course of numerical simulation and reservoir engineering, so that the students will have a macroscopic understanding of the whole knowledge system of petroleum engineering.

PDETOOL; fluid flow in porous media; partial differential equations; numerical simulation; visualization; teaching

2016-11-22

國家級大學生創新創業訓練計劃支持項目(201611414019);中國石油大學基礎研究基金(2462015YQ0604)

蘇關東(1995-)，男，廣東湛江人，本科生。

Tel.：18801322091;E-mail：su_gd@qq.com

趙蘭苓(1977-)，女，山東泰安人，博士，講師，研究方向：計算數學、復雜流體建模與計算。

Tel.：010-89731993；E-mail：zllblue@163.com

TE 312; O 242.1

A

1006-7167(2017)08-0137-06