油藏数值模拟毕业设计

本科毕业设计〔论文〕基于有限差分的油水两相渗流方程求解学生姓名：学号：专业班级：指导教师：2012年摘要为了保证油藏的稳定产油量以及最终采收率，以获得最大的经济效益，必须对油藏的压力以及饱和度等参数进行监控，因而需要对油藏进行模拟，以确定适宜的开采时间、注水量、开采速度等问题，因而油藏数值模拟对于提高油田效益是至关紧要的。油水两相渗流方程描述了原油开发过程中的油水两相流动过程，是一对耦合的偏微分方程，针对这组方程采用有限差分方法实现其数值求解，以对油藏的压力、饱和度进行即时监控。本文应用隐式压力——显式饱和度解法〔IMPES方法〕来对油藏油水两相一维渗流进行模型建立与求解，这一方法来自于Sheldon等人，以及Stone和Gardner的著作，根本思想是合并流体方程得到一个只含有压力的方程，某一时间步的压力求解出来后，饱和度采用显式更新。经过Matlab编写程序，实现对一维油水两相渗流方程的数值求解。关键词：油水两相；有限差分；压力；含水饱和度；油藏模拟ABSTRACTDisplacingoilbywaterisalwaysusedintheexploitationoftheoilreservoir.Inordertoobtainthemaximumeconomicbenefits,theparametersofreservoirpressureandsaturationshouldbemonitoredtoensurethestabilityofoilproductionandtheultimaterecoveryofthereservoir.Thusthereservoirmustbesimulatedsoastodeterminetheappropriaterecoverytime,waterinjectionrateandrecoveryrate.Itisnecessarytoestablishtheoil-watertwo-phaseflowmodelfortheoil-watertwo-phaseflowequation,coupledpartialdifferentialequations,candescribetheoil-watertwo-phaseflowintheprocessofcrudeoildevelopment.Achievingthenumericalsolutionaccordingtothisgroupofequationsandusingthefinite-differencemethodcandorealtimemonitoringforreservoirpressureandsaturationsothatappropriatemeasurescouldbetakentoensuretheminingofoilproductioninhighefficiencyandsustainable.Inthisthesis,implicitpressure-explicitsaturationsolution(IMPESmethod)whichisinventedbySheldon,etc.andfromworksofStone’sandGardner’sisappliedforthebuildingandsolvingofreservoir’sone-dimensionalwater-oiltwo-phaseseepagemodel.TThebasicideaiscombiningthefluidequationstoachieveaequationwhichcontainsonlythepressure.Aftersolvingthepressureattheverytimeandstep,thesaturationwillbedisplayedupdated.WritingprogramsthroughMatlab,one-dimensionalwater-oiltwo-phaseseepageequationswouldbesolvedKeywords:water-oiltwophase;finitedifference;pressure;watersaturation;reservoirsimulation目录TOC\o"1-3"\h\u12840第1章前言1287301.1油藏数值模拟在油田开发中的重要意义158141.2国内外研究方法现状1147441.3本文章节安排29376第2章相关根底知识3257112.1油藏渗流力学相关知识3282902.2有限差分法820279第3章油水两相渗流机理和求解11240583.1数学模型的建立11220123.2数学模型的求解的方法及参数处理1272323.3差分方程组得建立及求解15115313.4有关单位换算18175593.5计算程序框图205323.6计算实例2122380第4章结论238043致谢2423305参考文献2527215附录26第1章前言1.1油藏数值模拟在油田开发中的重要意义油气藏是在单一圈闭中具有统一压力系统的油气聚集单元。在原始条件下，油气藏处于平衡状态；当受到干扰〔如打井、生产)时，原来的不平衡状态被打破，油气藏处于动态变化中。油气藏从投入生产到最后废弃就是一个不断变化的动态过程。描述或实现油气藏动态变化的过程叫做模拟，油藏数值模拟是用数值的方法来求解描述油藏中流体渗流特征的数学模型，是一门将计算机、应用数学、油藏工程等结合起来的综合性工程应用学科，在油气田开发方案设计和动态分析中有十分重要的作用，为了保证油藏的稳定产油量以及最终采收率，以获得最大的经济效益，必须对油藏的压力以及饱和度等参数进行监控，因而需要对油藏进行模拟，以确定适宜的开采时间、注水量、开采速度等问题，需建立油水两相渗流模型，以对油藏的压力、饱和度进行即时监控，及时采取适宜的开采措施来保障油田的高效、高产、持久生产,它能够从油田实际出发对油藏进行精确地描述，以最少的投资、最适宜的速度去获得最高的最终采收率，也就是要获得最大的效益。1.2国内外研究方法现状油藏模型中有压力的、和饱和度、两组未知量，目前根本上有两类求解方法。一类是顺序求解法，即先求压力项，后求饱和度；另一类是联立求解法，即同时求解压力项和饱和度项。同时，由于方程中含非线性系数，他们依赖于压力和饱和度的变化，在求解数学模型时有如下几种处理方法：(1)隐式压力显式饱和度方法〔IMPES〕即隐式求解压力方程，显式求解饱和度方程。它属于顺序求解法的一种，是数值模拟中最常用的、最简单的一种方法。具有占内存小、计算工作量小、方法简便等优点。但该方法存在两个问题：第一，达西项的系数处理是显式的，因此对如锥进的问题，由于井底周围流速高，压差变化大，而存在较大误差，对于强非线性问题的适应性也差；第二，饱和度的计算是显式的，当时间步长较大时，会出现解的不稳定性。因此，IMPES方法只适用于一般的弱非线性渗流问题，对于某些非线性问题如注气、气锥或水锥等问题，IMPES方法无能为力，即使时间步长取得很小，仍会出现解的振荡或算出的压力和饱和度为负值的情况，以至于模拟计算无法正常进行。(2)半隐式方法半隐式方法属于联立求解的方法的一种，也是数值模拟中常用的一种方法。其根本思路是：联立求解油相方程和水相方程，同时求出压力和饱和度，因此压力和饱和度都是隐式求解。在计算过程中，半隐式方法对方程右端项的处理与IMPES方法完全相同，不同之处在于对方程左端项的处理。它需要对方程左端的达西系数、产量项及毛管压力进行泰勒级数展开，并忽略二阶小量，一阶导数项用n时刻的值。由于系数是处理近似的，并未真正用n+1时刻的值，所以这种方法叫做半隐式方法。(3)隐式压力隐式饱和度法〔IMPIMS〕隐式压力隐式饱和度法实际上是IMPES方法和半隐式方法的混合变种，它也是属于顺序求解方法的一种。它既有半隐式方法求解饱和度的特点，又保存了IMPES方法省内存、省工作量的特点。IMPIMS方法〔TheImplicitPressureImplicitSaturationMethod〕实际上是IMPIMS方法和半隐式方法的混合和变种，它也是属于顺序求解方法的一种。它既有半隐式方法求解饱和度的特点，又保存了IMPES方法省内存、省工作量的特点。IMPIMS方法的求解思路是：压力和饱和度分开顺序求解，求解压力时可直接利用IMPES方法的压力求解方法，然后将求出的压力值代入半隐式方法的水相差分方程中，将该方程化为只有含水饱和度的一个变量的差分方程，再用隐式计算格式求解即可。(4)全隐式方法对于某些强非线性渗流问题如高速渗流等，即使使用半隐式方法也会引起计算结果的波动，或者时间步长只能取到很小。为此提出了全隐式方法,这里简单介绍一种与半隐式方法相类似的全隐式方法，即达西系数项也用泰勒级数展开并忽略二阶小量，但一阶导数项不用时刻的值而用时刻的值，由此构成一个非线性的代数方程组，可用牛顿迭代法或其他非线性代数方程组的解法进行求解。1.3本文章节安排论文根据质量守恒原理建立了一维油水两相渗流的数学模型，经过逐步简化得到实验室进行单管模型一维水驱油实验时的数学模型，利用IMPES方法对此模型进行处理。第一局部首先列述了油藏渗流力学的相关知识，对油藏的油水渗流机理进行了详细的介绍以及油水两相渗流方程的推导过程，然后对有限差分法的根底知识进行了详细的介绍；第二局部主要是油水两相渗流模型的建立过程以及应用此方法对实例油藏的模拟。最后综述了油水两相渗流方程的建立过程中所遇到的问题、模型建立的理想条件及其准确性、针对实例模拟得出的结论。第2章相关根底知识2.1油藏渗流力学相关知识(一)、渗流过程中的力学分析油气水之所以能在岩石空隙中渗流，是各种力共同作用的结果。流体的重力和重力势能流体具有质量，在中立场中受到地球的吸引力而具有重力。在渗流过程中，流体的重力和它的相对位置联系起来，就表现为重力势能。这种重力势能用压力表示为：〔2-1〕式中表示重力势能的压力，Pa流体密度，重力加速度，相对位置高差，重力有时是动力，有时为阻力。流体的质量和惯性力惯性是物体所固有的一种物理特征，其大小取决于质量。流体质量的大小一般用密度表示。流体由于具有质量，因此也具有惯性。当流体运动时，惯性使其总要维持原状，因而惯性力在渗流过程中多表现为阻力。根据牛顿第二定律，惯性力用压力表示为：〔2-2〕式中，惯性力质量v速度加速度由于渗流时渗流速度通常很小，因此常忽略惯性力。流体的粘度及粘滞力粘滞力是流体的一种特殊属性。在流动的流体中，如果各层的流速不同，将有一对作用力及反作用力，使原来快的流层减速，而慢的加速。这一对等值而又反向的力阻碍着流层的相对运动。流体的这种属性叫粘滞力。度量粘滞性大小的参数叫粘度。由牛顿内摩擦定律表示为：〔2-3〕式中，A两流层的接触面积，沿流层法线方向的流速梯度，F内摩擦力，N粘滞系数，粘度的单位通常用表示：1粘度的单位以表示时成为泊〔P〕：1P=100cP〔厘泊〕在渗流中，粘滞力为阻力，且动力消耗主要用于克服渗流时流体的粘滞阻力。岩石及流体的压缩性及弹性力岩石和其中饱和的流体均具有压缩性〔或弹性〕，因此使得在油气渗流过程中产生弹性力。油气层岩石埋藏于地下几百米甚至几千米，油气层上面覆盖的岩柱压力被油气层本身骨架和其中饱和的流体所承受。因此，储层岩石和其中的流体都处于受压缩状态。油层除承受上覆岩柱压力外，本身也受油层压力。油气层岩石和其他固体一样，在外力作用下，它的形状和体积都要发生变化；当消除外力时，它又能恢复到原来的形状和大小。我们把岩石能恢复原状的性质叫岩石的弹性，又称压缩性。油气田在开采以前，油层内岩石和流体都处于均衡受力状态，各种力是相互平衡的；当油气层投入开采之后，油气层的压力不断下降，上覆岩柱压力和油层内流体压力之间形成压力差，使之失去平衡而迫使岩石颗粒变形，排列更加紧密，结果导致岩层孔隙体积减小。岩石颗粒变形后，孔隙体积的缩小程度取决于岩柱压力和油气层流体压力的差值，同时还与岩石本身的压缩性有关。岩石的压缩性通常用压缩系数表示：式中，孔隙体积的变化量；岩石的外表体积；岩石的压缩系数。岩石的压缩系数表示油层压力每降低单位压力时，单位体积岩石中孔隙提及的缩小量。由于岩石颗粒的组分不同，且岩石的孔隙结构不同，孔隙形状和承受力的作用点均不一样，因此各类岩石的压缩系数时不一样的。岩层中的流体〔油、气、水〕也具有压缩性，当作用与流体上的外力增加时，其体积会缩小，繁殖会膨胀。在油气开采过程中，流体体积膨胀，产生弹性力，推动流体流入井底。液体的压缩性常用液体的压缩系数来描述，可表示为：〔2-4〕式中，液体的绝对体积；压力改变时，液体体积相应的变化量；液体的压缩系数，表示改变单位压力时，单位体积液体的体积改变量。需注意的时，对于固体而言，弹性变形时有范围的，当外力超过某一限度时，将会发生形变。流体在弹性介质和塑性介质中的渗流规律是不同的。5、毛管力多孔介质时由无数个微小的毛管连接而成的。在油田开发过程中，油气水的渗流时由一种流体驱动另一种流体，在两相界面上会产生压力跳跃，它的大小取决于分界面的弯曲程度和流体的性质。这个跳跃的力就称为毛管压力或毛管力，用表示。毛管力与流体的性质和弯曲程度之间的关系可以用拉普拉斯方程表示：〔2-5〕式中，R，r分界面的曲率半径。液液界面的外表张力。在渗流中，毛管力既可以表现为动力，也可表现为阻力。在驱替压力不大时，假设油藏岩石亲水，那么水驱油时毛管力为动力；假设油藏岩石亲油，那么水驱油时毛管力为阻力。〔二〕、与油藏有关的压力概念以上对油气渗流过程中的力学现象和作用机理进行了讨论。油藏内的流体所受的各种力往往以压力的形式表示。外力所做的功将一起地层内液体能量的变化，这种变化将通过压力的变化来反映。因此，从本质上说压力是用来表征油藏能量的物理量。下面讨论与油气储集层有关的几个压力的概念。1.原始地层压力岩石地层压力是指油藏开发前流体所受的压力。史册的岩石地层压力一般是在油田所钻的第一批探井中测得的。2.供应压力供应压力是指当油藏中存在液源供应区时，在供应边缘上的压力。在人工注水条件下，水井井底压力即供应压力。3.井底压力井底压力是指油井正常生产时，在生产井井底所测得的压力。井底压力也成为井底流压，简称流压。4.折算压力油藏在开发前，整个油田处于平衡状态，油藏各点的流体所具有的总能量时相等的。由水动力学可知，单位质量流体所具有的能量由比位能、比压能和比动能组成。用M表示油层中的任意点，用z表示M点的标高，P表示M点的实测压力值，表示油层条件下的液体密度，v表示液体的流速，那么M点的液体所具有的总能量即总水头H为：〔2-6〕〔三〕、单相不可压缩液体稳定渗流的根本微分方程推导：由于是不可压缩液体，不必考虑状态变化，所以不需要状态方程，那么根本方程有：运动方程：〔2-7〕连续性方程：〔2-8〕将运动方程带入连续性方程进行数学模型推导，可得：〔2-9〕由于为常数，所以得到：〔2-10〕上式极为单相不可压缩流体在均值地层中的稳定渗流的根本微分方程。他的适用条件为：〔1〕单相均质液体；〔2〕符合线性运动规律；〔3〕不考虑多孔介质及液体的压缩性；〔4〕稳定渗流；〔5〕渗流是等温的；〔四〕、油水两相渗流的根本微分方程：1〕运动方程：对于油相：〔2-11〕对于水相：〔2-12〕2)油水两相渗流的连续性方程：对于油相：〔2-13〕对于水相：〔2-14〕将运动方程带入连续性方程得：〔2-15〕〔2-16〕上面两式也可以写成：〔2-17〕〔2-18〕此外，还有饱和度方程：〔2-19〕他的使用条件是：1〕彼此互不相溶，不起化学反响的油水两相同时流动；2〕不考虑毛管力及重力作用；3〕岩石和流体均不可压缩；4〕渗流符合线性渗流定律；5〕渗流过程是等温的。如果研究的是油水两相稳定渗流过程，那么液体饱和度将不随时间变化，即，得：〔2-20〕〔2-21〕以上两式即为油水两相稳定渗流的综合微分方程。2.2有限差分法有限差分法：数值求解常微分方程或偏微分方程的方法。物理学和其他学科领域的许多问题再被分析研究之后，往往可以归结为常微分方程和偏微分方程的求解问题。一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散值后对应的函数值。但是从原那么上说，这种方法仍然可以到达任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值差值计算来近似得到。这种方法是指随着计算机的诞生和应用而开展起来的。其计算格式和程序的设计都比拟直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成局部。有限差分法的具体操作可分为两个局部：用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；求解差分方程组。在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。通常采用的是规那么的分割方式，这样可以便于计算机自动实现和减少计算的复杂性。网路线划分的交点成为节点。假设与某个节点P相邻的节点都是定义在场域内的节点，那么P点称为非正那么节点。在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。有限差分法的差分格式：一个函数在x点上的一阶和和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。如对一个单变量函数f(x)，x为定义在区间[a,b]的连续变量。以步长将[a,b]区间离散化，我们得到一系列节点=a,,，然后求出f(x)在这些点上的近似值。显然步长h越小，近似解的精度就越好。与节点相邻的节点有和，因此在点可以构造如下形式的差值：节点的一阶向前差分节点的一阶向后差分节点的一阶中心差分与点相邻两代男的泰勒展开式可以写为〔2-22〕〔2-23〕〔2-22〕-〔2-23〕，并忽略h的平方和更高阶的项得到一阶微分的中心差商表示：。〔2-24〕利用〔2-22〕和〔2-23〕式我们还可以得到一阶微分的向前，向后一阶差商表示：，〔2-25〕。〔2-26〕将〔2-22〕和〔2-23〕式相加，忽略h的立方及更高阶的项得到二阶微分的中心差商表示：。〔2-27〕利用〔2-24〕~〔2-27〕式，我们就可以构造出微分方程的差分格式。这里要指出的是：在构造差分格式时，究竟应该选择向前，向后还是中间差分或差商来代替微分方程中的微分或微商，应当根据由此得到的差分方程解的稳定性和收敛性来考虑。同时兼顾到差分格式的简单和求解的方便。上述差分步骤应用于偏微分：例如，对于的情况，拉普拉斯算符在0点作用在此函数上的值，也可以用临近的点上的函数值来表示出来。〔见图2.1.1，且时〕〔2-28〕图2-1节点0及其附近节点对微分方程数值求解的误差来源：方法误差〔或截断误差〕。这时由于采用的计算方法所引起的误差。例如上面我们介绍的差商表示中，采用的泰勒展开式展开到第n+1项时的截断误差。具体方法的误差阶数取决于在离散化时的近似阶数。因此假设改良算法就可以减少截断误差。舍入误差〔或计算误差〕。这是由于计算机的有限字长而造成数据在计算机中的表示出现误差。在计算机运算的过程中，随着运算次数的增加舍入误差会积累的很大。如果在屡次运算后，舍入误差的精度影响是有限的，那么这个算法是稳定的，否那么是不稳定的。不稳定的算法是不能用的。第3章油水两相渗流机理和求解3.1数学模型的建立假设条件如下：油藏中仅存在油水两相渗流，油水互不溶解，且各自符合达西定律；岩石、流体均可压缩；考虑掩饰的非均质性及各向异性；不考虑油水之间毛管力的影响；忽略重力；这里根据数学模型的一般式，经逐步简化，得到实验室进行单管模型的一维水驱油实验室的数学模型。当考虑三维非均质油藏，油水互不相容，可压缩流体和岩石，考虑毛管力和重力时，数学模型的一般式为：〔〕〔3-1〕简化到一维，并忽略重力项：〔〕(3-2〕假设：①不考虑掩饰的压缩性〔即=常数〕，不考虑流体的体积变化〔即〕；②油水粘度为常数。于是得：水相：〔3-3〕油相：〔3-4〕式中，为地面标准状况下单位时间内单元体中注入〔或采出〕的体积流量。上述两个偏微分方程中的未知量有4个，即、、、，因此还需要写出两个辅助方程，即：〔3-5〕〔3-6〕初始条件为：〔3-7〕边界条件为：〔3-8〕上述边界条件中，注入量产出量均为，说明该水驱油试验为稳定驱替。上述〔3-1〕—〔3-8〕构成了该问题的完整的数学模型。利用数值方法进行求解后，可得到在不同的注入速率下，模型中任意一点的压力、饱和度随时间的分布和变化。3.2数学模型的求解的方法及参数处理对以上数学模型进行差分求解之前，这里首先对未知量的求解方法及有关参数的处理进行说明。1、数学模型的求解方法上述数学模型中有压力、和饱和度、两组未知量，本文应用隐式压力显式饱和度〔IMPES〕求解法进行求解。IMPES方法的根本思路：〔1〕通过乘以适当的系数，合并油方程和水方程，以消去微分方程组中的、，得到一个只含有、的方程。〔2〕由毛管力方程，可得，带入上面合并后的方程，得到一个只含有的方程，成为压力方程。〔3〕方程左端达西项系数上一时间段的值，同时毛管力也用上一时间阶段的值，即显示处理系数。于是可形成一个高阶现行代数方程组，用迭代法可以进行求解，先求出，然后得。将带入水相方程，用显式方法求出，然后得。IMPES方法具有所占内存小、计算工作量小、方法简便等优点。但该方法存在两个问题：第一、达西项的系数处理是显式的，因此对如锥进的问题，由于井底周围流速较高，压差变化大，而存在较大的误差，对于强非线性问题的适应性也差；第二、饱和度的计算是隐式的，当时间步长较大时，会出现解的不稳定性。因此，IMPES方法只适用于一般的弱非线性渗流的问题，对于某些非线性渗流的问题如注气、气锥或水锥等问题，IMPES方法无能为力，即使时间步长取得很小，仍会出现解的震荡或算出的压力和饱和度为负值的情况，以致模拟计算无法正常进行。2、参数处理在用有限差分对数学模型进行求解时，首先要将连续的油藏问题离散化为网格单元，然后对每一个网格单元，读入包括深度、有效厚度、孔隙度、渗透率、饱和度等根本参数。所有给定的这些参数都是网格节点处的值，在两个网格节点中间处的参数值是未知的，因此需要进行相应的处理。〔1〕渗透率的取值渗透率K是空间函数，其取值由以下几种方法：算数平均：〔3-9〕加权平均：〔3-10〕调和平均：〔3-11〕几何平均：〔3-12〕〔2〕相对渗透率的取值相对渗透率的取之原那么上是去流动方向上的上游节点值，通常称为上游权法。如下列图所示，取值方法是：当时，即由流向时〔3-13〕当时，即由流向时(3-14)式中，为第个节点的势。图3-1井点分布图采用上游权处理主要是对饱和度沿流向变化滞后进行修正。由于流动系数关于时间的显式处理造成流动的滞后，如流动系数在空间上按相邻节点的算数平均值取值，那么会加重滞后现象，影响解的精度。因此，上游权的处理实质上是将显式处理造成的时间上的滞后用空间上的向前来进行弥补。〔3〕毛管力曲线的处理由于油藏岩石的非均质性，即使用同一油层的岩心所测得的毛管压力曲线也有所不同，用J函数方法可以从实验室提供的大量毛管力资料中选择适宜的毛管力曲线。J函数是把流体界面张力、岩石润湿性及渗透率和孔隙度等的影响综合在一起来表征油层的毛管压力曲线特征的一个无因次函数。实践说明，它是处理毛管压力曲线的一种有效方法。J函数定义如下：〔3-15〕式中——油、水界面张力，；——水的润湿角，〔〕；——岩石毛管压力，；——渗透率，；——孔隙度，小数；——与相对应的无因次量。将众多的毛管压力资料按式〔3-15)计算，做出关系图，如图3-2所示，然后拟合出一条关系曲线。当油层岩性较为接近时，关系并根据油层，值分布的峰值区域可求出一组有代表性的关系。在进行油藏数值模拟计算时，一般需要将油、水相对渗透率曲线和毛管压力曲线离散化，然后在模拟计算过程中应用各种插值手段求出各网格在不同饱和度下的相对渗透率值和毛管压力值。图3-2曲线图3.3差分方程组得建立及求解下面用方法来建立差分方程组并进行求解。隐式求压力上面的数学模型式(3-3)至(3-6)中，有，，，4个未知量。首先我们利用及方程(3-3)和(3-4)的关系消去，，使其成为只含有压力，的方程。为此，式(3-3)＋式(3-4)得：〔3-16〕根据，可得，代入上式可得只含有压力的方程。这里为简化起见，假设毛管压力为0，于是可得，因此上式简化为：(3-17)令，分别表示油、水两相的流动系数：，表示总的流动系数：表示油、水两相的总流量：那么方程〔3-17)可简化为：(3-18)方程〔3-18)即为IMPES方法所化简得压力方程。下面写出该方程的隐式差分格式。假设采用块中心网格，且网格大小相等，均为，并假设网格为水注入处，网格为油和油、水产处处。式〔3-18〕的隐式差分方程为：(3-19)下面分3种情况来讨论方程〔3-19)。(1)对于第2个至第个网格，既无注入也无采出，因此方程〔3-18)中，于是可简化为：(3-20)假设系数按上游权原那么取值，并用n时刻的值〔即显示处理〕，那么上式方程可简化为：(3-21)即：(3-22)(2)对于第1个网格，单位体积中注入的体积流量为，式(3-18)中第二项由于取上游权，，于是可简化为：(3-23)整理上式，可得：(3-24)两端同乘以〔网格单元的体积〕，并令〔该网格的注入量〕，那么上式可化简为：(3-25)(3)对于第N个网格，其单位体积中采出的体积流量为，式(3-18)中第一项由于取上游权，，于是可简化为：(3-26)两端同乘以，令，得：式〔3-22),式〔3-25)和式〔3-27〕构成了从第一个网格到最后一个网格〔第个〕的线性代数方程组，矩阵方程如下：=〔3-28〕方程组的系数矩阵为三角矩阵，假设用定流量来进行求解，那么上述压力是与时间无关的量，不符合要求。可用方法进行求解，从而得到，，，在时刻的压力值。显式求饱和度数学模型中的水相方程〔3-3〕的差分方程为：〔3-29〕上式中的未知量为，也可分为以下三种情况来讨论：对于第2个至第个网格，既无注入也无采出，因此方程〔3-28〕中的，于是可得：〔3-30〕由于已经通过隐式求压力得到，故可用上式顺序〔〕求出每个网格节点的饱和度值。(2)对于第一个网格，〔即单位体积的注水速率〕，而且式〔3-29〕中的左面查分的第二项去上游权，，于是可化简为：〔3-31〕两边同乘以并化简，且，那么有：于是：〔3-32〕〔3〕对于最后一个〔第N个〕网格，其单位体积的产出量，实验中为定值。其中，产水量为，而且是跟随时间变化的，因此首先要求出的值。最后一个〔第N个〕网格的压力为，第N-1个网格的压力为，那么对于油和水来说均是相等的，故流出油、水提及的流量取决于流动系数、的大小。〔3-33〕〔3-34〕所以：〔3-35〕于是〔3-36〕方程〔3-28〕中的左面差分的第一项由于取上游权，，于是可化简为：〔3-37〕两边同乘以，并令，那么：〔3-38〕于是：〔3-39〕利用式〔3-30〕、式〔3-32〕和式〔3-39〕可求得网格的值。以上所表达的隐式求压力、显式求饱和度的计算过程即为IMPES方法从时刻到时刻的求解过程。这样从初始时刻开始，一步步依次求解下去，直到求得所要求的时间的压力和饱和度为止，油水两相的渗流方程就可得解，所建立的油藏模型就得以模拟出来。3.4有关单位换算以上的计算过程中需用到达西公式，但达西公式中所用的单位为水力学单位，而实际油田中常用的单位为工程单位，二者之间的换算关系见表3-1表3-1换算关系水力单位工程单位QK0.001uAm0.1MPaMPa〔1〕达西公式为，其中各变量均采用水力学单位。〔2〕当采用工程单位时，达西公式变为：。〔3〕因此，由于水力学单位换算到工程单位，达西公式前加了一个系数，即：3.5计算程序框图程序开始读入读入读入读入读入显式求计算值计算不同下的、读入显式求计算值计算不同下的、读入计算压力方程的系数矩阵和右边项计算压力方程的系数矩阵和右边项用追赶法用追赶法否否是是打印打印结束结束图3-3一维油水两相程序框图3.6计算实例下面把给定的油藏区块按一维方向均分成7个网格井点，以第一个网格井为注水井，第七个网格井为采油井，一注一采模型进行求解：注水井生产井图3-4井网分布图：孔隙度，渗透率，束缚水饱和度，原始含油饱和度，地层油粘度，地层水粘度，注水井井底压力，生产井井底流压，有谁渗透率曲线如表1所示.油层长70m，宽，油层厚度，综合压缩系数，取，时间步长取一天，模拟时间为，油水相对渗透率曲线见附录。模拟曲线如下：图3-5最后一个网格点〔生产井〕关系曲线总液量中水所占有的分量即为含水率,此曲线反映了油井的含水率随时间的变化趋势，开始后一段时间内含水率接近0，因为水驱前沿未到达产油井，采出的完全是油，水驱前沿到达产油井后，含水率开始上升，经过相当长时间后含水率近似可以到达1，此时油井开采完毕。根据此曲线可以综合评估油井的本钱与效益，以便选择适宜的开采时间。图3-61-7网格井点在T时刻的压力曲线此图反映了此模型从第一个井点到第七个井点的压力变化，可以看出是等差变化的。这是因为不考虑油藏的毛管力、重力以及岩石和流体的压缩性所致，模型为刚性水驱开采，故从注水井至开采井的压力是等差变化的。图3-7最后一个网格井点的曲线此曲线反映了最后一口网格井的随时间的变化趋势，随着水驱前沿的推进，水侵占油藏的孔隙以排除原油，所以油藏的含水饱和度会随时间而增大，最后到达稳定值，因为此时油藏孔隙中的剩余油不能被排出，仍然占据着油藏孔隙，所以随着时间的推移，含水饱和度不再变化。第4章结论1、本文应用隐式压力——显式饱和度解法〔IMPES方法〕来对油藏油水两相一维渗流进行模型建立与求解，通过合并流体方程得到一个只含有压力的方程，某一时间步的压力求解出来后，饱和度采用显式更新对油藏的模拟及求解。2、模型的求解有定流量和定压力两种方法，本文建立的模型不能用定流量来进行求解，因为假设用定流量进行求解那么此模型的压力方程与时间无关，会有无穷多解，这是因为只要保持注水井和采油井的压力恒定，就可以到达流量恒定，所以这样的压力解是无穷的。3、本文所用的IMPES方法具有占内存小、计算工作量小、方法简便等优点。4、本文建立的模型是在不考虑岩石和流体的压缩性以及油水的粘度、忽略重力，在这种刚性模型下对油藏的模拟准确度可能有一定的影响。IMPES方法达西项的系数处理是显式的，因此对如锥进的问题，由于井底周围流速高，压差变化大，而存在较大误差，对于强非线性问题的适应性也差；饱和度的计算是显式的，当时间步长较大时，会出现解的不稳定性。5、通过含水率曲线说明，相当长时间后，含水率接近100%，此时此油带的水驱油开采完毕。6、本文对一维油水两相渗流进行了数值模拟，对一维各个以划分的网格在个时刻的压力、含水饱和度、含水率等方面都有了一个全面的显示，能够对油藏较为准确的模拟，按此方法可以对二维、三维油水两相以及二维、三维油气水三相渗流进行拓展。致谢非常感谢张晓东老师在我大学的最后学习阶段——毕业设计阶段给自己的指导，从最初的定题，到资料收集，到写作、修改，到论文定稿，他给了我耐心的指导和无私的帮助。为了指导我们的毕业论文，他放弃了自己的休息时间，他的这种无私奉献的敬业精神令人钦佩，在此我向他表示我诚挚的谢意。同时，感谢所有任课老师和所有同学在这四年来给自己的指导和帮助，是他们教会了我专业知识，教会了我如何学习，教会了我如何做人。正是由于他们，我才能在各方面取得显著的进步，在此向他们表示我由衷的谢意，并祝所有的老师培养出越来越多的优秀人才，桃李满天下！通过这一阶段的努力，我的毕业论文《基于有限差分法的油水两相渗流方程求解》终于完成了，这意味着大学生活即将结束。在大学阶段，我在学习上和思想上都受益非浅，这除了自身的努力外，与各位老师、同学和朋友的关心、支持和鼓励是分不开的。在本论文的写作过程中，我的导师张晓东老师倾注了大量的心血，从选题到开题报告，从写作提纲，到一遍又一遍地指出每稿中的具体问题，严格把关，循循善诱，在此我表示衷心感谢。同时我还要感谢在我学习期间给我极大关心和支持的各位老师.写作毕业论文是一次再系统学习的过程，毕业论文的完成，同样也意味着新的学习生活的开始。我将铭记我曾是一名中国石油大学学子，在今后的工作中把惟真惟实的优良传统发扬光大。感谢各位专家的批评指导参考文献[1]哈利德.阿齐兹，安东尼.塞特瑞.油藏数值模拟第一版.北京：石油工业出版社，2004.[2]HenrikLof，MargotGerritsen,MarcoThiele.Parallelstreamlinesimulation.SPE113543,2008.[3]谢海兵，桓冠仁，郭尚平，等.PEBI网格二维两相流数值模拟.石油学报，1999，20〔20〕：57-61.[4]韩大匡，陈钦雷，阎存章．油藏数值模拟根底.北京：石油出版社，1993．[5]谢俊.剩余油饱和度平面分布方法研究及应用.西安石油学院学报，1998,13〔4〕:40-42.[6]张建国，杜殿发，侯健.油气层渗流力学．第二版．东营：石油大学出版社，2009.[7]李淑霞，谷建伟.油藏数值模拟根底东营：石油大学出版社2009.[8][9]邓建中，刘之行.计算方法.第二版.西安交通大学出版社2001.[10]胡良剑，邓晓君.MATLAB实验高等教育出版社2006.[11]卓金武.MATLAB在数学建模中的应用.北京航空航天大学出版社，2011.附录附录A表1油水相对渗透率曲线数据表序号10.200.780.0020.300.610.0330.400.460.0740.500.310.1150.600.190.1660.700.100.2370.800.000.30附录B程序代码：%******模型求解clearallclc%************赋初值*************N=7;%井数T=1800;%天数S=100;%面积dx=10;WI=5;Pwf=10;p0=[15;15];fai=0.25;p=zeros(N,1);q=zeros(N,1);A=zeros(N,N);K=1;Uo=5;Uw=1;Kro=zeros(N,T+2);Krw=zeros(N,T+2);tt=linspace(0,T,T+1);dert=1;Sw=0.2+p;%sw0Swtab=[0.20.30.40.50