1. 进入油藏数值模拟的世界#

1.1. 为什么要用数值模拟#

数值模拟（英文：Numerical Simulation）也叫计算机模拟。通俗点来说，就是建立与真实问题近似的数学模型，或者说建立能足够合理模拟实际的问题的数学模型，然后利用电子计算机来求解这个数学模型，实现对工程问题和物理问题乃至自然界各类问题研究的目的。这些数学模型往往是基于基础的理论模型，但是最后都会需要求解大量的未知数，少则数万，多则上百万上千万个，有多少未知数就有多少个方程。如此庞大的问题是人们手工计算不可能完成的，必须借助计算机来实现。用计算机来编写程序来自动求解问题。这样的程序就被称之为模拟器。数值模拟技术最早诞生于1953年Bruce G.H和Peaceman D.W模拟了一维气相不稳定径向和线形流。受当时计算机能力及数值解法限制，数值模拟技术只是初步应用于解一维一相问题。两相流动模拟诞生于1954年，West W J和Garvin W.W模拟了油藏不稳定两相流。早期的时候，受计算机硬件和编程语言的限制，模拟器都是针对特定的问题编写的，解决一个问题可能就需要编写一次程序。使用模拟器也需要经过专业的培训才能熟练操作。随着技术的逐步成熟和软件的不断完善，慢慢地产生了比较通用的模拟器。由于所使用的理论模型逐渐的统一，这些模拟器的大部分功能都趋于相同，仅仅可能是输入格式有些差异，有利于学习和研究油藏数值模拟应用。

在数值模拟出现之前，人们已经尝试过用建立物理模型来进行模拟。物理模型又可以分为比例模型和基本模型。比例模型是将实际问题的各种尺寸数量都按一定比例的缩小给定。通过研究比例模型，将其结果应用到解决实际问题上。比例模型是一种近似模型，很多参数要按比例来给定非常困难，建立完整的比例模型几乎是不可能的。基本模型是把实际问题的进行分解，对分解后的子问题选取具有一定特征问题，然后建立物理模型进行研究。基本模型很难解决整个的实际问题，但是它能解决机理上的一些问题，帮助科学家建立和完善理论基础。比如，数值模拟所使用的数学模型就是基于这些理论基础的。

数学模型是数值模拟的基础模型，但是数学模型有可能是复杂的偏微分方程，还需要选择合适的数值解法来近似求解，这个过程得到了可以用计算机求解的数值模型或者也可以称为数值方程。最后才是用计算语言进行编程，得到解决这个问题的模拟器。

1.2. 什么是油藏数值模拟#

既然我们的目的是要解决实际问题，那么就需要了解需要解决哪些问题，首先需要明确解决这些问题的目的。油藏所面临的问题就是油藏不确定性，即地下储层的构造不确定、油气资源储量不确定、生产动态不确定，等等这些不确定性问题，造成了开发生产的不确定。那么什么样的开发手段和途径能得到更高的最终采收率，就成为了油藏研究最终的需要解决的问题。而这个问题是经典的手段例如油藏工程方法不能解决的。因此需要借助数值模拟的方法来研究。 油藏数值模拟就是以数学方法建立模拟地下流体流动得数学模型，然后选择合适的离散化方法将数学模型变成可以求解的非线性方程组，再利用计算机来求解这个非线性方程组，得到可以反映油藏动态的参数。同时，可以反复利用这个模型，调整参数，对油藏的动态进行全面的预测和评估。油藏数值模拟是目前最先进地进行油藏动态研究的工具。

油藏数值模拟的数学模型也分为不同的类型。可以从需要研究的油藏特征来划分。对于常规的油气藏，使用黑油模型就足够了。如果要研究油气藏中的组分变化，或者油气藏中的组分变化对生产的影响较大不能忽略，例如凝析气藏、注二氧化碳开发等，就需要使用组分模型。而对于哪些稠油油藏，原油粘度非常高，需要靠注入高温蒸汽才能驱动原油流动的油藏，就需要使用热采模型。另外，如果在稠油油藏的开发中需要注入化学剂，考虑化学反应，就需要使用热采-组分模型。对于这些模型的选择，只要符合研究目的即可。没有必要用热采模型去模拟常规的油气藏，虽然在技术上也许是可行的，但是从参数设置和效率上来说，都远远没有使用黑油模型方便。

可以这样说，黑油模型是应用最广泛、最基础的油藏数值模拟数学模型。本文介绍中我们使用的是最基本的油藏数值模型——黑油模型来进行介绍。为了进一步了解更复杂的油藏数值模拟模型和使用方法，必须全面地理解和运用黑油模型。黑油模型是油藏数值模拟的基础模型。万丈高楼平地起，没有人能够在完全不懂黑油模型的情况下，把组分模型或者热采模型用的非常好。只有熟练的掌握黑油模型的理论和使用才能更好的研究和应用其他复杂的模型。

经典的黑油模型有一些假设条件，区别于其他的数学模型：

  • （1）油藏中有明显区别的相：油、水、气。也成为油相、水相、气相。

  • （2）其中油和水是不混相的，即它们之间不存在质量传递作用（传质），也不存在互相的相态变化。

  • （3）油藏储层中的油是由两种组分组成的混合物。这两种组分分别是地面油和标准条件下的气。即在不同压力下，油和气会存在混相的，即有传质和相态变化。

  • （4）油藏储层和流体的温度都是不变的，处于热力学平衡。储层和流体的物性主要取决于压力和体积的关系。

这里我们先介绍黑油模型的假设条件和简要介绍，详细的物性参数将会在后面章节中陆续介绍。黑油模型和组分模型的区别主要是在PVT（压力-体积-温度）物性参数的表达上。

刚才说的那些问题都比较笼统，下面通过一个例子我们来具体的说明一下油藏数值模拟如何来解决油藏的动态问题。假设现在我们发现了下图所示的一个油藏。以现有资料来看没有气层，也就是说不存在气顶，但是这个不确定。这个油藏含油面积约为10km2，油层的厚度50m，这些都是初步估计得数据，包括储层的物性参数都需要在开发中进一步的确定。储层是一个典型的背斜构造，东部（右侧）有一条断层，断层右侧的储层还没有完全的分析，有一定的可能性存在小油藏。地质分析发现该油藏具有边底水，但是目前水体的大小和能量还不确定。

_images/rs001.png

我们尽量把这个油藏的开发思路简化一下，假设需要确定的目的只有一个，这个油藏如何开采才能让10年后累积采出的油量最大。虽然问题只有一个，但是解决的过程是很复杂的。因为很多的参数具有不确定性，无法从一开始就确定一个准确的方案。如果没有油藏数值模拟作为手段，那么只能是逐步钻井，靠井的数据来进一步的深入研究。下图中列出了一些不确定的问题，这些都将在油藏实际生产时产生影响。

_images/rs002.png

钻井的成本每口井少则数百万，储层深一些的井的成本都能上千万，还需要修建管线、集输站等地面设施，所需要的总成本都是数以亿计。时间是不能倒退的，如果钻井失败，如果储层物性并不像我们预期的那种，又或者储层中的水体能量非常大，我们一开始所设计的方案没有考虑这个因素，导致生产井水淹，后期所产生的影响简直是不堪设想的。一旦油藏发生水淹，要想再挽救回来，这后期的投入是无法估计得。

如果油藏的上面地表有不可移动的大范围的物体如城市、村庄、山地等等，就不能靠从地面向储层钻直井，而是需要采用定向井技术，从这些范围之外向储层钻井。如下图所示：

_images/rs003.png

因此制定油藏的开发方案、评估影响开发的不确定性因素的影响、为决策团队提供规避风险的依据，这些都迫切地要求有这样一种手段，能比现有的方法都更合理更接近真实情况，能重复实验反复验证，能节约大量成本，包括研究的时间。油藏数值模拟技术就应运而生了。

对于上述的问题都可以用油藏数值模拟进行研究，例如井位和产量的设置，可以调整油藏数值模拟所使用的模型数据参数，设置不同的方案，计算完成后对产量和动态进行对比即可优选出方案。对于不确定性大的参数的分析也是如此，例如水层（水体）的大小分析、断层封堵性的分析，都可以设置不同的参数进行模拟研究，通过对结果参数分析，了解这些不确定性对油藏动态的影响，尽量控制和减小风险，制定开发方案。此外，在开发方案实施一段时间后，可以将实际的生产数据和数值模拟计算的结果进行对比，调整油藏模型的参数，尽可能让计算结果与实际的生产数据匹配，这样修正后的油藏模型所计算的结果就更合理。这个过程也是油藏数值模拟中很重要的一个环节，被称为历史拟合。

通过历史拟合，我们可以随时根据实际的油藏动态修正油藏模型，保证随后的预测更符合实际油藏的生产，增强了油藏模型的可用性。这个过程是不断打磨油藏模型的过程。之前我们讲过，在研究初期，油藏的地质认识是有非常大的不确定性的，因为我们只能使用间接的方法了解地下的情况。而使用数值模型，我们可以创建大量的油藏模型，每个模型反映这些不同的地质认识，使用这些模型进行历史拟合，哪些模型更能容易让计算结果拟合上实际生产数据，那么这些模型的地质情况就越接近地下的实际情况。

1.3. 油藏数值模拟的数学模型#

建立油藏数值模拟首先是需要建立数学模型，基于质量守恒、动量守恒和能量守恒这些基础定律来建立相应的数学模型。油藏数值模拟的流动问题是属于多孔介质中的流动问题，实际的应用中，我们不是在研究流体的流动问题，而更像是在研究扩散问题（方程的形式和属于扩散方程），目前还用不到动量守恒定律。不过如果是研究热采模型的数值模拟，就要用到能量守恒定律。

1.3.1. 质量守恒方程#

1.3.1.1. 单相流的方程#

质量守恒定律是建立油藏数值模拟相关数学模型的基础。由于地下储层中有油水气不同的流体，我们先解决单相流体的问题。先考虑一个基础模型问题，我们用岩芯来研究流动特征问题（这个是紧扣上面的一点）。取一段水平放置圆柱型岩芯，横截面积是A，物性均质，孔隙度为 \(\phi\) ，两端有压差，内部有单相流体（可以是混合物，不过物性是均匀的）的流动，流体得密度为ρ。如下图所示。假设其单位时间内单位面积上通过的质量流量是 \(m^{x}\) ，这里的上标x表示水平方向（X方向）。我们取这个岩性模型中的任意一个单元，由于是均质的一维模型，我们不考虑高度的变化，这个单元用∆x表示。

_images/rs004.png

注意在建立的数学模型中，我们取得任意单元∆x范围内的体积称为控制体积。控制体积是一个数据概念，可以认为是包含原模型物性的一个微单元。它没有实际体积的概念，因为如果控制体积很小，小于岩芯孔隙的级别，那么孔隙度对它就没有意义了。这里的孔隙度等属性对于这个控制体积仍然是有意义的。我们可以把岩性的物理性质参数都赋予它。在∆t时间内通过控制体积表面，即长度为x处流入的质量流量为：

\[A\cdot \bigtriangleup t\cdot m_{x}\]

在∆t时间内通过控制体积表面流出，即从长度为x+∆x流出的质量流量为：

\[A\cdot \bigtriangleup t\cdot m_{x+\bigtriangleup x}\]

同时，这个流入和流出的差值等于控制体积内质量累积的总量（称为累积量）。这个累积量分为两个部分。首先是∆t内由于压缩性引起的质量增量为如下形式。

\[\frac{\partial \rho \phi \bigtriangleup V}{\partial t}\]

上述表达式在整个公式中又称为累积项。如果控制体上存在强度为q ̃的源汇点。这个源汇点类似于岩芯这个部分有一个点，在流出流体。源汇点是用来模拟油藏中的井、水体等的模型。如果是注入称为源，而产出则称为汇。一般来说，油藏数值模拟中对于从油藏中流出是的方向是正向，而注入或流入是负向。因此q是一个向量，且强度的定义为单位时间单位体积的流体的质量。由此源汇点引起的质量减少为：

\[\bigtriangleup t \cdot \bigtriangleup V \cdot q\]

因此我们可以得到如下的一个等式：

\[A\cdot \bigtriangleup t\cdot m_{x+\bigtriangleup x} - A\cdot \bigtriangleup t\cdot m_{x} = \frac{\partial \rho \phi \bigtriangleup V}{\partial t} + \bigtriangleup t \cdot \bigtriangleup V \cdot q\]

这里再强调一下公式中的正负号。流出为正、流入为负，即汇为正、源为负。为了便于后面的推导，我们把上述方程中引入速度项，替换掉质量流量。令：

\[m_{x}=\rho u_{x}\]

上式中的 \(u_{x}\) 为X方向速度（也就是达西速度或表观速度），同时对方程取取∆x \(\longrightarrow\) 0的极限，就得到了该问题的质量守恒方程。于是质量守恒方程的形式变为：

\[-\frac{\partial (\rho u_x)}{\partial x}=\frac{\partial (\rho \phi )}{\partial t} + q\]

基于一维方向的推导思路，我们可以向二维和三维将方程进行扩展。二维坐标系和三维坐标系下的控制体积示意图见下图：

_images/rs005.png

下面给出三维直角坐标系下，任意多孔介质中单相流动的质量守恒方程：

\[-(\frac{\partial (\rho u_x)}{\partial x}+\frac{\partial (\rho u_y)}{\partial y}+\frac{\partial (\rho u_z)}{\partial z})=\frac{\partial (\rho \phi )}{\partial t} + q\]

可以将上述的方程可以写成散度算子的形式，即：

1.3.1.2. 辅助方程#

如果需要将多相流扩展为多相流，需要借助一些辅助方程。这里引入的辅助方程除了一些基本方程外，大多数都是黑油模型中的辅助方程。对于多相流来说，孔隙空间并不被某一相所填满，因此为油水气各相引入饱和度这个变量，饱和度就是这个相在孔隙中所占的体积跟总孔隙体积之比。对于油气水三相来说，其饱和度满足如下方程：

\[S_{o}+S_{w}+S_{g}=1\]

其中S_{o}、S_{w}、S_{g}分别是油水气的饱和度。

除了饱和度以外，我们还需要引入体积系数这个参数。体积系数也称为地层体积系数，它是黑油模型中最简单也最重要的状态方程，主要是为了实现地面体积与地下体积的转换。其定义为：

\[B_{l}=\frac{V_{l}^{r} }{V_{l}^{std}}\]

其中： \(V_{l}^{r}\) 是相 \(l\) 在地下的体积， \(l=o,w,g\) ；\(V_{l}^{std}\) 是相 \(l\) 在地面标况下的体积。

使用体积系数可以把流体的地下密度转换为由地面密度的表示形式。根据黑油模型的假设，水是相对独立的相。可以直接给出水的密度转换公式如下：

\[\rho_{w}=\frac{\rho_{w}^{std} }{B_{w}}\]

其中： \(\rho_{w}\) ，地下的水密度； \(\rho_w^{std}\) ，地面标况下的水密度；

油和气的情况要复杂一些，因为油和气存在传质和相态变化。如果考虑最完全的情况，在地下同时存在油气两相。此时油中溶解了一部分气，这些气在地下是存在于油中，但是到地面以后，油和气会分离。此时分离出的气和油的体积之比为溶解气油比。同时，在地下也有一部分油溶解进了地下的自由气中。当这些包含油的气到达地面后，气和油会分离。此时分离出的油和气的体积之比为挥发油气比，其描述公式如下：

\[\rho_{o}=\frac{\rho_{o}^{std}+R_{s}\cdot\rho_{g}^{std}}{B_{o}}\]
\[\rho_{g}=\frac{\rho_{g}^{std}+R_{v}\cdot\rho_{o}^{std}}{B_{g}}\]

其中： \(R_s\) 和 \(R_v\) 分别表示溶解气油比和挥发油气比。

1.3.1.3. 多相流的情况#

将辅助方程和单相流的质量守恒方程相结合，就可以推导出多相流的质量守恒方程。以上述各式为基础，可以得到不同相 \(l\) 的质量守恒方程如下：

\[-\nabla\cdot{\dot{\mathbin{m}}}_\mathbin{l}=\frac{\partial\left(m_l\right)}{\partial t}+{\widetilde{q}}_l\]

其中： \(m_l\) ，是多孔介质中单位体积下的 \(l\) 相的质量； ：math:{dot{mathbin{m}}}_mathbin{l} ，是 \(l\) 相的的质量通量（单位时间内单位面积上的质量流量）。此外，定义体积流量向量 \(\mathbin{u}_l\)，即单位时间没单位面积上的体积流量，可得：

\[{\dot{\mathbin{m}}}_\mathbin{l}=\mathbin{u}_l\cdot\rho_l\]

对于水相可以得到下述公式：

\[-\nabla\rho_w\mathbin{u}_\mathbin{w}=\frac{\partial\left(\rho_w\phi S_w\right)}{\partial t}+{\widetilde{q}}_w\]

两边同时除以 \(\rho_w^{std}\) ，即该相在地面标况时的体积，公式变为：

\[-\nabla\frac{1}{B_w}\mathbin{u}_\mathbin{w}=\frac{\partial}{\partial t}\left[\frac{1}{B_w}\phi S_w\right]+\frac{{\widetilde{q}}_w}{\rho_w^{std}}\]

令 \(\frac{{\widetilde{q}}_w}{\rho_w^{std}}=q_w\) ，把质量流量都转换为体积流量。可以得到油水气三相条件下的水相质量守恒方程：

\[-\nabla\frac{1}{B_w}\mathbin{u}_\mathbin{w}=\frac{\partial}{\partial t}\left[\frac{1}{B_w}\phi S_w\right]+q_w\]

对于油相，其质量通量可以表示如下：

\[{\dot{\mathbin{m}}}_\mathbin{o}=\mathbin{u}_o\cdot\rho_o+\mathbin{u}_g\cdot\rho_o^{std}\cdot\frac{R_v}{B_g}\]

跟水相类似，可以推出油水气三相条件下的油相质量守恒方程：

\[-\nabla\cdot\left[\frac{1}{B_o}\mathbin{u}_o+\frac{R_v}{B_g}\mathbin{u}_g\right]=\frac{\partial}{\partial t}\left[\phi\left(\frac{1}{B_o}S_o+\frac{R_v}{B_g}S_g\right)\right]+q_{fo}+q_g\cdot R_v\]

同理可得油水气三相条件下气相的质量守恒方程：

\[-\nabla\cdot\left[\frac{R_s}{B_o}\mathbin{u}_o+\frac{1}{B_g}\mathbin{u}_g\right]=\frac{\partial}{\partial t}\left[\phi\left(\frac{R_s}{B_o}S_o+\frac{1}{B_g}S_g\right)\right]+q_{fg}+q_o\cdot R_s\]

其中： \(q_{fo}\) 和 \(q_{fg}\) 分别表示自由由相和自由气相的体积流量； \(q_g\cdot R_v\) ，表示地下的气相到地面后分离出的油； \(q_o\cdot R_s\) ，表示地下的油相到地面后分离出的气。

1.3.2. 达西定律#

在建立质量守恒方程的时候我们使用了一个变量，叫达西速度。达西速度来自达西定律。这是科学家达西于1856年研究水动力学的时候发现的，对于单相流动达西定律可以表示为微分形式如下：

\[\mathbin{u}=-\frac{k}{\mu}\left(\nabla p+\rho\vec{g}\right)\]

上式中：

  • \(\mathbin{u}\) ，达西速度，即单相的流量；

  • \(k\) ，多孔介质中的渗透率张量，此渗透率为绝对渗透率；

  • \(\mu\) ，流体粘度；

  • \(\nabla p\) ，压力的梯度；

  • \(\vec{g}\) ，重力加速度向量。

在三维直角坐标系中，我们可以主要考虑张量的主对角元素，也就是XYZ方向的正方向。因此k可以用如下形式表示：

\[\begin{split}k=\begin{bmatrix} k_{xx} & & \\ & k_{yy} & \\ & & k_{zz}\end{bmatrix}\end{split}\]

将上述公式扩展为多相流，需要考虑每相的有效渗透率、粘度、压力以及密度等等特有的参数。推导过程比较复杂，我们这里不再赘述，直接给出最终的达西定律对于多相流的描述形式：

\[\mathbin{u}_l=-\frac{kk_{rl}}{\mu_l}\left(\nabla p_l+\rho_l\vec{g}\right)\]

即单相流的公式中的参数都加了一个下标 \(l，l=o,w,g\) ，分别代表油水气各相。 \(k_{rl}\) 为某相的相对渗透率。它跟绝对渗透率的乘积 \(kk_{rl}\) 就是该相的有效渗透率。

1.3.3. 基本流动方程#

将达西定律和质量守恒方程结合起来，可以得到多孔介质中流体的流动方程。单相流体的基本流动方程的表达形式为：

\[\nabla\cdot\rho\frac{k}{\mu}\left(\nabla p+\rho\vec{g}\right)=\frac{\partial\left(\rho\phi\right)}{\partial t}+\widetilde{q}\]

将多相流的达西定律公式代入多相流的物质守恒方程，可以得到油水气三相的流动方程如下：

水相方程：

\[\nabla\cdot\left[\frac{kk_{rw}}{{B_w\mu}_w}\left(\nabla p_w+\rho_w\vec{g}\right)\right]=\frac{\partial}{\partial t}\left[\frac{1}{B_w}\phi S_w\right]+q_w\]

油相方程：

\[\nabla\cdot\left[\frac{kk_{ro}}{B_o\mu_o}\left(\nabla p_o+\rho_o\vec{g}\right)+\frac{R_vkk_{rg}}{B_g\mu_g}\left(\nabla p_g+\rho_g\vec{g}\right)\right]=\frac{\partial}{\partial t}\left[\phi\left(\frac{1}{B_o}S_o+\frac{R_v}{B_g}S_g\right)\right]+q_{fo}+q_g\cdot R_v\]

气相方程：

\[\nabla\cdot\left[\frac{R_skk_{ro}}{B_o\mu_o}\left(\nabla p_o+\rho_o\vec{g}\right)+\frac{kk_{rg}}{B_g\mu_g}\left(\nabla p_g+\rho_g\vec{g}\right)\right]=\frac{\partial}{\partial t}\left[\phi\left(\frac{R_s}{B_o}S_o+\frac{1}{B_g}S_g\right)\right]+q_{fg}+q_o\cdot R_s\]

我们可以定义方程中的传导率项 \(\lambda_l\) 为如下公式：

\[\lambda_l=\frac{kk_{rl}}{{B_l\mu}_l}\]

代入上面三个公式，可以得到如下简化形式的方程：

\[\nabla\cdot\left[\lambda_w\left(\nabla p_w+\rho_w\vec{g}\right)\right]=\frac{\partial}{\partial t}\left[\frac{1}{B_w}\phi S_w\right]+q_w\]
\[\nabla\cdot\left[\lambda_o\left(\nabla p_o+\rho_o\vec{g}\right)+R_v\lambda_g\left(\nabla p_g+\rho_g\vec{g}\right)\right]=\frac{\partial}{\partial t}\left[\phi\left(\frac{1}{B_o}S_o+\frac{R_v}{B_g}S_g\right)\right]+q_{fo}+q_g\cdot R_v\]
\[\nabla\cdot\left[R_s\lambda_o\left(\nabla p_o+\rho_o\vec{g}\right)+\lambda_g\left(\nabla p_g+\rho_g\vec{g}\right)\right]=\frac{\partial}{\partial t}\left[\phi\left(\frac{R_s}{B_o}S_o+\frac{1}{B_g}S_g\right)\right]+q_{fg}+q_o\cdot R_s\]

上述推导中并没有考虑源汇点的类型，而是直接引入源汇对系统的影响，即源汇点产生的流量变化。实际的油藏数值模拟方程还需要将源汇点分成不同的类型，有井、水体以及其他可以供应流量的模型。为了能方便的说明后面继续推导的思路，不至于把问题弄得太复杂，这里我们会保留源汇点直接以产量表示的方式。

重要

在专业油藏数值模拟软件中，对井的处理是非常重要的核心之一。可以说井处理的好坏很大程度上决定了算法的稳定性和计算的速度快慢。