（油气田开发工程专业论文）时变非线性渗流理论与数值计算.pdf

t h e t h e o r yo f n o n l i n e a r f l u i df l o ww i t ht i m e l y c h a n g i n g a n dn u m e r i c a lc a l c u l a t i o n b y s o n g y o n g ( o i l g a sf i e l dd e v e l o p m e n te n g i n e e r i n g ，u n i v e r s i t yo fp e t r o l e u m ，c h i n a ) s u p e r v i s e r ：l iy a n g ( p r o f e s s o r ) a b s t r a c t i ti si n d i c a t e df r o ml a b o r a t o r ye x p e r i m e n t sa n dl o c a lp r a c t i c et h a tl o n g t e r m s c o u r i n gt h r o u g hi n j e c t i n gw a t e rh a sc h a n g e ds o l n ep a r a m e t e r si nt h er e s e r v o i r s ， a n dt h i s c h a n g eh a s e f f e c to nt h em e c h a n i s mo fw a t e rd i s p l a c i n go i la n dt h e r e m a i n i n g o i ld i s t r i b u t i o n i nt h i sp a p e r , o nt h eb a s eo ft h er u l eo ft h ep a r a m e t e r sc h a n g i n gn o n l i n e a r l y a n dt i m e l yi nt h er e s e r v o i r sw h i c hi se s t a b l i s h e db yl a b o r a t o r ye x p e r i m e n t sa n dt h e s t a t i s t i c d a t a ，t h e m e t h o do fn u m e r i c a ls i m u l a t i o ni s u s e d ，a t w o p h a s e s m a t h e m a t i c a lm o d e li s e s t a b l i s h e di nw h i c h p a r a m e t e r sc h a n g i n gt i m e l y a r e c o n s i d e r e d t h ef u l l yi m p l i c i tm e t h o di ss t u d i e d ，b e s i d e s ，ap r a c t i c a b l en u m e r i c a l s i m u l a t o ri sd e v e l o p e d ， t h i sp a p e ra n a l y z e st h ei n f l u e n c eo ft h ep a r a m e t e r sc h a n g i n gt i m e l ys u c ha s o nt h ea b s o l u t ep e r m e a b i l i t y 、p o r o s i t y 、t h eo i ld e n s i t ya n dt h eo i lv i s c o s i t yo nt h e d i s t r i b u t i o no ft h er e m a i n i n go i l i ti ss h o w e dt h a tt h et i m e l yc h a n g i n go fa b s o l u t e p e r m e a b i l i t y h a st h ee f f e c to nt h ed i s t r i b u t i o no ft h e r e m a i n i n g a n dt h e w a t e r f l o o d i n ge f f e c t i tw i l l i n c r e a s et h ew a t e r - c u ta n dd e c r e a s et h eo i lr e c o v e r y ， t h ed i f f e r e n tc h a n g i n gr u l e sa l s oh a v et h ed i f f e r e n te f f e c t t h er e s u l ti nt h i sp a p e r c o m p a r e sw i t ht h ec o m m e r c i a ls i m u l a t o rc m g a n di ti n d i c a t e st h a tt h er e s u l ti s l e l i a b l e k e y w o r d s ：n h i n e l j c a 】s i m u l a t i o n p a r a m e t e r sc h a n g i n gt i m e l y , t h ef l l l l yi m p l i c i tm e t h o d 独创性声明 本人声明所呈交的论文是我个人在导m 口指导f 进行的研究工作及取得的 研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含 其他人已经发表或撰写过的研究成果，也不包含为获得石油大学或其它教育机 构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任饲f 献均已在论文中作了明确的说明并表示了谢意。 签名 彝压 苷铝l d 恤年。月p 日 关于论文使用授权的说明 本人完全了解石油大学有关保留、使用学位论文的规定，即：学校有权保 留送交论文的复e o q t 及电子版，允许论文被查阅和借阅；学校可以公布论文的 全部或部分内容，可以采用影印、缩印或其他复制手段保存论文。 f 保密论文在解密后应遵守此规定) 学生签名 导师签名 盟 9 i 拶 。铲年，o 月加日 石油大学( 华东) 硕士论文第1 章前言 第1 章前言 目前，我国东部大部分油田都已经进入了中高含水后期，而高含水期油田 稳产的关键取决于对剩余油分布的了解。当油田一旦进入中后期，了解和掌握 油藏中剩余油饱和度的空间分布，确定其剩余储量，剩余的可采量及其品位， 是油藏经营管理决策的重要依据。从已开发的油田增加产量和提高采收率是当 务之急，而确定剩余油分布是重中之重。 研究意义及目标 截至1 9 9 5 年底，全国油田综合含水率已高达8 1 3 6 ，而仅采出可采储量 的6 5 3 3 ，仍还有约3 5 的可采储量需在高含水期开采。目前我国油田平均 每采1 吨油将产水5 6 吨，随着含水率的增高，采1 吨油产出水量越来越大 【”，因此，认清高含水期剩余油分布特征，是十分迫切和亟待解决的重大问题。 油田由原始油层状态，经长期注水进入高含水期，主要发生了五个方面的 变化【1 】= 1 宏观上地下原油由集中连续分布变为油水交错分布； 2 油层参数主要表现出“三高、两低、四改变”一渗透率、孔隙度、粒 度中值增高；泥质含量、束缚水含量降低：地层水矿化度、孔隙结构、润湿性、 相渗透率曲线改变，并由三维空间变化变为四维空间变化，增加了随时间发生 变化； 3 ，以电法为基础的测井系列及解释图版，由于油层参数及地层水矿化度 改变，已不能准确判断油水层及水淹级别； 4 开采对象由集中分布的主力油层转向差油层、薄油层；主力油层进入 层内流动单元控制油水分布： 石油大学( 华东) 硕士论文第l 章前言 5 开发前期以井点和一定厚度层平均参数为基础的油田开发工作，已不 能把握油田开发主动权，要求进入精细定量划分油水层和井间不同环境沉积砂 体的分布及其参数的预测。 本文的研究目标就是在考虑地层参数随着时间发生变化的情况下研究高 含水期剩余油的分布。 研究现状 油藏数值模拟技术从5 0 年代开始研究至今，已发展成为一项较为成熟的 技术。在油田开发方案的编制和确定，油田开采中生产措施的调整和优化，以 及提高油藏采收率方面，已逐渐成为一种不可欠缺的主要研究手段 ”。 一般情况来说，油藏数值模拟研究可以概括为以下几个步骤p 1 明确此项模拟研究的任务、目的和具体要求。 2 收集各项必需的原始资料和数据，分析其完整程度和可靠性，必要时 提出需要补充的资料清单及获取它们的要求和方法。 3 整理所获得的各项资料，进行油藏描述，建立地质模型。 4 进行数值模拟方法的设计。 5 进行模拟计算及历史拟合。 6 预测各种不同开发方案或调整方案的动态及技术经济指标。 7 整理和分析模拟结果，提出建议，形成报告。 其中，数值模拟方法设计的基本过程可见下图： 2 石油大学( 华东) 硕士论文第l 章前言 l郁昆差分方程组 l i线i 妨程组 l 求解压力、饱和度等 图卜1 数值模拟方法设计的基本过程 1 9 9 1 年吕平发表了“毛管数对天然岩心渗流特征的影响”1 4 1 ，运用理论 分析和室内实验相结合的方法，对天然岩心驱油效率的影响因素进行了研究， 提出毛管数是影响驱油效率的主要因素之一。实验研究给出了作为饱和度和毛 管数的二元函数的一组相对渗透率曲线，使得贝克莱一列维尔特的两相渗流理 论得以发展，并为三次采油物理化学渗流理论和数值模拟研究提供了实验基 础。 1 9 9 9 年刘爱武和陈其荣发表了“谈谈油水渗流过程中的几个问题” 5 1 ，在 大量室内实验的基础上，探讨了油水单相流的特点；油水两相渗流的贾敏效应； 有机物沉积对油水两相渗流的影响等几个问题。其结论是：1 由于水膜的存在， 油水的单相流体在渗流过程中，表现出不同的特点；2 在油水两相渗流过程中， 存在着贾敏效应，要采取措施减小这种损害；3 原油组分中沥青质沉积吸附有 可能对油水两相渗流产生影响，应加强该方面的研究工作。 1 9 9 9 年文光耀，刘红等发表了“水驱油藏高含水期精细模拟技术研究与应 用”【6 ，在剖析高含水期开发过程中所表现出来的渗流现状和规律的基础上， 以解决开发阶段油藏工程问题为目的，系统地提出了油藏精细模拟的原理、思 路和方法，提高了模拟的精度，拓宽了油藏模拟的应用领域，已经在z w 油 田高含水开发单元剩余油研究和综合调整实际应用中发挥了重要作用，取得了 较满意的应用效果。 陈晓楼和刘爱玲发表了“水渗流模拟装置的研制”1 7 1 ，大庆油田由于长期 石油大学( 华东) 硕士论文第1 章前言 注水开发，地层中流体的渗流严重地影响固井、二界面的胶结质量。通过对 大庆油田的现场调研，确定了高含水期影响水泥环胶结的主要参数，并研制成 功了水渗流模拟装置。该装置通过建立井下水泥环微单元实现了对地层的渗透 率、孔隙压力、井下温度和地层中水的渗流状态模拟，模拟地层中水的渗流对 水泥环的冲刷、侵蚀等。该装置的研制成功，为进一步研究压差和水渗流对固 井质量的影响提供了依据。 2 0 0 0 年冯秀芳翻译的“高含水油藏的开采特点”1 8 中指出，在具有高原始 含水饱和度( o 4 o 6 ) 的储层中，在特定的孔道几何形态特征和热力条件下， 从油藏一开始投入开采就有可能出现油和水混合渗流。同时，随着储层孔隙中 含水量的升高，原油渗流条件变差，并且比最初达到临界含油饱和度的油藏更 强烈。这个特点在低渗透储层中更明显，因为，对低渗透储层中原油的渗流来 说需要更大的压力梯度。 2 0 0 1 年同登科等发表了“压力依赖于地层渗透率的分形油藏的数值研究” 【9 1 ，在文章中考虑了地层渗透率对于压力是敏感的，认为渗透率要发生变化， 并且是压力的函数。建立了应力敏感地层分形介质的数学模型，并采用 d o u g l a s - j o n e s 预估一校正法获得了无限大地层有限半径井定产量生产时的数 值解。 文光耀在“水驱油藏高含水期精细模拟技术研究与应用”【1 0 1 中考虑到了绝 对渗透率的变化，将其由传统的静态参数改为动态参数，动态变化为含水率的 函数。 目前，对于高含水油藏普遍认为储层物性参数与原始的参数有很大的不 同，特别是对于绝对渗透率，同登科考虑的是单相渗流，文光耀在考虑水驱两 相时只是认为绝对渗透率是含水率的函数。实际上还要考虑饱和度的因素，这 时绝对渗透率将还会成为饱和度的函数，此外，注入倍数也会影响绝对渗透率 的变化。 研究内容 石油大学( 华东) 硕士论文第1 章前言 本论文主要研究时变非线性渗流的理论以及数值计算。即建立能够准确描 述水驱渗流过程的时变非线性渗流模型、研究考虑时变规律的渗流场特征以及 与剩余油分布的关系。 1 时变非线性渗流数学模型的研究 在时变规律研究的基础上，考虑岩石物性变化、流体物性变化，建立参数 时变非线性渗流数学模型。 2 时变非线性渗流模型的数值求解方法研究 运用现代非线性动力学和数学方法研究时变非线性渗流模型的数值求解 方法，主要研究时交渗流模型数值求解方法。 3 研究考虑时变规律的不同含水阶段渗流场特征以及与剩余油分布的关 系。 4 编制时变非线性数值模拟软件。 5 ，现场应用。 技术路线 本论文的研究方法是采用数值模拟方法进行求解。技术路线如下： 】建立数学模型 采用的模型为简化的黑油模型。 2 建立数值模型 采用有限差分方法进行离散，对时间的离散采用向前差分，对空间的离散 采用中心差分。 3 ，建立线性方程组 石油大学( 华东) 硕士论文第1 罩前言 离散化后，得到的是一个非线性差分方程组，需要将其线性化，同时由于 本文考虑了绝对渗透率等参数随时间的变化，所以在数值求解时没有采用较常 用的i m p e s s 方法( 隐式压力显式饱和度) ，而是采用了全隐式方法。 4 线性方程组求解 采用预处理共轭梯度法来求解线性方程组。 5 计算机实现 采用v i s u a lc + + 编写了较为实用的变参数油水两相数值模拟器。 本文中公式、列表以及插图中所用参数的单位，除了特别加以注明之外均 采用达西单位制。 6 石油大学( 华东) 硕士论文第2 章变参数油水两相数学模型建立 第2 章变参数油水两相数学模型建立 应用数值模拟方法来研究油田开发问题，首先必须根据油藏的实际渗流情 况建立数学模型，即建立基本渗流方程式及相应的定解条件，形成一个完整的 数学方程组，如图所示【1 1 。 旦医霜 圈一圜 流动辅助方程 参数辅助方程 化学辅助方程 物理辅助方程 ，偏微分方程( 组) 图2 - i 建立渗流方程关系图 2 1 方程推导 本文选取的数学模型为简化的黑油模型，在模型的维数选择上采取了三维 模型。并有如下假设条件： ( 1 ) 油藏中存在油水两相流体渗流； ( 2 ) 油藏中岩石和流体均可压缩； ( 3 ) 油藏流体渗流符合达西定律； ( 4 ) 油藏岩石具有各向异性和非均质性； 7 石油大学( 华东) 硕士论文第2 章变参数油水两相数学模型建立 ( 5 ) 不考虑毛管力和重力的影响 ( 6 ) 考虑渗透率的变化。 1 _ 运动方程： 油相：g o ：一k k ，- 。v p 。 水相：吒：一k k 2 v p 。 2 质量守恒方程： 油组分：一v 如z o ) + q 。= 昙。咒) 优 水组分：一v 。v w ) + q 。= 昙如。s 。) o t ( 2 1 1 ) ( 2 1 2 ) ( 2 1 3 ) ( 2 1 4 ) 口 警v p 圯= 昙帅。) ( 2 ) 口 警叫圯= 昙帅。) ( 2 ) 式( 2 1 5 ) 和( 2 1 6 ) q b f f _ 端可分别展开为： v ( 警v p = 昙 警蝌冰警矧+ 兰 ( 警瑚 ( 2 1 7 ) v 降即) = 却警) ；讣球警) 斟昙 降瑚 ( 2 1 8 ) 8 石油大学( 华东) 硕十论文第2 章变参数油水两相数学模型建立 将 昙小舻。詈鹕鲁魄鲁 昙咖 ) = p w s w 扰城等慨等 等讯c 。詈 s 。= 1 一s 。 o s w ：一堡 o to t 代入( 2 1 9 ) 和( 2 1 1 0 ) 上式，得： 昙。兄) = 见s 。c 。+ 筘。成c o ) 瓦o p 一百o s w 昙小几s 鼬 ) 詈+ 驴爿 即： 去水 陋d r “o1 k 到 j + 兰夺 ( 2 1 ，9 ) ( 2 1 1 0 ) f 2 1 ，1 2 ) 牺+ 纰g ) 等一矧 r 2 1 1 3 ) 差永永扣枇+ 碳+ 矧 9 鲨研 望西 厶 呜 叫 型=芑监西 石油大学( 华东) 硕士论文第2 章变参数油水两相数学模型建立 3 辅助方程 4 初始条件 5 外边界条件 6 内边界条件 s 。+ s ，= 1 p k = 只。 【s 。b = s 。 竺1 一：0 o n l ( 2 1 1 4 ) ( 2 1 1 5 ) ( 2 1 ，1 6 ) f 2 1 1 7 ) q i ( x 彤郅) 。：。，= q f t ) ( 2 1 ，1 8 ) 2 2 变参数的处理 2 2 1 绝对渗透率的处理 本文考虑了绝对渗透率的变化，对它的处理采取了如下几种方法： 1 假设绝对渗透率是压力和含水饱和度的函数，将绝对渗透率转换成压 力和含水饱和度的二元函数参与计算。函数表达式如下： 世：( 2 一p 一耖肌去训)( 2 肌1 ) 式中： k 。一初始渗透率，a n 2 ； 只一初始压力，m p a ： s 。一初始含水饱和度。 2 假设绝对渗透率与压力存在函数关系，作如下两种假设： 1 ) 世= k 。+ 世。p 【o 1 ( p r ) 2 + o 2 ( p r ) j ( 2 2 1 2 ) j o 石油大学( 华东) 硕士论文第2 章变参数油水两相数学模型建立 式中： k 。初始渗透率，f a n 2 ； 只初始压力，m p a ； 此式中，渗透率与压差存在二次幂的函数关系，随着压差的增大，渗透率 也增大。 ：) k ：= k o 。+ 缸k + o ：e x 。b 2 ( p ( p 一- r p ) o + ) 2 ，j ：j 。2 p 只) jp一-rpok k l n ( 5 p ：( z z ，s ) 叫i = 。缸+ 2 ( p r ) + 1 ) 】一r 4 ”。 式中： 石油大学( 华东) 硕士论文第2 章变参数油水两相数学模型建立 k 。一初始渗透率，a n 2 ； 只一初始压力，m p a ； 此式中，渗透率的表达式为分段函数，当压差小于某一值时，渗透率与压 差存在二次幂函数关系，渗透率变化较快；当压差大于这一值后，渗透率与压 差存在对数函数关系，渗透率变化较慢。 3 假设绝对渗透率与注入水的孔隙体积倍数存在函数关系 算完毕后，对绝对渗透率进行修正，再计算第n + 1 时步。 1 ) k = k o ( 1 + 0 0 1 却2 + 0 0 1 - l 印) 式中： k 。一初始渗透率，舯2 ； i n p 一注入倍数： 大。 第h 时步计 ( 2 2 1 4 ) 此式中，渗透率为注入倍数的二次函数，随着注入倍数的增大，渗透率增 f k = k 。o + l n ( i n p + 1 ) )i n p 5 2 ) k = k 。g 刚5 + o 5 15 却 o ，一。 。 + o 。，用差分方程 l u ( x ，o ) = 鲈( x ) ，一o o x l 时，差分格式是不稳定的。 目前比较常用的稳定性分析方法有矩阵分析法和f o u r i e r 级数法。矩阵分 析法是一种比较严格的方法，但由于使用这种方法需要用到矩阵的特征向量并 要求出它的特征值，难度较大，因此并不是所有的问题都可以用矩阵分析法来 进行分析，特别是对油藏数值模拟中的一些强非线性问题。f o u r i e r 分析法理 论上没有矩阵分析法严格，例如没有考虑边界条件的影响等。但方法简单，使 用方便，因而在有油藏数值模拟中得到了广泛的应用。 1 f o u r i e r 分析法 定义在区间( o ，上) 上的任何一个函数，g ) ，只要它满足d i r i c h l e t 条件( 即 处处连续，或在区间上只有有限个间断点，并且间断点的跃度是有限的； 区闻内只有有限个极值) ，就可以展开成f o u r i e r 级数 ，g ) ：艺k 。s 等+ 删n 竽1 扣0 l l 根据e u l e r 公式 c 。s z = 昙g 矗+ e - j z ) s i n z = 昙g 豇- e - 血) 式中：j = 4 7 - 1 。 则上式写成复数形式为： 厂g ) = c ；e 。了 + 瓜 式中c t 为f o u r i e r 系数。 2 0 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 由于误差占。在区间( o ，三) 上是以离散形式存在的，即只有在节点上才有值。 为了使误差在区间( o ，) 上满足d i r i c h l e t 条件，需要补充误差在节点间的定义。 例如在点i 与点i + 1 之问，定义 心 i = 1 , 2 ， 一出 这样，就可以将误差在整个求解区间上展开为f o u r i e r 级数了。展开后， 在第t ”时刻点i 上的误差可以表示为： s 譬”z c ：e i 导 现在只拿出级数中的某一项 e 7 = c ；e 。 来进行分析。如果对应于式( 3 1 3 1 ) 0 0 的每一项 绝对值都不增大，那么，也将不增加。 n t n 式( 3 1 3 1 ) 既要能表示成e 随时间的变化 等于初始误差。假定c 具有形式 q = e “ ( 3 1 - 3 1 ) 在时间增加的过程中其 又要在t = 0 时保证误差 式中a 为某一常数( 实数或复数) ，一为时间步长指标，纽为时间步长。 代入( 3 1 3 1 ) 得 腻 e ? 2 e o n a t e 。 相邻两个时间步长的e 之比为 2 1 ( 3 ，1 3 2 ) 占 占 ，tl = 、jb s 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 i e 7 + 1 ：尝r ：善( 3 ) 霹 。一m 。卑 7。 式中善= e “称为误差放大因子。显然，若蚓1 ，则误差只能随着n 值的增加 而减少或至少不被扩大，从而保证了差分方程的稳定性。因此 蚓1 就成为差分方程稳定的条件，称为v o nn e u m a n n 条件。 为了使用方便，可将式( 3 1 3 2 ) 写成如下形式 e ? ：e i p 0 “ 式中羼：孥 l 用f o u r i e r 分析法对差分方程进行稳定性分析的一般步骤是：先求出相应 于某一种差分格式的误差方程，将误差方程中的占以式( 3 1 3 3 ) 形式的单项代 入，运算后求出放大因子占。若吲1 能够得到满足，则差分格式就是稳定的。 这时所求出的使蚓1 得以满足的条件，称为稳定性条件。若旧1 不能满足， 则差分格式是不稳定的。 2 矩阵分析法 考虑 三“= 詈一昙( a ( 彬) 塞) + d ( 州m = ，( 彬) 0 x 1 0 t t u ( x ，0 ) = 9 ( x ) ，0 x 1 u ( o ，f ) = u ( 1 ，f ) = 0 ，0 t t 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 的差分格式，可以写为a u “1 = 占。“+ 矿，设n 一1 阶矩阵a 的逆存在 则相应的齐次格式可写为l ，“1 = ( a m ) 。1 曰v 。 记h o = ( a ) 1 b ，并称之为过渡矩阵，则有 v “= h ( 。) v 当系数与k 无关时，a = a ，b = b ，h = a b ，则得到 v “1 = h v 。 通过对日( 1 或日进行直接估计来给出稳定性条件，这就是矩阵方法。 除了上述较常用的f o u r i e r 分析法和矩阵分析法，还有一种方法，称为h i r t 启示性方法，这种方法是一种近似分析方法。主要是把差分格式在某确定点上 作t a y l o r 基数近似展开，把高阶误差略去，只留下最低阶的误差项。如果差分 格式是相容的，那么这样得到的新的微分方程( 称之为第一微分方程或修正微 分方程) 与原来的微分方程相比只增加了一些含有小参数的较高阶导数的附加 项。h i r t 方法就是利用第一微分近似的适定性来研究微分格式的稳定性。h i r t 方法的判别准则是这样的：如果第一微分近似是适定的，那么原来微分方程的 差分格式是稳定的，否则是不稳定的。 设u ( x ，f ) 是偏微分方程的解，“? 是逼近这个偏微分方程的差分格式的“真 解”。这里所指的“真解”是指在求解差分格式过程中，忽略了各种类型的误 差，比如舍入误差等。即是说求解查分割似的过程是严格精确的。差分方程的 收敛性，是指对于一个差分格式，如果在所考虑的区域内( 包括空间坐标与时 间坐标) 的任一点，当时间步长f 和空间步长 趋向0 时，e ：= u ( x ，f 。) 一u ：哼0 ， 即差分方程的真解趋于原微分方程的真解，则称该差分格式时收敛的：否则称 差分格式是不收敛的。显然，收敛性对于一个差分格式来说是必不可少的，一 个不收敛的差分格式是没有任何意义的。 对差分方程的收敛性直接进行理论分析，通常要比稳定性分析困难得多。 然而，由于收敛性、稳定性和相容性之间有比较密切的关系，可以根据相容性 石油大学( 华东) 硕十论文第3 章变参数油水两相数值模拟研究 和稳定性的分析来推断其收敛性。差分方程分析中的l a x 等价理论就阐明了 这种关系：对一个相容逼近于原微分方程的差分方程来说，稳定性是收敛性的 必要和充分条件，或者说收敛性和稳定性是等价的。也就是说，只要证明了差 分方程的稳定性和相容性，也就证明了它的收敛性。 3 - 2 线陛方程组求解方法研究 描述多维多相流动的偏微方程组经过差分格式导出了差分方程，这些非线 性方程最终都是通过线性化后，变成一个联立的线性方程组来求解。求解方法 主要分为两种，一种为直接法，一种为迭代法。直接法所需的计算量随着网格 数的增加而快速的增长，因此对于网块数大的油藏模型，直接法可能会因为花 费太高而在实际中无法使用。迭代法包括线性超松弛法( l s o r ) 、强隐式方 法( s i p ) 、预处理正交极小化方法( o r t h o m i n ) 。 对于大型线性方程组，常用的线松弛法和强隐式方法也存在一些弱点，如 收敛速度强烈依赖于迭代因子的选取，选择不当，则收敛很慢，而迄今为止， 迭代因子的选取仍无成熟有效的办法。共轭梯度法早在5 0 年代提出，已经经 历了较长的发展阶段，7 0 年代以来，在共轭梯度法的基础上经过不断的改进， 并将其与矩阵的不完全因子分解法相结合，形成了一类新的算法，统称为预处 理共轭梯度方法，具体算法很多，但归纳起来都是把一种矩阵的不完全分解方 法和一种在共轭梯度法基础上发展出来的迭代加速法结合而成的算法，迭代求 解的基本原理如下【1 4 】： 设有方程组 丘石= b 假定m 为非奇异矩阵，方程组变换为： m x = 似一4 弦+ b 构造下列迭代公式： ( 3 2 1 ) ( 3 2 2 ) 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 m x “1 = 一a ) x + 6 ( 3 2 3 ) 记x 。州= x 。州一x ，r 。= b a z 则迭代公式变化为： f r + 1 = 尺。 ( 3 2 4 ) 若迭代是收敛的，则当k 足够大时，x “1 * x 。构造迭代方法的关键是 如何选取矩阵m ，使之以最少的迭代次数得到满足要求的解。显然，矩阵m 越接近系数矩阵彳，则达到收敛标准所需要的迭代工作量越少，但计算m 的 逆矩阵所需要时间也要增加。在直接解法中，矩阵m 就是a 本身。 把矩阵a 直接分解成和u 的乘积而进行求解所需的计算工作量很大， 而且由于在油藏模拟中求解的都是含有大量零元素的稀疏矩阵，这种分解常使 矩阵a 中大量本来为零的元素在上+ 【，中变为非零元素，增加了内存占用量。 因此，当线性方程组的阶数很大时，一般采用不完全的l u 分解。 为了求得一个矩阵的不完全l u 分解，首先要进行矩阵的符号消元。所谓 符号消元并不是根据矩阵的元素的具体数值进行消元计算，而是根据矩阵中非 零元素的所在位置求出在g a u s s 消元过程中每一个可充填位置的充填级次。其 含义是：由于在进行g a u s s 消元法时，一方面消去了一部分的非零元素，同时 另一方面在带宽范围内也把一部分位置上的零元素变成了非零元素，这也可以 称为在这个位置上充填了新的非零元素，而消元过程中充填的先后次序则称为 充填级次。一般来说，新产生( 充填) 的非零元素值要显著的小于原有的非零 元素，后充填的( 充填级次大的) 非零元素值n d , 于先充填的( 级次小的) 非 零元素。所以为了要通过不完全l u 分解获得较简单的近似矩阵，必须设法减 少后充填的非零元素，即保留充填级次较低的非零元素，去掉级次较高的非零 元素而把它们赋为零值。不过，保留的元素太少了，矩阵求解时迭代的次数会 增多，总起来也不一定合算；而保留的非零元素太多了，迭代次数虽可减少， 但消元过程就和常规的g a u s s 消去法差不多了，又失去了不完全l u 分解的意 义。所以，在符号消元时要确定一个较优的保留级次。 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 符号消元时，充填级次是这样规定的：首先，在矩阵的每一个可充填位置 上赋一个初始充填级次，( z f ) 。 刮a ， + 刚。 其它 其次，设符号消元已进行到第m 步，即用第i n 行的主对角元素消去m 行以 下第m 列上的各非零元素，则第m 行阻下各充填位置的充填级次重新定义为 f ( x # ) = m i n p ( 石f ) ，z ( 工。) + ，( z 柳) + 1 j ，= m + 1 ， 矩阵的f 级不完全分解，就是在l u 分解过程中，指保留充填级次低于或 等于某一预先确定的值，的位置上的非零元素，形成原矩阵的一个近似分解。 显然，分解级次，定的越大，分解后的矩阵就越接近于原矩阵。对于足够大 的f ，这种不完全l u 分解实际上可以等价于矩阵的完全l u 分解。 至于分解后的近似矩阵l u 中各元素的值的计算，可按照完全l u 分解时 的计算方法进行。当遇有充填级次大于f 的位置时，可以直接置其元素为零； 也可以将它的值经过适当的处理后，加到相邻的充填级次小于或等于，的位 置上，然后再置该位置的元素为零。 有了矩阵的不完全的上u 分解后，迭代格式可写为： l u a x “1 = r ( 3 2 5 ) 式中r = b 一爿z u = m = 彳+ e 其中e 为误差矩阵。因此。方程组的求解过程为 y = l - 1 r a r “= u 一y ( 3 2 6 ) ( 3 2 7 ) k ；p 1 | ) v x ( f 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 的。 x “：x + a y “ ( 3 2 8 ) 因为三、u 分别为下三角、上三角矩阵，所以上述公式的求解是非常容易 但是，如无特殊的加速收敛措施，上述迭代过程的收敛速度非常慢。如 果稀疏矩阵是对称的，共轭梯度法是一种很有效的加速方法。但是，对于非对 称矩阵，常规的共轭梯度法就不能使用了，而在这种方法基础上发展起来的 o r t h o m i n 方法却能非常有效的解决非对称矩阵的问题。o r t h o m i n 方法 采用的加速措施有两个：一是正交化，二是极小化。所谓正交化是指：如果把 a x “1 看成维空间的一个向量，则它在维空间中确定了一个方向，在大多 数迭代方法中，对z 的修正都是沿该方向进行的。而o r t h o m i n 方法不是 采用a x “1 的方向作为修正方向，而是以a x “1 与本次迭代以前的各次迭代的 修正向量q 构造一个新的向量g “1 ，使m q “1 与以前各次迭代的m q 正交，并 以q “1 的方向为本次迭代对z 2 的修正方向。而极小化是指在确定了对z 的修 正方向口“1 后，其修正值要再乘上一个“1 丽成为国“1 q “1 ，其中国“1 称为极 小化因子，它的选取使新的余量恤一( d k + l m q “1 队为最小，以保证新的余量不 l z 大于以前的所有余量而加速收敛。 o r t h o m n 的整个迭代过程可写为 x 制= f 一1 丑 k a y q i 。l x “1 = x + + 1 q + 式中：一正交化系数： 口”= ( 岣7 ，脚“1 ) ( 坳，岣1 ) ( 3 2 9 ) r 3 2 1 0 ) ( 3 2 1 1 ) f 3 2 1 2 ) 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 伍。，峋“1 ) 巧m。qi+i,mqtm) f 3 2 1 3 ) 实践表明，实际计算时，并不需要使m q “1 与以前所有迭代m q 的正交， 而只要使其与前面的有限个( 称为n o r t h 个) g q l 正交就可阻了。这样，式 ( 3 3 1 0 ) 可以重新写为： q “。( 3 ，2 1 4 ) n o r t h 值的选取随问题的规模、难易程度、所需要的分解方法而异。例 如对于黑油模型，n o r t h 值约为5 1 0 为宜，而对热采问题约在1 0 1 5 左右。 还有一种处理方法则不是每次迭代都作n o r t h 次正交化，而是采取所谓重新 启动的方式，即在进行了n o r t h + 1 次迭代后，再重新从头开始进行正交化。 这样，在进行了n o r t h 次后每次迭代所进行的正交化次数平均只有 n o r t h 2 次，从而减少了工作量。实际计算表明，前述两种正交化处理方法 在收敛速度上并无明显差异；因此，目前大多数采用后一种方法。为了将这两 种处理方法加以区别，一般称前一种方法为非重新启动型o r t h o m i n 法，而 稼后一种处理方法为重薪启动型0 r t h o m 矾法。 3 3 全隐式方法研究 求解多相渗流问题的方法从原理上大体可分为两类：一类是对多个独立变 量一个一个的逐次求解的方法，一般是先求油层压力，然后再求饱和度或泡点 压力等其他变量。如果是先隐式的求出压力后再显式的求饱和度，称为隐式压 力显式饱和度方法，或称i m p e s ( i m p l i c i t p r e s s u r ee x p l i c i ts a t u r a t i o n ) 方法； 如果在踌式求出压力后再用隐式来求其他参数如饱和度等，则称为隐式交营僻 法，也称隐式序贯解法，简称s e q ( s e q u e n t i a ls o l u t i o n ) 方法。另一类是联立 求解方程组，从而同时解出各独立变量的方法，称为联立解法，简称s s ( s i m u l t a n e o u ss o l u t i o n ) 方法，在联立方法中，根据对方程中各有关系数线 性化方法的不同，又可以分为半隐式方法( s e m i i m p l i c i t m e t h o d ) 和全隐式方 始一 曲 一簖 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 法( f u l l yi m p l i c i tm e t h o d ) 。 i m p e s 方法的主要思路如下【1 5 】 1 ) 通过乘以适当的系数，合并油方程和水方程，以消去差分方程组中的 饱和度变量，从而得到一个只含压力的方程； 2 ) 由毛管压力公式可将p w 表示为p o 和p c 的形式，得到一个只含变量 p o 的压力差分方程： 3 ) 方程左端达西项系数用上一时间段的值，毛管力也用上一时间段的值， 即显式处理系数； 4 ) 在解出压力差分方程式之后，将求解出的压力值作为已知条件代入油 方程或水方程，显式求出饱和度值； 5 ) 井点所在网格的产量项均作显式处理，由上一时间段的饱和度值计算 井点网格的油水产量。 本文对方程组的差分采用了隐式差分，这样方程组无论是左端或右端都 是非线性的，线性化处理采用了全隐式方法。 用全隐式方法解方程时，每一时间步都要经过一系列的迭代过程。它的 基本做法是在每一个n + 1 时间步的开始，先按第n 时间步末所得的求解变量的 值( 在每一个时间步开始时要选择一组迭代初值) ，求出方程组内各系数的值， 接着解方程组，开始迭代。每一次迭代都求出求解变量的一组新值，以此求出 新的系数来更新原来的近似系数值，然后再解方程组，这样逐次迭代下去，直 至求出一组满足精度要求的值为止，这最后的一组迭代值便可作为n + 1 时间 步的终值。这就是用全隐式方法求解方程组一个时间步的迭代过程。然后再转 入下一个时间步的迭代，这样一步一步迭代下去，就是全隐式方法求解的全过 程。 对于任意一个从t ”n t ”1 的时间步长，在用全隐式方法迭代求解时，设对 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 任意变量彳，定义任意两次迭代j 和j + 1 之间z 的差值为 f i x = x ”1 一x s 即爿：x5 + g x 而且，当s = 0 时工o = ” 经过多次迭代，当防一石5l 占时，x “= x ”+ 其中占为所给的精度要求。 3 4 差分方程的建立 1 方程左端项的差分化 由( 2 1 5 ) ： v ( 訾即h ( 警料球警料却警阍o p 在离散时，( p i k k ，- tv p ) “用s + 1 次迭代解( 旦! 望鱼v p ) m 来代替，即： ( 警= 降寸( 引”1 v 限川 + ( 甜店凇+ ( 玎”芷5 如f v p ) 5 4 ， + k s + 1 嗡h 酬5 上式可近似为 如砜卯忖“一( 1 1 5 互垫查堂! 望查! 堡主笙塞 笙! 至壅童堑迪查堕塑塑堕堡垫翌壅 ( 等卯卜( 鲁寸降卜“圳l，l f ， + ( 钭弦3 + ( 鲁 。k “卅b 聊 a 固 埘时嗡h m 卵p 悯”1 俐 p t k k ，- j ：旯，堕鱼：，生：屈，一i ：n ，整理得 | “ i i弘 舭翳北。0 雌p 、1 s + l 。一 二 ( 3 4 3 ) b 瓤f k 掣缝掣 + 等犁吼。) + 。b 丁咖訾剖慨。一。，信l 。+ 嗽肚一) j 吼。萼剖斯划 + 盯掣斛， 3 1 ，lj k j l 2 石油大学( 华东) 硕士论文第3 章变参数油水两相数值模拟研究 令： b 瓤肿屯裂 + * 助 瑞。) 一眩。 等犁斛嘛。一龟船)毡+ 趣1l 咧，、2 删 4 啦船， 趣+ h ( 3 4 4 ) + 么。嘞。裂敞。略一，+ 匮i 胁嘛。一锄_ 巩掣孰b - y 2 圹s 私) + 巩掣戡吟， 蛳t2 碉翻9 ( 3 4 5 ) 吒彤。 霸乙。厂训 + ：。x 。 毒 ：k 。忆k m 一只：，t 【筹 ：。必 ； 州叫。( “厂t ) ( 3 4 6 ) 石油人学( 华东) 硕士论文 第3 章变参数油水两相数值模拟研究 2 血i ( a x 。+ z ) 9 r 嘣儿2 碌翻 嗤。加。厂t 彤 t 喘。( 锐。僻t l , j , k - - i , j , k - - 叫。 ( 3 4 ，7 ) 叫叫。叫。广t ：( 孰坳，。 w 咧s 。o 啪。厂t 助，。 吒2 ( + 。tp j + 必。t r m + 。e s “+ 必。t ) 瓯2 ( 一胁e 胁一一, j , k s o , i - 胁) 于是，( 3 4 4 ) 和( 3 4 5 ) 可化简为 ( 3 4 8 ) ( 3 4 9 ) ( 3 4 1 0 ) b 瓤k 瓷鲁掣。+ s 瓤划s a b 瓤。至簧享犁+ 气咖嚷气咖s 勘厂瓯c s a - z ， 则，( 3 4 3 ) 可化简为： 互垄奎堂! 兰堡! 堡丝塞 塑! 童壅叁塑迪查亘塑垫堡堡型堕壅 昙r 务i 1 同理，可以得到 h 。訾掣+ p 叫z + 1 。广蹴肚一_ + h 。訾掣+ 诹揽工。峨渺嘞。一瓯 ar 。a p l1 万卜刮2 瓦 孤of ，丑。o 沈p f l = 面1 ( 3 4 1 3 ) ( 芑，+ 。! ! ! ! i ! ：；丢害；塑+ ，咒，。只篁必。+ r 枷，必，。s ，艺：m ，。一气 + f 芑“。三! ! i i i ! ：；云掣十。一巧。e z “，。+ r 埘，见，。s 。毛，。一g 一 f 吒。叫! ! a ! z ： j 与a 二z 掣+ 乇。+ 鼻了：+ + 。+ s o i l 。叫t 十 k 十1 垭“4 + 以“抖儿删小舯，2 ”“2 + f 芑。! a 墨z ：；! 二+ ! a 二z ；塑+ 乙，；- y ：，只：一必+ r 。，。一s ，艺：。一 c 9 。1 嘎 k 倒小2 u 4 一悲。扩 4 叫 f 3 4 1 4 ) 一_ 。) ( 3 4 。1 5 ) 2 方程右端项的差分化 篝墅= 古脚删