流线数值模拟非规则网格处理方法

1、流线数值模拟中的无效网格处理方法王宇1,李治平2（1陕西延长石油（集团）有限责任公司研究院，陕西 西安710075 ; 2中国地质大学（北京）， 北京 100083）摘要：对于砂体横向连续性差，为呈透镜状，局部区域为局部断层影响的水动力滞留区、或 泥岩的油区进行数值模拟时必然涉及到不可流动区域的处理，不可流动区域网格化后变为不 可流动网格，又称无效网格或死网格。不可流动区域的存在会影响流线数值模拟的压力场计 算以及流线追踪过程，导致流线数值模拟失败，因此死网格的研究对流线数值模拟计算具有 中重要意义。对于差分过程不可流动网格的处理方式，给出了三种处理方法，分别是直接计 算法，镜像反映法以及网格

2、设置法，并介绍流线追踪、流线饱和度模型中涉及到的死网格处 理方法。通过例子计算分析，证明了不可流动网格处理方法的有效性。网格设置方法具有良 好的规律性，计算简便，直接计算法的方程系数矩阵规律不好掌握，规律复杂，而镜像反映 法要不断的对死网格重新赋值，影响了计算速度。当模拟区存在不可流动区域时，无效网格 处理方法能够有效的解决流线数值模拟过程出现的问题。关键词：流线数值模拟；无效网格；流线追踪；数值模拟文献标识码：A(3)0引言流线模拟技术能将三维模拟模型 转化为一系列的一维流线模型，实践表 明流线方法是一种快速有效的数值模 拟法计算方法，流线模拟技术能更好地 认识地下流体的分布、运移，改善油田

3、 开发效果和提高油田最终采收率1。在 实际油田开采过程中，油田内往往存在 无油非连接区域，有些局部区域为局部 断层影响的水动力滞留区，或为泥岩， 此时在油藏数值模拟过程中对应为非 流动区域，网格化后即为死网格。死 网格的存在会影响影响流线数值模拟 的压力场计算以及流线追踪过程，导致 流线数值模拟失败。另外，油藏形状各式各样，进行网 格划分后往往是非规则形状，形成统一 的差分过程及对角矩阵过程复杂，而将 油藏划分为规则的矩形，对其内部流体 不可流动区域按死网格处理，形成流程 化计算，其计算过程更为高效，简化了计算过程。因此，研究死网格在数值模 拟区域的处理方法对于数值模拟计算 具有重要意思。1压

4、力数学模型与无效网格的处理1.1 压力方程在二维两相的前提下，考虑流体与 岩石的压缩性，忽略重力及毛管压力， 考虑维数因子和产量项，根据达西定 律、质量守恒定律以及流体状态方程可 得到二维两相流数学模型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 差分过程无效网格的处理方法对存在无效网格的压力方程差分 时，有三种方法处理，第一种是直接计 算法，对研究区域网格化处理后，只对 有效网格节点差分，而与无效网格相连 的有效网格得到差分方程后，求解五对 角矩阵方程，此方法会因无效网格数及 所处模拟区域位置的变化而

6、需重进行 差分，重新计算五对角系数矩阵，这种 方法无法形成统一、有规律的计算方 法，尤其不适应程序化、软件化数值模 拟。第二种方法是采用“镜像反映法” 3 O将死网格作为周边与其连接的可流 动网格的镜像，使其压力等于周边网格 的压力，此时镜像网格在计算过程中无 流体发生渗流。如图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、格即为死网格。此方法尽管不需要随着无效网格 的改变而重新计算方程的系数矩阵， 能 形成一个统一的算法，但需要每一个时 间步长内设置无效网格的压力，使其与 周边网格压力相等，其计算过程繁杂， 模拟过程的不断赋值，大大的降低了计 算的效率。第三种方法为网格设置法，可直接 设置死网格的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。此方法可 直接将不规则形状油藏划为矩形油藏进行计算，并形成了统一的系数矩阵计 算方法，其系数矩阵计算规律简单，计 算速度快。2流线数学模型与无效网格的处理流线数值模拟方法作为新

9、兴的数 值模拟技术方法，与传统的数值模拟方 法相比有很多的优良性，能为生产实践 提供指导和技术支持。当模拟区块外 部、内部存在不可流动区域时，无效网 格处理方法能够有效的解决数值模拟 过程中的压力求解、流线追踪以及流线 饱和度计算等因存在无效网格而出现 的问题。2.1 流线追踪过程与无效网格的处理若存在死网格的流线数值模拟中， 流线应避开死网格区域。Pollock提出 了追踪流线轨迹的方法是在网格系统 中的压力场已知的情况下，应用达西方 程建立流体真实流动速度场，然后再此 基础上追踪流线。考虑二维情况，根据达西定律计算 每个网格界面上的速度分量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 井区长8油藏内的四口采油井，一口注 水井作为应用对象，根据所研究区域储 层主要参数及储层流体特性资料作为 流线数值模拟的几本参数。网格系统采

13、用块中心网格系统；纵 向不划分网格，平面上不考虑储层的非 均值性，将平面划分为19M 19父1网格系 统。无效网格区域分布在三个区域。模 拟区域中央为注水井，4 口采油井分布 在储层四边，无效网格分布在下面三个 部分，分别为左下边的单个无效网格、 两个连接无效网格以及右下边6个相连 的无效网格。根据前面的数学模型及第三种无 效网格处理方法，编写了流线数值模拟 计算程序，模拟3个月后，得到了模拟 区域一注四采井网形式的流线分布特 征，如图3所示。I 2 I 4 5 I T B S IU H 12 0 M IT II 194V图3存在死网格流线分布从图3中可以看出流线从中间一口 注水井出发，流向周

14、边的四口生产，流 线碰到死网格区域会明显的出现变化， 而三部分不可流区域无流线通过，流线 绕过不可流动区域，流向了注水井，流 线明显改变了流线原有的形状，证明了 可不流动网格处理的有效性及合理性。由于不可流动区域的存在，靠近不 可流动区域，流线分布较密，而流线分 布密集地区反映了注水流量大，驱油效 果好。4结论1）文中给出了压力方程的三种死 网格处理方法中，第三种方法计算最为 简便，方程系数矩阵算法最容易得出， 第一种方法无法形成统一、有规律的计 算方法，尤其不适应程序化、软件化数 值模拟。而第二种方法计算过程繁杂， 模拟过程的不断赋值，大大的降低了计 算的效率。而第三种方法弥补了前两种 方法

15、的不足，即形成了统一的系数矩阵 计算方法，又有相对较快的计算速度。2 ）直接设置法可将不规则形状油 藏按规则的矩形油藏进行数值模拟计 算，并形成了统一的系数矩阵计算方 法，简化了模拟计算过程。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韩大匡，陈钦雷，闫存章.油藏数值模拟 基础.石油工业出版社，19934王洪宝，苏振阁，陈忠云.油藏水驱开发三 维流线模型J.石油勘探与开发，2004, 31(2): 99-1035R.C.M.Portella,T.

18、A.Hewett. Fast 3-D Reservoir Simulation and Applications Using StreamlinesJ.SPE39061 6刘慧卿.油藏数值模拟方法专题.东营： 石油大学出版社，20017姚军，吴明录，戴卫华等.流线数值试井解 释模型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期作者简介：王宇，男，陕西省西安市人，助 理工程师，硕士，主要从事油气藏开发研究 工作，邮编 710075,联系方式：陕西省延安市甘泉县下寺

20、湾镇下 寺湾采油厂， 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