三维三相达西流法及其在库车坳陷的应用

石广仁，马进山，常军华

（中国石油天然气股份有限公司 石油勘探开发研究院，北京 100083）

三维三相达西流法及其在库车坳陷的应用

石广仁，马进山，常军华

（中国石油天然气股份有限公司 石油勘探开发研究院，北京 100083）

多相达西流法是油气二次运移模拟技术中最复杂的方法，存在三大技术难点：①模型十分复杂，求解比油藏模拟的难度大得多；②收敛性、稳定性不易保证，导致计算结果精度低、计算中途异常停止；③计算机耗时量巨大，在微机上算完一个大型盆地要花几天时间。为此，采用了基于PEBI（直交平分）网的有限体积法、基于牛顿迭代的全隐式法及ORTHMIN（正交极小化）法，摸索出了“正确地选取一个模拟工区的最大允许时间步长”是解决上述三大技术难点的最大关键点，从而攻克了这三大技术难点。研发了三维三相达西流软件，并应用于库车坳陷，获得了较好的油气运聚史模拟结果。

有限体积法；PEBI网；最大允许时间步长；多相达西流法；盆地模拟；二次运移；库车坳陷

烃类运移聚集史模拟是盆地模拟系统中最重要的部分，也是迄今为止技术上最薄弱的环节。当今世界上有4种烃类运聚定量模拟技术［1，2］：多相达西流法、流径法、混合法及侵入逾渗法。其中，多相达西流法是技术最复杂也是采用最早的方法，曾被誉为“主力技术”，但其存在三大技术难点：①模型十分复杂，求解比油藏模拟的难度大得多；②收敛性、稳定性不易保证，导致计算结果精度低、计算中途异常停止；③计算机耗时量巨大，在微机上算完一个大型盆地要花几天时间。为此，笔者攻克了这三大技术难点，研发了三维三相达西流软件，并将此应用于库车坳陷，获得了较好的油气运聚史模拟结果。

1 三维三相达西流软件

目前，世界上有5个多相达西流软件（表1），其中PetroFlow 3D和Temis 3D应用较多。多相达西流软件的模拟器多数为改造的三维三相黑油模型，只是采用的网格及数值解法不同而已。

表1 本三维三相达西流软件与世界上同类软件的主要特性比较Table 1 Com parison ofmajor characteristics of this 3-D 3-phase Darcy flow software w ith other sim ilar softwares in the world

表1中提到了3种类型的网格：①矩形网，是指固定矩形面组成的平行六面体；②三角网，是指可变三角面组成的四面体；③PEBI网，是指具有PEBI（perpendicular and bisection，译为直交平分或垂直平分）性质的可变多角面组成的棱柱体。三维地质体的剖分按棱柱体进行，能确保每个网格块的中心点与任一相邻网格块的中心点之间的连线垂直于这两个网格块的交面且被交面平分。

黑油模型在盆地模拟与油藏模拟中的模拟条件是不同的（表2），它在盆地模拟中的求解远比它在油藏模拟中的求解困难，尤其是在求解的稳定性和收敛性方面，这主要是由盆地模拟可变的模拟范围所引起的。

1.1 三维地质体的网格化

在盆地演化期间，盆地的几何体积因沉降及沉积物填充而大量地全局性地增大，因沉积压实而小量地全局性地变小，因剥蚀而大量或小量地局部性地变小，因断层及裂缝而大量或小量地局部性地增大，甚至因沉积间断而基本不变。因此，地质家所研究的对象实为一个随时间变化的地质体。为了定量地模拟一个盆地，其随时间变化的地质体必须剖分为相当数量的小块（即称网格块）。本软件选用了PEBI网。

1.2 基于PEBI网的有限体积法

三维三相黑油模型的流动方程组［6］太复杂，以致无法对一个实际的地质模型给出解析解，应对之进行离散化而变成线性代数方程组。离散化的方法有有限差分法、有限元法和有限体积法。本软件采用基于PEBI网的有限体积法，从而获得：

1）水流动方程

表2 盆地模拟与油藏模拟的不同模拟条件Table 2 Different conditions between basin modeling and reservoir simulation

2）油流动方程

3）气流动方程

上三式中：v——任一PEBI网格块的符号；

Vi——任一PEBI网格块（由x，y，z表达）；

∰——关于v的三重积分；

div——向量场的散度；

al——等于［K］Krl/（μlBl）（l=w，o，g），因为［K］是绝对渗透率张量（即有方向性的绝对渗透率），故al也是张量；

Krw，Kro，Krg——水、油、气的相对渗透率，可由油层物理试验及改进的Stone公式确定；

μw，μo，μg——水、油、气的粘度，可由PVT函数确定；

Bw，Bo，Bg——水、油、气的地层体积因子，可由PVT函数确定；

grad——两个相邻PEBI网格块之间的梯度；

Φw，Φo，Φg——水、油、气的势，分别是水压力Pw、油压力Po、气压力Pg的函数；

Qw，Qo，Qg——水、油、气的源或汇，若源取正值，若汇取负值；

Rs——溶解气油比，可由PVT函数确定；

t——地质时间；

φ——地层孔隙度；

Sw，So，Sg——水、油、气的饱和度。

由于PEBI网的优点（相邻两网格块中心点的连线恰好就是这两网格块的公共面的法线方向），使方程组（1），（2），（3）的离散化变得比较容易，最终化成为水、油、气三相离散方程组。然而，这个离散方程组只有3个方程，而未知数却有6个（水、油、气压力Pw，Po，Pg和水、油、气饱和度Sw，So，Sg）。所以，还需要其他3个独立的方程。这3个独立的方程是：饱和度平衡方程、油水系统中的毛细管力方程及油气系统中的毛细管力方程。它们分别表示如下：

上三式中：Pcow——油水系统中的毛细管力；

Pcog——油气系统中的毛细管力。

Pcow和Pcog可由油层物理试验确定。

当给出了一个三维地质体的边界条件、初始条件［6］时，就可以通过式（1）—（6）这6个方程，求解出该三维地质体的Pw，Po，Pg及Sw，So，Sg6个未知数。这6个未知数表示了该三维地质体油气二次运移的基本状况。

1.3 基于牛顿迭代的全隐式法

由于式（1）—（3）及其离散方程组都是非线性方程组，因此，①采用全隐式法，而不采用半隐式法，因为半隐式法可能因采用过时的饱和度而得到不合理的计算结果；②采用牛顿迭代法，而不采用其他迭代法，因为牛顿迭代法既有超线性收敛速度，又能保持矩阵的稀疏性。

牛顿算法分两步进行：①生成Jacobi矩阵；②求解问题变为大型线性方程组的Jacobi方程。

Jacobi方程是一个大型、稀疏、非对称的线性代数方程组，即：

式中：A——二维（3I，3I）的系数矩阵，即Jacobi矩阵，其中I为网格块个数；

X——一维（3I）的未知数矩阵，含有3I个未知数（δSw1，δSg1，δPo1，δSw2，δSg2，δPo2，…，δSwI，δSgI，δPoI）；

B——一维（3I）的常数矩阵，含有3I个余量。

1.4 大型线性方程组的求解

求解Jacobi方程，选取存储空间尽可能小的存储方式以及速度尽可能快的迭代算法是十分重要的。对于这类线性代数方程组，迄今尚缺少一致公认最优的迭代算法。正交极小化（ORTHMIN，即orthogonalminimum）法经油藏模拟的实算表明是行之有效的方法，故选择了正交极小化法。正交极小化法的具体步骤为：首先对矩阵A进行近似的不完全因子分解，这比其他方法更接近真解；再经过正交极小化法对AX=B进行求解。正交极小化法的优点是：快速收敛，对迭代参数不敏感，对矩阵带宽不敏感，对矩阵是否对称不敏感。

1.5 解决三大技术难点的关键点

如何解决前述的三大技术难点？除了采用上述的基于PEBI网的有限体积法、基于牛顿迭代的全隐式法及正交极小化法外，还摸索出“正确地选取一个模拟工区的最大允许时间步长”是解决上述三大技术难点的最大关键点。

在理论上，全隐式法具有无条件的收敛性和稳定性。然而，实际应用时发现并非如此。为了确保收敛性和稳定性，合适地选择正交极小化法的时间步长是至关重要的。一般来说，时间步长越大，计算机耗时量就越小，但收敛性和稳定性就越差；时间步长越小，收敛性和稳定性就越好，但计算机耗时量就越大。

在下面的库车坳陷实例中，对时间步长的选取进行了试验，结果发现：①若选取的时间步长太大，会影响模拟的稳定性（即正交极小化发生浮点溢出而导致运行失败）；②若选取的时间步长不足够小，虽然没有出现稳定性问题，但再采用比它小一点的时间步长，会出现两者的模拟结果存在不允许的差异（即收敛性问题）；③若选取的时间步长太小，会造成过长的计算时间而浪费机时。因此，应选取在确保稳定性和收敛性条件下的最大时间步长，这个步长就称之谓最大允许时间步长（0.025 Ma）。这里说明一下如何获得0.025 Ma。令时间步长为Δt。试算表明：当Δt≥0.5 Ma时（例如Δt=5，2.5，1，0.5 Ma），计算发生浮点溢出而导致运行失败；当0.5 Ma＞Δt≥0.025 Ma时（例如Δt=0.25，0.1，0.05，0.025 Ma），虽然计算没有发生浮点溢出而导致运行失败，但使用这4个时间步长算出的计算结果的差异超过了允许的误差；当Δt=0.01 Ma时，计算结果与Δt= 0.025 Ma时的差异没有超过允许的误差；当Δt= 0.005 Ma时，计算结果与Δt=0.025，0.01 Ma时的差异都没有超过允许的误差。这表示Δt=0.025，0.01，0.005Ma都是允许的。但Δt=0.025Ma时的计算机耗时量为9 h 38 min；而Δt=0.01，0.005 Ma时的计算机耗时量分别为19 h 39 min，37 h 8 min。显然，Δt应选取0.025 Ma，它正是这3个Δt中的最大者。这是在平面网格块最大个数取为900时的试算结果。实际上，由于最大允许时间步长基本上不依赖于平面网格块最大个数，所以这种试算可在平面网格块最大个数较小的条件下进行。例如，平面网格块最大个数取为100，也可以获得最大允许时间步长为0.025 Ma的结论，这样可大大节省试算的机时（计算机耗时量为59 min）。

2 在库车坳陷的应用

2.1 地质背景

库车坳陷位于塔里木盆地的北部，北倚南天山造山带，南界塔北隆起，东西长约380 km，南北宽约30～145 km，西宽东窄，面积约30 000 km2。截至2003底，该坳陷内已钻探井102口，其中获工业油流井2口，获工业气流井56口。油气藏类型有凝析油气藏、干气气藏和油藏3种。在该坳陷中已经发现2个油田及13个气田。库车坳陷已成为中国最富的天然气聚集区之一（图1；表3）。

图1 库车坳陷构造及油气田分布［8］Fig.1 Structures and oil/gas fields in the Kuqa Depression［8］

表3 库车坳陷地层年龄Table 3 Geological ages of formations in the Kuqa Depression

库车坳陷是在海西运动基础上于晚二叠世开始发育起来的中、新生代叠合坳陷。坳陷演化具有多阶段性和复杂性，中、新生代陆相地层发育齐全，后期褶皱-冲断变形改造强烈。主力烃源岩为中、上三叠统湖相泥岩和中、下侏罗统湖泊沼泽相煤系地层，平均厚度达600 m。暗色泥岩有机碳含量多在1.5％～5.0％，煤层的残余有机碳均值为58％。源岩干酪根类型以腐殖型（Ⅲ型）为主，约占90％；其次是腐泥-腐殖型（Ⅱ）型，约为10％。源岩现今主体成熟度Ro在0.7％～1.8％，处于成熟-高成熟阶段，以生气为主。第三系巨厚（100～1 000 m）的盐-膏-泥优质区域盖层封盖了整个坳陷，而这个区域盖层之下分布着古近系和白垩系的良好砂岩储集层。优越的源、盖条件形成了克拉2、大北等大型气藏以及牙哈、英买7、羊塔克等中型凝析油气田［7～9］。

2.2 模拟条件

石广仁和张庆春［10］曾使用简易的方法对库车坳陷做过油气运聚史模拟，介绍并使用了两类输入参数：①坳陷内71口已钻井的数据以及从地震剖面上所取的577口人工模拟井的数据，共计90 000余个；②地史、热史、成岩史、生烃史及排烃史的模拟结果。该简易方法视地层为一个具有厚度的平面，故为拟三维，而不是真三维；只算油、气二相，不算油、气、水三相；使用简易的流体流动方法，故不属于多相达西流法。本文讨论的才是三维三相达西流法，但也同样使用了上述的两类输入参数。

运移聚集史模型是正演模型，故历史模拟是从古到今进行的。一般来说（假设盆地的最大埋深是在现今），三维地质体的网格块个数随之由小到大变化着，则现今网格块个数为最大，这一个数是很大的。被模拟的库车坳陷的面积约为30 000 km2，东西向最长为383 km，南北向最长为145 km，最大历史埋深为10.4 km。在三维地质体的剖分中，网格块垂向长度取为100 m，平面网格块最大个数取为900，则实际的最大网格块个数I=31 660（节点）。因此，式（7）中的二维系数矩阵A（3I，3I）最大时为A（94 980，94 980）。

计算机耗时量与最大网格块个数I有关，还与地层个数、烃源层个数，尤其是模拟历史、模拟时间步长有关。被模拟的库车坳陷有8个地层（表3）、2个烃源层（T，J），模拟历史为208～0 Ma，模拟时间步长取最大允许时间步长（0.025 Ma），计算机耗时量为9 h 38 min（这是用户可以接受的）。该实例表明，目前微机可承受30 000节点的三维三相达西流模拟。

2.3 模拟结果

使用模拟结果，可以用三维立体图及其他形式的图件来显示各地层的油、气聚集量史。但是这里由于篇幅所限，①为了简明起见，改用曲线形式表示各地层的油、气聚集量史（图2）；②因为白垩系是最大的油气储集层（图2），用三维显示其现今（0 Ma）的油聚集量（图3a）以及气聚集量（图3b）；③为了分析古近系和白垩系的储气层，给出了平面等值线图（图4）。

图2 库车坳陷各相关地层油、气聚集量史Fig.2 Oil and gas accumulation histories of each related formation in the Kuqa Depression

2.4 地质分析

2.4.1 坳陷的模拟总量与勘探成果吻合

表4列出了库车坳陷生、排烃的模拟总量，其中排油系数=排油量/生油量，排气系数=排气量/生气量［10］。在排烃史模拟结果的基础上，这次采用三维三相达西流软件继续进行运聚史的计算，获得了聚烃的模拟总量，其中聚油系数=聚油量/排油量，聚气系数=聚气量/排气量（图2；表4）。在以前库车坳陷的圈闭评价中，曾获得了该坳陷的油气地质储量（含预测、控制、探明储量）（表4）［10］。从表4可见，模拟的聚油量（409.5× 106t）稍大于油的地质储量（380×106t），相对误差为7.8％；模拟的聚气量（2 277×109m3）也稍大于气的地质储量（2 100×109m3），相对误差为8.4％。这表明了模拟结果的合理性，也表明库车坳陷还有深入勘探的前景。此外，模拟的聚气量为模拟的聚油量的5.56倍，气的地质储量为油的地质储量的5.53倍，两者十分接近，均说明库车坳陷是富气的。

图3 白垩系在0 Ma时的累计油、气聚集量Fig.3 Cumulative oil and gas accumulations in the Cretaceous（K）at0 Ma

图4 在0 Ma时a）古近系累计油聚集量和b）白垩系累计气聚集量等值线Fig.4 Cumulative oil accumulation in the Paleogene（a）and cumulative gas accumulation in the Cretaceous（b）at 0 Ma

表4 库车坳陷烃类生、排、聚总量（模拟）及地质储量Table 4 Simulated total amount of hydrocarbons generated，expelled and accumulated in the Kuqa Depression and its reserves in p lace

由图2可见，坳陷的油、气聚集时期主要是新近纪（24～2 Ma）。烃类排、聚时期与主要构造的形成时期基本是吻合的。这个模拟结果与王庭斌教授通过多方面的分析所得的结论（即“库车坳陷的气藏主要形成、定型于新近纪”）［11］基本符合。

2.4.2 各相关地层的模拟总量与勘探成果吻合

从图2a可见，新近系（N1j，N1-2k，N2k）无油；古近系（E）只有极少量的油；而白垩系（K）才是最大的储油层（图3a），约占坳陷储油量的78％。从图2b可见，新近系（N1j，N1-2k）只有极少量的气；而白垩系（K）又是最大的储气层（图3b），储气量约占坳陷的74％；其次是古近系（E），储气量约占坳陷的15％。这些分层的模拟结果与勘探成果吻合。

2.4.3 模拟的油气田分布与勘探成果吻合

观察一下运聚模拟结果（图4）与已发现的油气田（图1）的符合程度。对比表明：①图4a上最大的油聚集带恰好是牙哈（YH）油田（图1），其余较小的油聚集带也有深入勘探的前景；②图4b上较大的气聚集带恰好是大北（DB）、克拉2（KL2）、羊塔克（YT）、玉东（YD）、英买7（YM7）、红旗（HQ）、提尔根（T）及牙哈（YH）等气田（图1）。十分明显，模拟结果与目前勘探成果吻合。

2.5 结论和讨论

由库车坳陷的应用实例可以得出如下结论：①采用本三维三相达西流来研究油气二次运移是有效可行的；②库车坳陷油气运移的模拟结果不仅与实际情况符合，而且表明还有深入勘探的前景。

然而，由于油气运移机理的复杂性及地质因素的不确定性，要想建立一个较完美的全定量模型并得到较好的应用效果，对地质家来说仍是一个严重的挑战。要想使用上述方法得到一个较好的应用，应注意两点：①尽可能把地史、热史、成岩史、生烃史和排烃史模拟得精确，这是油气运移定量模拟的基础；②尽量取准油气运移的敏感性参数。

1 Hantschel T，Kauerauf A I.Fundamentals of basin modeling and petroleum systemsmodeling［M］.Berlin：Springer-Verlag，2009

2 石广仁.油气运聚定量模拟技术现状、问题及设想［J］.石油与天然气地质，2009，30（1）：1～10

3 袁益让，韩玉笈，赵卫东，等.多层油资源运移聚集的数值模拟和实际应用［J］.应用数学和力学，2002，23（8）：827～836

4 李长峰，袁益让.三维两相渗流驱动问题迎风区域分裂显隐差分法［J］.计算数学，2007，29（2）：113～136

5 Mello U T，Rodrigues JR P，Rossa A L.A control-volume finiteelementmethod for three-dimensional multiphase basin modeling［J］.Marine＆amp；Petroleum Geology，2009，26（4）：504-518

6 石广仁.油气盆地数值模拟方法（第三版）［M］.北京：石油工业出版社，2004

7 周兴熙.库车油气系统成藏作用与成藏模式［J］.石油勘探与开发，2001，28（2）：8～10

8 赵靖舟，戴金星.库车前陆逆冲带天然气成藏期与成藏史［J］.石油学报，2002，23（2）：6～10

9 王庭斌.中国气田的成藏特征分析［J］.石油与天然气地质，2003，24（2）：103～110

10 石广仁，张庆春.库车坳陷的油气运移全定量模拟［J］.地球科学——中国地质大学学报，2004，29（4）：391～396

11 王庭斌.中国气藏主要形成、定型于新近纪以来的构造运动［J］.石油与天然气地质，2004，25（2）：126～132

（编辑 李 军）

3-D 3-phase Darcy flow method and its application to the Kuqa Depression

Shi Guangren，Ma Jinshan and Chang Junhua
（PetroChina Research Institute of Petroleum Exploration＆amp；Development，Beijing 100083，China）

Themulti-phase Darcy flow method is the most complicated among all the modeling methods for secondary hydrocarbonmigration.It has threemajor technical problems：（1）itsmodel is too complicated to easily generate a solution，at leastmuch difficult than that in reservoir simulation；（2）it is not easy to ensure the convergence and stability of themodel and thus rendering a calculation of low accuracy or causing abnormal shutdown；and（3）it is very time-consuming as the calculation of a larger basinmay take several days to complete.To tackle these problems，we tried the finite volumemethod based on PEBI（perpendicular and bisection）cells，the fully implicit solution based on the Newton iteration，and the ORTHMIN（orthogonalminimum）algorithm，and came to a conclusion that the key to deal with the problems is to correctly select the allowable maximum time step for a simulated block.Finally，we developed a program of3D 3-phase Darcy Flow.An application of the software to the Kuqa Depression had gained positive results.

finite volume method，PEBI gridding，allowable maximum time-step，multi-phase Darcy flow，basin modeling，secondarymigration，Kuqa Depression

TE19 < class="emphasis_bold">文献标识码：A

A

0253-9985（2010）04-0403-07

2010-05-28。

石广仁（1940—），男，教授级高级工程师，地学定量。

中石油科技部项目（2008A-0602）。