基于有限差分的油水两相渗流方程求解油藏数值模拟毕业设计.doc

本本 科科 毕毕 业业 设设 计 论文 计 论文 基于有限差分的油水两相渗流方程求解 中国石油大学 华东 本科毕业设计 论文 摘 要 为了保证油藏的稳定产油量以及最终采收率 以获得最大的经济效益 必须对油 藏的压力以及饱和度等参数进行监控 因而需要对油藏进行模拟 以确定合适的开采 时间 注水量 开采速度等问题 因而油藏数值模拟对于提高油田效益是至关紧要的 油水两相渗流方程描述了原油开发过程中的油水两相流动过程 是一对耦合的偏微分 方程 针对这组方程采用有限差分方法实现其数值求解 以对油藏的压力 饱和度进 行即时监控 本文应用隐式压力 显式饱和度解法 IMPES 方法 来对油藏油水两相 一维渗流进行模型建立与求解 这一方法来自于 Sheldon 等人 以及 Stone 和 Gardner 的著作 基本思想是合并流体方程得到一个只含有压力的方程 某一时间步的压力求 解出来后 饱和度采用显式更新 经过 Matlab 编写程序 实现对一维油水两相渗流方 程的数值求解 关键词 关键词 油水两相 有限差分 压力 含水饱和度 油藏模拟 中国石油大学 华东 本科毕业设计 论文 ABSTRACT Displacing oil by water is always used in the exploitation of the oil reservoir In order to obtain the maximum economic benefits the parameters of reservoir pressure and saturation should be monitored to ensure the stability of oil production and the ultimate recovery of the reservoir Thus the reservoir must be simulated so as to determine the appropriate recovery time water injection rate and recovery rate It is necessary to establish the oil water two phase flow model for the oil water two phase flow equation coupled partial differential equations can describe the oil water two phase flow in the process of crude oil development Achieving the numerical solution according to this group of equations and using the finite difference method can do real time monitoring for reservoir pressure and saturation so that appropriate measures could be taken to ensure the mining of oil production in high efficiency and sustainable In this thesis implicit pressure explicit saturation solution IMPES method which is invented by Sheldon etc and from works of Stone s and Gardner s is applied for the building and solving of reservoir s one dimensional water oil two phase seepage model TThe basic idea is combining the fluid equations to achieve a equation which contains only the pressure After solving the pressure at the very time and step the saturation will be displayed updated Writing programs through Matlab one dimensional water oil two phase seepage equations would be solved Keywords water oil two phase finite difference pressure water saturation reservoir simulation 中国石油大学 华东 本科毕业设计 论文 目目 录录 第 1 章 前 言 1 1 1 油藏数值模拟在油田开发中的重要意义 1 1 2 国内外研究方法现状 1 1 3 本文章节安排 2 第 2 章 相关基础知识 3 2 1 油藏渗流力学相关知识 3 2 2 有限差分法 8 第 3 章 油水两相渗流机理和求解 11 3 1 数学模型的建立 11 3 2 数学模型的求解的方法及参数处理 12 3 3 差分方程组得建立及求解 15 3 4 有关单位换算 18 3 5 计算程序框图 20 3 6 计算实例 21 第 4 章 结 论 23 致 谢 24 参考文献 25 附 录 26 中国石油大学 华东 本科毕业设计 论文 1 第 1 章 前 言 1 1 油藏数值模拟在油田开发中的重要意义 油气藏是在单一圈闭中具有统一压力系统的油气聚集单元 在原始条件下 油气 藏处于平衡状态 当受到干扰 如打井 生产 时 原来的不平衡状态被打破 油气藏 处于动态变化中 油气藏从投入生产到最后废弃就是一个不断变化的动态过程 描述或实现油气藏动态变化的过程叫做模拟 油藏数值模拟是用数值的方法来求 解描述油藏中流体渗流特征的数学模型 是一门将计算机 应用数学 油藏工程等结 合起来的综合性工程应用学科 在油气田开发方案设计和动态分析中有十分重要的作 用 为了保证油藏的稳定产油量以及最终采收率 以获得最大的经济效益 必须对油 藏的压力以及饱和度等参数进行监控 因而需要对油藏进行模拟 以确定合适的开采 时间 注水量 开采速度等问题 需建立油水两相渗流模型 以对油藏的压力 饱和 度进行即时监控 及时采取合适的开采措施来保障油田的高效 高产 持久生产 它能 够从油田实际出发对油藏进行精确地描述 以最少的投资 最合适的速度去获得最高 的最终采收率 也就是要获得最大的效益 1 2 国内外研究方法现状 油藏模型中有压力的 和饱和度 两组未知量 目前基本上有两类求解 w P o P w S o S 方法 一类是顺序求解法 即先求压力项 后求饱和度 另一类是联立求解法 即同 时求解压力项和饱和度项 同时 由于方程中含非线性系数 他们依赖于压力和饱和 度的变化 在求解数学模型时有如下几种处理方法 1 隐式压力显式饱和度方法 IMPES 即隐式求解压力方程 显式求解饱和度方程 它属于顺序求解法的一种 是数值 模拟中最常用的 最简单的一种方法 具有占内存小 计算工作量小 方法简便等优 点 但该方法存在两个问题 第一 达西项的系数处理是显式的 因此对如锥进的问 题 由于井底周围流速高 压差变化大 而存在较大误差 对于强非线性问题的适应 性也差 第二 饱和度的计算是显式的 当时间步长较大时 会出现解的不稳定性 因此 IMPES 方法只适用于一般的弱非线性渗流问题 对于某些非线性问题如注气 气 锥或水锥等问题 IMPES 方法无能为力 即使时间步长取得很小 仍会出现解的振荡或 算出的压力和饱和度为负值的情况 以至于模拟计算无法正常进行 2 半隐式方法 中国石油大学 华东 本科毕业设计 论文 2 半隐式方法属于联立求解的方法的一种 也是数值模拟中常用的一种方法 其基 本思路是 联立求解油相方程和水相方程 同时求出压力和饱和度 因此压力和饱和 度都是隐式求解 在计算过程中 半隐式方法对方程右端项的处理与 IMPES 方法完全 相同 不同之处在于对方程左端项的处理 它需要对方程左端的达西系数 产量项及 毛管压力进行泰勒级数展开 并忽略二阶小量 一阶导数项用 n 时刻的值 由于系数 是处理近似的 并未真正用 n 1 时刻的值 所以这种方法叫做半隐式方法 3 隐式压力隐式饱和度法 IMPIMS 隐式压力隐式饱和度法实际上是 IMPES 方法和半隐式方法的混合变种 它也是属 于顺序求解方法的一种 它既有半隐式方法求解饱和度的特点 又保留了 IMPES 方法 省内存 省工作量的特点 IMPIMS 方法 The Implicit Pressure Implicit Saturation Method 实际上是 IMPIMS 方法和半隐式方法的混合和变种 它也是属于顺序求解方法的一种 它既有半 隐式方法求解饱和度的特点 又保留了 IMPES 方法省内存 省工作量的特点 IMPIMS 方法的求解思路是 压力和饱和度分开顺序求解 求解压力时可直接利用 IMPES 方法的压力求解方法 然后将求出的压力值代入半隐式方法的水相差分方程 1n 0 p 中 将该方程化为只有含水饱和度的一个变量的差分方程 再用隐式计算格式求解即 可 4 全隐式方法 对于某些强非线性渗流问题如高速渗流等 即使使用半隐式方法也会引起计算结 果的波动 或者时间步长只能取到很小 为此提出了全隐式方法 这里简单介绍一种与 半隐式方法相类似的全隐式方法 即达西系数项也用泰勒级数展开并忽略二阶小量 但一阶导数项不用时刻的值而用时刻的值 由此构成一个非线性的代数方程组 n1n 可用牛顿迭代法或其他非线性代数方程组的解法进行求解 1 3 本文章节安排 论文根据质量守恒原理建立了一维油水两相渗流的数学模型 经过逐步简化得到 实验室进行单管模型一维水驱油实验时的数学模型 利用 IMPES 方法对此模型进行处 理 第一部分首先列述了油藏渗流力学的相关知识 对油藏的油水渗流机理进行了详 细的介绍以及油水两相渗流方程的推导过程 然后对有限差分法的基础知识进行了详 细的介绍 第二部分主要是油水两相渗流模型的建立过程以及应用此方法对实例油藏 的模拟 最后综述了油水两相渗流方程的建立过程中所遇到的问题 模型建立的理想 中国石油大学 华东 本科毕业设计 论文 3 条件及其准确性 针对实例模拟得出的结论 第 2 章 相关基础知识 2 1 油藏渗流力学相关知识 一 渗流过程中的力学分析 油气水之所以能在岩石空隙中渗流 是各种力共同作用的结果 1 流体的重力和重力势能 流体具有质量 在中立场中受到地球的吸引力而具有重力 在渗流过程中 流体的重 力和它的相对位置联系起来 就表现为重力势能 这种重力势能用压力表示为 2 gzPz 1 式中 表示重力势能的压力 Pa z P 流体密度 3 mkg 重力加速度 g 2 sm 相对位置高差 zm 重力有时是动力 有时为阻力 2 流体的质量和惯性力 惯性是物体所固有的一种物理特征 其大小取决于质量 流体质量的大小一般用密 度表示 流体由于具有质量 因此也具有惯性 当流体运动时 惯性使其总要维持原 状 因而惯性力在渗流过程中多表现为阻力 根据牛顿第二定律 惯性力用压力表示为 2 ma dt v P d m a 2 式中 惯性力 a P 质量m v 速度 加速度a 由于渗流时渗流速度通常很小 因此常忽略惯性力 3 流体的粘度及粘滞力 中国石油大学 华东 本科毕业设计 论文 4 粘滞力是流体的一种特殊属性 在流动的流体中 如果各层的流速不同 将有一 对作用力及反作用力 使原来快的流层减速 而慢的加速 这一对等值而又反向的力 阻碍着流层的相对运动 流体的这种属性叫粘滞力 度量粘滞性大小的参数叫粘度 由牛顿内摩擦定律表示为 中国石油大学 华东 本科毕业设计 论文 4 2 dy dv uAF 3 式中 A 两流层的接触面积 2 m 沿流层法线方向的流速梯度 dy dv msm 1 F 内摩擦力 N 粘滞系数 粘度的单位通常用表示 uspa smPa 1smPasPa 10 3 粘度的单位以表示时成为泊 P scmg 1P 100cP 厘泊 在渗流中 粘滞力为阻力 且动力消耗主要用于克服渗流时流体的粘滞阻力 4 岩石及流体的压缩性及弹性力 岩石和其中饱和的流体均具有压缩性 或弹性 因此使得在油气渗流过程中产生 弹性力 油气层岩石埋藏于地下几百米甚至几千米 油气层上面覆盖的岩柱压力被油 气层本身骨架和其中饱和的流体所承受 因此 储层岩石和其中的流体都处于受压缩 状态 油层除承受上覆岩柱压力外 本身也受油层压力 油气层岩石和其他固体一样 在外力作用下 它的形状和体积都要发生变化 当 消除外力时 它又能恢复到原来的形状和大小 我们把岩石能恢复原状的性质叫岩石 的弹性 又称压缩性 油气田在开采以前 油层内岩石和流体都处于均衡受力状态 各种力是相互平衡 的 当油气层投入开采之后 油气层的压力不断下降 上覆岩柱压力和油层内流体压力之 间形成压力差 使之失去平衡而迫使岩石颗粒变形 排列更加紧密 结果导致岩层孔 隙体积减小 岩石颗粒变形后 孔隙体积的缩小程度取决于岩柱压力和油气层流体压力的差值 同时还与岩石本身的压缩性有关 岩石的压缩性通常用压缩系数表示 式中 孔隙体积的变化量 V 岩石的外表体积 f V 岩石的压缩系数 f C 岩石的压缩系数表示油层压力每降低单位压力时 单位体积岩石中孔隙提及的缩 中国石油大学 华东 本科毕业设计 论文 6 小量 中国石油大学 华东 本科毕业设计 论文 5 由于岩石颗粒的组分不同 且岩石的孔隙结构不同 孔隙形状和承受力的作用点 均不一样 因此各类岩石的压缩系数时不一样的 岩层中的流体 油 气 水 也具有压缩性 当作用与流体上的外力增加时 其 体积会缩小 繁殖会膨胀 在油气开采过程中 流体体积膨胀 产生弹性力 推动流 体流入井底 液体的压缩性常用液体的压缩系数来描述 可表示为 2 PV V C l l 1 1 4 式中 液体的绝对体积 l V 压力改变时 液体体积相应的变化量 l V P 液体的压缩系数 表示改变单位压力时 单位体积液体的体积改变 f C 量 需注意的时 对于固体而言 弹性变形时有范围的 当外力超过某一限度时 将 会发生形变 流体在弹性介质和塑性介质中的渗流规律是不同的 5 毛管力 多孔介质时由无数个微小的毛管连接而成的 在油田开发过程中 油气水的渗流 时由一种流体驱动另一种流体 在两相界面上会产生压力跳跃 它的大小取决于分界 面的弯曲程度和流体的性质 这个跳跃的力就称为毛管压力或毛管力 用表示 c P 毛管力与流体的性质和弯曲程度之间的关系可以用拉普拉斯方程表示 2 11 rR Pc 5 式中 R r 分界面的曲率半径 液液界面的表面张力 在渗流中 毛管力既可以表现为动力 也可表现为阻力 在驱替压力不大时 若油藏 岩石亲水 则水驱油时毛管力为动力 若油藏岩石亲油 则水驱油时毛管力为阻力 二 与油藏有关的压力概念 以上对油气渗流过程中的力学现象和作用机理进行了讨论 油藏内的流体所受的 各种力往往以压力的形式表示 外力所做的功将一起地层内液体能量的变化 这种变化将通过压力的变化来反映 因此 从本质上说压力是用来表征油藏能量的物理量 下面讨论与油气储集层有关的 中国石油大学 华东 本科毕业设计 论文 8 几个压力的概念 1 原始地层压力 0 P 岩石地层压力是指油藏开发前流体所受的压力 史册的岩石地层压力一般是在油 田 中国石油大学 华东 本科毕业设计 论文 6 所钻的第一批探井中测得的 2 供给压力 c P 供给压力是指当油藏中存在液源供给区时 在供给边缘上的压力 在人工注水条 件下 水井井底压力即供给压力 3 井底压力 w P 井底压力是指油井正常生产时 在生产井井底所测得的压力 井底压力也成为井 底流压 简称流压 4 折算压力 r P 油藏在开发前 整个油田处于平衡状态 油藏各点的流体所具有的总能量时相等 的 由水动力学可知 单位质量流体所具有的能量由比位能 比压能和比动能组成 用 M 表示油层中的任意点 用 z 表示 M 点的标高 P 表示 M 点的实测压力值 表 示油层条件下的液体密度 v 表示液体的流速 则 M 点的液体所具有的总能量即总水头 H 为 2 g v g P zH 2 2 6 三 单相不可压缩液体稳定渗流的基本微分方程推导 由于是不可压缩液体 不必考虑状态变化 所以不需要状态方程 则基本方程有 1 运动方程 2 P grad K v 7 2 连续性方程 2 0v div 8 将运动方程带入连续性方程进行数学模型推导 可得 2 0P grad K div 9 由于为常数 所以得到 u k 2 0 z P y P x P 2 2 2 2 2 2 10 中国石油大学 华东 本科毕业设计 论文 10 上式极为单相不可压缩流体在均值地层中的稳定渗流的基本微分方程 他的适 用条件为 1 单相均质液体 2 符合线性运动规律 中国石油大学 华东 本科毕业设计 论文 11 3 不考虑多孔介质及液体的压缩性 4 稳定渗流 5 渗流是等温的 四 油水两相渗流的基本微分方程 1 运动方程 对于油相 2 P grad K v o o o 11 对于水相 2 P grad K v w w w 12 2 油水两相渗流的连续性方程 对于油相 2 t S z v y v x v ozo yo ox 13 对于水相 2 t S z v y v x v wzw yw xw 14 将运动方程带入连续性方程得 2 t S z PK zy PK yx PK x o o o o o o o 15 2 t S z PK zy PK yx PK x w w w w w w w 16 上面两式也可以写成 2 t S P K o o o 17 2 t S P K w w w 18 中国石油大学 华东 本科毕业设计 论文 12 此外 还有饱和度方程 2 1SS wo 19 他的使用条件是 1 彼此互不相溶 不起化学反应的油水两相同时流动 2 不考虑毛管力及重力作用 3 岩石和流体均不可压缩 4 渗流符合线性渗流定律 5 渗流过程是等温的 如果研究的是油水两相稳定渗流过程 那么液体饱和度将不随时间变化 即 得 0 t S t S wo 2 0P K o o 20 2 0P K w w 21 以上两式即为油水两相稳定渗流的综合微分方程 2 2 有限差分法 有限差分法 数值求解常微分方程或偏微分方程的方法 物理学和其他学科领域 的许多问题再被分析研究之后 往往可以归结为常微分方程和偏微分方程的求解问题 一般说来 处理一个特定的物理问题 除了需要知道它满足的数学方程外 还应当同 时知道这个问题的定解条件 然后才能设计出行之有效的计算方法来求解 有限差分 法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值 在有限差分方法中 我们放弃了微分方程中独立变量可以取连续值的特征 而关注独 立变量离散值后对应的函数值 但是从原则上说 这种方法仍然可以达到任意满意的 计算精度 因为方程的连续数值解可以通过减小独立变量离散取值的间格 或者通过 离散点上的函数值差值计算来近似得到 这种方法是指随着计算机的诞生和应用而发 展起来的 其计算格式和程序的设计都比较直观和简单 因而 它的实际应用已经构 成了计算数学和计算物理的重要组成部分 有限差分法的具体操作可分为两个部分 中国石油大学 华东 本科毕业设计 论文 13 1 用差分代替微分方程中的微分 将连续变化的变量离散化 从而得到差分方 程组的数学形式 2 求解差分方程组 在第一步中 我们通过所谓的网络分割法 将函数定义域分成大量相邻而不重合 的子区域 通常采用的是规则的分割方式 这样可以便于计算机自动实现和减少计算 的复杂性 网路线划分的交点成为节点 若与某个节点 P 相邻的节点都是定义在场域 内的节点 则 P 点称为非正则节点 在第二步中 数值求解的关键就是要应用适当的 计算方法 求得特定问题在所有这些节点上的离散近似值 有限差分法的差分格式 一个函数在 x 点上的一阶和和二阶微商 可以近似地用它所临近的两点上的函数 值的差分来表示 如对一个单变量函数 f x x 为定义在区间 a b 的连续变量 以步 长将 a b 区间离散化 我们得到一系列节点 a xh 1 x 然后求出 f x 在这些点上的近x2axx xx 2312 hhbh nn xx 1 似值 显然步长 h 越小 近似解的精度就越好 与节点相邻的节点有和 i xh i xh i x 因此在点可以构造如下形式的差值 i x 节点的一阶向前差分 ii xfhxf i x 节点的一阶向后差分 hxfxf ii i x 节点的一阶中心差分 hxfhxf ii i x 与点相邻两代男的泰勒展开式可以写为 i x 2 4 3 2 4 3 2 iiiiii xf h xf h xf h xhfxfhxf 22 2 4 3 2 4 3 2 iiiiii xf h xf h xf h xhfxfhxf 23 2 22 2 23 并忽略 h 的平方和更高阶的项得到一阶微分的中心差商表示 2 h hxfhxf xf ii i 2 24 利用 2 22 和 2 23 式我们还可以得到一阶微分的向前 向后一阶差商表示 2 h xfhxf xf ii i 25 2 h hxfxf xf ii i 中国石油大学 华东 本科毕业设计 论文 14 26 将 2 22 和 2 23 式相加 忽略 h 的立方及更高阶的项得到二阶微分的中心差商 表示 2 2 2 h hxfxfhxf xf iii i 27 利用 2 24 2 27 式 我们就可以构造出微分方程的差分格式 这里要指出 的是 在构造差分格式时 究竟应该选择向前 向后还是中间差分或差商来代替微分 方程中的微分或微商 应当根据由此得到的差分方程解的稳定性和收敛性来考虑 同 时兼顾到差分格式的简单和求解的方便 上述差分步骤应用于偏微分 例如 对于的情况 拉普拉斯算符在 0 点作用在此函数上的值 yxff 也可以用临近的点上的函数值来表示出来 见图 2 1 1 且 2 2 2 2 2 y f x f f 时 hhhhh 4321 2 4 4 4 42 2 04321 2 4 24 y f x fh h fffff f 28 图 2 1 节点 0 及其附近节点 对微分方程数值求解的误差来源 1 方法误差 或截断误差 这时由于采用的计算方法所引起的误差 例如上 中国石油大学 华东 本科毕业设计 论文 15 面我们介绍的差商表示中 采用的泰勒展开式展开到第 n 1 项时的截断误差 1 n hO 具体方法的误差阶数取决于在离散化时的近似阶数 因此若改进算法就可以减少截断 误差 2 舍入误差 或计算误差 这是由于计算机的有限字长而造成数据在计算机 中的表示出现误差 在计算机运算的过程中 随着运算次数的增加舍入误差会积累的 很大 如果在多次运算后 舍入误差的精度影响是有限的 那么这个算法是稳定的 否则是不稳定的 不稳定的算法是不能用的 第 3 章 油水两相渗流机理和求解 3 1 数学模型的建立 假设条件如下 1 油藏中仅存在油水两相渗流 油水互不溶解 且各自符合达西定律 2 岩石 流体均可压缩 3 考虑掩饰的非均质性及各向异性 4 不考虑油水之间毛管力的影响 5 忽略重力 这里根据数学模型的一般式 经逐步简化 得到实验室进行单管模型的一维水驱油实 验室的数学模型 1 当考虑三维非均质油藏 油水互不相容 可压缩流体和岩石 考虑毛管力和重 力时 数学模型的一般式为 3 LLLLL L rL L S t qDg P KK w oL 1 简化到一维 并忽略重力项 3 LLL L L rL L S t q xu KK x w oL 中国石油大学 华东 本科毕业设计 论文 16 2 2 假设 不考虑掩饰的压缩性 即 常数 不考虑流体的体积变化 即 11Bo w B 油水粘度为常数 于是得 水相 3 t S q xu KK x o ov o o or 3 油相 3 t S q xu KK x o ov o o or 4 式中 为地面标准状况下单位时间内单元体中注入 或采出 的体积流 scl L lv q q 量 上述两个偏微分方程中的未知量有 4 个 即 因此还需要写出两个 w o w S o S 辅助方程 即 3 1 ow SS 5 3 wocow PPP 6 初始条件为 3 lx SS P 0 0 x 00 x wcw 7 边界条件为 3 8 0t qqq q qq q vovwvlxv vwv0 xv 中国石油大学 华东 本科毕业设计 论文 17 上述边界条件中 注入量产出量均为 表明该水驱油试验为稳定驱替 上述 3 1 v q 3 8 构成了该问题的完整的数学模型 利用数值方法进行求解后 可得到在不同的注入速率下 模型中任意一点的压力 饱和度随时间的分布和变化 3 2 数学模型的求解的方法及参数处理 对以上数学模型进行差分求解之前 这里首先对未知量的求解方法及有关参数的 处理进行说明 1 数学模型的求解方法 上述数学模型中有压力 和饱和度 两组未知量 本文应用隐式压力 o P w P w S o S 显式饱和度 IMPES 求解法进行求解 IMPES 方法的基本思路 1 通过乘以适当的系数 合并油方程和水方程 以消去微分方程组中的 w S 得到一个只含有 的方程 o S o S w P 2 由毛管力方程 可得 带入上面合并后的方程 得 wocow PPP cowow PPP 到一个只含有的方程 成为压力方程 o P 3 方程左端达西项系数上一时间段的值 同时毛管力也用上一时间阶段的值 即显示处理系数 于是可形成一个高阶现行代数方程组 用迭代法可以进行求解 先 求出 然后得 1n o P n c 1n o 1n w PPP 4 将带入水相方程 用显式方法求出 然后得 1n w P 1n w S 1n w 1n o 1 SS IMPES 方法具有所占内存小 计算工作量小 方法简便等优点 但该方法存在两个问题 第一 达西项的系数处理是显式的 因此对如锥进的问题 由于井底周围流速较高 压差变化大 而存在较大的误差 对于强非线性问题的适应性也差 第二 饱和度的 计算是隐式的 当时间步长较大时 会出现解的不稳定性 因此 IMPES 方法只适t 用于一般的弱非线性渗流的问题 对于某些非线性渗流的问题如注气 气锥或水锥等 问题 IMPES 方法无能为力 即使时间步长取得很小 仍会出现解的震荡或算出的压力 和饱和度为负值的情况 以致模拟计算无法正常进行 2 参数处理 在用有限差分对数学模型进行求解时 首先要将连续的油藏问题离散化为网格单 元 然后对每一个网格单元 读入包括深度 有效厚度 孔隙度 渗透率 饱和度等 中国石油大学 华东 本科毕业设计 论文 18 基本参数 所有给定的这些参数都是网格节点处的值 在两个网格节点中间处的参数 值是未知的 因此需要进行相应的处理 1 渗透率的取值K 渗透率 K 是空间函数 其取值由以下几种方法 算数平均 3 2 1 2 1 ii i KK K 9 加权平均 3 1ii 1i1iii 2 1i xx xKxK K 10 调和平均 3 11i 1i 2 1 x x xx iii i i KK K 11 几何平均 3 1ii2 1i KKK 12 2 相对渗透率的取值 相对渗透率的取之原则上是去流动方向上的上游节点值 通常称为上游权法 如 下图所示 取值方法是 当时 即由流向 时 i1 i 1i i 3 1rl2 1 irl iw SKK 13 当时 即由 流向时 1 ii i1i 3 14 iw SKK rl2 1 irl 式中 为第 个节点的势 i i 中国石油大学 华东 本科毕业设计 论文 19 图 3 1 井点分布图 采用上游权处理主要是对饱和度沿流向变化滞后进行修正 由于流动系数关于时间 的显式处理造成流动的滞后 如流动系数在空间上按相邻节点的算数平均值取值 则 会加重滞后现象 影响解的精度 因此 上游权的处理实质上是将显式处理造成的时 间上的滞后用空间上的向前来进行弥补 3 毛管力曲线的处理 由于油藏岩石的非均质性 即使用同一油层的岩心所测得的毛管压力曲线也有所 不同 用 J 函数方法可以从实验室提供的大量毛管力资料中选择合适的毛管力曲线 J 函数是把流体界面张力 岩石润湿性及渗透率和孔隙度等的影响综合在一起来表征油 层的毛管压力曲线特征的一个无因次函数 实践表明 它是处理毛管压力曲线的一种 有效方法 J 函数定义如下 3 2 1 cos 1000 KSP SJ owow wc w 15 式中 油 水界面张力 ow mmN 水的润湿角 ow 岩石毛管压力 wc SPMPa 渗透率 K 2 m 孔隙度 小数 与相对应的无因次量 w SJ w S 将众多的毛管压力资料按式 3 15 计算 做出关系图 如图 3 2 所示 然后拟 w SJ 合出一条关系曲线 当油层岩性较为接近时 关系并根据油层 值分 w SJ w SJ k 布的峰值区域可求出一组有代表性的关系 wwc SsP 在进行油藏数值模拟计算时 一般需要将油 水相对渗透率曲线和毛管压力曲线 离散化 然后在模拟计算过程中应用各种插值手段求出各网格在不同饱和度下的相对 中国石油大学 华东 本科毕业设计 论文 20 渗透率值和毛管压力值 图 3 2 曲线图 w SJ 3 3 差分方程组得建立及求解 下面用方法来建立差分方程组并进行求解 IMPES 隐式求压力 上面的数学模型式 3 3 至 3 6 中 有 4 个未知量 首先我们利 w S 0 S o P w P 用及方程 3 3 和 3 4 的关系消去 使其成为只含有压力 的1SS wo o S w S 0 P w P 方程 为此 式 3 3 式 3 4 得 3 0 ovwv o o row w rw qq x PKK xx PKK x 16 根据 可得 代入上式可得只含有压力的方程 这 wocow ppp cowow ppp 0 p 里为简化起见 假设毛管压力为 0 于是可得 因此上式简化为 ppp wo 3 17 0qq x p kk xx p KK x ovwv o ro w rw 令 分别表示油 水两相的流动系数 o w o ro o KK w rw w KK 表示总的流动系数 ow 表示油 水两相的总流量 v q vovwv qqq 则方程 3 17 可简化为 3 18 0q x P x v 方程 3 18 即为 IMPES 方法所化简得压力方程 下面写出该方程的隐式差分格式 中国石油大学 华东 本科毕业设计 论文 21 假设采用块中心网格 且网格大小相等 均为 并假设网格为水注入处 x 1 i 网格为油和油 水产处处 ni 式 3 18 的隐式差分方程为 3 19 0q x x PP x PP vi 1n 1i 1n i i 1n i 1n 1i i 2 1 2 1 下面分 3 种情况来讨论方程 3 19 1 对于第 2 个至第个网格 既无注入也无采出 因此方程 3 18 中 1N 0 vi q 于是可简化为 3 20 0 PP PP 1n 1i 1n i i 1n i 1n 1i i 2 1 2 1 若系数按上游权原则取值 并用 n 时刻的值 即显示处理 则上式方程可简化为 3 21 0 PP PP 1n 1i 1n i n 1i 1n i 1n 1i n i 即 3 22 0PP P 1n 1i n i 1n i n i n 1i 1n 1i n 1i 2 对于第 1 个网格 单位体积中注入的体积流量为 式 3 18 中第二项由于取 vi q 上游权 于是可简化为 0 0 n 3 23 0q x PP v 2 1n 1 1n 2 n 1 整理上式 可得 3 24 n 1 v 2 1n 2 1n 1 qx PP 两端同乘以 网格单元的体积 并令 该网格的注入量 xA xAqQ vv 则上式可化简为 3 25 A xQ PP n 1 v 1n 2 1n 1 3 对于第 N 个网格 其单位体积中采出的体积流量为 式 3 18 中第一项由于 v q 取上游权 于是可简化为 0 n i 3 26 0q x PP v 2 1n 1N 1n N n 1N 两端同乘以 令 得 xA xAqQ vv A xQ PP n 1N v 1n N 1n 1N 中国石油大学 华东 本科毕业设计 论文 22 式 3 22 式 3 25 和式 3 27 构成了从第一个网格到最后一个网格 第个 N 的线性代数方程组 矩阵方程如下 1 1 1 1 1N12N2N 3322 2211 i A xWI N 1N 3 2 1 P P P P P A xWIP 0 0 0 A xWIP wf wf 3 28 方程组的系数矩阵为三角矩阵 若用定流量来进行求解 则上述压力是与时间无 关的量 不符合要求 可用方法进行求解 从而得到 在Thomas 1 P 2 P N P 时刻的压力值 1n 2 显式求饱和度 数学模型中的水相方程 3 3 的差分方程为 3 t S S q x x P P x P P n wi 1n wi wvi 1n 1 i 1n i n 1 wi 1n i 1n 1i n wi 29 上式中的未知量为 也可分为以下三种情况来讨论 1n wi S 1 对于第 2 个至第个网格 既无注入也无采出 因此方程 3 28 中的1 N 0qwvi 于是可得 3 1n 1 i 1n i n 1 wi 1n i 1n 1i n wi 2 n 1w 1n 1w PPPP x t SS 30 由于已经通过隐式求压力得到 故可用上式顺序 求出每个网 1n i P 1 n 2 i 格节点的饱和度值 2 对于第一个网格 即单位体积的注水速率 而且式 3 29 中的左面查 vwvi qq 分的第二项去上游权 于是可化简为 0 n 0w 中国石油大学 华东 本科毕业设计 论文 23 3 t SS q x PP n 1w 1n 1w v 2 1n 1 1n 2 n 1w 31 两边同乘以并化简 且 则有 xA xAqQ vv t SSx Q x PP n 1w 1n 1w v 1n 1 1n 2 n 1w 于是 3 x PP A Q x t SS 1n 2 1n 1 n 1w v n 1w 1n 1w 32 3 对于最后一个 第 N 个 网格 其单位体积的产出量 实验中为 ovwvv qqq v q 定值 其中 产水量为 而且 是跟随时间变化的 因此首先要求出的值 wv q wv q wv q 最后一个 第 N 个 网格的压力为 第 N 1 个网格的压力为 则对于 N P 1 N P N1 N P P 油和水来说均是相等的 故流出油 水提及的流量取决于流动系数 的大小 o w 3 1n N 1n 1Noov PPq 33 3 1n N 1n 1Nwwv PPq 34 所以 3 1n N 1n 1Nwov PPq 35 于是 3 wo w v wv q q 36 方程 3 28 中的左面差分的第一项由于取上游权 于是可化简为 0 n wn 3 t SS q x PP n wN 1n wN wv 2 1n 1N 1n N n 1 wN 37 两边同乘以 并令 则 xA xAqQ wvwv 3 t SSx Q x PP n wN 1n wN wv 1n N 1n 1N n 1 wN 38 于是 中国石油大学 华东 本科毕业设计 论文 24 3 n 1N ow wv 1n N 1n 1 N n 1 wN n Nw 1n Nw A Q x PP x t SS 39 利用式 3 30 式 3 32 和式 3 39 可求得网格的 值 N 1 2 i 1n wi S 以上所叙述的隐式求压力 显式求饱和度的计算过程即为 IMPES 方法从时刻到n 时刻的求解过程 这样从初始时刻开始 一步步依次求解下去 直到求得所要求1n 的时间的压力和饱和度为止 油水两相的渗流方程就可得解 所建立的油藏模型就得 以模拟出来 3 4 有关单位换算 以上的计算过程中需用到达西公式 但达西公式中所用的单位为水力学单位 而 实际油田中常用的单位为工程单位 二者之间的换算关系见表 3 1 表 3 1 换算关系 水力单位工程单位 Q s cm3d cm3 K Dum20 001 Dum2 usmPa smPa A 2 cm 2 m x cm m P 0 1MPa MPa 1 达西公式为 其中各变量均采用水力学单位 x P A K Q 2 当采用工程单位时 达西公式变为 2 4 3 6 10 x 10P 10A 10K 86400 10 Q 3 因此 由于水力学单位换算到工程单位 达西公式前加了一个系数 即 x P A K 0 0864Q 中国石油大学 华东 本科毕业设计 论文 25 3 5 计算程序框图 程序开始 读入 xLK wo 读入N2 1iiKiS rww 读入 iwc PS 0t 读入 vmax Qxt 计算不同下的 w S rw K ro K 计算值 ow 计算压力方程的系数矩阵和右边项 idibiaic 用追赶法 iP 显式求 iS iS ow 中国石油大学 华东 本科毕业设计 论文 26 图 3 3 一维油水两相程序框图 是 打印 iSiSiPt ow 结束 max tt 否 中国石油大学 华东 本科毕业设计 论文 21 3 6 计算实例 下面把给定的油藏区块按一维方向均分成 7 个网格井点 以第一个网格井为注水 井 第七个网格井为采油井 一注一采模型进行求解 注水井 生产井 图 3 4 井网分布图 已知 孔隙度 渗透率 束缚水饱和度 原始含油饱0 25 2 m1K 2 0Swc 和度 地层油粘度 地层水粘度 注水井井底压力8 0Soi smPa5 o smPa1 w 生产井井底流压 有谁渗透率曲线如表 1 所示 油层长MPa15Piwf MPa10Pwf 70m 宽 油层厚度 综合压缩系数 取 时间m7010mh 14MPa 101C x 10m 步长取一天 模拟时间为 油水相对渗透率曲线见附录 模拟曲线如下 1360d 0500100015002000 0 0 2 0 4 0 6 0 8 1 t d fw 图 3 5 最后一个网格点 生产井 关系曲线tfw 总液量中水所占有的分量即为含水率 此曲线反映了油井的含水率随时间的变化 w f 趋势 开始后一段时间内含水率接近 0 因为水驱前沿未到达产油井 采出的完全是油 水驱前沿到达产油井后 含水率开始上升 经过相当长时间后含水率近似可以达到 1 此时油井开采完毕 根据此曲线可以综合评估油井的成本与效益 以便选择合适的开 采时间 1 P 2 P 3 P 4 P 5 P6 P 7 P 中国石油大学 华东 本科毕业设计 论文 22 1234567 11 11 5 12 12 5 13 13 5 14 14 5 i P Mpa 图 3 6 1 7 网格井点在 T 时刻的压力曲线 此图反映了此模型从第一个井点到第七个井点的压力变化 可以看出是等差变化 的 这是因为不考虑油藏的毛管力 重力以及岩石和流体的压缩性所致 模型为刚性 水驱开采 故从注水井至开采井的压力是等差变化的 020040060080010001200140016001800 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 t d Sw 图 3 7 最后一个网格井点的曲线tSw 此曲线反映了最后一口网格井的随时间的变化趋势 随着水驱前沿的推进 水 w S 侵占油藏的孔隙以排除原油 所以油藏的含水饱和度会随时间而增大 最后达到稳定 值 因为此时油藏孔隙中的残余油不能被排出 仍然占据着油藏孔隙 所以随着时间 的推移 含水饱和度不再变化 中国石油大学 华东 本科毕业设计 论文 23 第 4 章 结 论 1 本文应用隐式压力 显式饱和度解法 IMPES 方法 来对油藏油水两相一维 渗流进行模型建立与求解 通过合并流体方程得到一个只含有压力的方程 某一时间 步的压力求解出来后 饱和度采用显式更新对油藏的模拟及求解 2 模型的求解有定流量和定压力两种方法 本文建立的模型不能用定流量来进行 求解 因为若用定流量进行求解则此模型的压力方程与时间无关 会有无穷多解 这 是因为只要保持注水井和采油井的压力恒定 就可以达到流量恒定 所以这样的压力 解是无穷的 3 本文所用的 IMPES 方法具有占内存小 计算工作量小 方法简便等优点 4 本文建立的模型是在不考虑岩石和流体的压缩性以及油水的粘度 忽略重力 在这种刚性模型下对油藏的模拟准确度可能有一定的影响 IMPES 方法达西项的系数处 理是显式的 因此对如锥进的问题 由于井底周围流速高 压差变化大 而存在较大 误差 对于强非线性问题的适应性也差 饱和度的计算是显式的 当时间步长较大时 会出现解的不稳定性 5 通过含水率曲线表明 相当长时间后 含水率接近 100 此时此油带的水驱油 开采完毕 6 本文对一维油水两相渗流进行了数值模拟 对一维各个以划分的网格在个时刻 的压力 含水饱和度 含水率等方面都有了一个全面的显示 能够对油藏较为准确的 模拟 按此方法可以对二维 三维油水两相以及二维 三维油气水三相渗流进行拓展 中国石油大学 华东 本科毕业设计 论文 24 致 谢 非常感谢张晓东老师在我大学的最后学习阶段 毕业设计阶段给自己的指导 从最初的定题 到资料收集 到写作 修改 到论文定稿 他给了我耐心的指导和无 私的帮助 为了指导我们的毕业论文 他放弃了自己的休息时间 他的这种无私奉献 的敬业精神令人钦佩 在此我向他表示我诚挚的谢意 同时 感谢所有任课老师和所 有同学在这四年来给自己的指导和帮助 是他们教会了我专业知识 教会了我如何学 习 教会了我如何做人 正是由于他们 我才能在各方面取得显著的进步 在此向他 们表示我由衷的谢意 并祝所有的老师培养出越来越多的优秀人才 桃李满天下 通过这一阶段的努力 我的毕业论文 基于有限差分法的油水两相渗流方程求解 终 于完成了 这意味着大学生活即将结束 在大学阶段 我在学习上和思想上都受益非 浅 这除了自身的努力外 与各位老师 同学和朋友的关心 支持和鼓励是分不开的 在本论文的写作过程中 我的导师张晓东老师倾注了大量的心血 从选题到开题报告 从写作提纲 到一遍又一遍地指出每稿中的具体问题 严格把关 循循善诱 在此我 表示衷心感谢 同时我还要感谢在我学习期间给我极大关心和支持的各位老师 写作毕业论文是一次再系统学习的过程 毕业论文的完成 同样也意味着新的学习生 活的开始 我将铭记我曾是一名中国石油大学学子 在今后的工作中把惟真惟实的优 良传统发扬光大 感谢各位专家的批评指导 中国石油大学 华东 本科毕业设计 论文 25 参考文献 1 哈利德 阿齐兹 安东尼 塞特瑞 油藏数值模拟 第一版 北京 石油工业 出版社 2004 2 Henrik Lof