油藏数值模拟原理黑油模型(for培训处)

1、- 1 - 2006年年4月月13日日 李建芳李建芳 油藏数值模拟原理油藏数值模拟原理黑油模型- 2 -黑油模型是指油质较重性质的油藏类型。黑油模型是指油质较重性质的油藏类型。 黑油模型是最完善、最成熟，也是应用最黑油模型是最完善、最成熟，也是应用最为广泛的模型。为广泛的模型。是油藏数值模拟的基础，其它模型大都是黑油模是油藏数值模拟的基础，其它模型大都是黑油模型的扩展。型的扩展。目的：软件这样一个目的：软件这样一个“黑箱黑箱”是如何操作的，以是如何操作的，以便指导我们更好地去应用软件。便指导我们更好地去应用软件。- 3 - 黑油模型的基本假设黑油模型的基本假设 单相渗流黑油模型基本微分方程单相

2、渗流黑油模型基本微分方程 三维三相黑油模型基本微分方程三维三相黑油模型基本微分方程 初边值条件初边值条件5. 网格系统网格系统6. 黑油模型的差分方程黑油模型的差分方程7. 差分方程的线性化差分方程的线性化8. 线性代数方程组的求解线性代数方程组的求解9. 黑油模型的主要数据流黑油模型的主要数据流- 4 -（1） 油藏中的渗流是等温渗流。油藏中的渗流是等温渗流。（2） 油藏中最多只有油、气、水三相，每一相均遵守达油藏中最多只有油、气、水三相，每一相均遵守达西定律。西定律。（3）油藏烃类只含有油、气两个组分。在油藏状态下，）油藏烃类只含有油、气两个组分。在油藏状态下，油气两组分可能形成油气两相，

3、油组分完全存在于油油气两组分可能形成油气两相，油组分完全存在于油相内，气组分则可以以自由气的方式存在于气相中，相内，气组分则可以以自由气的方式存在于气相中，也可以以溶解气的方式存在于油相中，所以地层内油也可以以溶解气的方式存在于油相中，所以地层内油相为油组分和气组分的某种组合。在常规油田中，一相为油组分和气组分的某种组合。在常规油田中，一般不考虑油组分向气组分挥发的现象。般不考虑油组分向气组分挥发的现象。 （4）油藏中气体的溶解和逸出是瞬间完成的，即认为油）油藏中气体的溶解和逸出是瞬间完成的，即认为油藏中油气两相瞬时达到相平衡状态。藏中油气两相瞬时达到相平衡状态。 （5）油水之间不互溶；天然气

4、也假定不溶于水。）油水之间不互溶；天然气也假定不溶于水。- 5 -单相渗流黑油模型基本微分方程单相渗流黑油模型基本微分方程黑油模型基本微分方程主要是由黑油模型基本微分方程主要是由达西公式达西公式(达西达西定定律律)和和连续性方程连续性方程(物质守恒定律物质守恒定律)所构成所构成 一维单相渗流的一维单相渗流的Darcy公式为：公式为： 三维单相渗流的三维单相渗流的Darcy公式为：公式为：QK dPU = -Adx()KPgZ U- 6 -渗流速度是一个空间向量，用分量表示为：渗流速度是一个空间向量，用分量表示为： 渗流速度矢量；渗流速度矢量； 压力；压力； 势函数势函数 绝对渗透率；绝对渗透率

5、； 重力加速度常数；重力加速度常数； 粘度；粘度； 垂向坐标；垂向坐标； 密度；密度；UK()()()xyzKPgZKxxKPgZKyyKPgZKzz UUUPgZ- 7 -Gringarten, 2002Data from all sources- 8 -根据物质守恒原理，在单位时间内根据物质守恒原理，在单位时间内流入单元内的流体质量流入单元内的流体质量-流出单元的流体质量流出单元的流体质量=单元内流体质量的变化单元内流体质量的变化取渗流场中一个微小六面体体积单元来研究取渗流场中一个微小六面体体积单元来研究- 9 -212121212121jjjiiiyyyxxxjjjjjjiiiiiijj

6、jiiiyyyyyyxxxxxxyyyxxx12112112112121212121- 10 -该六面体单元的中心点坐标为（该六面体单元的中心点坐标为（x，y，z），），其其每一侧面的质量流速均以其侧面的中心点的质量每一侧面的质量流速均以其侧面的中心点的质量来代替。来代替。 方向流入质量方向流入质量= 方向流出质量方向流出质量= 方向流入质量方向流入质量= 方向流出质量方向流出质量= 方向流入质量方向流入质量= 方向流出质量方向流出质量= 单元体内流体质量的变化单元体内流体质量的变化= x(, , )2()xxxy zy zu x(,)2()xxxy zyzu yy(, )2()yyx yzx

7、zu (, )2()yyx yzxzu zz(,)2()zzx y zxyu ( , ,)2()zzx y zx yu ()xyzt - 11 -根据根据物质守恒定律物质守恒定律：将上式两边同除以将上式两边同除以 ，(, , )(, , )22( , )( , )22( , ,)( , ,)22()()()()()()()xxxxxy zxy zyyyyx yzx yzzzzzx y zx y zy zuy zux zux zux yux yux y zt x y z (,)(,)22(,)(,)22(,)(,)22()()()()()()()xxxxxy zxy zyyyyxyzxyzzzz

8、zxy zxy zuuxuuyuuzt- 12 -并令并令 ， ，把源汇项（井）加入到方程中，得到把源汇项（井）加入到方程中，得到单相渗流单相渗流连续性方程连续性方程 ：0 x 0y0z ()()()()yxzuuuxyzt()()()()yxzuuuqxyzt- 13 -用用Hamilton算子算子 表示，并将达西公式代入，得表示，并将达西公式代入，得到，到，单相渗流微分方程（流动方程）：单相渗流微分方程（流动方程）： 生产井生产井 注入井注入井 关关 井井()()Kqt 0q 0q 0q 流动项源汇项累积项- 14 -这就是地下流体所遵循的渗流规律。这就是地下流体所遵循的渗流规律。总之，渗

9、流基本微分方程包含三部分内容总之，渗流基本微分方程包含三部分内容: 流动项流动项、源汇项源汇项以及以及累积项累积项 渗流基本微分方程由两个定律控制：渗流基本微分方程由两个定律控制： 达西定律与物质守恒定律达西定律与物质守恒定律压力场和饱和度场的分布是流体流动的结果，流压力场和饱和度场的分布是流体流动的结果，流动规律由达西定律来确定，流量由物质守恒来确动规律由达西定律来确定，流量由物质守恒来确定，这两个规律是数值模拟的定，这两个规律是数值模拟的关键关键所在。所在。- 15 -将单相渗流达西定律推广到多相渗流中，得到将单相渗流达西定律推广到多相渗流中，得到三相流动达西公式：三相流动达西公式： 式中

10、，式中，w，o，g水、油、气诸相下标。水、油、气诸相下标。 相对渗透率；相对渗透率；()()()rwrwwwwwwwrorooooooorgrgggggggK KK KPgZK KK KPgZK KK KPgZ UUUrK- 16 - 对油、水、气三相，每一相运用物质守恒定律对油、水、气三相，每一相运用物质守恒定律，可以建立以下，可以建立以下连续性方程连续性方程 ()()()()()()()()wwwwwoooooooogggggooooggogSqtSqtSSqqttUUUU- 17 - 把达西公式代入连续性方程，得到把达西公式代入连续性方程，得到三维三相渗流基本微分方程三维三相渗流基本微分

11、方程 ()()()()()()()()rwwwwwwwooorooooooorgggroooggogoggggooK KSqtK KSqtK KK KqqSStt相对渗透率、毛管力是饱和度函数相对渗透率、毛管力是饱和度函数,而粘度、体积系数而粘度、体积系数孔隙度与密度均是压力的函数孔隙度与密度均是压力的函数- 18 -令令 rwwwwoorooooggroooorggggK KTK KTK KTK KT- 19 -上述上述T的物理意义是反映相邻网格间的流动性，的物理意义是反映相邻网格间的流动性，又叫传导率（又叫传导率（transmissibility) - 20 - oP()()()()()(

12、)()()wwwwwoooooooogggggooooggogSTqtSTqtSSTTqqttoSwSoP求解变量求解变量wPgPgS- 21 -还要补充一些方程还要补充一些方程饱和度饱和度 ：油水毛管力：油水毛管力：油气毛管力：油气毛管力：相对渗透率和毛管力是饱相对渗透率和毛管力是饱和度函数：和度函数：cowowPPPcoggoPPP1wogSSS()rwrwwKKS()rgrggKKS(,)rorowgKKSS()cowcowwPPS()cogcoggPPS- 22 -Stone公式：公式：式中，式中， 束缚水饱和度束缚水饱和度 油水两相流动时油相相对渗透率油水两相流动时油相相对渗透率 油

13、气两相流动时油相相对渗透率油气两相流动时油相相对渗透率 ()() ()()()()()rowwrwwroggrggrorwwrggrowwirKSKSKSKSKKSKSKSwirSrowKrogK- 23 -流体流体PVT性质性质水相：水相： 气相：气相： （Z) 油相：随压力变化，相态会发生变化。油相：随压力变化，相态会发生变化。 三相状态：三相状态： 两相状态：两相状态：岩石的高压物性岩石的高压物性相渗曲线以及流体相渗曲线以及流体PVT性质，以表格形式给出，性质，以表格形式给出，称为静态数据。称为静态数据。wC( )ooooP( )ggooP( )ooP( )wPw( ,)oooosP P

14、( ,)ggoosP P( ,)oosP P( )gPgwBgBgC( )ooCCP( ,)oosCCP P( )rrCC P( )P- 24 -实际求解变量实际求解变量初始化输入数据初始化输入数据ljfinit.txtoPwSgS油相流体性质与压力和饱和压力的关系油相流体性质与压力和饱和压力的关系- 25 -产量项（作源汇项处理）（动态数据）产量项（作源汇项处理）（动态数据）()()()()rwwwBHwBHwooroooBHoBHoggroooBHoBHorgggBHgBHgKqPIPPgZZKqPIPPgZZKqPIPPgZZKqPIPPgZZ- 26 -其中，其中， 在井基准面深度在井

15、基准面深度 上的井底压力。上的井底压力。 节点等效泄油半径。节点等效泄油半径。 井筒半径。井筒半径。X X，Y Y，Z Z 网格节点网格节点x x，y y和和z z方向步长方向步长 PI叫采油指数，叫采油指数， 通常写成通常写成KH,可以直接输入可以直接输入2lnewK ZPIrrerwr0.2erxy BHPBHzK Z井数据又叫动态数据井数据又叫动态数据ljfmodel.txt- 27 -油气藏外边界条件油气藏外边界条件: :第一类边界第一类边界 条件条件 P=f(x,y,z,t)第二类边界第二类边界 条件条件 第三类边界条件第三类边界条件: - 28 -4 外边界外边界 封闭边界封闭边界

16、:尖灭、断层、圈闭尖灭、断层、圈闭 流动边界：边、底水流动边界：边、底水4 内边界，即内边界，即井的生产条件井的生产条件 定产边界条件（定产边界条件（rate constraint):给定产液量、产油量给定产液量、产油量、产气量、注水量、注气量、产气量、注水量、注气量 定压边界条件定压边界条件(pressure constraint):给定井底流压给定井底流压0),(ntzyxpCtzyxp),(- 29 -油藏初始条件油藏初始条件: :( , , ,0)( , , )wwSx y zSx y z( , , ,0)( , , )ggSx y zSx y z- 30 - 31 -在一般情况下在一

17、般情况下流动方程流动方程求不出解析解，只能用离散求不出解析解，只能用离散化方法求数值解。化方法求数值解。 离散化方法，首先要把求解区域按一定的网离散化方法，首先要把求解区域按一定的网格系统进行剖分。网格剖分方法有很多，像直角格系统进行剖分。网格剖分方法有很多，像直角坐标网格、径向网格、角点网格、坐标网格、径向网格、角点网格、PEBI网格等。网格等。目前，商用软件中大都采用目前，商用软件中大都采用矩形网格系统矩形网格系统。离散。离散化方法主要有：有限差分化方法主要有：有限差分（Finite Difference)、有有限元限元(Finite Element)以及流线法以及流线法(Stream L

18、ine)。比较成熟的是比较成熟的是有限差分方法。有限差分方法。有限差分方法大多有限差分方法大多空间变量采用空间变量采用网格块中心（网格块中心（Block Center)差分差分, 时时间变量按间变量按向前差分方法向前差分方法。- 32 -Producer 2Water InjectorProducer 1- 33 -212121212121jjjiiiyyyxxxjjjjjjiiiiiijjjiiiyyyyyyxxxxxxyyyxxx12112112112121212121- 34 -读取网格数据时，R方向循环最快，然后是 方向，最后是Z方向径向网格数据读取规则- 35 -Corner Poi

19、nt- 36 -块中心网格表示的含有断层的倾斜构造块中心网格表示的含有断层的倾斜构造 角点网格表示的含有断层的倾斜构造角点网格表示的含有断层的倾斜构造- 37 -l模型被分为模型被分为20510个网格区域。个网格区域。l网格从（网格从（1，1，1）开始向下倾斜。）开始向下倾斜。l各网格块在各网格块在X方向上是方向上是300ft，在，在Y方向上是方向上是1000ft。l层厚度从上至下分别是层厚度从上至下分别是32，22，20，4，32，4，26，26，4，28ft- 38 -前前20个个TOPS定义网格定义网格 (1, 1, 1) 到到 (20, 1, 1)6855.000 6865.00068

20、75.0006885.0006895.0006905.0006915.0006925.0006935.0006945.0007005.0007015.0007025.0007035.0007045.0007055.0007065.0007075.0007085.0007095.000-下面下面20个个TOPS定义网格定义网格 (1, 2, 1) 到到 (20, 2, 1)6930.0006940.0006950.0006960.0006970.0006980.0006990.0007000.0007010.0007020.0007080.0007090.0007100.0007110.0007

21、120.0007130.0007140.0007150.0007160.0007170.000-下面下面20个个TOPS定义网格定义网格 (1, 3, 1) 到到 (20, 3, 1)7030.0007040.0007050.0007060.0007070.0007080.0007090.0007100.0007110.0007120.0007180.0007190.0007200.0007210.0007220.0007230.0007240.0007250.0007260.0007270.000-下面下面20个个TOPS定义网格定义网格 (1, 4, 1) 到到 (20, 4, 1)713

22、0.0007140.0007150.0007160.0007170.0007180.0007190.0007200.0007210.0007220.0007280.0007290.0007300.0007310.0007320.0007330.0007340.0007350.0007360.0007370.000-下面下面20个个TOPS定义网格定义网格 (1, 5, 1) 到到 (20, 5, 1)7205.0007215.0007225.0007235.0007245.0007255.0007265.0007275.0007285.0007295.0007355.0007365.00073

23、75.0007385.0007395.0007405.0007415.0007425.0007435.0007445.000/ 这完成了第一层这完成了第一层TOPS的定义的定义- 39 -COORD-此关键字定义坐标线此关键字定义坐标线-X1Y1Z1X2Y2Z2 0.0.6825.0000.0.7023.000 300.0000 0.6835.000300.00000.7033.000 600.0000 0.6845.000600.00000.7043.000 900.0000 0.6855.000900.00000.7053.000 1200.000 0.6865.0001200.0000.

24、7063.000./对于对于20 * 10 的网格模型的网格模型, 需要需要21*11 条坐标线，即条坐标线，即 231条。条。-每一条坐标线都用每一条坐标线都用6个数来定义，所以关键字个数来定义，所以关键字COORD包含有包含有1386个数。个数。ZCORN-这个关键字定义网格拐点的深度，数据按照这个关键字定义网格拐点的深度，数据按照X (或或 R) 方向循环最快，然后是方向循环最快，然后是Y (或或 THETA) 方向，最后是方向，最后是Z方向方向的顺序排列。的顺序排列。6825.0006835.0006835.0006845.0006845.0006855.0006855.0006865

25、.0006865.0006875.0006875.0006885.0006885.0006895.0006895.0006905.0006905.0006915.0006915.0006925.0006975.0006985.0006985.0006995.0006995.0007005.0007005.0007015.0007015.0007025.0007025.0007035.0007035.0007045.0007045.0007055.0007055.0007065.0007065.0007075.000./ 对于对于1000个网格，需要个网格，需要 8000 个个ZCORN值值-这

26、些只是角点网格模型的数据体中很少的一部分这些只是角点网格模型的数据体中很少的一部分- 40 -BC（块中心网格）CP（角点网格）网格描述简单网格描述复杂前处理软件不是必须的必须使用前处理软件与其它许多模拟器兼容与很少的其它模拟器兼容很难为非常规构造建立模型能够准确的为非常规构造建立模型网格数据量小网格数据量大无法区分倾斜面和断层区分倾斜面和断层很难如实的建立尖灭和剥蚀面的模型如实的建立尖灭和剥蚀面的模型沿着断层面建立的网格连接不正确，需要手工进行修改。能够正确的模拟沿着断层面层之间的接触关系容易建立径向网格模型若没有前处理软件，则很难建立径向网格模型。- 41 -x1ii正常连接的网格关系正常

27、连接的网格关系粗细网格交界处的粗细网格交界处的流量计算流量计算非正常连接技术1, 11/2riii iiik kyzQBx 两点流量公式：两点流量公式：- 42 -4 传统的结构网格传统的结构网格-即规则矩形剖即规则矩形剖分的网格分的网格,具有如下不足具有如下不足:4 存在较严重的网格取向效应存在较严重的网格取向效应4 不能保证每口井都在网格中心不能保证每口井都在网格中心4 不能很好地描述油藏形态不能很好地描述油藏形态4 粗细网格间的过渡衔接性差粗细网格间的过渡衔接性差4 研究新的网格技术势在必行研究新的网格技术势在必行.结构网格的不足xy好好差差井井好好开发软件- 43 -Gurpinar,

28、 2001Castellini, 2001Wolfsteiner et al., 2002Prevost 2003- 44 -4 常用的非结构网格常用的非结构网格4 各网格中心之间连线构成三角各网格中心之间连线构成三角网网4 PEBI网格（又叫垂直平分网格网格（又叫垂直平分网格Perpendicular Bisectors):各三各三角形边的垂直平分线构成网格角形边的垂直平分线构成网格边界边界4 CVFE网格网格-控制体有限元网格控制体有限元网格(Control volume finite element):各三角形边的中点与重心的连各三角形边的中点与重心的连线构成网格边界线构成网格边界4 待

29、解决的问题待解决的问题4 非结构网格条件下代数方程组非结构网格条件下代数方程组的快速求解的快速求解CVFE网格PEBI网格开发软件- 45 -PEBI- 46 -Locally different grids- 47 -New Faulted GridHeavily Faulted ReservoirNew Fault Traces- 48 - 49 -Geological Model of a Heavily Faulted ReservoirSimulation Model- 50 -1. 1. 概述概述 局部网格加密局部网格加密(Local Grid Refinement(Local G

30、rid Refinement，简称，简称LGR)LGR)技术主技术主要分：要分：静态局部网格加密静态局部网格加密(Static Local Grid Refinement)(Static Local Grid Refinement)局部加局部加密区域在整个计算过程中固定不变的，常用于断层、尖灭密区域在整个计算过程中固定不变的，常用于断层、尖灭、裂缝以及井附近的加密。该技术已经成熟，在多加软件、裂缝以及井附近的加密。该技术已经成熟，在多加软件公司的软件产品中都有这一功能。公司的软件产品中都有这一功能。动态局部网格加密动态局部网格加密(Dynamic Local Grid Refinement)(

31、Dynamic Local Grid Refinement)局部加局部加密区域在整个计算过程中随流体的运动位置而变化，主要密区域在整个计算过程中随流体的运动位置而变化，主要用于追踪驱替前缘。该方法在研究之中，最早是用于追踪驱替前缘。该方法在研究之中，最早是HeinemannHeinemann将动态局部网格用于水驱油藏数值模拟将动态局部网格用于水驱油藏数值模拟 。- 51 -Productive Area 1Productive Area 2Sealing InterbeddingAquiferAquifer- 52 -AquiferProductive AreaVertical Refinem

32、ent Zone- 53 -随驱替前缘移动的动态局部网格加密系统随驱替前缘移动的动态局部网格加密系统- 54 -开发软件- 55 -Single Well BlockHorizontal Well GridComputed Results- 56 -Independent time stepping At any locationSeveral different gridding options are available within one model in the same position- 57 -Window GridBasic GridWindow polygonW in d

33、o w G rid B a s ic G rid- 58 -Automatically and implicitly calculated in/out-flux will be used as a von Neumann boundary condition- 59 -Good match required windowing- 60 - 差分方程差分方程 11111111111()()() ()()() ()()()() nnnnwwwwwwwnnooononooooooonnnnggooggoggngnooggooggVTqSStVTqSStTTqqVSSSStxxxyyyzzzTTTT

34、 - 61 - 111122()()xxxiiiiiiTTT 111122()()yyyjjjjjjTTT 111122()()zzzkkkkkkTTT - 62 - 该差分方程组无论是左端还是右端都是非线该差分方程组无论是左端还是右端都是非线性的，是一组性的，是一组非线性代数方程组非线性代数方程组。非线性最强的。非线性最强的参数是相对渗透率和毛管力。这两个参数在时间参数是相对渗透率和毛管力。这两个参数在时间方面的处理不同，对计算的稳定性以及计算结果方面的处理不同，对计算的稳定性以及计算结果影响最大，也是形成形形色色解法的关键所在。影响最大，也是形成形形色色解法的关键所在。 根据隐式程度的不同

35、，数值模拟的算法有根据隐式程度的不同，数值模拟的算法有IMPES (Implicit Pressure Explicit Saturation) 方方法、法、半隐式半隐式(Semi-Implicit Method)方法方法和和全隐式全隐式方法方法(Fully Implicit Method )。IMPES：隐式求解压力显式求解饱和度。对流度隐式求解压力显式求解饱和度。对流度中与时间有关的量均取上一时间步的值。中与时间有关的量均取上一时间步的值。- 63 -IMPES方法的优点是速度快，缺点是稳定性差方法的优点是速度快，缺点是稳定性差,尤其是锥进问题以及过泡点问题不能很好解决。尤其是锥进问题以及

36、过泡点问题不能很好解决。因此时间步长不能太大。目前，因此时间步长不能太大。目前，IMPES方法在化方法在化学驱模型，尤其是流线法中发挥着作用。学驱模型，尤其是流线法中发挥着作用。全隐式方法全隐式方法主要采用了主要采用了Newton-Raphson迭代技迭代技术术, 对方程中各项系数，包括井系数全部隐式处对方程中各项系数，包括井系数全部隐式处理。全隐式方法较之半隐式方法，具有更好的稳理。全隐式方法较之半隐式方法，具有更好的稳定性，时间步长更长，彻底解决了过泡点问题，定性，时间步长更长，彻底解决了过泡点问题，是八十年代数值模拟技术走向成熟的关键技术，是八十年代数值模拟技术走向成熟的关键技术，也是大

37、型油藏数值模拟的主流方法。也是大型油藏数值模拟的主流方法。 - 64 - 全隐式方法和半隐式方法的根本差别在于：全隐式方法和半隐式方法的根本差别在于： 半隐式方法半隐式方法是将非线性方程直接线性化形成是将非线性方程直接线性化形成线性代数方程组，来求解线性代数方程组，来求解n+1时刻的值。时刻的值。 全隐式方法全隐式方法采用采用Newton-Raphson方法，迭方法，迭代地求解非线性方程组，得到代地求解非线性方程组，得到n+1时刻的变量。时刻的变量。 用一维非线性方程表示：用一维非线性方程表示： 给定初始近似给定初始近似 对函数对函数 在点处在点处 进行一阶泰勒展开进行一阶泰勒展开, 可近可近

38、似地表示为似地表示为000f( ) f( ) f( )() 0 xxx x xf( ) 0 x 0 xf( ) 0 x 0 x- 65 -半隐式：半隐式：全隐式方法：全隐式方法：10()()nnnnf xxxf xx为初值000( )( )f xxxf x- 66 -下面介绍全隐式方法。为了减小计算机舍入误下面介绍全隐式方法。为了减小计算机舍入误差的影响，差的影响， 模型中不直接求压力和饱和度，而模型中不直接求压力和饱和度，而是求其增量。是求其增量。其中：其中：x代表求解变量代表求解变量 n时间迭代步时间迭代步 牛顿迭代步牛顿迭代步1nnxxx1llxxx11nllxxxxl- 67 -求解变

39、量：求解变量： ， ，对于二相状态：对于二相状态： 对于三相状态：对于三相状态： 1111111111111()()()()()()()()()()()()(llllnwwwwwwwwwlolololonooooooooollllglglooggogglgngooggooggoVTqSSStVTqSSStTTqqVSSSSt )oggSSPSwSxSxSoSxPs- 68 -对方程组进行线性化处理。对所有变量进行对方程组进行线性化处理。对所有变量进行Taylor展开，保留线性部分。变量进行运算后展开，保留线性部分。变量进行运算后也保留线性部分，舍去高阶项。形成线性代数也保留线性部分，舍去高阶项

40、。形成线性代数方程组。方程组。1llTTTTTPSwSxPSwSx1llwPSwSxPSwSx - 69 -全隐式方法形成的线性代数方程组的结构全隐式方法形成的线性代数方程组的结构 xbA矩阵结构示意图矩阵结构示意图- 70 - 油藏数值模拟的线性代数方程组是大型稀疏矩阵油藏数值模拟的线性代数方程组是大型稀疏矩阵，其求解效率在油藏数值模拟计算中举足轻重。，其求解效率在油藏数值模拟计算中举足轻重。因为，其计算量占整个计算量约因为，其计算量占整个计算量约70%80%。因。因此，研制快速求解大型稀疏矩阵的新方法成为油此，研制快速求解大型稀疏矩阵的新方法成为油藏数值模拟发展的一个重要方向。藏数值模拟发展的一个重要方向。 线性代数方程组的求解方法主要有三类：线性代数方程组的求解方法主要有三类： 1）直接法：）直接法：Gauss消元、消元、D4排序的排序的Gauss消元消元 2）迭代法：）迭代法：SOR(松弛迭代法）松弛迭代法） 3）共轭梯度类法、预处理共轭梯度类法（）共轭梯度类法、预处理共轭梯度类法（PCG方法）方法） - 71 - 目前最先进的方法是预处理共轭梯度类方目前最先进的方法是预处理共轭梯度类方法。法。BO软件中采用软件中采用RSVP方法。该方法计算快方法。该方法计算快，精度高。，精度高。- 72 -黑油模型油藏数