油藏数值模拟原理黑油模型.ppt

- 1 -,2006年4月13日,李建芳,油藏数值模拟原理,黑油模型,- 2 -,黑油模型简介,黑油模型是指油质较重性质的油藏类型。 黑油模型是最完善、最成熟，也是应用最为广泛的模型。 是油藏数值模拟的基础，其它模型大都是黑油模型的扩展。 目的：软件这样一个“黑箱”是如何操作的，以便指导我们更好地去应用软件。,- 3 -,黑油模型简介,黑油模型的基本假设 单相渗流黑油模型基本微分方程 三维三相黑油模型基本微分方程 初边值条件 5. 网格系统 6. 黑油模型的差分方程 7. 差分方程的线性化 8. 线性代数方程组的求解 9. 黑油模型的主要数据流,- 4 -,1 、黑油模型基本假设,（1） 油藏中的渗流是等温渗流。 （2） 油藏中最多只有油、气、水三相，每一相均遵守达西定律。 （3）油藏烃类只含有油、气两个组分。在油藏状态下，油气两组分可能形成油气两相，油组分完全存在于油相内，气组分则可以以自由气的方式存在于气相中，也可以以溶解气的方式存在于油相中，所以地层内油相为油组分和气组分的某种组合。在常规油田中，一般不考虑油组分向气组分挥发的现象。 （4）油藏中气体的溶解和逸出是瞬间完成的，即认为油藏中油气两相瞬时达到相平衡状态。 （5）油水之间不互溶；天然气也假定不溶于水。,- 5 -,2.单相渗流基本微分方程,单相渗流黑油模型基本微分方程 黑油模型基本微分方程主要是由达西公式(达西定 律)和连续性方程(物质守恒定律)所构成 一维单相渗流的Darcy公式为： 三维单相渗流的Darcy公式为：,- 6 -,渗流速度是一个空间向量，用分量表示为： 渗流速度矢量； 压力； 势函数 绝对渗透率； 重力加速度常数； 粘度； 垂向坐标； 密度；,2. 单相渗流基本微分方程,- 7 -,Reservoir Simulation III,Gringarten, 2002,Data from all sources,- 8 -,2.单相渗流基本微分方程,根据物质守恒原理，在单位时间内 流入单元内的流体质量-流出单元的流体质量 =单元内流体质量的变化 取渗流场中一个微小六面体体积单元来研究,- 9 -,- 10 -,2.单相渗流基本微分方程,该六面体单元的中心点坐标为（x，y，z），其 每一侧面的质量流速均以其侧面的中心点的质量 来代替。 方向流入质量= 方向流出质量= 方向流入质量= 方向流出质量= 方向流入质量= 方向流出质量= 单元体内流体质量的变化=,- 11 -,根据物质守恒定律： 将上式两边同除以 ，,2.单相渗流基本微分方程,- 12 -,并令 ， ， 把源汇项（井）加入到方程中，得到 单相渗流连续性方程 ：,2.单相渗流基本微分方程,- 13 -,用Hamilton算子 表示，并将达西公式代入，得 到，单相渗流微分方程（流动方程）： 生产井 注入井 关 井,2.单相渗流基本微分方程,流动项,源汇项,累积项,- 14 -,这就是地下流体所遵循的渗流规律。 总之，渗流基本微分方程包含三部分内容: 流动项、源汇项以及累积项 渗流基本微分方程由两个定律控制： 达西定律与物质守恒定律 压力场和饱和度场的分布是流体流动的结果，流 动规律由达西定律来确定，流量由物质守恒来确 定，这两个规律是数值模拟的关键所在。,2.单相渗流基本微分方程,- 15 -,将单相渗流达西定律推广到多相渗流中，得到 三相流动达西公式： 式中，w，o，g水、油、气诸相下标。 相对渗透率；,3. 三维三相渗流基本微分方程,- 16 -,对油、水、气三相，每一相运用物质守恒定律，可以建立以下连续性方程,3. 三维三相渗流基本微分方程,- 17 -,把达西公式代入连续性方程，得到 三维三相渗流基本微分方程,3. 三维三相渗流基本微分方程,相对渗透率、毛管力是饱和度函数,而粘度、体积系数孔隙度与密度均是压力的函数,- 18 -,3. 三维三相渗流基本微分方程,令,- 19 -,上述T的物理意义是反映相邻网格间的流动性，又叫传导率（transmissibility),3. 三维三相渗流基本微分方程,- 20 -,3. 三维三相渗流基本微分方程,求解变量,- 21 -,3. 三维三相渗流基本微分方程,还要补充一些方程 饱和度 ： 油水毛管力： 油气毛管力： 相对渗透率和毛管力是饱和度函数：,- 22 -,Stone公式： 式中， 束缚水饱和度 油水两相流动时油相相对渗透率 油气两相流动时油相相对渗透率,3. 三维三相渗流基本微分方程,- 23 -,3. 三维三相渗流基本微分方程,流体PVT性质 水相： 气相： （Z) 油相：随压力变化，相态会发生变化。 三相状态： 两相状态： 岩石的高压物性 相渗曲线以及流体PVT性质，以表格形式给出， 称为静态数据。,- 24 -,3. 三维三相渗流基本微分方程,实际求解变量 初始化输入数据ljfinit.txt,油相流体性质与压力和饱和压力的关系,- 25 -,产量项（作源汇项处理）（动态数据）,3. 三维三相渗流基本微分方程,- 26 -,3. 三维三相渗流基本微分方程,其中， 在井基准面深度 上的井底压力。 节点等效泄油半径。 井筒半径。 X，Y，Z 网格节点x，y和z方向步长 PI叫采油指数， 通常写成KH,可以直接输入,井数据又叫动态数据ljfmodel.txt,- 27 -,4 .初边值条件,油气藏外边界条件:,第一类边界 条件 P=f(x,y,z,t) 第二类边界 条件 第三类边界条件:,- 28 -,4 .初边值条件,外边界 封闭边界:尖灭、断层、圈闭 流动边界：边、底水 内边界，即井的生产条件 定产边界条件（rate constraint):给定产液量、产油量、产气量、注水量、注气量 定压边界条件(pressure constraint):给定井底流压,- 29 -,4 .初边值条件,油藏初始条件:,- 30 -,5. 网格系统,- 31 -,在一般情况下流动方程求不出解析解，只能用离散化方法求数值解。 离散化方法，首先要把求解区域按一定的网格系统进行剖分。网格剖分方法有很多，像直角坐标网格、径向网格、角点网格、PEBI网格等。目前，商用软件中大都采用矩形网格系统。离散化方法主要有：有限差分（Finite Difference)、有限元(Finite Element)以及流线法(Stream Line)。比较成熟的是有限差分方法。有限差分方法大多空间变量采用网格块中心（Block Center)差分, 时间变量按向前差分方法。,5. 网格系统,- 32 -,5.网格系统:矩形网格,- 33 -,- 34 -,5.网格系统:径向网格系统,读取网格数据时，R方向循环最快，然后是 方向，最后是Z方向,径向网格数据读取规则,- 35 -,5.网格系统:角点网格,Corner Point,- 36 -,块中心网格与角点网格的比较,块中心网格表示的含有断层的倾斜构造,角点网格表示的含有断层的倾斜构造,- 37 -,5.网格系统:块中心网格的数据结构,模型被分为20510个网格区域。 网格从（1，1，1）开始向下倾斜。 各网格块在X方向上是300ft，在Y方向上是1000ft。 层厚度从上至下分别是32，22，20，4，32，4，26，26，4，28ft,- 38 -,5.网格系统:块中心网格的数据结构,-前20个TOPS定义网格 (1, 1, 1) 到 (20, 1, 1) 6855.000 6865.000 6875.000 6885.000 6895.000 6905.000 6915.000 6925.000 6935.000 6945.000 7005.000 7015.000 7025.000 7035.000 7045.000 7055.000 7065.000 7075.000 7085.000 7095.000 -下面20个TOPS定义网格 (1, 2, 1) 到 (20, 2, 1) 6930.000 6940.000 6950.000 6960.000 6970.000 6980.000 6990.000 7000.000 7010.000 7020.000 7080.000 7090.000 7100.000 7110.000 7120.000 7130.000 7140.000 7150.000 7160.000 7170.000 -下面20个TOPS定义网格 (1, 3, 1) 到 (20, 3, 1) 7030.000 7040.000 7050.000 7060.000 7070.000 7080.000 7090.000 7100.000 7110.000 7120.000 7180.000 7190.000 7200.000 7210.000 7220.000 7230.000 7240.000 7250.000 7260.000 7270.000 -下面20个TOPS定义网格 (1, 4, 1) 到 (20, 4, 1) 7130.000 7140.000 7150.000 7160.000 7170.000 7180.000 7190.000 7200.000 7210.000 7220.000 7280.000 7290.000 7300.000 7310.000 7320.000 7330.000 7340.000 7350.000 7360.000 7370.000 -下面20个TOPS定义网格 (1, 5, 1) 到 (20, 5, 1) 7205.000 7215.000 7225.000 7235.000 7245.000 7255.000 7265.000 7275.000 7285.000 7295.000 7355.000 7365.000 7375.000 7385.000 7395.000 7405.000 7415.000 7425.000 7435.000 7445.000 / 这完成了第一层TOPS的定义,- 39 -,角点网格的数据结构,-COORD -此关键字定义坐标线 -X1 Y1 Z1 X2 Y2 Z2 0. 0. 6825.000 0. 0. 7023.000 300.0000 0. 6835.000 300.0000 0. 7033.000 600.0000 0. 6845.000 600.0000 0. 7043.000 900.0000 0. 6855.000 900.0000 0. 7053.000 1200.000 0. 6865.000 1200.000 0. 7063.000 /对于20 * 10 的网格模型, 需要21*11 条坐标线，即 231条。-每一条坐标线都用6个数来定义，所以关键字COORD包含有1386个数。 ZCORN -这个关键字定义网格拐点的深度，数据按照X (或 R) 方向循环最快，然后是Y (或 THETA) 方向，最后是Z方向的顺序排列。 6825.000 6835.000 6835.000 6845.000 6845.000 6855.000 6855.000 6865.000 6865.000 6875.000 6875.000 6885.000 6885.000 6895.000 6895.000 6905.000 6905.000 6915.000 6915.000 6925.000 6975.000 6985.000 6985.000 6995.000 6995.000 7005.000 7005.000 7015.000 7015.000 7025.000 7025.000 7035.000 7035.000 7045.000 7045.000 7055.000 7055.000 7065.000 7065.000 7075.000 . / 对于1000个网格，需要 8000 个ZCORN值 -这些只是角点网格模型的数据体中很少的一部分,- 40 -,块中心网格与角点网格的比较,- 41 -,正常连接的网格关系,粗细网格交界处的流量计算,非正常连接技术,两点流量公式：,- 42 -,非结构网格技术,传统的结构网格-即规则矩形剖分的网格,具有如下不足: 存在较严重的网格取向效应 不能保证每口井都在网格中心 不能很好地描述油藏形态 粗细网格间的过渡衔接性差 研究新的网格技术势在必行.,结构网格的不足,x,y,好,差,井,好,开发软件,- 43 -,非结构网格技术,Gurpinar, 2001,Castellini, 2001,Wolfsteiner et al., 2002,Prevost 2003,- 44 -,非结构网格技术(Unstructured Grids),常用的非结构网格 各网格中心之间连线构成三角网 PEBI网格（又叫垂直平分网格Perpendicular Bisectors):各三角形边的垂直平分线构成网格边界 CVFE网格-控制体有限元网格(Control volume finite element):各三角形边的中点与重心的连线构成网格边界 待解决的问题 非结构网格条件下代数方程组的快速求解,CVFE网格,PEBI网格,开发软件,- 45 -,5.网格系统: PEBI网格,PEBI,- 46 -,unstructured PEBI grids,Locally different grids,- 47 -,Advanced Fault Modeling,Heavily Faulted Reservoir,- 48 -,非结构网格技术 Channels can be modelled,- 49 -,2. 非结构网格技术 (complex) faults,Geological Model of a Heavily Faulted Reservoir,Simulation Model,- 50 -,1. 概述 局部网格加密(Local Grid Refinement，简称LGR)技术主要分： 静态局部网格加密(Static Local Grid Refinement)局部加密区域在整个计算过程中固定不变的，常用于断层、尖灭、裂缝以及井附近的加密。该技术已经成熟，在多加软件公司的软件产品中都有这一功能。 动态局部网格加密(Dynamic Local Grid Refinement)局部加密区域在整个计算过程中随流体的运动位置而变化，主要用于追踪驱替前缘。该方法在研究之中，最早是Heinemann将动态局部网格用于水驱油藏数值模拟 。,局部网格加密(Local Grid Refinement),- 51 -,Productive Area 1,Productive Area 2,Sealing Interbedding,Aquifer,Aquifer,局部网格加密技术,- 52 -,局部网格加密技术,Aquifer,Productive Area,Vertical Refinement Zone,- 53 -,随驱替前缘移动的动态局部网格加密系统,化学驱驱替前缘动态追踪数值模拟研究,- 54 -,1.3 动态局部网格加密技术 (注化学剂时的So）,开发软件,- 55 -,杂交网格技术(Hybrid Grid Refinement),Horizontal Well Grid,Computed Results,- 56 -,开窗技术(Windowing Technique),Independent time stepping At any location Several different gridding options are available within one model in the same position,- 57 -,开窗技术,- 58 -,开窗技术,Automatically and implicitly calculated in/out-flux will be used as a von Neumann boundary condition,- 59 -,开窗技术,Good match required windowing,- 60 -,差分方程,6. 差分方程(Difference Equations),- 61 -,6. 差分方程,- 62 -,该差分方程组无论是左端还是右端都是非线性的，是一组非线性代数方程组。非线性最强的参数是相对渗透率和毛管力。这两个参数在时间方面的处理不同，对计算的稳定性以及计算结果影响最大，也是形成形形色色解法的关键所在。 根据隐式程度的不同，数值模拟的算法有IMPES (Implicit Pressure Explicit Saturation) 方法、半隐式(Semi-Implicit Method)方法和全隐式方法(Fully Implicit Method )。 IMPES：隐式求解压力显式求解饱和度。对流度中与时间有关的量均取上一时间步的值。,6. 差分方程,- 63 -,IMPES方法的优点是速度快，缺点是稳定性差,尤其是锥进问题以及过泡点问题不能很好解决。因此时间步长不能太大。目前，IMPES方法在化学驱模型，尤其是流线法中发挥着作用。 全隐式方法主要采用了Newton-Raphson迭代技术, 对方程中各项系数，包括井系数全部隐式处理。全隐式方法较之半隐式方法，具有更好的稳定性，时间步长更长，彻底解决了过泡点问题，是八十年代数值模拟技术走向成熟的关键技术，也是大型油藏数值模拟的主流方法。,6. 差分方程,- 64 -,全隐式方法和半隐式方法的根本差别在于： 半隐式方法是将非线性方程直接线性化形成 线性代数方程组，来求解n+1时刻的值。 全隐式方法采用Newton-Raphson方法，迭 代地求解非线性方程组，得到n+1时刻的变量。 用一维非线性方程表示： 给定初始近似 对函数 在点处 进行一阶泰勒展开, 可近 似地表示为,6. 差分方程,- 65 -,半隐式： 全隐式方法：,6. 差分方程,- 66 -,下面介绍全隐式方法。为了减小计算机舍入误 差的影响， 模型中不直接求压力和饱和度，而是求其增量。 其中：x代表求解变量 n时间迭代步 牛顿迭代步,6. 差分方程,- 67 -,求解变量： ， ， 对于二相状态： 对于三相状态：,6. 差分方程,- 68 -,对方程组进行线性化处理。对所有变量进行Taylor展开，保留线性部分。变量进行运算后也保留线性部分，舍去高阶项。形成线性代数方程组。,7. 差分方程的线性化,- 69 -,全隐式方法形成的线性代数方程组的结构,7. 差分方程的线性化,矩阵结构示意图,- 70 -,油藏数值模拟的线性代数方程组是大型稀疏矩阵，其求解效率在油藏数值模拟计算中举足轻重。因为，其计算量占整个计算量约70%80%。因此，研制快速求解大型稀疏矩阵的新方法成为油藏数值模拟发展的一个重要方向。 线性代数方程组的求解方法主要有三类： 1）直接法：Gauss消元、D4排序的Gauss消元 2）迭代法：SOR(松弛迭代法） 3）共轭梯度类法、预处理共轭梯度类法（PCG方法）,8. 线性代数方程组的求解,- 71 -,目前最先进的方法是预处理共轭梯度类方法。BO软件中采用RSVP方法。该方法计算快，精度高。,8. 线性代数方程组的求解,- 72 -,黑油模型油藏数值模拟器