1 黑油模型理论基础1.1 基本假设（1）油藏中渗流是等温渗流；（2）油藏中有油、气、水三相，各相流体的渗流均符合达西定律； （3）模型考虑油组分、气组分、水组分三组分； （4）气组分在油气相、水气相之间发生质量交换； （5）相平衡瞬间完成；（6）水组分只存在于水相中，与油气相之间没有质量交换； （7）油藏岩石微可压缩，各向异性；（8）油藏流体可压缩，且考虑渗流过程中重力、毛管力的影响。

1.2 数学模型()()()()()rw w w w vw w w w ro o o o vo o o o rg so ro g g o o g g o o g sw rw so o sw w w w vg w w go w kk s p gD q B t B kk s p gD q B t B kk R kk p gD p gD B B s R kk R s R s p gD q B t B B B φρμφρμρρμμρφμ⎡⎤⎛⎫∂∇∇-+= ⎪⎢⎥∂⎣⎦⎝⎭⎡⎤⎛⎫∂∇∇-+= ⎪⎢⎥∂⎣⎦⎝⎭⎡⎤⎡⎤∇∇-+∇∇-+⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎛⎡⎤∂∇∇-+=++ ⎢⎥∂⎣⎦⎝⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎛⎫⎫⎪ ⎪⎪ ⎪ ⎪⎪⎭⎝⎭⎩(1)辅助方程：1o w g cow o w cgo g o s s s p p p p p p ++=⎫⎪=-⎬⎪=-⎭(2)初始条件：()()()()()()000000,,,,,,,,,,,,,,,t w t w o t o p x y z t p x y z s x y z t s x y z s x y z t s x y z ===⎧=⎪⎪=⎨⎪=⎪⎩(3)边界条件：()()()()0,,,()(,,),,,,,L v v wf wf pn Q x y z t Q t x y z p x y z t p t x y z δδ∂⎧=⎪∂⎪=⎨⎪=⎪⎩(4)2 黑油模型程序整体结构图3 组员及分工4 主程序4.1 主程序主要功能（1）定义运算所需数组；（2）需要调用和生成的文件的打开和关闭；（3）通过调用子程序给模型赋基础数据和初始数据；（4）通过调用子程序给模型的运行确定各种控制；（5）在运行过程中反复读入井点数据（包括产量和注入量、井底流压、流动指数等）以及打印输出控制码；（6）分层计算油气水地质储量；（7）进行井点产量项处理和形成压力矩阵；（8）通过调用子程序求解压力方程；（9）显式计算饱和度；（10）进行过泡点处理；（11）根据压力和饱和度增量控制，自动调节计算时间步长；（12）未饱和网格块饱和度计算；（13）变泡点处理；（14）在每一运算时间步末进行物质平衡检验；（15）打印油井、水井、气井的分层报告和分井报告；（16）通过调用子程序进行计算结果打印输出和形成文件；（17）重启动运行方式的选择和重启动文件的生成和调用；（18）运算终止的控制及错误信息的提示。

4.2 主程序流程图5 子程序CGIN5.1 功能介绍用不完全LU 分解预处理共轭梯度法求解压力方程隐式差分得到的大型稀疏线性方程组。

5.2 理论基础实践表明，预处理共轭梯度型方法无论是适应性还是计算速度，都远远超过了前面的方法。

这类算法的优点是迭代收敛速度不依赖于迭代因子的选取，收敛速度极快，应用范围极广，包括化学驱、混相驱、热力驱等复杂模型和各种难以求解的系数矩阵。

例如，可以解决传统的共轭梯度法解决不了的非对称矩阵问题及一般算法难以求解的病态系数矩阵。

因此，目前这种类型的算法已成为油藏数值模拟中最先进的求解大型线性方程组的方法。

这里，黑油模型的子程序CGIN 中运用的是将系数矩阵的不完全LU 分解和ORTHOMIN 加速法结合的预处理共轭梯度法。

5.2.1 迭代求解的基本原理设有方程组：AX B =(5)并设M 为非奇异矩阵，则可以构造迭代公式： ()()()1k k MX M A X B +=-+(6)或写成：()()1k k M X B AX +∆=-(7)式中，()()()11k k k X X X ++∆=- ，表示两次迭代之间的增量。

式(7)就是通常所用的迭代公式。

如果迭代是收敛的，则当k 足够大时，()()()11,0k k k X X X ++≈∆≈，式(7)就成为(5)。

构造迭代方法的关键问题之一是如何选取矩阵M ，使之能以最少的迭代次数得到满足要求的解。

显然，矩阵M 越接近于系数矩阵A ，则达到收敛标准所需要的迭代次数越少，但相应的，求解方程(7)所需要的时间就要增加。

在前面所讲的简单迭代法中，M 为对角矩阵，其元素为矩阵A 主对角线上的相应元素，这是构造矩阵M 的最简单的方法。

高斯-赛德尔迭代法中，矩阵M 为矩阵A 的下三角部分。

在直接解法中，矩阵M 为矩阵A 本身，这时式(7)完全等价于(5)，只要一次迭代就可以求得式(5)的解。

在黑油模型的子程序CGIN 运用的预处理共轭梯度法中，矩阵M 为矩阵A 进行不完全LU 分解后得到的近似矩阵。

5.2.2 矩阵的不完全LU 分解把大型矩阵A 直接分解成L 与U 的乘积而进行求解所需要的计算工作量大，而且由于在油藏数值模拟中求解的大多是含有大量零元素的稀疏矩阵，这种分解常使矩阵A 中大量本来为零的元素在L U +中变为非零元素，增加了内存占用量。

即使是带状矩阵，如果矩阵的阶数很大，在带状区域内，这种由零变为非零元素的数量也仍然很大。

因此，在油藏数值模拟中，当线性方程组的阶数很大时，实际上一般不采用高斯消元法或相应的LU 分解方法。

矩阵的不完全LU 分解就是为了解决这一问题而提出的。

这种方法可以尽量保持矩阵A 原有的稀疏性质，从而节约了内存，减少了计算工作量。

充填级次越大，分解后的矩阵就越接近于原矩阵，对于足够大的充填级次，这种不完全LU 分解实际上可以等价于矩阵的完全LU 分解。

而矩阵的零级不完全LU 分解的非零元素全部位于原矩阵中非零元素所在的位置上。

这里，黑油模型的子程序CGIN 中运用的就是系数矩阵A 的零级不完全LU 分解。

5.2.3 ORTHOMIN 加速法有了矩阵的不完全LU 分解之后，令()()1M D E D D F -=--(8)式中，E -和F -是矩阵A 的严格下三角矩阵和严格上三角矩阵。

D 是待求解的对角阵。

则迭代公式(7)可写成：()()()()11k k D E D D F X R +---∆=(9)式中，()()kkR B AX =-，因此，式(5)的求解过程变为：()()1k D E X R -=(10)121D X X -=(11) ()()12k D F X X +-∆=(12)()()()11k k k X X X ++=+∆(13)因为1,,D E D D F ---分别为下三角矩阵、对角阵、上三角矩阵，所以上述公式的求解是非常容易的。

但是如果没有特殊的加速收敛措施，上述迭代过程的收敛速度非常慢。

如果系数矩阵是对称的，则共轭梯度法是一种很有效的加速方法。

但是，对于非对称矩阵，常规的共轭梯度法就不能使用了，而在这种方法的基础上发展起来的ORTHOMIN 方法却能非常有效的解决非对称矩阵的问题。

ORTHOMIN 方法所采用的加速措施有两个：一是正交化，二是极小化。

所谓正交化是指如果把()1k X +∆看成N 维空间的一个向量，则它在N 维空间中就确定了一个方向，在大多数迭代方法中，对()kX 的修正都是沿着该方向进行的。

而ORTHOMIN 方法不是采用()1k X +∆的方向作为修正方向，而是以()1k X +∆与本次迭代以前的各次迭代修正量()i q 构造一个新的向量()1k q +，使()1k Mq +与以前各次迭代的()i Mq 正交，并以()1k q +的方向为本次迭代对()k X 的修正方向。

而极小化是指在确定了对()kX 的修正方向()1k q +后，其修正值要再乘上一个()1k ω+而成为()()11k k q ω++，其中()1k ω+成为极小化因子，它的选取使新的余量()()()112kk k R Mq ω++-为最小，以保证新的余量不大于（一般是小于）以前的所有余量而加速收敛。

整个ORTHOMIN 方法的迭代过程为： ()()11k k X M R +-∆=(14)()()()()1111kk k k ii i qXa q +++==∆-∑(15) ()()()()111k kk k X X q ω+++=+(16) ()()()()111k k k k R R Aq ω+++=-(17)其中：()()()()()11,,i k k i i i Mq M X a Mq Mq ++⎡⎤∆⎣⎦=⎡⎤⎣⎦(18)()()()()()1111,,k k k k k R Mq Mq Mq ω++++⎡⎤⎣⎦=⎡⎤⎣⎦(19)式中，()1k i a +为正交化系数，其目的是让所有的()i Mq 都正交。

实际应用上述方法时，并不需要使()1k Mq +与以前所有迭代的()i Mq 都正交。

根据式(15)，当迭代次数k 增加时，每次迭代所需计算的求和项增加，不但增加了工作量，而且增加了计算中的舍入误差。

实践表明，只要使()1k Mq +与其前面的有限个（成为NORTH 个）()i Mq 正交就可以了。

这样，式(15)可以重新写为：()()()()111kk k k ii i k NORTHqXa q +++=-=∆-∑(20)其中，NORTH 值的选取随研究问题的规模、难易程度及所需要的分解方法而异。

例如对于黑油模型，NORTH 值以5~10为宜，而对热采问题，在10~15左右为宜。

还有一种处理方法不是每次迭代都做NORTH 次正交化，而是采取所谓重启动的方式，即在进行了NORTH+1次迭代后，再重新从头开始正交化。

这样，在进行了NORTH 次迭代后，每次迭代所进行的正交化次数平均只有2NORTH次，从而减少了工作量。

实际计算表明，前述两种正交化处理方法在收敛速度上并无明显差异，因此目前大多采用后一种方法。

为了将这两种方法加以区别，一般称前一种方法为非重新启动型ORTHOMIN 方法，而称后一种处理方法为重新启动型ORTHOMIN 方法。

最后，将前面介绍的不完全LU 分解与ORTHOMIN 加速法相结合，就是黑油模型的子程序CGIN 中运用的预处理共轭梯度法。

这种方法可用于求解常规迭代法难以求解的各种复杂的系数矩阵方程，甚至是病态的系数矩阵方程等。

另外，这种方法具有适应性强、计算速度快、收敛速度快等优点，已成为目前油藏数值模拟中求解大型线性方程组的最有效的方法，在许多数值模拟软件中都被广泛应用。

5.3 主要变量说明A ：压力方程七对角系数矩阵； PVEC ：网格块压力； QVEC ：压力方程右端向量；UL1：系数矩阵下三角最右边的一条对角线，与主对角线紧邻，是与(I,J,K-1)网格块压力有关的系数；UL2：系数矩阵下三角中间的一条对角线，是与(I,J-1,K)网格块压力有关的系数； UL3：系数矩阵下三角最左边的一条对角线，是与(I-1,J,K)网格块压力有关的系数； U1：系数矩阵上三角最左边的一条对角线，与主对角线紧邻，是与(I,J,K+1)网格块压力有关的系数；U2：系数矩阵上三角中间的一条对角线，是与(I,J+1,K)网格块压力有关的系数； U3：系数矩阵上三角最右边的一条对角线，是与(I+1,J,K)网格块压力有关的系数； DE ：主对角线；DEE ：系数矩阵不完全LU 分解中D -1对角元； KNK ：共轭梯度法中正交化次数的限定； NMAX ：总网格数；NN1：Y 方向和Z 方向网格数的乘积； NZ ：Z 方向网格数； EFC ：误差限。