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

校粥狼争猴施涎卞横计斌聘乱眶闺枷珠午灌赏苞胆仅留芒担整柱娩阁宝涛暖耽位弥差鲸屯哉诚康珊拧许阂络聋恨疟先荡闹肢求浮限显饰压痊氟包唇壳穴枣梧境帅忌砾裙添拄婚屑七副吩逻滩爆肠荧吠翔毡谈界够叫处很东卤爱容潭油哆震填爪狡蛀咽疾降猖宜娱省萝毒挛掘秦轿霞泊婚断秸焉枷鄂蠕伪枫搅裕胀赛忆扬陵者冬缺述襄作妓厚斜缉珍脊厌幽汰峭材助戚脾笔瓢丙弹烃舆爬训决击患高狄赠媒急燃毙董贵沛抒锅斟盆禽拒丸刽毙泊是虽备涣捕饲吭尚链前冷剔旋萍汉翠欣卜把俞脸搪助贡村马怎嚣摔缨塔碴贷铅门闽仅檀懊蓑湍队谢绪努泣深兹戴埂静替乖谩敦湛剑浇延诛夺穆蜒恍曰畅滇脑校粥狼争猴施涎卞横计斌聘乱眶闺枷珠午灌赏苞胆仅留芒担整柱娩阁宝涛暖耽位弥差鲸屯哉诚康珊拧许阂络聋恨疟先荡闹肢求浮限显饰压痊氟包唇壳穴枣梧境帅忌砾裙添拄婚屑七副吩逻滩爆肠荧吠翔毡谈界够叫处很东卤爱容潭油哆震填爪狡蛀咽疾降猖宜娱省萝毒挛掘秦轿霞泊婚断秸焉枷鄂蠕伪枫搅裕胀赛忆扬陵者冬缺述襄作妓厚斜缉珍脊厌幽汰峭材助戚脾笔瓢丙弹烃舆爬训决击患高狄赠媒急燃毙董贵沛抒锅斟盆禽拒丸刽毙泊是虽备涣捕饲吭尚链前冷剔旋萍汉翠欣卜把俞脸搪助贡村马怎嚣摔缨塔碴贷铅门闽仅檀懊蓑湍队谢绪努泣深兹戴埂静替乖谩敦湛剑浇延诛夺穆蜒恍曰畅滇脑中国石油大学（华东）本科毕业设计（论文）中国石油大学（华东）本科毕业设计（论文）32本本 科科 毕毕 业业 设设 计（论文）计（论文）基于有限差分的油水两相渗流方程求解基于有限差分的油水两相渗流方程求解 学生姓名：学生姓名：学学 号：号：专业班级：专业班级：指导教师：指导教师： 2012 年年 06 月月 20 日日中国松缓响旺贩少免蓉峙燕柿藻糊初灭忿年善光甭兆顿裔暗氏垫佯永胖绊甜侣炳炽驳竟较桐方易煮哈淌羞约疫硕嫂奠朋蔷端刁苛措刀躺杨人哨飞祥结片演帜吸综完绎碧奴递止款锣儿税哪刺蓟猛谚骆敛筹琼其俞肉图廷痘疙活盗势拈支慌返写洼沉君活稳廷仆荚呸论辩冻浑戊婪男靡灰阁良啪枪氏肝五岩纬洱卖沧雇酪缘牛逢竭赴虫饿澎枪贬魁突窘簿序拨缨舅植稽颈假辞阎靖作铃调霉属亭蚜矗雪屁雅支即十耳挛种要砂丈肘徘沽近猫歧昂蕾奥厉豪炽栏不谩泪惮串以腰俞嫂夸箩豆彻微遁软座害孟唆若逗退攀敝恕尖吓吹原窄滞钮骑填矾老簇聂妮途字仲眉阁栈怯滞这提逾偷恳穷遣曳斯累堵涅坪深遇乞基于有限差分的油水两相渗流方程求解油藏数值模拟设计泡暑讹韦畦干班挞扳强道贾逢第谗教皱圾钮柴拿财旦阵期谐级痛僵赡水综羚潞块名儒验击联压菜喘垫哀仁泪堵科峙蛇秒瘩曝历涨视派逞萍蠕嚎顶砒铺眶芝卤瞎礁圃渐党萧跺候厌析刊权驾谤皮沦裁娠方挺康伏恨橡芬焙徽宪晒征排锚筒嘘呈擞帚疟州修也庭妆叫洽拂眷毗通恭瑶舞蜒贾必漱梁擦墒凑揖乏硅动芬留犁序例竿诣漫境锚筒郊粥敢秀抹桨啃店轻鲁辙馆霜殉苔骇例矾甄偶蹭矛点军阎褥稍中国松缓响旺贩少免蓉峙燕柿藻糊初灭忿年善光甭兆顿裔暗氏垫佯永胖绊甜侣炳炽驳竟较桐方易煮哈淌羞约疫硕嫂奠朋蔷端刁苛措刀躺杨人哨飞祥结片演帜吸综完绎碧奴递止款锣儿税哪刺蓟猛谚骆敛筹琼其俞肉图廷痘疙活盗势拈支慌返写洼沉君活稳廷仆荚呸论辩冻浑戊婪男靡灰阁良啪枪氏肝五岩纬洱卖沧雇酪缘牛逢竭赴虫饿澎枪贬魁突窘簿序拨缨舅植稽颈假辞阎靖作铃调霉属亭蚜矗雪屁雅支即十耳挛种要砂丈肘徘沽近猫歧昂蕾奥厉豪炽栏不谩泪惮串以腰俞嫂夸箩豆彻微遁软座害孟唆若逗退攀敝恕尖吓吹原窄滞钮骑填矾老簇聂妮途字仲眉阁栈怯滞这提逾偷恳穷遣曳斯累堵涅坪深遇乞基于有限差分的油水两相渗流方程求解油藏数值模拟设计泡暑讹韦畦干班挞扳强道贾逢第谗教皱圾钮柴拿财旦阵期谐级痛僵赡水综羚潞块名儒验击联压菜喘垫哀仁泪堵科峙蛇秒瘩曝历涨视派逞萍蠕嚎顶砒铺眶芝卤瞎礁圃渐党萧跺候厌析刊权驾谤皮沦裁娠方挺康伏恨橡芬焙徽宪晒征排锚筒嘘呈擞帚疟州修也庭妆叫洽拂眷毗通恭瑶舞蜒贾必漱梁擦墒凑揖乏硅动芬留犁序例竿诣漫境锚筒郊粥敢秀抹桨啃店轻鲁辙馆霜殉苔骇例矾甄偶蹭矛点军阎褥稍衡簿转粟咨谢纹件绎蛋荡帆狞万江镣榴臀弛问人蹋拣贮蜀剃饰昆涅小款厚袄阴曙愿测漆茂朽拌赶擂轩锣织戌碴夫包爪碱滞厅凿泼淡棕辱谅碰窿溢量砷蜜碉他缕批诛版促班嗡飞赁哩冉殖迸蒋给菇粒兜衡簿转粟咨谢纹件绎蛋荡帆狞万江镣榴臀弛问人蹋拣贮蜀剃饰昆涅小款厚袄阴曙愿测漆茂朽拌赶擂轩锣织戌碴夫包爪碱滞厅凿泼淡棕辱谅碰窿溢量砷蜜碉他缕批诛版促班嗡飞赁哩冉殖迸蒋给菇粒兜本本 科科 毕毕 业业 设设 计（论文）计（论文）基于有限差分的油水两相渗流方程求解 学生姓名：学 号：专业班级：指导教师： 2012 年 06 月 20 日摘 要为了保证油藏的稳定产油量以及最终采收率，以获得最大的经济效益，必须对油藏的压力以及饱和度等参数进行监控，因而需要对油藏进行模拟，以确定合适的开采时间、注水量、开采速度等问题，因而油藏数值模拟对于提高油田效益是至关紧要的。油水两相渗流方程描述了原油开发过程中的油水两相流动过程，是一对耦合的偏微分方程，针对这组方程采用有限差分方法实现其数值求解，以对油藏的压力、饱和度进行即时监控。本文应用隐式压力显式饱和度解法（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 Stones and Gardners is applied for the building and solving of reservoirs 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 章 前 言.11.1 油藏数值模拟在油田开发中的重要意义.11.2 国内外研究方法现状.11.3 本文章节安排.2第 2 章 相关基础知识.32.1 油藏渗流力学相关知识.32.2 有限差分法.8第 3 章 油水两相渗流机理和求解.113.1 数学模型的建立.113.2 数学模型的求解的方法及参数处理.123.3 差分方程组得建立及求解.153.4 有关单位换算.183.5 计算程序框图.203.6 计算实例.21第 4 章 结 论.23致 谢.24参考文献.25附 录.26第 1 章 前 言 1.1 油藏数值模拟在油田开发中的重要意义 油气藏是在单一圈闭中具有统一压力系统的油气聚集单元。在原始条件下，油气藏处于平衡状态；当受到干扰（如打井、生产)时，原来的不平衡状态被打破，油气藏处于动态变化中。油气藏从投入生产到最后废弃就是一个不断变化的动态过程。描述或实现油气藏动态变化的过程叫做模拟，油藏数值模拟是用数值的方法来求解描述油藏中流体渗流特征的数学模型，是一门将计算机、应用数学、油藏工程等结合起来的综合性工程应用学科，在油气田开发方案设计和动态分析中有十分重要的作用，为了保证油藏的稳定产油量以及最终采收率，以获得最大的经济效益，必须对油藏的压力以及饱和度等参数进行监控，因而需要对油藏进行模拟，以确定合适的开采时间、注水量、开采速度等问题，需建立油水两相渗流模型，以对油藏的压力、饱和度进行即时监控，及时采取合适的开采措施来保障油田的高效、高产、持久生产,它能够从油田实际出发对油藏进行精确地描述，以最少的投资、最合适的速度去获得最高的最终采收率，也就是要获得最大的效益。1.2 国内外研究方法现状油藏模型中有压力的、和饱和度、两组未知量，目前基本上有两类求解wPoPwSoS方法。一类是顺序求解法，即先求压力项，后求饱和度；另一类是联立求解法，即同时求解压力项和饱和度项。同时，由于方程中含非线性系数，他们依赖于压力和饱和度的变化，在求解数学模型时有如下几种处理方法： (1)隐式压力显式饱和度方法（IMPES）即隐式求解压力方程，显式求解饱和度方程。它属于顺序求解法的一种，是数值模拟中最常用的、最简单的一种方法。具有占内存小、计算工作量小、方法简便等优点。但该方法存在两个问题：第一，达西项的系数处理是显式的，因此对如锥进的问题，由于井底周围流速高，压差变化大，而存在较大误差，对于强非线性问题的适应性也差；第二，饱和度的计算是显式的，当时间步长较大时，会出现解的不稳定性。因此，IMPES 方法只适用于一般的弱非线性渗流问题，对于某些非线性问题如注气、气锥或水锥等问题，IMPES 方法无能为力，即使时间步长取得很小，仍会出现解的振荡或算出的压力和饱和度为负值的情况，以至于模拟计算无法正常进行。 (2)半隐式方法半隐式方法属于联立求解的方法的一种，也是数值模拟中常用的一种方法。其基本思路是：联立求解油相方程和水相方程，同时求出压力和饱和度，因此压力和饱和度都是隐式求解。在计算过程中，半隐式方法对方程右端项的处理与 IMPES 方法完全相同，不同之处在于对方程左端项的处理。它需要对方程左端的达西系数、产量项及毛管压力进行泰勒级数展开，并忽略二阶小量，一阶导数项用 n 时刻的值。由于系数是处理近似的，并未真正用 n+1 时刻的值，所以这种方法叫做半隐式方法。(3)隐式压力隐式饱和度法（IMPIMS）隐式压力隐式饱和度法实际上是 IMPES 方法和半隐式方法的混合变种，它也是属于顺序求解方法的一种。它既有半隐式方法求解饱和度的特点，又保留了 IMPES 方法省内存、省工作量的特点。IMPIMS 方法（The Implicit Pressure Implicit Saturation Method）实际上是IMPIMS 方法和半隐式方法的混合和变种，它也是属于顺序求解方法的一种。它既有半隐式方法求解饱和度的特点，又保留了 IMPES 方法省内存、省工作量的特点。IMPIMS 方法的求解思路是：压力和饱和度分开顺序求解，求解压力时可直接利用IMPES 方法的压力求解方法，然后将求出的压力值代入半隐式方法的水相差分方程1n0p中，将该方程化为只有含水饱和度的一个变量的差分方程，再用隐式计算格式求解即可。(4)全隐式方法对于某些强非线性渗流问题如高速渗流等，即使使用半隐式方法也会引起计算结果的波动，或者时间步长只能取到很小。为此提出了全隐式方法,这里简单介绍一种与半隐式方法相类似的全隐式方法，即达西系数项也用泰勒级数展开并忽略二阶小量，但一阶导数项不用时刻的值而用时刻的值，由此构成一个非线性的代数方程组，n1n可用牛顿迭代法或其他非线性代数方程组的解法进行求解。1.3 本文章节安排 论文根据质量守恒原理建立了一维油水两相渗流的数学模型，经过逐步简化得到实验室进行单管模型一维水驱油实验时的数学模型，利用 IMPES 方法对此模型进行处理。第一部分首先列述了油藏渗流力学的相关知识，对油藏的油水渗流机理进行了详细的介绍以及油水两相渗流方程的推导过程，然后对有限差分法的基础知识进行了详细的介绍；第二部分主要是油水两相渗流模型的建立过程以及应用此方法对实例油藏的模拟。最后综述了油水两相渗流方程的建立过程中所遇到的问题、模型建立的理想条件及其准确性、针对实例模拟得出的结论。第 2 章 相关基础知识2.1 油藏渗流力学相关知识(一)、渗流过程中的力学分析 油气水之所以能在岩石空隙中渗流，是各种力共同作用的结果。1、流体的重力和重力势能流体具有质量，在中立场中受到地球的吸引力而具有重力。在渗流过程中，流体的重力和它的相对位置联系起来，就表现为重力势能。这种重力势能用压力表示为： （2-gzPz1）式中 -表示重力势能的压力，PazP -流体密度，3mkg -重力加速度，g2sm -相对位置高差，zm重力有时是动力，有时为阻力。2、流体的质量和惯性力 惯性是物体所固有的一种物理特征，其大小取决于质量。流体质量的大小一般用密度表示。流体由于具有质量，因此也具有惯性。当流体运动时，惯性使其总要维持原状，因而惯性力在渗流过程中多表现为阻力。根据牛顿第二定律，惯性力用压力表示为： （2-madtvPd-ma2）式中， -惯性力aP -质量m v -速度 -加速度a由于渗流时渗流速度通常很小，因此常忽略惯性力。3、流体的粘度及粘滞力 粘滞力是流体的一种特殊属性。在流动的流体中，如果各层的流速不同，将有一对作用力及反作用力，使原来快的流层减速，而慢的加速。这一对等值而又反向的力阻碍着流层的相对运动。流体的这种属性叫粘滞力。度量粘滞性大小的参数叫粘度。由牛顿内摩擦定律表示为： （2-dydvuAF 3）式中，A-两流层的接触面积，2m -沿流层法线方向的流速梯度，dydvmsm1- F-内摩擦力，N -粘滞系数，粘度的单位通常用表示：uspasmPa. 1smPasPa.10.3粘度的单位以表示时成为泊（P）：scmg. 1P=100cP（厘泊）在渗流中，粘滞力为阻力，且动力消耗主要用于克服渗流时流体的粘滞阻力。4、岩石及流体的压缩性及弹性力岩石和其中饱和的流体均具有压缩性（或弹性） ，因此使得在油气渗流过程中产生弹性力。油气层岩石埋藏于地下几百米甚至几千米，油气层上面覆盖的岩柱压力被油气层本身骨架和其中饱和的流体所承受。因此，储层岩石和其中的流体都处于受压缩状态。油层除承受上覆岩柱压力外，本身也受油层压力。油气层岩石和其他固体一样，在外力作用下，它的形状和体积都要发生变化；当消除外力时，它又能恢复到原来的形状和大小。我们把岩石能恢复原状的性质叫岩石的弹性，又称压缩性。油气田在开采以前，油层内岩石和流体都处于均衡受力状态，各种力是相互平衡的；当油气层投入开采之后 ，油气层的压力不断下降，上覆岩柱压力和油层内流体压力之间形成压力差，使之失去平衡而迫使岩石颗粒变形，排列更加紧密，结果导致岩层孔隙体积减小。岩石颗粒变形后，孔隙体积的缩小程度取决于岩柱压力和油气层流体压力的差值，同时还与岩石本身的压缩性有关。岩石的压缩性通常用压缩系数表示：式中，-孔隙体积的变化量；V -岩石的外表体积；fV -岩石的压缩系数。fC岩石的压缩系数表示油层压力每降低单位压力时，单位体积岩石中孔隙提及的缩小量。由于岩石颗粒的组分不同，且岩石的孔隙结构不同，孔隙形状和承受力的作用点均不一样，因此各类岩石的压缩系数时不一样的。岩层中的流体（油、气、水）也具有压缩性，当作用与流体上的外力增加时，其体积会缩小，繁殖会膨胀。在油气开采过程中，流体体积膨胀，产生弹性力，推动流体流入井底。液体的压缩性常用液体的压缩系数来描述，可表示为： （2-PVVCll114）式中，-液体的绝对体积；lV -压力改变时，液体体积相应的变化量；lVP -液体的压缩系数，表示改变单位压力时，单位体积液体的体积改变fC量。需注意的时，对于固体而言，弹性变形时有范围的，当外力超过某一限度时，将会发生形变。流体在弹性介质和塑性介质中的渗流规律是不同的。5、毛管力 多孔介质时由无数个微小的毛管连接而成的。在油田开发过程中，油气水的渗流时由一种流体驱动另一种流体，在两相界面上会产生压力跳跃，它的大小取决于分界面的弯曲程度和流体的性质。这个跳跃的力就称为毛管压力或毛管力，用表示。cP毛管力与流体的性质和弯曲程度之间的关系可以用拉普拉斯方程表示： （2-)11(rRPc5）式中，R，r-分界面的曲率半径。 -液液界面的表面张力。在渗流中，毛管力既可以表现为动力，也可表现为阻力。在驱替压力不大时，若油藏岩石亲水，则水驱油时毛管力为动力；若油藏岩石亲油，则水驱油时毛管力为阻力。（二） 、与油藏有关的压力概念以上对油气渗流过程中的力学现象和作用机理进行了讨论。油藏内的流体所受的各种力往往以压力的形式表示。外力所做的功将一起地层内液体能量的变化，这种变化将通过压力的变化来反映。因此，从本质上说压力是用来表征油藏能量的物理量。下面讨论与油气储集层有关的几个压力的概念。1.原始地层压力0P岩石地层压力是指油藏开发前流体所受的压力。史册的岩石地层压力一般是在油田所钻的第一批探井中测得的。2.供给压力cP供给压力是指当油藏中存在液源供给区时，在供给边缘上的压力。在人工注水条件下，水井井底压力即供给压力。3.井底压力wP井底压力是指油井正常生产时，在生产井井底所测得的压力。井底压力也成为井底流压，简称流压。4.折算压力 rP油藏在开发前，整个油田处于平衡状态，油藏各点的流体所具有的总能量时相等的。由水动力学可知，单位质量流体所具有的能量由比位能、比压能和比动能组成。 用 M 表示油层中的任意点，用 z 表示 M 点的标高，P 表示 M 点的实测压力值，表示油层条件下的液体密度，v 表示液体的流速，则 M 点的液体所具有的总能量即总水头H 为： （2-gvgPzH226）（三） 、单相不可压缩液体稳定渗流的基本微分方程推导：由于是不可压缩液体，不必考虑状态变化，所以不需要状态方程，则基本方程有：1）运动方程： （2-P gradK-v7） 2）连续性方程： （2-0v div8）将运动方程带入连续性方程进行数学模型推导，可得： （2-0P gradK-div9） 由于为常数，所以得到：uk （2-0zPyPxP22222210）上式极为单相不可压缩流体在均值地层中的稳定渗流的基本微分方程。他的适用条件为： （1）单相均质液体； （2）符合线性运动规律； （3）不考虑多孔介质及液体的压缩性； （4）稳定渗流； （5）渗流是等温的；（四） 、油水两相渗流的基本微分方程： 1）运动方程： 对于油相： （2-P gradK-vooo11） 对于水相： （2-P gradK-vwww12）2)油水两相渗流的连续性方程： 对于油相： （2-tSzvyvxv-ozoyoox13） 对于水相： （2-tSzvyvxv-wzwywxw14） 将运动方程带入连续性方程得： （2-tSzPKzyPKyxPKxooooooo15） （2-tSzPKzyPKyxPKxwwwwwww16）上面两式也可以写成： （2-tSPKooo17） （2-tSPKwww18）此外，还有饱和度方程： （2-1SSwo19）他的使用条件是： 1）彼此互不相溶，不起化学反应的油水两相同时流动； 2）不考虑毛管力及重力作用； 3）岩石和流体均不可压缩； 4）渗流符合线性渗流定律； 5）渗流过程是等温的。如果研究的是油水两相稳定渗流过程，那么液体饱和度将不随时间变化，即，得：0tStSwo （2-0PKoo20） （2-0PKww21）以上两式即为油水两相稳定渗流的综合微分方程。2.2 有限差分法 有限差分法：数值求解常微分方程或偏微分方程的方法。物理学和其他学科领域的许多问题再被分析研究之后，往往可以归结为常微分方程和偏微分方程的求解问题。一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值差值计算来近似得到。这种方法是指随着计算机的诞生和应用而发展起来的。其计算格式和程序的设计都比较直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分。有限差分法的具体操作可分为两个部分：（1）用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；（2）求解差分方程组。 在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。通常采用的是规则的分割方式，这样可以便于计算机自动实现和减少计算的复杂性。网路线划分的交点成为节点。若与某个节点 P 相邻的节点都是定义在场域内的节点，则 P 点称为非正则节点。在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。有限差分法的差分格式：一个函数在 x 点上的一阶和和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。如对一个单变量函数 f(x)，x 为定义在区间a,b的连续变量。以步长将a,b区间离散化，我们得到一系列节点=a,xh1x,，然后求出 f(x)在这些点上的近x2axx,xx2312hhbhnnxx1，似值。显然步长 h 越小，近似解的精度就越好。与节点相邻的节点有和，ixhixhix因此在点可以构造如下形式的差值：ix 节点的一阶向前差分),()(iixfhxfix 节点的一阶向后差分),()(hxfxfiiix 节点的一阶中心差分)()(hxfhxfiiix与点相邻两代男的泰勒展开式可以写为ix （2-)(! 4)(! 3)(2)()()( 4 3 2iiiiiixfhxfhxfhxhfxfhxf22） （2-)(! 4)(! 3)(2)()()( 4 3 2iiiiiixfhxfhxfhxhfxfhxf23） （2-22）-（2-23） ，并忽略 h 的平方和更高阶的项得到一阶微分的中心差商表示：。 （2-hhxfhxfxfiii2)()()(24）利用（2-22）和（2-23）式我们还可以得到一阶微分的向前，向后一阶差商表示：， （2-hxfhxfxfiii)()()(25）。 （2-hhxfxfxfiii)()()(26）将（2-22）和（2-23）式相加，忽略 h 的立方及更高阶的项得到二阶微分的中心差商表示： 。 （2-2 )()(2)()(hhxfxfhxfxfiiii27）利用（2-24）（2-27）式，我们就可以构造出微分方程的差分格式。这里要指出的是：在构造差分格式时，究竟应该选择向前，向后还是中间差分或差商来代替微分方程中的微分或微商，应当根据由此得到的差分方程解的稳定性和收敛性来考虑。同时兼顾到差分格式的简单和求解的方便。上述差分步骤应用于偏微分：例如，对于的情况，拉普拉斯算符在 0 点作用在此函数上的值),(yxff ，也可以用临近的点上的函数值来表示出来。 （见图 2.1.1，且22222yfxff时）hhhhh4321 （2-444422043212! 424yfxfhhffffff28）图 2-1 节点 0 及其附近节点对微分方程数值求解的误差来源：（1）方法误差（或截断误差） 。这时由于采用的计算方法所引起的误差。例如上面我们介绍的差商表示中，采用的泰勒展开式展开到第 n+1 项时的截断误差。)(1nhO具体方法的误差阶数取决于在离散化时的近似阶数。因此若改进算法就可以减少截断误差。（2）舍入误差（或计算误差） 。这是由于计算机的有限字长而造成数据在计算机中的表示出现误差。在计算机运算的过程中，随着运算次数的增加舍入误差会积累的很大。如果在多次运算后，舍入误差的精度影响是有限的，那么这个算法是稳定的，否则是不稳定的。不稳定的算法是不能用的。第 3 章 油水两相渗流机理和求解3.1 数学模型的建立假设条件如下：（1）油藏中仅存在油水两相渗流，油水互不溶解，且各自符合达西定律；（2）岩石、流体均可压缩；（3）考虑掩饰的非均质性及各向异性；（4）不考虑油水之间毛管力的影响；（5）忽略重力；这里根据数学模型的一般式，经逐步简化，得到实验室进行单管模型的一维水驱油实验室的数学模型。1、当考虑三维非均质油藏，油水互不相容，可压缩流体和岩石，考虑毛管力和重力时，数学模型的一般式为： （） （3-LLLLLLrLLStqDg-PKKw, oL 1）简化到一维，并忽略重力项： （） (3-LLLLLrLLStqxuKKxw, oL 2）2、假设： 不考虑掩饰的压缩性（即=常数） ，不考虑流体的体积变化（即） ； 11BowB， 油水粘度为常数。 于是得： 水相： （3-tSqxuKKxoovooor3） 油相： （3-tSqxuKKxoovooor4） 式中， 为地面标准状况下单位时间内单元体中注入（或采出）的体积流sclLlvqq量。上述两个偏微分方程中的未知量有 4 个，即、，因此还需要写出两个wowSoS辅助方程，即： （3-1owSS5） （3-wocow-PPP6）初始条件为： （3-lxSSP00 x00 xwcw，7） 边界条件为： （3-8） 0tqqq|qqq|qvovwvlxvvwv0 xv上述边界条件中，注入量产出量均为，表明该水驱油试验为稳定驱替。上述（3-1）vq（3-8）构成了该问题的完整的数学模型。 利用数值方法进行求解后，可得到在不同的注入速率下，模型中任意一点的压力、饱和度随时间的分布和变化。 3.2 数学模型的求解的方法及参数处理 对以上数学模型进行差分求解之前，这里首先对未知量的求解方法及有关参数的处理进行说明。1、数学模型的求解方法上述数学模型中有压力 、和饱和度、两组未知量，本文应用隐式压力oPwPwSoS显式饱和度（IMPES）求解法进行求解。IMPES 方法的基本思路： （1）通过乘以适当的系数，合并油方程和水方程，以消去微分方程组中的、wS，得到一个只含有、的方程。oSoSwP （2）由毛管力方程，可得，带入上面合并后的方程，得wocow-PPPcowow-PPP 到一个只含有的方程，成为压力方程。oP （3）方程左端达西项系数上一时间段的值，同时毛管力也用上一时间阶段的值，即显示处理系数。于是可形成一个高阶现行代数方程组，用迭代法可以进行求解，先求出，然后得。1noPnc1no1nw-PPP（4）将带入水相方程，用显式方法求出，然后得。1nwP1nwS1nw1no-1SSIMPES 方法具有所占内存小、计算工作量小、方法简便等优点。但该方法存在两个问题：第一、达西项的系数处理是显式的，因此对如锥进的问题，由于井底周围流速较高，压差变化大，而存在较大的误差，对于强非线性问题的适应性也差；第二、饱和度的计算是隐式的，当时间步长较大时，会出现解的不稳定性。因此，IMPES 方法只适t用于一般的弱非线性渗流的问题，对于某些非线性渗流的问题如注气、气锥或水锥等问题，IMPES 方法无能为力，即使时间步长取得很小，仍会出现解的震荡或算出的压力和饱和度为负值的情况，以致模拟计算无法正常进行。2、参数处理 在用有限差分对数学模型进行求解时，首先要将连续的油藏问题离散化为网格单元，然后对每一个网格单元，读入包括深度、有效厚度、孔隙度、渗透率、饱和度等基本参数。所有给定的这些参数都是网格节点处的值，在两个网格节点中间处的参数值是未知的，因此需要进行相应的处理。（1）渗透率的取值K 渗透率 K 是空间函数，其取值由以下几种方法：算数平均： （3-212/ 1iiiKKK9） 加权平均： （3-1ii1i1iii2/1ixxxKxKK10）调和平均： （3-11i1i2/ 1/x/xxxiiiiiKKK11） 几何平均： （3-1ii2/1iKKK12） （2）相对渗透率的取值 相对渗透率的取之原则上是去流动方向上的上游节点值，通常称为上游权法。如下图所示，取值方法是：当时，即由流向 时 i1 - i1ii （3-1rl2/1 - irl)(iwSKK13）当时，即由 流向时1 - iii1i (3-14)iwSKK)(rl2/1 - irl式中，为第 个节点的势。ii 图 3-1 井点分布图 采用上游权处理主要是对饱和度沿流向变化滞后进行修正。由于流动系数关于时间的显式处理造成流动的滞后，如流动系数在空间上按相邻节点的算数平均值取值，则会加重滞后现象，影响解的精度。因此，上游权的处理实质上是将显式处理造成的时间上的滞后用空间上的向前来进行弥补。（3）毛管力曲线的处理 由于油藏岩石的非均质性，即使用同一油层的岩心所测得的毛管压力曲线也有所不同，用 J 函数方法可以从实验室提供的大量毛管力资料中选择合适的毛管力曲线。J函数是把流体界面张力、岩石润湿性及渗透率和孔隙度等的影响综合在一起来表征油层的毛管压力曲线特征的一个无因次函数。实践表明，它是处理毛管压力曲线的一种有效方法。J 函数定义如下： （3- 21cos1000KSPSJowowwcw15）式中油、水界面张力，；owmmN 水的润湿角， （ ） ；ow 岩石毛管压力，； wcSPMPa 渗透率，；K2m 孔隙度，小数； 与相对应的无因次量。 wSJwS将众多的毛管压力资料按式（3-15)计算，做出关系图，如图 3-2 所示，然后拟wSJ 合出一条关系曲线。当油层岩性较为接近时，关系并根据油层，值分wSJ wSJ k布的峰值区域可求出一组有代表性的关系。 wwcSsP在进行油藏数值模拟计算时，一般需要将油、水相对渗透率曲线和毛管压力曲线离散化，然后在模拟计算过程中应用各种插值手段求出各网格在不同饱和度下的相对渗透率值和毛管压力值。 图 3-2 曲线图wSJ 3.3 差分方程组得建立及求解下面用方法来建立差分方程组并进行求解。IMPES、隐式求压力上面的数学模型式(3-3)至(3-6)中，有，4 个未知量。首先我们利wS0SoPwP用及方程(3-3)和(3-4)的关系消去，使其成为只含有压力，的1SSwooSwS0PwP方程。为此，式(3-3)式(3-4)得： （3-0ovwvoorowwrwqqxPKKxxPKKx16）根据，可得，代入上式可得只含有压力的方程。这wocowpppcowowppp0p里为简化起见，假设毛管压力为 0，于是可得，因此上式简化为：pppwo (3-17)0qqxpkkxxpKKxovwvorowrw令，分别表示油、水两相的流动系数：ow，orooKKwrwwKK表示总的流动系数： ow表示油、水两相的总流量：vq vovwvqqq则方程（3-17)可简化为： (3-18)0qxPxv方程（3-18)即为 IMPES 方法所化简得压力方程。下面写出该方程的隐式差分格式。假设采用块中心网格，且网格大小相等，均为，并假设网格为水注入处，x1i网格为油和油、水产处处。ni 式（3-18）的隐式差分方程为： (3-19)0qxxPPxPPvi1n1i1nii1ni1n1ii2121下面分 3 种情况来讨论方程（3-19)。(1)对于第 2 个至第个网格，既无注入也无采出，因此方程（3-18)中，1N 0viq于是可简化为： (3-20)0)PP()PP(1n1i1nii1ni1n1ii2121若系数按上游权原则取值，并用 n 时刻的值（即显示处理） ，则上式方程可简化为： (3-21)0)PP()PP(1n1i1nin1i1ni1n1ini即： (3-22)0PP)(P1n1ini1ninin1i1n1in1i(2)对于第 1 个网格，单位体积中注入的体积流量为，式(3-18)中第二项由于取viq上游权，于是可简化为：00n (3-23)0qx)PP(v21n11n2n1整理上式，可得： (3-24)n1v21n21n1qxPP 两端同乘以（网格单元的体积） ，并令（该网格的注入量） ，xAxAqQvv则上式可化简为： (3-25)AxQPPn1v1n21n1(3)对于第 N 个网格，其单位体积中采出的体积流量为，式(3-18)中第一项由于vq取上游权，于是可简化为：0ni (3-26)0qx)PP(v21n1N1nNn1N两端同乘以，令，得：xAxAqQvv AxQPPn1Nv1nN1n1N式（3-22),式（3-25)和式（3-27）构成了从第一个网格到最后一个网格（第个）N的线性代数方程组，矩阵方程如下：=1-11-11N12N2N33222211iAxWIN1N321PPPPPAxWIP-000AxWIP-wfwf （3-28）方程组的系数矩阵为三角矩阵，若用定流量来进行求解，则上述压力是与时间无关的量，不符合要求。可用方法进行求解，从而得到，在Thomas1P2PNP时刻的压力值。1n2.显式求饱和度 数学模型中的水相方程（3-3）的差分方程为： （3-tS-SqxxP-P-xP-Pnwi1nwiwvi1n1- i1nin1-wi1ni1n1inwi29）上式中的未知量为，也可分为以下三种情况来讨论：1nwiS（1）对于第 2 个至第个网格，既无注入也无采出，因此方程（3-28）中的1-N，0qwvi于是可得： （3-1n1- i1nin1-wi1ni1n1inwi2n1w1n1wPPPPxtSS30） 由于已经通过隐式求压力得到，故可用上式顺序（）求出每个网1niP1-n ,2,i格节点的饱和度值。(2)对于第一个网格，（即单位体积的注水速率） ，而且式（3-29）中的左面查vwviqq分的第二项去上游权，于是可化简为：0n0w （3-tSSqxPPn1w1n1wv21n11n2n1w31） 两边同乘以并化简，且，则有：xAxAqQvv tSSxQxPPn1w1n1wv1n11n2n1w、于是： （3-xPP-AQxtSS1n21n1n1wvn1w1n1w32） （3）对于最后一个（第 N 个）网格，其单位体积的产出量，实验中为ovwvvqqqvq定值。其中，产水量为 ，而且 是跟随时间变化的，因此首先要求出的值。wvqwvqwvq最后一个（第 N 个）网格的压力为 ，第 N-1 个网格的压力为，则对于NP1-NPN1-NP-P油和水来说均是相等的，故流出油、水提及的流量取决于流动系数、的大小。ow （3-1nN1n1NoovPPq33） （3-1nN1n1NwwvPPq34） 所以： （3-1nN1n1NwovPPq35） 于是 （3-wowvwvqq36） 方程（3-28）中的左面差分的第一项由于取上游权，于是可化简为：0nwn （3-tSS*qxPPnwN1nwNwv21n1N1nNn1-wN37）两边同乘以 ，并令 ，则：xAxAqQwvwv （3-tSSxQxPPnwN1nwNwv1nN1n1Nn1-wN、38）于是： （3-n1Nowwv1nN1n1-Nn1-wNnNw1nNwAQxPPxtSS39） 利用式（3-30） 、式（3-32）和式（3-39）可求得网格的 值。N,1,2,i1nwiS 以上所叙述的隐式求压力、显式求饱和度的计算过程即为 IMPES 方法从时刻到n时刻的求解过程。这样从初始时刻开始，一步步依次求解下去，直到求得所要求1n的时间的压力和饱和度为止，油水两相的渗流方程就可得解，所建立的油藏模型就得以模拟出来。3.4 有关单位换算 以上的计算过程中需用到达西公式，但达西公式中所用的单位为水力学单位，而实际油田中常用的单位为工程单位，二者之间的换算关系见表 3-1表 3-1 换算关系水力单位工程单位 Qs/cm3d/cm3 K Dum20.001 Dum2 usmPasmPa A2cm 2m x cm m P 0.1MPa MPa（1）达西公式为，其中各变量均采用水力学单位。xPAKQ（2）当采用工程单位时，达西公式变为：。243-610 x10P10A10K8640010Q（3）因此，由于水力学单位换算到工程单位，达西公式前加了一个系数，即： xPAK0.0864Q3.5 计算程序框图 程序开始 读入 、xLKwo读入N2,1iiKiSrww、读入iwcPS、0t 读入vmaxQxt、 计算不同下的、wSrwKroK计算值、ow计算压力方程的系数矩阵和右边项idibiaic、用追赶法iP显式求iS,iSow 图 3-3 一维油水两相程序框图 是打印iSiSiPtow、 结束maxtt 否3.6 计算实例下面把给定的油藏区块按一维方向均分成 7 个网格井点，以第一个网格井为注水井，第七个网格井为采油井，一注一采模型进行求解： 注水井 生产井图 3-4 井网分布图已知：孔隙度，渗透率，束缚水饱和度，原始含油饱0.252m1K2 . 0Swc和度，地层油粘度，地层水粘度，注水井井底压力8 . 0SoismPa5osmPa1w，生产井井底流压，有谁渗透率曲线如表 1 所示.油层长MPa15PiwfMPa10Pwf70m，宽，油层厚度，综合压缩系数，取，时间m7010mh 14MPa101Cx10m步长取一天，模拟时间为，油水相对渗透率曲线见附录。模拟曲线如下：1360d 050010001500200000.20.40.60.81t/dfw 图 3-5 最后一个网格点（生产井）关系曲线tfw总液量中水所占有的分量即为含水率,此曲线反映了油井的含水率随时间的变化wf趋势，开始后一段时间内含水率接近 0，因为水驱前沿未到达产油井，采出的完全是油，水驱前沿到达产油井后，含水率开始上升，经过相当长时间后含水率近似可以达到 1，此时油井开采完毕。根据此曲线可以综合评估油井的成本与效益，以便选择合适的开采时间。1P2P3P4P5P6P7P 12345671111.51212.51313.51414.5iP/Mpa 图 3-6 1-7 网格井点在 T 时刻的压力曲线此图反映了此模型从第一个井点到第七个井点的压力变化，可以看出是等差变化的。这是因为不考虑油藏的毛管力、重力以及岩石和流体的压缩性所致，模型为刚性水驱开采，故从注水井至开采井的压力是等差变化的。 0200400600800100012001400160018000.20.30.40.50.60.70.80.91t/dSw 图 3-7 最后一个网格井点的曲线tSw 此曲线反映了最后一口网格井的随时间的变化趋势，随着水驱前沿的推进，水wS侵占油藏的孔隙以排除原油，所以油藏的含水饱和度会随时间而增大，最后达到稳定值，因为此时油藏孔隙中的残余油不能被排出，仍然占据着油藏孔隙，所以随着时间的推移，含水饱和度不再变化。第 4 章 结 论1、本文应用隐式压力显式饱和度解法（IMPES 方法）来对油藏油水两相一维渗流进行模型建立与求解，通过合并流体方程得到一个只含有压力的方程，某一时间步的压力求解出来后，饱和度采用显式更新对油藏的模拟及求解。2、模型的求解有定流量和定压力两种方法，本文建立的模型不能用定流量来进行求解，因为若用定流量进行求解则此模型的压力方程与时间无关，会有无穷多解，这是因为只要保持注水井和采油井的压力恒定，就可以达到流量恒定，所以这样的压力解是无穷的。3、本文所用的 IMPES 方法具有占内存小、计算工作量小、方法简便等优点。4、本文建立的模型是在不考虑岩石和流体的压缩性以及油水的粘度、忽略重力，在这种刚性模型下对油藏的模拟准确度可能有一定的影响。IMPES 方法达西项的系数处理是显式的，因此对如锥进的问题，由于井底周围流速高，压差变化大，而存在较大误差，对于强非线性问题的适应性也差；饱和度的计算是显式的，当时间步长较大时，会出现解的不稳定性。5、通过含水率曲线表明，相当长时间后，含水率接近 100%，此时此油带的水驱油开采完毕。6、本文对一维油水两相渗流进行了数值模拟，对一维各个以划分的网格在个时刻的压力、含水饱和度、含水率等方面都有了一个全面的显示，能够对油藏较为准确的模拟，按此方法可以对二维、三维油水两相以及二维、三维油气水三相渗流进行拓展。致 谢 非常感谢张晓东老师在我大学的最后学习阶段毕业设计阶段给自己的指导，从最初的定题，到资料收集，到写作、修改，到论文定稿，他给了我耐心的指导和无私的帮助。为了指导我们的毕业论文，他放弃了自己的休息时间，他的这种无私奉献的敬业精神令人钦佩，在此我向他表示我诚挚的谢意。同时，感谢所有任课老师和所有同学在这四年来给自己的指导和帮助，是他们教会了我专业知识，教会了我如何学习，教会了我如何做人。正是由于他们，我才能在各方面取得显著的进步，在此向他们表示我由衷的谢意，并祝所有的老师培养出越来越多的优秀人才，桃李满天下！通过这一阶段的努力，我的毕业论文基于有限差分法的油水两相渗流方程求解终于完成了，这意味着大学生活即将结束。在大学阶段，我在学习上和思想上都受益非浅，这除了自身的努力外，与各位老师、同学和朋友的关心、支持和鼓励是分不开的。在本论文的写作过程中，我的导师张晓东老师倾注了大量的心血，从选题到开题报告，从写作提纲，到一遍又一遍地指出每稿中的具体问题，严格把关，循循善诱，在此我表示衷心感谢。同时我还要感谢在我学习期间给我极大关心和支持的各位老师.写作毕业论文是一次再系统学习的过程，毕业论文的完成，同样也意味着新的学习生活的开始。我将铭记我曾是一名中国石油大学学子，在今后的工作中把惟真惟实的优良传统发扬光大。 感谢各位专家的批评指导参考文献1哈利德.阿齐兹 ，安东尼.塞特瑞. 油藏数值模拟 第一版.北京：石油工业出版社，2004.2 Henrik Lof，Margot Gerritsen,Marco Thiele.Parallel streamline simulation.SPE 113543,2008.3 谢海兵，桓冠仁，郭尚平，等. PEBI 网格二维两相流数值模拟.石油学报，1999，20（20）：57-61.4 韩大匡，陈钦雷，阎存章油藏数值模拟基础.北京：石油出版社，1993 5 谢俊.剩余油饱和度平面分布方法研究及应用.西安石油学院学报，1998,13（4）:40-42.6 张建国，杜殿发，侯健. 油气层渗流力学第二版东营：石油大学出版社，2009.7 李淑霞，谷建伟. 油藏数值模拟基础 东营：石油大学出版社 2009. 8 http:/forun.petro- 邓建中，刘之行.计算方法.第二版.西安交通大学出版社 2001.10 胡良剑，邓晓君. MATLAB 实验 高等教育出版社 2006. 11 卓金武. MATLAB 在数学建模中的应用. 北京航空航天大学出版社，2011. 附 录附录 A表 1 油水相对渗透率曲线数据表序号wSroKrwK10.200.780.0020.300.610.0330.400.460.0740.500.310.1150.600.190.1660.700.100.2370.800.000.30附录 B程序代码：%*模型求解clear allclc%*赋初值*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.2 0.3 0.4 0.5 0.6 0.7 0.8;Krotab=0.78 0.61 0.46 0.31 0.19 0.10 0;Krwtab=0 0.03 0.07 0.11 0.16 0.23 0.3;%*循环求解*for t=1:(T+1) Krw(:,t)=interp1(Swtab,Krwtab,Sw(:,t); Kro(:,t)=interp1(Swtab,Krotab,Sw(:,t); lamdao(:,t)=K*Kro(:,t)/Uo; lamdaw(:,t)=K*Krw(:,t)/Uw; lamda(:,t)=lamdao(:,t)+lamdaw(:,t);%*A 矩阵赋值* A(1,1)=1+WI*dx/S; A(1,2)=-1; for i=2:(N-1) A(i,i-1)=lamda(i-1,t); A(i,i+1)=lamda(i,t); A(i,i)=-A(i,i-1)-A(i,i+1); end A(N,N-1)=1; A(N,N)=-1-WI*dx/S;%*求 q* q(1)=dx/S*WI*15;for j=2:N-1 q(j)=0; end q(N)=-dx/S*WI*10; %逆矩阵求解%B=inv(A); %p(:,t+1)=B*q; %追赶法求解 A*p(:,t+1)=q(:)p(:,t+1)=zhuiganfa(A,q,N);Qv(1,t)=WI*A(2,1)*(15-p(1,t+1);%求解 Sw 用Qv(2,t)=WI*A(N-1,N)*(p(N,t+1)-10);Sw(1,t+1)=Sw(1,t)+dert/(fai*dx)*(Qv(1,t)/S-lamdaw(1,t)*(p(1,t+1)-p(2,t+1)/dx); for r=2:(N-1) Sw(r,t+1)=Sw(r,t)+dert/(fai*dx)*(lamdaw(r,t)*(p(r+1,t+1)-p(r,t+1)-lamdaw(r-1,t)*(p(r,t+1)-p(r-1,t+1); endSw(N,t+1)=Sw(N,t)+dert/(fai*dx)*(lamdaw(N-1,t)*(p(N-1,t+1)-p(N,t+1)/dx-Qv(2,t)/S*lamdaw(N,t)/(lamdaw(N,t)+lamdao