1、1 黑油模型理论基础1.1 基本假设（1）油藏中渗流是等温渗流；（2）油藏中有油、气、水三相，各相流体的渗流均符合达西定律；（3）模型考虑油组分、气组分、水组分三组分；（4）气组分在油气相、水气相之间发生质量交换；（5）相平衡瞬间完成；（6）水组分只存在于水相中，与油气相之间没有质量交换；（7）油藏岩石微可压缩，各向异性；（8）油藏流体可压缩，且考虑渗流过程中重力、毛管力的影响。1.2 数学模型* MERGEFORMAT (1)rw wwvro oovorg srg oogswr sswwvgokspgDqBtBtkRkppDBBR RDqtBAAA 辅助方程：* MERGEFORMAT (2
2、)1owgcgoosp初始条件：* MERGEFORMAT (3)00,twtwotopxyzpxyzss边界条件：* MERGEFORMAT (4)0,(),)Lvvwf wfpnQxyzttxyzp2 黑油模型程序整体结构图3 组员及分工姓名 负责程序张敏 主程序，CGIN, MF1, INTERP, INTPVT高博 主程序，TRAN1, SOLMAT, MATBAL, PRTPS郝延征 主程序，TABLE, UINIT1, CODES, NODES徐杨 主程序， GRID1, PARM1, QRATE4 主程序4.1 主程序主要功能（1）定义运算所需数组；（2）需要调用和生成的文件的打
3、开和关闭；（3）通过调用子程序给模型赋基础数据和初始数据；（4）通过调用子程序给模型的运行确定各种控制；（5）在运行过程中反复读入井点数据（包括产量和注入量、井底流压、流动指数等）以及打印输出控制码；（6）分层计算油气水地质储量；（7）进行井点产量项处理和形成压力矩阵；（8）通过调用子程序求解压力方程；（9）显式计算饱和度；（10）进行过泡点处理；（11）根据压力和饱和度增量控制，自动调节计算时间步长；（12）未饱和网格块饱和度计算；（13）变泡点处理；（14）在每一运算时间步末进行物质平衡检验；（15）打印油井、水井、气井的分层报告和分井报告；（16）通过调用子程序进行计算结果打印输出和形成
4、文件；（17）重启动运行方式的选择和重启动文件的生成和调用；（18）运算终止的控制及错误信息的提示。4.2 主程序流程图5 子程序 CGIN5.1 功能介绍用不完全 LU 分解预处理共轭梯度法求解压力方程隐式差分得到的大型稀疏线性方程组。5.2 理论基础实践表明，预处理共轭梯度型方法无论是适应性还是计算速度，都远远超过了前面的方法。这类算法的优点是迭代收敛速度不依赖于迭代因子的选取，收敛速度极快，应用范围极广，包括化学驱、混相驱、热力驱等复杂模型和各种难以求解的系数矩阵。例如，可以解决传统的共轭梯度法解决不了的非对称矩阵问题及一般算法难以求解的病态系数矩阵。因此，目前这种类型的算法已成为油藏数
5、值模拟中最先进的求解大型线性方程组的方法。这里，黑油模型的子程序 CGIN 中运用的是将系数矩阵的不完全 LU 分解和 ORTHOMIN 加速法结合的预处理共轭梯度法。5.2.1 迭代求解的基本原理设有方程组：* MERGEFORMAT (5)AXB并设 为非奇异矩阵，则可以构造迭代公式：M* MERGEFORMAT (6)1kkM或写成：* MERGEFORMAT (7)1kkXBA式中， ，表示两次迭代之间的增量。11kkkX式* MERGEFORMAT (7)就是通常所用的迭代公式。如果迭代是收敛的，则当 足够大k时， ，式* MERGEFORMAT (7)就成为* MERGEFORMA
6、T (5)。构11,0kk造迭代方法的关键问题之一是如何选取矩阵 ，使之能以最少的迭代次数得到满足要求的M解。显然，矩阵 越接近于系数矩阵 ，则达到收敛标准所需要的迭代次数越少，但相应MA的，求解方程* MERGEFORMAT (7)所需要的时间就要增加。在前面所讲的简单迭代法中，为对角矩阵，其元素为矩阵 主对角线上的相应元素，这是构造矩阵 的最简单的方法。M高斯-赛德尔迭代法中，矩阵 为矩阵 的下三角部分。在直接解法中，矩阵 为矩阵 本A身，这时式* MERGEFORMAT (7)完全等价于* MERGEFORMAT (5)，只要一次迭代就可以求得式* MERGEFORMAT (5)的解。在
7、黑油模型的子程序 CGIN 运用的预处理共轭梯度法中，矩阵 为矩阵 进行不完全 LU 分解后得到的近似矩阵。MA5.2.2 矩阵的不完全 LU 分解把大型矩阵 A直接分解成 与 的乘积而进行求解所需要的计算工作量大，而且由于在LU油藏数值模拟中求解的大多是含有大量零元素的稀疏矩阵，这种分解常使矩阵 A中大量本来为零的元素在 中变为非零元素，增加了内存占用量。即使是带状矩阵，如果矩阵的阶L数很大，在带状区域内，这种由零变为非零元素的数量也仍然很大。因此，在油藏数值模拟中，当线性方程组的阶数很大时，实际上一般不采用高斯消元法或相应的 LU 分解方法。矩阵的不完全 LU 分解就是为了解决这一问题而提
8、出的。这种方法可以尽量保持矩阵 A原有的稀疏性质，从而节约了内存，减少了计算工作量。充填级次越大，分解后的矩阵就越接近于原矩阵，对于足够大的充填级次，这种不完全LU 分解实际上可以等价于矩阵的完全 LU 分解。而矩阵的零级不完全 LU 分解的非零元素全部位于原矩阵中非零元素所在的位置上。这里，黑油模型的子程序 CGIN 中运用的就是系数矩阵 A的零级不完全 LU 分解。5.2.3 ORTHOMIN 加速法有了矩阵的不完全 LU 分解之后，令* MERGEFORMAT (8)1MDEF式中， 和 是矩阵 A的严格下三角矩阵和严格上三角矩阵。 是待求解的对角阵。EF D则迭代公式* MERGEFO
9、RMAT (7)可写成：* MERGEFORMAT (9)11kkEFXR式中， ，因此，式* MERGEFORMAT (5)的求解过程变为：kkRBAX* MERGEFORMAT (10)1kDEXR* MERGEFORMAT (11)2* MERGEFORMAT (12)12kF* MERGEFORMAT (13)1kkkX因为 分别为下三角矩阵、对角阵、上三角矩阵，所以上述公式的求解1,DEF是非常容易的。但是如果没有特殊的加速收敛措施，上述迭代过程的收敛速度非常慢。如果系数矩阵是对称的，则共轭梯度法是一种很有效的加速方法。但是，对于非对称矩阵，常规的共轭梯度法就不能使用了，而在这种方法
10、的基础上发展起来的 ORTHOMIN 方法却能非常有效的解决非对称矩阵的问题。ORTHOMIN 方法所采用的加速措施有两个：一是正交化，二是极小化。所谓正交化是指如果把 看成 N 维空间的一个向量，则它在 N 维空间中就确定了一个1kX方向，在大多数迭代方法中，对 的修正都是沿着该方向进行的。而 ORTHOMIN 方法不k是采用 1kX的方向作为修正方向，而是以 1kX与本次迭代以前的各次迭代修正量 构iq造一个新的向量 ，使 与以前各次迭代的 正交，并以 1kq的方向为本次迭代1kq1kM iMq对 k的修正方向。而极小化是指在确定了对 k的修正方向 后，其修正值要再乘上一个 而成为 ，其中
11、 成为极小化因子，它的选取使新的余量11k1k为最小，以保证新的余量不大于（一般是小于）以前的所有余量而加速2kRq收敛。整个 ORTHOMIN 方法的迭代过程为：* MERGEFORMAT (14)1kkXMR* MERGEFORMAT (15)1kkkiiqaq* MERGEFORMAT (16)11kkk* MERGEFORMAT (17)RAq其中：* MERGEFORMAT (18)11,ikki iiMqXa* MERGEFORMAT (19)11,kkkRq式中， 为正交化系数，其目的是让所有的 i都正交。1kia实际应用上述方法时，并不需要使 1kM与以前所有迭代的 iMq都正
12、交。根据式* MERGEFORMAT (15)，当迭代次数 增加时，每次迭代所需计算的求和项增加，不但增加了工作量，而且增加了计算中的舍入误差。实践表明，只要使 1k与其前面的有限个（成为 NORTH 个） iMq正交就可以了。这样，式* MERGEFORMAT (15)可以重新写为：* MERGEFORMAT (20)111kkkiiiNORTHqXaq其中，NORTH 值的选取随研究问题的规模、难易程度及所需要的分解方法而异。例如对于黑油模型，NORTH 值以 510 为宜，而对热采问题，在 1015 左右为宜。还有一种处理方法不是每次迭代都做 NORTH 次正交化，而是采取所谓重启动的方
13、式，即在进行了NORTH+1 次迭代后，再重新从头开始正交化。这样，在进行了 NORTH 次迭代后，每次迭代所进行的正交化次数平均只有 次，从而减少了工作量。实际计算表明，前述两2NORTH种正交化处理方法在收敛速度上并无明显差异，因此目前大多采用后一种方法。为了将这两种方法加以区别，一般称前一种方法为非重新启动型 ORTHOMIN 方法，而称后一种处理方法为重新启动型 ORTHOMIN 方法。最后，将前面介绍的不完全 LU 分解与 ORTHOMIN 加速法相结合，就是黑油模型的子程序 CGIN 中运用的预处理共轭梯度法。这种方法可用于求解常规迭代法难以求解的各种复杂的系数矩阵方程，甚至是病态
14、的系数矩阵方程等。另外，这种方法具有适应性强、计算速度快、收敛速度快等优点，已成为目前油藏数值模拟中求解大型线性方程组的最有效的方法，在许多数值模拟软件中都被广泛应用。5.3 主要变量说明A：压力方程七对角系数矩阵；PVEC：网格块压力；QVEC：压力方程右端向量；UL1：系数矩阵下三角最右边的一条对角线，与主对角线紧邻，是与(I,J,K-1)网格块压力有关的系数；UL2：系数矩阵下三角中间的一条对角线，是与(I,J-1,K)网格块压力有关的系数；UL3：系数矩阵下三角最左边的一条对角线，是与(I-1,J,K)网格块压力有关的系数；U1：系数矩阵上三角最左边的一条对角线，与主对角线紧邻，是与
15、(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：误差限。5.4 不完全 LU 分解预处理共轭梯度法黑油模型的子程序 CGIN 中运用的不完全 LU 分解预处理共轭梯度法的迭代过程如下：* MERGEFORMAT (21)1kkXMR* MER
16、GEFORMAT (22)11kkkiiiKNqaq* MERGEFORMAT (23)11kkk* MERGEFORMAT (24)RAq其中，* MERGEFORMAT (25)211,ikkiiiMRa* MERGEFORMAT (26)11,kkkAQ5.5 程序流程图6 子程序 MF16.1 功能介绍利用系数矩阵的不完全 LU 分解，求出两次迭代之间的增量。6.2 主要变量说明X：两次迭代之间的增量；Y：上一次迭代的残余量；N：总网格数；NN1：Y 方向和 Z 方向网格数的乘积；NZ：Z 方向的网格数；DEE：系数矩阵不完全 LU 分解中 D-1 对角元。6.3 算法说明系数矩阵 A
17、 的不完全 LU 分解， 1MEF，目的是求 ，即1XMY，可分如下几步：MXY（1） , 是下三角矩阵；DER（2） , 是对角矩阵；1S1（3） , 是上三角矩阵；FX6.4 程序流程图7 子程序 INTERP7.1 功能介绍与泡点压力无关的参数插值。7.2 主要变量说明INTERP(X,Y,N,XO,YO)X：PVT 列表中的压力值；Y：已知的与泡点压力无关的参数；N：PVT 列表中压力值的数目；XO：各网格块的压力；YO：待求参数。7.3 程序流程图8 子程序 INTPVT8.1 功能介绍与泡点压力有关的参数插值。8.2 主要变量说明INTPVT(BPT,RM,X,Y,N,XO,YO)BPT：泡点压力；RM：斜率；X：PVT 列表中的压力值；Y：已知的与泡点压力有关的参数；N：PVT 列表中参数值数目；XO：各网格块的压力；YO：待求参数。8.3 程序流程图9 总结通过解剖黑油模型，加深了对黑油模型程序执行的总流程的认识。尤其，通过对 CGIN和 MF1 子程序的解剖，透彻理解了“目前求解大型稀疏线性方程组的不完全 LU 分解预处理共轭梯度法” ，收获颇丰。也锻炼了团结合作的能力。感谢李淑霞老师的指导！
