黑油模型解剖

1黑油模型理论基础1.1基本假设（1）油藏中渗流是等温渗流；（2）油藏中有油、气、水三相，各相流体的渗流均符合达西定律；（3）模型考虑油组分、气组分、水组分三组分；（4）气组分在油气相、水气相之间发生质量交换；（5）相平衡瞬间完成；（6）水组分只存在于水相中，与油气相之间没有质量交换；（7）油藏岩石微可压缩，各向异性；（8）油藏流体可压缩，且考虑渗流过程中重力、毛管力的影响。1.2数学模型r(1)辅助方程：1owgcowowcgogossspppppp(2)初始条件：000000,,,,,,,,,,,,,,,twtwotopxyztpxyzsxyztsxyzsxyztsxyz(3)边界条件：0,,,()(,,),,,,,LvvwfwfpnQxyztQtxyzpxyztptxyz(4)2黑油模型程序整体结构图3组员及分工姓名负责程序张敏主程序，CGIN,MF1,INTERP,INTPVT高博主程序，TRAN1,SOLMAT,MATBAL,PRTPS郝延征主程序，TABLE,UINIT1,CODES,NODES徐杨主程序，GRID1,PARM1,QRATE4主程序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迭代求解的基本原理设有方程组：AXB(5)并设M为非奇异矩阵，则可以构造迭代公式：1kkMXMAXB(6)或写成：1kkMXBAX(7)式中，11kkkXXX，表示两次迭代之间的增量。式(7)就是通常所用的迭代公式。如果迭代是收敛的，则当k足够大时，11,0kkkXXX，式(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中大量本来为零的元素在LU中变为非零元素，增加了内存占用量。即使是带状矩阵，如果矩阵的阶数很大，在带状区域内，这种由零变为非零元素的数量也仍然很大。因此，在油藏数值模拟中，当线性方程组的阶数很大时，实际上一般不采用高斯消元法或相应的LU分解方法。矩阵的不完全LU分解就是为了解决这一问题而提出的。这种方法可以尽量保持矩阵A原有的稀疏性质，从而节约了内存，减少了计算工作量。充填级次越大，分解后的矩阵就越接近于原矩阵，对于足够大的充填级次，这种不完全LU分解实际上可以等价于矩阵的完全LU分解。而矩阵的零级不完全LU分解的非零元素全部位于原矩阵中非零元素所在的位置上。这里，黑油模型的子程序CGIN中运用的就是系数矩阵A的零级不完全LU分解。5.2.3ORTHOMIN加速法有了矩阵的不完全LU分解之后，令1MDEDDF(8)式中，E和F是矩阵A的严格下三角矩阵和严格上三角矩阵。D是待求解的对角阵。则迭代公式(7)可写成：11kkDEDDFXR(9)式中，kkRBAX，因此，式(5)的求解过程变为：1kDEXR(10)121DXX(11)12kDFXX(12)11kkkXXX(13)因为1,,DEDDF分别为下三角矩阵、对角阵、上三角矩阵，所以上述公式的求解是非常容易的。但是如果没有特殊的加速收敛措施，上述迭代过程的收敛速度非常慢。如果系数矩阵是对称的，则共轭梯度法是一种很有效的加速方法。但是，对于非对称矩阵，常规的共轭梯度法就不能使用了，而在这种方法的基础上发展起来的ORTHOMIN方法却能非常有效的解决非对称矩阵的问题。ORTHOMIN方法所采用的加速措施有两个：一是正交化，二是极小化。所谓正交化是指如果把1kX看成N维空间的一个向量，则它在N维空间中就确定了一个方向，在大多数迭代方法中，对kX的修正都是沿着该方向进行的。而ORTHOMIN方法不是采用1kX的方向作为修正方向，而是以1kX与本次迭代以前的各次迭代修正量iq构造一个新的向量1kq，使1kMq与以前各次迭代的iMq正交，并以1kq的方向为本次迭代对kX的修正方向。而极小化是指在确定了对kX的修正方向1kq后，其修正值要再乘上一个1k而成为11kkq，其中1k成为极小化因子，它的选取使新的余量112kkkRMq为最小，以保证新的余量不大于（一般是小于）以前的所有余量而加速收敛。整个ORTHOMIN方法的迭代过程为：11kkXMR(14)1111kkkkiiiqXaq(15)111kkkkXXq(16)111kkkkRRAq(17)其中：11,,ikkiiiMqMXaMqMq(18)1111,,kkkkkRMqMqMq(19)式中，1kia为正交化系数，其目的是让所有的iMq都正交。实际应用上述方法时，并不需要使1kMq与以前所有迭代的iMq都正交。根据式(15)，当迭代次数k增加时，每次迭代所需计算的求和项增加，不但增加了工作量，而且增加了计算中的舍入误差。实践表明，只要使1kMq与其前面的有限个（成为NORTH个）iMq正交就可以了。这样，式(15)可以重新写为：111kkkkiiikNORTHqXaq(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：误差限。5.4不完全LU分解预处理共轭梯度法黑油模型的子程序CGIN中运用的不完全LU分解预处理共轭梯度法的迭代过程如下：11kkXMR(21)111kkkkiiikKNKqXaq(22)111kkkkXXq(23)111kkkkRRAq(24)其中，2111,,ikkiiiMRRaRMR(25)1111,,kkkkkRMRAQQ(26)5.5程序流程图6子程序MF16.1功能介绍利用系数矩阵的不完全LU分解，求出两次迭代之间的增量。6.2主要变量说明X：两次迭代之间的增量；Y：上一次迭代的残余量；N：总网格数；NN1：Y方向和Z方向网格数的乘积；NZ：Z方向的网格数；DEE：系数矩阵不完全LU分解中D-1对角元。6.3算法说明系数矩阵A的不完全LU分解，1MDEDDF，目的是求1XMY，即MXY，可分如下几步：（1）DERY,DE是下三角矩阵；（2）1DSR,1