1、流線數值模擬中的無效網格處理方法王宇1,李治平2（1陜西延長石油（集團）有限責任公司研究院，陜西 西安710075 ; 2中國地質大學（北京）， 北京 100083）摘要：對于砂體橫向連續(xù)性差，為呈透鏡狀，局部區(qū)域為局部斷層影響的水動力滯留區(qū)、或 泥巖的油區(qū)進行數值模擬時必然涉及到不可流動區(qū)域的處理，不可流動區(qū)域網格化后變?yōu)椴?可流動網格，又稱無效網格或死網格。不可流動區(qū)域的存在會影響流線數值模擬的壓力場計 算以及流線追蹤過程，導致流線數值模擬失敗，因此死網格的研究對流線數值模擬計算具有 中重要意義。對于差分過程不可流動網格的處理方式，給出了三種處理方法，分別是直接計 算法，鏡像反映法以及網格

2、設置法，并介紹流線追蹤、流線飽和度模型中涉及到的死網格處 理方法。通過例子計算分析，證明了不可流動網格處理方法的有效性。網格設置方法具有良 好的規(guī)律性，計算簡便，直接計算法的方程系數矩陣規(guī)律不好掌握，規(guī)律復雜，而鏡像反映 法要不斷的對死網格重新賦值，影響了計算速度。當模擬區(qū)存在不可流動區(qū)域時，無效網格 處理方法能夠有效的解決流線數值模擬過程出現的問題。關鍵詞：流線數值模擬；無效網格；流線追蹤；數值模擬文獻標識碼：A(3)0引言流線模擬技術能將三維模擬模型 轉化為一系列的一維流線模型，實踐表 明流線方法是一種快速有效的數值模 擬法計算方法，流線模擬技術能更好地 認識地下流體的分布、運移，改善油田

3、 開發(fā)效果和提高油田最終采收率1。在 實際油田開采過程中，油田內往往存在 無油非連接區(qū)域，有些局部區(qū)域為局部 斷層影響的水動力滯留區(qū)，或為泥巖， 此時在油藏數值模擬過程中對應為非 流動區(qū)域，網格化后即為死網格。死 網格的存在會影響影響流線數值模擬 的壓力場計算以及流線追蹤過程，導致 流線數值模擬失敗。另外，油藏形狀各式各樣，進行網 格劃分后往往是非規(guī)則形狀，形成統(tǒng)一 的差分過程及對角矩陣過程復雜，而將 油藏劃分為規(guī)則的矩形，對其內部流體 不可流動區(qū)域按死網格處理，形成流程 化計算，其計算過程更為高效，簡化了計算過程。因此，研究死網格在數值模 擬區(qū)域的處理方法對于數值模擬計算 具有重要意思。1壓

4、力數學模型與無效網格的處理1.1 壓力方程在二維兩相的前提下，考慮流體與 巖石的壓縮性，忽略重力及毛管壓力， 考慮維數因子和產量項，根據達西定 律、質量守恒定律以及流體狀態(tài)方程可 得到二維兩相流數學模型35: V（KLAvp）+qo =馬帆so）otKk P（p）+qw=（e PwSw）（2）wCt由方程（1）、（2）可得到壓力控制方 程（3）:P) (- -qw)二: 5三 .o . w t對方程（3）五點隱式差分，把引 起非線性的系數作顯式處理后得：(CO/ ：o q/ 1w)=.1如果把式子（4）寫成矩陣的形式,n 1 n 1n 1 n 1工(.2/21(2 3/2)1/2一 州改)n

5、1 n 1n 1 n 1,1/M、n、Pj f-Pj* 弗 n、Pj 1-Pj、(4)+(/T)j4/2- -Xi/J)jl/2 )xXj 1/2為 1/2即AU =g,其中A稱為該線性方程組的系數矩陣，則其矩陣結構為五對角形 式。如果在某一個網格處有井存在，則 把它作為點源或點匯來處理，在對網格 建立的差分方程中增加一個產量項；1.2 差分過程無效網格的處理方法對存在無效網格的壓力方程差分 時，有三種方法處理，第一種是直接計 算法，對研究區(qū)域網格化處理后，只對 有效網格節(jié)點差分，而與無效網格相連 的有效網格得到差分方程后，求解五對 角矩陣方程，此方法會因無效網格數及 所處模擬區(qū)域位置的變化而

6、需重進行 差分，重新計算五對角系數矩陣，這種 方法無法形成統(tǒng)一、有規(guī)律的計算方 法，尤其不適應程序化、軟件化數值模 擬。第二種方法是采用“鏡像反映法” 3 O將死網格作為周邊與其連接的可流 動網格的鏡像，使其壓力等于周邊網格 的壓力，此時鏡像網格在計算過程中無 流體發(fā)生滲流。如圖1所示，對死網格 i j點差值時后，采用鏡像反映法處理n 1 n 1 n 1 n 1 n 1 n 1 時，口仞 P+=P ,P=P , Pj4=Pj ，Pr = Pn + , P； + = Pn，此時在 i+1/2, iT/2,1/2, j-1/2網格邊界處，無流量通過, （i,j）網格內無流量變化，這時（i,j）網

7、格即為死網格。此方法盡管不需要隨著無效網格 的改變而重新計算方程的系數矩陣， 能 形成一個統(tǒng)一的算法，但需要每一個時 間步長內設置無效網格的壓力，使其與 周邊網格壓力相等，其計算過程繁雜， 模擬過程的不斷賦值，大大的降低了計 算的效率。第三種方法為網格設置法，可直接 設置死網格的4及Ki為零，即不可流 動網格，此時網格就成為死網格。因綜 合壓縮系數Cf =2+s,Cw +（1 -Sv）Co,此 時方程（3）右邊的綜合壓縮系數意義 已經不在適用，應以展開形式帶入方 程，以避免存在無意義情況。另將式（4） 中的人TT/2等系數，均按“上游權”的 取值原則。由于計算過程中網格的滲透 率取“上游權”原

8、則，但對于死網格則 不需要此判斷過程，只需要確定死網格 邊界滲透率的值為零；如圖1所示， （i, j）網格為死網格點，其邊界滲透率分別為 kri 書/2、kri -4/2、krj +/2 和 krj 4/2均為 零。由于不可流動網格的存在，對應網 格點的滲透率為零，差分后的五對角矩 陣行列式為零，無法求解方程組，因此 需要將矩陣中為零的行、列除去。壓力 方程求解完畢后，得到對應網格在 tn.1 時刻的壓力值，即壓力場I。此方法可 直接將不規(guī)則形狀油藏劃為矩形油藏進行計算，并形成了統(tǒng)一的系數矩陣計 算方法，其系數矩陣計算規(guī)律簡單，計 算速度快。2流線數學模型與無效網格的處理流線數值模擬方法作為新

9、興的數 值模擬技術方法，與傳統(tǒng)的數值模擬方 法相比有很多的優(yōu)良性，能為生產實踐 提供指導和技術支持。當模擬區(qū)塊外 部、內部存在不可流動區(qū)域時，無效網 格處理方法能夠有效的解決數值模擬 過程中的壓力求解、流線追蹤以及流線 飽和度計算等因存在無效網格而出現 的問題。2.1 流線追蹤過程與無效網格的處理若存在死網格的流線數值模擬中， 流線應避開死網格區(qū)域。Pollock提出 了追蹤流線軌跡的方法是在網格系統(tǒng) 中的壓力場已知的情況下，應用達西方 程建立流體真實流動速度場，然后再此 基礎上追蹤流線?？紤]二維情況，根據達西定律計算 每個網格界面上的速度分量7:(H/2,j) =- 卻2,j(p 為p,j)

10、/(x -x) (5)Vy(i,j 1/2) = -1i,j故(p-p)/(y-y) (6)如圖2所示，采用Pollock方法進其中穿過網格塊得速度梯度為：0，0,速度的定義為：Vx =dx/dt, .x對(4)式進行積分，可得到到達x方 向到達出口界面所需的時間% =1儼0 2-尹潮,同理可得y mxVx,o +m(x -xo)方向所需時間：&yin%m土印,myVyo4mXy-y)二者之中值小者即為流線到達出口界 面所需的時間,即：Ate = minte,x,&ey), 其所代表的方向即為流線出口的方向。在得到了流線出口方向后，出口的確切位置由 xe =x0Vx,i exp(mx&e)-V

11、x,0，mxye =VoVyi exp(my&e)-Vy0 瓶走。my對于無效網格，應沒有流體流入或 流出，網格邊界封閉，無效網格邊界不 存在流體速度，因此無效網格邊界速度 應取為零；為了避免速度為零而出現無意義 情況，可將死網格邊界速度在一定精度 下取值無窮小，這樣可模擬流線在遇到0.60.81行流線追蹤。x方向的速度定義為:Vx =Vx,0 mx(x-x0)死網格時的情況，此時流線流動到死網 格邊界的時間無窮大，即不可到達。完 整的流線追蹤過程可得到流線流動的完整路徑8。j + 1Vy,.y：2.2 流線飽和度模型與無效網格處理,0流線i-1要采用流線方法求解生產階段的 克數學模型，首先要

12、將基于網格的二 Lift模理轉化為沿流線的一維模i+1定義沿著流線的傳播時間(Time Of對滲透率對滲透率Vx,y入口Vy,0j-1Flight, TOF )為：i(s)=()djd j)|Vt( )|圖2流線穿過網格過程處理后最終得到： -0,即流線 :t飽和度方程9 0無效網格內不存在流體流動、變 化，因此流線不經過無效網格，飽和度 方程僅跟流線所經過的網格有關，因 此無效網格不影響其他網格飽和度的 計算。3應用分析本例中選取長慶定邊姬堀油田羅 1 井區(qū)長8油藏內的四口采油井，一口注 水井作為應用對象，根據所研究區(qū)域儲 層主要參數及儲層流體特性資料作為 流線數值模擬的幾本參數。網格系統(tǒng)采

13、用塊中心網格系統(tǒng)；縱 向不劃分網格，平面上不考慮儲層的非 均值性，將平面劃分為19M 19父1網格系 統(tǒng)。無效網格區(qū)域分布在三個區(qū)域。模 擬區(qū)域中央為注水井，4 口采油井分布 在儲層四邊，無效網格分布在下面三個 部分，分別為左下邊的單個無效網格、 兩個連接無效網格以及右下邊6個相連 的無效網格。根據前面的數學模型及第三種無 效網格處理方法，編寫了流線數值模擬 計算程序，模擬3個月后，得到了模擬 區(qū)域一注四采井網形式的流線分布特 征，如圖3所示。I 2 I 4 5 I T B S IU H 12 0 M IT II 194V圖3存在死網格流線分布從圖3中可以看出流線從中間一口 注水井出發(fā)，流向周

14、邊的四口生產，流 線碰到死網格區(qū)域會明顯的出現變化， 而三部分不可流區(qū)域無流線通過，流線 繞過不可流動區(qū)域，流向了注水井，流 線明顯改變了流線原有的形狀，證明了 可不流動網格處理的有效性及合理性。由于不可流動區(qū)域的存在，靠近不 可流動區(qū)域，流線分布較密，而流線分 布密集地區(qū)反映了注水流量大，驅油效 果好。4結論1）文中給出了壓力方程的三種死 網格處理方法中，第三種方法計算最為 簡便，方程系數矩陣算法最容易得出， 第一種方法無法形成統(tǒng)一、有規(guī)律的計 算方法，尤其不適應程序化、軟件化數 值模擬。而第二種方法計算過程繁雜， 模擬過程的不斷賦值，大大的降低了計 算的效率。而第三種方法彌補了前兩種 方法

15、的不足，即形成了統(tǒng)一的系數矩陣 計算方法，又有相對較快的計算速度。2 ）直接設置法可將不規(guī)則形狀油 藏按規(guī)則的矩形油藏進行數值模擬計 算，并形成了統(tǒng)一的系數矩陣計算方 法，簡化了模擬計算過程。3） Pollock方法進行追蹤流線軌跡 過程中，因死網格的存在，在處理網格 邊界，會影響流線變化，進而改變飽和 度的計算過程。4）另外應用無效格處理方法可解 決含啟動壓力梯度的情況下油藏動邊 界問題，因油水兩相方程動邊界很難給 出解析形式，所以需要應用數值方法求 取確定，通過小時間步長計算確定動邊 界范圍，再將動邊界外部確定為死網 格，這樣可解決因啟動壓力梯度的存 在，動邊界外部流體滲流速度為負的情 況

16、。5符號說明%=K（Kw+Ko）,為總流度， Nf2/（mPas）； K w o為油藏的絕對滲透率，o, w分別為油、水相粘度，Pa s； kr。，心分別為油、水相的相對滲透率，無因此；Krww 2,為水的分流量， ；,K rw4Cf =?Cr +SwCw +(1 _Sw)Co為油層綜合壓縮系數，1/MPa ； p為油藏壓力，MPa； qo,qw分別為油、水相在單位時間、單位體積巖石注入(或采出)的質量流量，kg/(m3 s)；P。，心分別為油、水相的密度,kg/m3；S。、Sw為油、水想飽和度，；金為孔隙度， ； t為時間，s。參考文獻1Batycky, R.P.A Three-dimens

17、ional Two-phase Scale Streamline Simulator , ph.D. Thesis, Department of Petroleum Engineering, School of Earth Science, Stanford University: 19972郭鳴黎，程東風，李大勇.文25東復雜斷 塊油藏剩余油分布研究.江漢石油學院學報， 20033韓大匡，陳欽雷，閆存章.油藏數值模擬 基礎.石油工業(yè)出版社，19934王洪寶，蘇振閣，陳忠云.油藏水驅開發(fā)三 維流線模型J.石油勘探與開發(fā)，2004, 31(2): 99-1035R.C.M.Portella,T.

18、A.Hewett. Fast 3-D Reservoir Simulation and Applications Using StreamlinesJ.SPE39061 6劉慧卿.油藏數值模擬方法專題.東營： 石油大學出版社，20017姚軍，吳明錄，戴衛(wèi)華等.流線數值試井解 釋模型J.石油學報,2006 , 27(3) : 96-99 8Batycky,Rob P. Martin J.Blunt,and Marco R.Thiele. A 3D Field Scale StreamlineSimulatorWith GravityandChanging Well Conditions.SPE3

19、6726, 1996 9RomanA.Berenblvum,AlexanderA.Shapiro,KristianJesson. Black OilStreamlineSimulator With CapillaryEffects. SPE 84037, 200310尹虎，王新海，劉洪等流考慮啟動壓力 梯度的頁巖氣藏數值模擬J.天然氣與石油,2012年04期11王平，姜瑞忠，王公昌.高含水期水驅 狀況影響因素數值模擬研究 J.天然氣與石 油,2012年04期作者簡介：王宇，男，陜西省西安市人，助 理工程師，碩士，主要從事油氣藏開發(fā)研究 工作，郵編 710075,聯系方式：陜西省延安市甘泉縣下寺

20、灣鎮(zhèn)下 寺灣采油廠， 716100,手, 身份證號碼,電郵 wangyu79667643.項目來源：國家重大科技專項課題“新一代 油藏數值模擬軟件” (2009ZX05009-006 )The Processing Method of Non-Flowable Grid in Streamline NumericalSimulationAbstract: It is inevitable involves the deal with of non-flowable region if there are poor continu

21、ity of the sand show as lentoid in lateral distribution , hydrodynamic hold-up area caused by fault or shale in local area. Non-flowable area can change to immobile grid that also know as dead grid after be grided. The presence of a non-flowable region can affect the calculation of pressure filed, streamline tracing of streamline numerical simulation and lead to failure. So it is important to study the non-flowable grid. There are three methods to handle non-flowable grid in differential process, they are direct calculation method, “mirror image” method