一种超压驱动下原油二次运移速率的计算方法

摘要

本发明涉及地质勘探技术领域，公开了一种超压驱动下原油二次运移速率的计算方法，包括：通过连续性方程、达西定律、流体势表达式和地层流体物性关系，得到超压驱动流体流动基本方程，考虑毛细管压力和油水饱和度，将基本方程转化为油-水两相流动方程；通过有限差分显式格式法对流动方程进行离散化，得到显式表达式；对显式表达式进行求解得到孔压，基于孔压计算得到超压驱动下原油的二次运移速率。本发明结合超压储层中流体流动特点和物性特点，经地质建模、数学推导和求解，计算得到超压储层内超压驱动原油二次运移速率，可进一步得到超压驱动油气运移的距离和原油分布范围，为油气地质评价和勘探以及新的油气储量的发现提供科学依据。

说明书

技术领域

本发明涉及地质勘探技术领域，主要涉及一种超压驱动下原油二次运移速率的计算方法。

背景技术

地层压力是影响地下油气运移和聚集的一个极为重要的因素，异常高压(超压)与油气运聚有着密切的关系。目前，国内外在对异常高压与油气运聚关系的研究中，大多是单纯从流体势特征出发给出油气运移指向和定性解释。另外，也有个别学者通过浮力驱动建模的方法计算油气运聚效率。由于发育有超压的油气储层的特殊性和研究的复杂性，目前还没有涉及超压驱动下原油二次运移速率的计算方法。

发明内容

本发明提供了一种超压驱动下原油二次运移速率的计算方法，解决了或部分解决了现有技术中无法计算出超压驱动下原油二次运移速率的技术问题。

为解决上述技术问题，本发明提供了一种超压驱动下原油二次运移速率的计算方法，包括：

通过连续性方程、达西定律、流体势表达式和地层流体物性关系，得到超压驱动流体流动基本方程$>▿·\left[{\rho }_l\frac{{\mathrm{kk}}_{\mathrm{rl}}}{{\mu }_l}({▿P}_l-{\rho }_{l}g▿D)\right]+w_l=\frac{\partial }{\partial t}\left({\rho }_l\phi S_l\right);$>其中，▽为关于空间矢量微分算子；l为o或w，o表示油相、w表示水相；φ为岩石的有效孔隙度；P_l为l相流体的压力；S_l为l相流体的饱和度；k_rl为l相流体的相对渗透率；k为岩石的绝对渗透率；ρ_l为l相流体的密度；μ_l为l相流体的粘度；D为距离基准面以下的距离，取正值；w_l为l相流体源汇项，源取正号、汇取负号；

考虑毛细管力和油水饱和度，将所基本方程转化为超压驱动下油-水两相流动方程：$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial y}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{\partial P_o}{\partial t}\end{array}\right);$>

其中，k_x为岩石在x方向的绝对渗透率，k_y为岩石在y方向的绝对渗透率；φ为岩石的有效孔隙度；P_o为地层原油的压力，P_w为地层水的压力；S_o为地层原油的饱和度，S_w为地层水的饱和度；k_ro为地层原油的相对渗透率，k_rw为地层水的相对渗透率；ρ_o为地层原油的密度，ρ_w为地层水的密度；μ_o为地层原油的动力粘度，μ_w为地层水的动力粘度；

通过有限差分显式格式法对所述油-水两相流动方程进行离散化，得到显式表达式；对所述显式表达式进行求解得到孔压，基于所述孔压计算得到超压驱动下原油的二次运移速率。

进一步地，所述对所述显式表达式进行求解得到孔压，包括：通过高斯-赛德尔迭代法对所述显式表达式进行求解得到所述孔压。

进一步地，所述考虑毛细管力和油水饱和度，将所述基本方程转化为油-水两相流动方程，包括：

将所述超压驱动流体流动基本方程改写为油的流动方程：

$>▿·\left[{\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}({▿P}_o-{\rho }_{o}g▿D)\right]=\frac{\partial }{\partial t}\left({\rho }_o\phi S_o\right);$>

设流体密度不随空间变化，将所述油的流动方程的等号左边化为：

$>▿·\left[{\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}({▿P}_o-{\rho }_{o}g▿D)\right]={\rho }_o▿·\left[\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}(▿P_o-{\rho }_{o}g▿D)\right];$>

将所述油的流动方程的等号右边展开，即

$>\left(\begin{array}{c}\frac{\partial }{\partial t}\left({\rho }_o\phi S_o\right)=\phi S_o\frac{\partial {\rho }_o}{\partial t}+{\rho }_o{S}_o\frac{\partial \phi }{\partial t}+{\rho }_o\phi \frac{\partial S_o}{\partial t}\\ ={\rho }_o\phi S_o\frac{1}{{\rho }_o}\frac{\partial {\rho }_o}{\partial P_o}\frac{\partial P_o}{\partial t}+{\rho }_o\phi S_o\frac{1}{\phi }\frac{\partial \phi }{\partial P_{\phi }}\frac{{\partial P}_{\phi }}{\partial t}+{\rho }_o\phi \frac{{\partial S}_o}{\partial t}\\ ={\rho }_o\phi S_o{C}_o\frac{\partial P_O}{\partial t}+{\rho }_o\phi S_o{C}_{\phi }\frac{{\partial P}_{\phi }}{\partial t}+{\rho }_o\phi \frac{{\partial S}_o}{\partial t}\end{array}\right);$>其中，$>C_o=\frac{1}{{\rho }_o}\frac{{\partial \rho }_o}{{\partial P}_p},C_{\phi }=\frac{1}{\phi }\frac{\partial \phi }{{\partial P}_{\phi }};$>考虑毛细管力，可知$>P_{\phi }=P_o-\frac{1}{2}{P}_{\mathrm{cow}}=P_w+\frac{1}{2}{P}_{\mathrm{cow}};$>

在一个模拟时间步长Δt内，饱和度是不变的，故设P_cow在Δt时间内不随时间变化，因而于是将所述油的流动方程的等号右边进一步转化为：$>\frac{\partial }{\partial t}\left({\rho }_o\phi S_o\right)={\rho }_o\phi S_o(C_o+C_{\phi })\frac{{\partial P}_o}{\partial t}+{\rho }_o\phi \frac{{\partial S}_o}{\partial t};$>

所述油的流动方程的等号左边等于等号右边，故所述油的流动方程转化为简化后：$>▿·\left[\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}({▿P}_o-{\rho }_{o}g▿D)\right]=\phi S_o(C_o+C_{\phi })\frac{\partial P_o}{\partial t}+\phi \frac{{\partial S}_o}{\partial t};$>

将所述超压驱动流体流动基本方程改写为水的流动方程

$>▿·\left[{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}(▿P_w-{\rho }_{w}g▿D)\right]=\frac{\partial }{\partial t}\left({\rho }_w\phi S_w\right);$>

基于所述油的流动方程的转化过程，将所述水的流动方程转化为：

$>▿·\left[\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}({▿P}_w-{\rho }_{w}g▿D)\right]=\phi S_w(C_w+C_{\phi })\frac{\partial P_w}{\partial t}+\phi \frac{{\partial S}_w}{\partial t};$>

其中，$>C_w=\frac{1}{{\rho }_w}\frac{{\partial \rho }_w}{{\partial P}_w};$>

考虑毛细管力，结合将所述转化后的水的流动方程进一步转化为：$>\left(\begin{array}{c}▿·\left[\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}({▿P}_o-{▿P}_{\mathrm{cow}}-{\rho }_{w}g▿D)\right]=\phi S_w(C_w+C_{\phi })(\frac{{\partial P}_o}{\partial t}-\frac{{\partial P}_{\mathrm{cow}}}{\partial t})+\phi \frac{{\partial S}_w}{\partial t}\\ =\phi S_w(C_w+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_w}{\partial t}\end{array}\right)$>

将所述转化后的油的流动方程加上所述进一步转化的水的流动方程，得到第一转化的油-水两相运动方程：

$>\left(\begin{array}{c}▿·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g▿D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿P}_{\mathrm{cow}}]\\ ={\mathrm{\phi S}}_o(C_o+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_o}{\partial t}+\phi S_w(C_w+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_w}{\partial t}\\ =\phi [S_o(C_o+C_{\phi })+S_w(C_w+C_{\phi })]\frac{{\partial P}_o}{\partial t}+\phi (\frac{{\partial S}_o}{\partial t}+\frac{{\partial S}_w}{\partial t})\end{array}\right);$>

考虑油水饱和度，可知代入所述第一转化的油-水两相运动方程化简得到第二转化的油-水两相运动方程：

$>\left(\begin{array}{c}▿·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g▿D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿P}_{\mathrm{cow}}]\\ ={\mathrm{\phi S}}_o(C_o+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_o}{\partial t}+\phi S_w(C_w+C_{\phi })\frac{{\partial P}_o}{\partial t}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{{\partial P}_o}{\partial t}\end{array}\right);$>

当空间微分算子▽取x方向时，将所述第二转化的油-水两相运动方程的等号左边化为：

$>\left(\begin{array}{c}{▿}_x·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿}_x{P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g{▿}_{x}D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿}_x{P}_{\mathrm{cow}}]\\ =(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\end{array}\right);$>

当空间微分算子▽取y方向时，将所述第二转化的油-水两相运动方程的等号左边化为：

$>\left(\begin{array}{c}{▿}_y·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿}_y{P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g{▿}_{y}D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿}_y{P}_{\mathrm{cow}}]\\ =(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\end{array}\right);$>

故将所述第二转化的油-水两相流动方程化简为最终的油-水两相流动方程：$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial y}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{\partial P_o}{\partial t}\end{array}\right).$>

进一步地，C_o＝β_o、C_w＝β_w和C_φ＝(1-φ)C_P；

其中，β_o原油体积压缩系数：地层中单位流体压降的原油体积变化率，表达式为：

$>{\beta }_o=\frac{1}{V_o}\frac{{\partial V}_o}{\partial P};$>

其中，V_o为单位质量地层原油的体积；

β_w为地层水体积压缩系数：地层中单位流体压降的水体积变化率，表达式为：

$>{\beta }_w=\frac{1}{V_w}\frac{{\partial V}_w}{\partial P};$>

其中，V_w为单位质量地层水的体积；

C_P为孔隙体积压缩系数：地层中单位流体压降的孔隙体积变化率，表达式为：

$>C_P=\frac{1}{V_{\phi }}\frac{{\partial V}_{\phi }}{\partial P};$>

其中，V_φ为一定体积的地层岩石中的孔隙体积。

进一步地，所述通过有限差分显式格式法对所述最终的油-水两相流动方程进行离散化，得到显式表达式，包括：

设时间步长：Δt＝t(n+1)-t(n)，空间步长：Δx＝x(i+1)-x(i)、Δy＝y(j+1)-y(j)，可知P_oⁿ(i,j)为网格点(i,j)处的初始孔压经过n个Δt时间后的孔压，P_oⁿ⁺¹(i,j)为网格点(i,j)处的孔压P_oⁿ(i,j)经过Δt时间后的孔压；P_cow(i,j)为网格点(i,j)处的毛细管压力；D(i,j)为目的层在网格点(i,j)处的深度；k_x(i,j)为网格点(i,j)处岩石在x方向的绝对渗透率；k_y(i,j)为网格点(i,j)处岩石在y方向的绝对渗透率；

结合微分方程、泰勒公式和有限差分解法原理，用差商代替微商，得到：

$>\frac{\partial {P_o}^n(i,j)}{\partial t}=\frac{{P_o}^{n+1}(i,j)-{P_o}^n(i,j)}{\mathrm{\Delta t}};$>

$>\frac{{\partial }^2{P_o}^n(i,j)}{{\partial x}^2}=\frac{{P_o}^n(i-1,j)-2{P_o}^n(i,j)+{P_o}^n(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2};$>

$>\frac{{\partial }^2{P_o}^n(i,j)}{{\partial y}^2}=\frac{{P_o}^n(i-1,j)-2{P_o}^n(i,j)+{P_o}^n(i+1,j)}{{\left(\mathrm{\Delta y}\right)}^2};$>

$>\frac{{\partial }^2{P}_{\mathrm{cow}}(i,j)}{{\partial x}^2}=\frac{P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2};$>

$>\frac{{\partial }^{2}D(i,j)}{{\partial x}^2}=\frac{D(i-1,j)-2D(i,j)+D(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2};$>

其中，为第n个时间步长时网格点(i,j)处孔压P_oⁿ(i,j)对时间的一阶偏导数；为第n个时间步长时网格点(i,j)处的孔压P_oⁿ(i,j)对x方向的二阶偏导数；为第n个时间步长时网格点(i,j)处的孔压P_oⁿ(i,j)对y方向的二阶偏导数；为网格点(i,j)处毛细管压力对x方向的二阶偏导数；为网格点(i,j)处毛细管压力对y方向的二阶偏导数；为网格点(i,j)处高程对x方向的二阶偏导数；为网格点(i,j)处高程对y方向的二阶偏导数；

将和代入所述最终的油-水两相运动方程$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial y}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{\partial P_o}{\partial t}\end{array}\right)$>中，得到再次转化的油-水两相流动方程：

$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{P_o}^n(i-1,j)-2{P_o}^n(i,j)+{P_o}^n(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}\\ -({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{D(i-1,j)-2D(i,j)+D(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}\\ -\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{P_o}^n(i,j-1)-2{P_o}^n(i,j)+{P_o}^n(i,j+1)}{{\left(\mathrm{\Delta y}\right)}^2}\\ -({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{D(i,j-1)-2D(i,j)+D(i,j+1)}{{\left(\mathrm{\Delta y}\right)}^2}\\ -\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{P_{\mathrm{cow}}(i,j-1)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i,j+1)}{{\left(\mathrm{\Delta y}\right)}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{{P_o}^{n+1}(i,j)-{P_o}^n(i,j)}{\mathrm{\Delta t}}\end{array}\right);$>

引入两个大于0的无量纲变量：

$>{\lambda }_x(i,j)=\frac{(\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w})·\mathrm{\Delta t}}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2};$>

$>{\lambda }_y(i,j)=\frac{(\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w})·\mathrm{\Delta t}}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2};$>

并且设置以下四个中间参数：

$>{\alpha }_x(i,j)=\frac{({\rho }_o\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w})g}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2};$>

$>{\alpha }_y(i,j)=\frac{({\rho }_o\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w})g}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2};$>

$>{\beta }_x(i,j)=\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2};$>

$>{\beta }_y(i,j)=\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2};$>

其中，λ_x(i,j)为地层流体在x方向的传导系数，λ_y(i,j)为地层流体在y方向的传导系数；

将所述两个无量纲变量λ_x(i,j)和λ_y(i,j)和四个中间参数α_x(i,j)、α_y(i,j)、β_x(i,j)、β_y(i,j)代入所述再次转简化的油-水两相流动方程中，得到第(i,j)个网格点在t(n+1)＝t(n)+Δt时孔压P_oⁿ⁺¹(i,j)关于上一时刻孔压P_oⁿ的显式表达式：$>\left(\begin{array}{c}{P_o}^{n+1}(i,j)=[1-2{\lambda }_x(i,j)-2{\lambda }_y(i,j)]{P_o}^n(i,j)\\ +{\lambda }_x(i-1,j)·{P_o}^n(i-1,j)+{\lambda }_x(i+1,j)·{P_o}^n(i+1,j)\\ +{\lambda }_y(i,j-1)·{P_o}^n(i,j-1)+{\lambda }_y(i,j+1)·{P_o}^n(i,j+1)\\ -{\alpha }_x(i,j)·[D(i-1,j)-2D(i,j)+D(i+1,j)]·\mathrm{\Delta t}\\ -{\alpha }_y(i,j)·[D(i-1,j)-2D(i,j)+D(i,j+1)]·\mathrm{\Delta t}\\ -{\beta }_x(i,j)·[P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i+1,j)]·\mathrm{\Delta t}\\ -{\beta }_y(i,j)·[P_{\mathrm{cow}}(i,j-1)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i,j+1)]·\mathrm{\Delta t}\end{array}\right).$>

进一步地，所述通过高斯-赛德尔迭代法对所述显式表达式进行求解得到所述孔压，包括：

通过所述高斯-赛德尔迭代法对所述显式表达式进行迭代收敛；具体包括：

在第n+1个时间步长(n+1)Δt时，在计算第(i,j)个网格点在第k+1次迭代后的孔压P_o^n+1,k+1(i,j)时，将所述显式表达式中第(i-1,j)和第(i,j-1)个网格点的压力值P_o^n,k(i-1,j)和P_o^n,k(i,j-1)改为当前次迭代已经产生的新值P_o^n,k+1(i-1,j)和P_o^n,k+1(i,j-1)，故在每一个时间步长Δt内，将所述显式表达式改写为最终显式差分格式表达式：

$>\left(\begin{array}{c}{P_o}^{n,k+1}(i,j)=[1-2{\lambda }_x(i,j)-2{\lambda }_y(i,j)]{P_o}^{n,k}(i,j)\\ +{\lambda }_x(i-1,j)·{P_o}^{n,k+1}(i-1,j)+{\lambda }_x(i+1,j)·{P_o}^{n,k}(i+1,j)\\ +{\lambda }_y(i,j-1)·{P_o}^{n,k+1}(i,j-1)+{\lambda }_y(i,j+1)·{P_o}^{n,k}(i,j+1)\\ -{\alpha }_x(i,j)·[D(i-1,j)-2D(i,j)+D(i+1,j)]·\mathrm{\Delta t}\\ -{\alpha }_y(i,j)·[D(i,j-1)-2D(i,j)+D(i,j+1)]·\mathrm{\Delta t}\\ -{\beta }_x(i,j)·[P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i+1,j)]·\mathrm{\Delta t}\\ -{\beta }_y(i,j)·[P_{\mathrm{cow}}(i,j-1)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i,j+1)]·\mathrm{\Delta t}\end{array}\right);$>

其中，P_o^n+1,k+1(i,j)为第(i,j)个网格点在第n到第n+1个时间步长内、第k+1次迭代后的孔压；

对所述最终显式差分格式表达式进行求解得到所述孔压。

进一步地，所述对所述最终显式差分格式表达式进行求解得到所述孔压，包括：

步骤(1)：获得初始参数，设定初始条件和边界条件，设置空间步长、时间步长和时间步长数；

步骤(2)：把t(n)时刻的孔压赋给指标，并将时间t(n)增加一个时间步长，进入下一个时间步长；

步骤(3)：计算2个无量纲的传导系数和四个中间参数；

步骤(4)：计算Δt时间后所有网格点第k+1次迭代后的孔压P_o^n+1,k+1(i,j)；

其中，在计算第n+1个时间步长时、第(i,j)个网格点在第k+1次迭代后的孔压P_o^n+1,k+1(i,j)时，将所述显式表达式中第(i-1,j)和第(i,j-1)个网格点的压力值P_o^n,k(i-1,j)和P_o^n,k(i,j-1)改为当前次迭代已经产生的新值P_o^n,k+1(i-1,j)和P_o^n,k+1(i,j-1)；因为在计算第(i,j)个网格点的孔压时，它前面的网格点(第(i-1,j)和第(i,j-1)个网格点)第k+1次迭代后的孔压P_o^n,k+1(i-1,j)和P_o^n,k+1(i,j-1)已经求出；

步骤(5)：判断所有网格点在该次迭代后得到的孔压P_o^n+1,k+1(i,j)与指标oldP_oⁿ(i,j)之间的相对误差的最大值是否在设定精度eps之内，即其是否小于eps；

步骤(6)：若否，将其赋值给指标，同时将迭代次数k加1，进入下一次迭代，然后转到步骤(4)，依次执行步骤(4)和步骤(5)；

步骤(7)：若是，则将其作为t(n+1)时刻的孔压P_oⁿ⁺¹(i,j)，并判断此时的时间步长数n+1是否大于或等于设定的时间步长数N，即判断n+1≥N是否成立；

步骤(8)：如果否，说明当前时间并未达到成藏期时间，并将t(n+1)时刻的孔压P_oⁿ⁺¹(i,j)覆盖上一个时间步长的孔压，然后转到步骤(2)，依次执行步骤(2)、(3)、(4)、(5)、(6)和步骤(7)；

步骤(9)：如果是，说明当前时间达到成藏期时间，则t(n+1)时刻的孔压P_oⁿ⁺¹(i,j)为所求的孔压。

进一步地，所述基于所述孔压计算得到超压驱动下原油的二次运移速率，包括：

基于所述孔压，结合达西定律和流体势方程，计算得到所述超压驱动下原油的二次运移速率。

进一步地，所述基于所述孔压，结合达西定律和流体势方程，计算得到所述超压驱动下原油的二次运移速率，包括：

基于所述孔压，结合达西定律和流体势方程，结合中心差商法得到超压驱动下原油的二次运移速率，得到

$>{v_{\mathrm{ox}}}^{n+1}(i,j)=-8.64\times {10}^{-8}{\rho }_o\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}\frac{{{\Phi }_o}^{n+1}(i+1,j)-{{\Phi }_o}^{n+1}(i-1,j)}{2\mathrm{\Delta x}}$>

$>{v_{\mathrm{oy}}}^{n+1}(i,j)=-8.64\times {10}^{-8}{\rho }_o\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}\frac{{{\Phi }_o}^{n+1}(i,j+1)-{{\Phi }_o}^{n+1}(i,j-1)}{2\mathrm{\Delta y}}$>

$>{v_{\mathrm{oxy}}}^{n+1}(i,j)=3.65\times {10}^5\sqrt{{\left[{v_{\mathrm{ox}}}^{n+1}(i,j)\right]}^2+{\left[{v_{\mathrm{oy}}}^{n+1}(i,j)\right]}^2}$>

其中，为第(i,j)个网格点在第n+1个时间步长时的油势；V_oxⁿ⁺¹(i,j)为第(i,j)个网格点在第n+1个时间步长时在x方向的运移速度；V_oyⁿ⁺¹(i,j)为第(i,j)个网格点在第n+1个时间步长时在y方向的运移速度；V_oxyⁿ⁺¹(i,j)为第(i,j)个网格点在第n+1个时间步长时在x、y两个方向叠加后的运移速率，即为所求的超压驱动原油二次运移的速率。

进一步地，所述显式表达式成立所要满足的条件为：λ_x(i,j)+λ_y(i,j)＜1/2且λ_x(i,j)和λ_y(i,j)均大于零。

本发明的有益效果在于：

本发明提供的超压驱动下原油二次运移速率的计算方法，结合超压储层中流体流动特点和物性特点，经地质建模、数学推导和求解，计算得出超压储层内超压释放过程中驱动原油二次运移速率。其中，本发明考虑了超压驱动下油气二次运移的达西流的动力学过程和油-水两相的毛细管阻力和重力作用，从而实现了超压的动态释放和原油二次运移速率的计算，并针对渤海湾盆地沾化凹陷大尺度的超压储层进行了石油运移速率的实际应用计算，从而可以进一步得到超压驱动油气运移的距离和原油分布范围，可为油气地质评价和勘探以及新的油气储量的发现提供科学依据。

附图说明

图1为本发明实施例提供的超压驱动下原油二次运移速率的计算方法的流程图；

图2为本发明实施例提供的超压驱动下原油二次运移速率的计算方法中利用高斯--赛德尔迭代法对计算模型进行求解的流程图；

图3为通过本发明实施例计算得到的渤海湾盆地沾化凹陷渤南洼陷Es3x亚段超压驱动原油二次运移速率的等值线图。

具体实施方式

为进一步阐述本发明为达成预定发明目的所采取的技术手段及功效，以下结合附图及较佳实施例，对依据本发明提出的超压驱动下原油二次运移速率的计算方法的具体实施方式及工作原理进行详细说明。

在对本发明实施例提供的方法进行说明之前，需要先对超压储层(多孔介质)中原油二次运移的驱动机制进行说明：

促使油气二次运移的动力和阻力有：浮力、水动力、重力和毛细管压力。根据已有研究，在超压层中，从烃源岩排出的烃类进入饱含水的储集层中，为油-水两相条件下，石油在浮力和水动力的驱使下克服毛细管阻力而发生油气的二次运移。在超压作用下，储集层中油水以分散状快速运移，超压储层中的原油很难达到克服毛细管阻力所需的连片油柱长度。因此，在超压储层内，原油从深洼区的超压带向四周做二次运移时，浮力的驱使动力作用很小，相对超压流体的动力作用可以忽略不计。

因此，在关键成藏期超压储层中的油气二次运移主要为超压驱动下的达西流的动力学的驱动机制。油水在受重力作用的同时，油-水两相将凹陷(或洼陷)区的异常高压作为运移动力、毛细管压力作为阻力进行运移，超压流体动力、储层物性和毛细管力控制着油气二次运移的方向和距离。

根据以上对多孔介质超压系统中油气二次运移的驱动机制的分析，目的层位超压被上下致密岩层封闭，从平面分布特征来看，原油的运移以层内侧向运移为主，将地层模型简化为多孔介质，利用达西渗流方程计算超压环境下原油的二次运移速率。在计算过程中，对原油二次运移时的地层模型和相态特征做了如下假设条件和简化：

1)、原油从超压烃源岩排进超压储层后，超压不断积累，而且会随时驱动石油运移。石油运移的时机应与关键成藏期保持一致，而且触发条件如构造作用、超压驱动力大于流体运移的阻力等。在这里，将超压作为石油运移的驱动力，将毛细管压力作为运移阻力，也考虑重力的作用；

2)、在超压作用下，原油在多孔介质中的二次运移符合达西渗流条件；

3)、在超压储层内，油以分散的形式与地层水进行油-水两相沿地层侧向运移；

4)、超压储层中多孔岩石的骨架不可压缩，地层水的密度不变，原油可与温度的变化建立关系；

5)、模拟期为超压油气系统中的某一关键成藏期，该时期内的超压分布与储层物性特征空间分布一定。

参见图1，本发明实施例提供的超压驱动下原油二次运移速率的计算方法，包括：

通过连续性方程、达西定律、流体势表达式和地层流体物性关系，得到超压驱动流体流动基本方程，如式(1)：

$>▿·\left[{\rho }_l\frac{{\mathrm{kk}}_{\mathrm{rl}}}{{\mu }_l}({▿P}_l-{\rho }_{l}g▿D)\right]+w_l=\frac{\partial }{\partial t}\left({\rho }_l\phi S_l\right)---\left(1\right)$>

其中，▽为关于空间矢量微分算子；l为o或w，o表示油相、w表示水相；φ为岩石的有效孔隙度；P_l为l相流体的压力；S_l为l相流体的饱和度；k_rl为l相流体的相对渗透率；k为岩石的绝对渗透率；ρ_l为l相流体的密度；μ_l为l相流体的粘度；D为距离基准面以下的距离，取正值；w_l为l相流体源汇项，源取正号、汇取负号，在储层中不考虑生烃作用，故w_l＝0；

其中，连续性方程的物理意义：当流体通过多孔介质时，在任一时间间隔内从该截面流出的质量等于向该截面流入的质量。表达式如下：

其中，为流体的质量率(源或汇)，源取正号、汇取负号。

流体势的物理意义：单位质量流体所具有的机械能之和，主要包括压能和位能(势能)，表达式如下：

$>\Phi ={\int }_{P_0}^P\frac{\mathrm{dP}}{\rho \left(P\right)}+\mathrm{gz}$>

其中，Φ为流体势；P为任意位置z处的压力；P₀为基准点(z＝0)处的压力；z为距离基准面以下的距离，z与D的数值相同，取号相反。

由于储层条件下石油流体的可压缩性较小，若将密度视为定值，则流体势表达式的压能项中的积分号可取消，故可将流体势表达式简化为：

$>\Phi =\frac{P}{\rho }+\mathrm{gz}---\left(3\right)$>

达西定律的物理意义：流体的流动速度与其渗透率、密度、流体势微分成正比，与其动力粘度成反比。如下式所示：

$>v=-\frac{k}{\mu }\rho ·▿\Phi ---\left(4\right)$>

其中，ν为流体流动速度；

地层流体物性关系的物理意义：流体的密度随压力的变化而变化。其表达式如下：

ρ＝ρ(P)；

其中，ρ(P)为流体密度关于压力的函数；

考虑毛细管力和油水饱和度，将超压驱动流体流动基本方程(式(1))转化为油-水两相流动方程：

$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial y}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{\partial P_o}{\partial t}\end{array}\right)$>

其中，k_x为岩石在x方向的绝对渗透率，k_y为岩石在y方向的绝对渗透率；φ为岩石的有效孔隙度；P_o为地层原油的压力，P_w为地层水的压力；S_o为地层原油的饱和度，S_w为地层水的饱和度；k_ro为地层原油的相对渗透率，k_rw为地层水的相对渗透率；ρ_o为地层原油的密度，ρ_w为地层水的密度；μ_o为地层原油的动力粘度，μ_w为地层水的动力粘度。

对本步骤进行说明，本步骤具体包括：

将超压驱动流体流动基本方程(式(1))改写为油和水的流动方程。其中，油的流动方程为式(5)，水的流动方程为式(6)：

$>▿·\left[{\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}(▿P_o-{\rho }_{o}g▿D)\right]=\frac{\partial }{\partial t}\left({\rho }_o\phi S_o\right)---\left(5\right)$>

$>▿·\left[{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}(▿P_w-{\rho }_{w}g▿D)\right]=\frac{\partial }{\partial t}\left({\rho }_w\phi S_w\right)---\left(6\right)$>

设流体密度不随空间变化，即同一储层内油的密度和水的密度均不变，将油的流动方程(式(5))的等号左边化为：

$>▿·\left[{\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}({▿P}_o-{\rho }_{o}g▿D)\right]={\rho }_o▿·\left[\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}(▿P_o-{\rho }_{o}g▿D)\right]$>

将油的流动方程(式(5))的等号右边展开，即

$>\left(\begin{array}{c}\frac{\partial }{\partial t}\left({\rho }_o\phi S_o\right)=\phi S_o\frac{\partial {\rho }_o}{\partial t}+{\rho }_o{S}_o\frac{\partial \phi }{\partial t}+{\rho }_o\phi \frac{\partial S_o}{\partial t}\\ ={\rho }_o\phi S_o\frac{1}{{\rho }_o}\frac{\partial {\rho }_o}{\partial P_o}\frac{\partial P_o}{\partial t}+{\rho }_o\phi S_o\frac{1}{\phi }\frac{\partial \phi }{\partial P_{\phi }}\frac{{\partial P}_{\phi }}{\partial t}+{\rho }_o\phi \frac{{\partial S}_o}{\partial t}\\ ={\rho }_o\phi S_o{C}_o\frac{\partial P_o}{\partial t}+{\rho }_o\phi S_o{C}_{\phi }\frac{{\partial P}_{\phi }}{\partial t}+{\rho }_o\phi \frac{{\partial S}_o}{\partial t}\end{array}\right)---\left(7\right)$>

其中，$>C_o=\frac{1}{{\rho }_o}\frac{{\partial \rho }_o}{{\partial P}_o},C_{\phi }=\frac{1}{\phi }\frac{\partial \phi }{{\partial P}_{\phi }};P_{\phi }=\frac{1}{2}(P_o+P_w),$>为孔隙压力(简称孔压)。考虑毛细管压力(式(8))，可知

P_w＝P_o-P_cow   (8)

其中，P_cow为油水两相间的毛细管力。

在一个模拟时间步长Δt内，饱和度是不变的，故设P_cow在Δt时间内不随时间变化，因而于是将油的流动方程的等号右边进一步转化为：

$>\frac{\partial }{\partial t}\left({\rho }_o\phi S_o\right)={\rho }_o\phi S_o(C_o+C_{\phi })\frac{{\partial P}_o}{\partial t}+{\rho }_o\phi \frac{{\partial S}_o}{\partial t}---\left(9\right)$>

油的流动方程的等号左边等于等号右边，故将油的流动方程转化为：

$>▿·\left[\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}({▿P}_o-{\rho }_{o}g▿D)\right]=\phi S_o(C_o+C_{\phi })\frac{\partial P_o}{\partial t}+\phi \frac{{\partial S}_o}{\partial t}---\left(10\right)$>

基于油的流动方程的转化过程，将水的流动方程(式(6))转化为：

$>▿·\left[\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}({▿P}_w-{\rho }_{w}g▿D)\right]=\phi S_w(C_w+C_{\phi })\frac{\partial P_w}{\partial t}+\phi \frac{{\partial S}_w}{\partial t}.$>

其中，$>C_w=\frac{1}{{\rho }_w}\frac{{\partial \rho }_w}{{\partial P}_w};$>

考虑毛细管力(式(8))，结合将转化后的水的流动方程进一步转化为：

$>\left(\begin{array}{c}▿·\left[\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}({▿P}_o-{▿P}_{\mathrm{cow}}-{\rho }_{w}g▿D)\right]=\phi S_w(C_w+C_{\phi })(\frac{{\partial P}_o}{\partial t}-\frac{{\partial P}_{\mathrm{cow}}}{\partial t})+\phi \frac{{\partial S}_w}{\partial t}\\ =\phi S_w(C_w+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_w}{\partial t}\end{array}\right)---\left(11\right)$>

将所述转化后的油的流动方程(式(10))加上进一步转化的水的流动方程(式(11))，得到第一转化的油-水两相运动方程：

$>\left(\begin{array}{c}▿·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g▿D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿P}_{\mathrm{cow}}]\\ ={\mathrm{\phi S}}_o(C_o+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_o}{\partial t}+\phi S_w(C_w+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_w}{\partial t}\\ =\phi [S_o(C_o+C_{\phi })+S_w(C_w+C_{\phi })]\frac{{\partial P}_o}{\partial t}+\phi (\frac{{\partial S}_o}{\partial t}+\frac{{\partial S}_w}{\partial t})\end{array}\right)---\left(12\right)$>

考虑油水饱和度S_o+S_w＝1，可知，代入第一转化的油-水两相流动方程(式(12))得到第二转化的油-水两相流动方程：

$>\left(\begin{array}{c}▿·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g▿D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿P}_{\mathrm{cow}}]\\ ={\mathrm{\phi S}}_o(C_o+C_{\phi })\frac{{\partial P}_o}{\partial t}+\phi \frac{{\partial S}_o}{\partial t}+\phi S_w(C_w+C_{\phi })\frac{{\partial P}_o}{\partial t}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{{\partial P}_o}{\partial t}\end{array}\right)---\left(13\right)$>

当空间微分算子▽取x方向时，将第二转化的油-水两相流动方程(式(13))的等号左边化为：

$>\left(\begin{array}{c}{▿}_x·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿}_x{P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g{▿}_{x}D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿}_x{P}_{\mathrm{cow}}]\\ =(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\end{array}\right)---\left(14\right)$>

当空间微分算子▽取y方向时，第二转化的油-水两相流动方程(式(13))等号左边化为：

$>\left(\begin{array}{c}{▿}_y·[(\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}){▿}_y{P}_o-({\rho }_o\frac{{\mathrm{kk}}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w})g{▿}_{y}D-\frac{{\mathrm{kk}}_{\mathrm{rw}}}{{\mu }_w}{▿}_y{P}_{\mathrm{cow}}]\\ =(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\end{array}\right)---\left(15\right)$>

故将第二转化的油-水两相运动方程转化为最终的油-水两相流动方程。具体地，将式(14)和式(15)相加，去掉空间微分算子符号，可得到最终的油-水两相流动方程(式(16))：

$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial x}^2}-({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial x}^2}-\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial x}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{\partial }^2{P}_o}{{\partial y}^2}-({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{{\partial }^{2}D}{{\partial y}^2}-\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{{\partial }^2{P}_{\mathrm{cow}}}{{\partial y}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{\partial P_o}{\partial t}\end{array}\right)---\left(16\right)$>

在此，对C_o、C_w和C_φ进行说明：

根据岩石物理学和油层物理学的相关知识可知：在地层条件下，地层岩石骨架一般不具有压缩性，且地层原油体积压缩系数、地层水体积压缩系数和储层孔隙体积压缩系数概念和表达式分别为：

原油体积压缩系数：地层中单位流体压降的原油体积变化率，表达式为；

$>{\beta }_o=\frac{1}{V_o}\frac{{\partial V}_o}{\partial P}$>

其中，V_o为单位质量地层原油的体积。

地层水体积压缩系数：地层中单位流体压降的水体积变化率，表达式为；

$>{\beta }_w=\frac{1}{V_w}\frac{{\partial V}_w}{\partial P}$>

其中，V_w为单位质量地层水的体积。

孔隙体积压缩系数：地层中单位流体压降的孔隙体积变化率，表达式为：

$>C_P=\frac{1}{V_{\phi }}\frac{{\partial V}_{\phi }}{\partial P}$>

其中，V_φ为一定体积的地层岩石中的孔隙体积。

经过基本的物理公式和数学微分变换，可知C_o＝β_o、C_w＝β_w和C_φ＝(1-φ)C_P。β_o、β_w和C_P可根据国家相关标准和石油与天然气行业相关标准在油田现场和室内通过测试得到。鉴于此，在本发明实施例中，就根据C_o＝β_o、C_w＝β_w和C_φ＝(1-φ)C_P进行参数替换，且C_o、β_o、C_w、β_w、C_φ和C_P的单位都为MPa^-1。

通过有限差分显式格式法对最终的转化方程进行离散化，得到显式表达式；对显式表达式进行求解得到孔压，基于孔压计算得到超压驱动下原油的二次运移速率。

其中，通过有限差分显式格式法对简化方程进行离散化，得到显式表达式，包括：

设时间步长：Δt＝t(n+1)-t(n)，空间步长：Δx＝x(i+1)-x(i)、Δy＝y(j+1)-y(j)，可知P_oⁿ(i,j)为网格点(i,j)处的初始孔压经过n个Δt时间后的孔压，P_oⁿ⁺¹(i,j)为网格点(i,j)处的孔压P_oⁿ(i,j)经过Δt时间后的孔压；P_cow(i,j)为网格点(i,j)处的毛细管压力；D(i,j)为目的层在网格点(i,j)处的深度(一般取海平面为基准面)；k_x(i,j)为网格点(i,j)处岩石在x方向的绝对渗透率；k_y(i,j)为网格点(i,j)处岩石在y方向的绝对渗透率。

结合微分方程、泰勒公式和有限差分解法原理，用差商代替微商，得到：

$>\frac{\partial {P_o}^n(i,j)}{\partial t}=\frac{{P_o}^{n+1}(i,j)-{P_o}^n(i,j)}{\mathrm{\Delta t}}---(17-a)$>

$>\frac{{\partial }^2{P_o}^n(i,j)}{{\partial x}^2}=\frac{{P_o}^n(i-1,j)-2{P_o}^n(i,j)+{P_o}^n(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}---(17-b)$>

$>\frac{{\partial }^2{P_o}^n(i,j)}{{\partial y}^2}=\frac{{P_o}^n(i-1,j)-2{P_o}^n(i,j)+{P_o}^n(i+1,j)}{{\left(\mathrm{\Delta y}\right)}^2}---(17-c)$>

$>\frac{{\partial }^2{P}_{\mathrm{cow}}(i,j)}{{\partial x}^2}=\frac{P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}---(17-d)$>

$>\frac{{\partial }^{2}D(i,j)}{{\partial x}^2}=\frac{D(i-1,j)-2D(i,j)+D(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}---(17-e)$>

其中，为第n个时间步长时网格点(i,j)处孔压P_oⁿ(i,j)对时间的一阶偏导数；为第n个时间步长时网格点(i,j)处的孔压P_oⁿ(i,j)对x方向的二阶偏导数；为第n个时间步长时网格点(i,j)处的孔压P_oⁿ(i,j)对y方向的二阶偏导数；为网格点(i,j)处毛细管压力对x方向的二阶偏导数；为网格点(i,j)处毛细管压力对y方向的二阶偏导数；为网格点(i,j)处高程对x方向的二阶偏导数；为网格点(i,j)处高程对y方向的二阶偏导数。

将式(17-a)到(17-e)变形代入最终的油-水两相运动方程(式(16))中，得到再次转化的油-水两相流动方程：

$>\left(\begin{array}{c}(\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{P_o}^n(i-1,j)-2{P_o}^n(i,j)+{P_o}^n(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}\\ -({\rho }_o\frac{k_x{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{D(i-1,j)-2D(i,j)+D(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}\\ -\frac{k_x{k}_{\mathrm{rw}}}{{\mu }_w}\frac{P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i+1,j)}{{\left(\mathrm{\Delta x}\right)}^2}\\ +(\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})\frac{{P_o}^n(i,j-1)-2{P_o}^n(i,j)+{P_o}^n(i,j+1)}{{\left(\mathrm{\Delta y}\right)}^2}\\ -({\rho }_o\frac{k_y{k}_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w})g\frac{D(i,j-1)-2D(i,j)+D(i,j+1)}{{\left(\mathrm{\Delta y}\right)}^2}\\ -\frac{k_y{k}_{\mathrm{rw}}}{{\mu }_w}\frac{P_{\mathrm{cow}}(i,j-1)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i,j+1)}{{\left(\mathrm{\Delta y}\right)}^2}\\ =\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]\frac{{P_o}^{n+1}(i,j)-{P_o}^n(i,j)}{\mathrm{\Delta t}}\end{array}\right)---\left(18\right)$>

在这里，引入两个大于0的无量纲变量：

$>{\lambda }_x(i,j)=\frac{(\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w})·\mathrm{\Delta t}}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2}---(19-a)$>

$>{\lambda }_y(i,j)=\frac{(\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w})·\mathrm{\Delta t}}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2}(19-b)$>

并且设置以下四个中间参数：

$>{\alpha }_x(i,j)=\frac{({\rho }_o\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w})g}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2}---(20-a)$>

$>{\alpha }_y(i,j)=\frac{({\rho }_o\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+{\rho }_w\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w})g}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2}---(20-b)$>

$>{\beta }_x(i,j)=\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2}---(20-c)$>

$>{\beta }_y(i,j)=\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2}---(20-d)$>

其中，λ_x(i,j)为地层流体在x方向的传导系数，λ_y(i,j)为地层流体在y方向的传导系数。

将两个无量纲变量λ_x(i,j)和λ_y(i,j)和四个中间参数α_x(i,j)、α_y(i,j)、β_x(i,j)、β_y(i,j)代入再次转化的油-水两相流动方程(式(18))中，得到第(i,j)个网格点在t(n+1)＝t(n)+Δt时(第n+1个时间步长后)孔压P_oⁿ⁺¹(i,j)关于上一时刻(第n个时间步长后)孔压P_oⁿ的显式表达式：

$>\left(\begin{array}{c}{P_o}^{n+1}(i,j)=[1-2{\lambda }_x(i,j)-2{\lambda }_y(i,j)]{P_o}^n(i,j)\\ +{\lambda }_x(i-1,j)·{P_o}^n(i-1,j)+{\lambda }_x(i+1,j)·{P_o}^n(i+1,j)\\ +{\lambda }_y(i,j-1)·{P_o}^n(i,j-1)+{\lambda }_y(i,j+1)·{P_o}^n(i,j+1)\\ -{\alpha }_x(i,j)·[D(i-1,j)-2D(i,j)+D(i+1,j)]·\mathrm{\Delta t}\\ -{\alpha }_y(i,j)·[D(i-1,j)-2D(i,j)+D(i,j+1)]·\mathrm{\Delta t}\\ -{\beta }_x(i,j)·[P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i+1,j)]·\mathrm{\Delta t}\\ -{\beta }_y(i,j)·[P_{\mathrm{cow}}(i,j-1)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i,j+1)]·\mathrm{\Delta t}\end{array}\right)---\left(21\right)$>

式(21)即为孔压显式差分格式的表达式。

对显式表达式进行求解得到孔压，包括：

通过高斯-赛德尔迭代法对显式表达式进行求解得到孔压。

具体地，为了加快迭代收敛到设定精度的速度，本发明实施例先通过高斯-赛德尔迭代法对显式表达式进行迭代收敛，即在每一个时间步长Δt内，在计算当前网格点在第k+1次迭代后的孔压时，采用前面网格点在第k+1次迭代已经求出的孔压新值，以达到加快收敛的目的。在这里，就是在第n+1个时间步长(n+1)Δt时，在计算第(i,j)个网格点在第k+1次迭代后的孔压P_o^n+1,k+1(i,j)时，将显式表达式(式(21))中第(i-1,j)和第(i,j-1)个网格点的压力值P_o^n,k(i-1,j)和P_o^n,k(i,j-1)改为当前次(第k+1次)迭代已经产生的新值P_o^n,k+1(i-1,j)和P_o^n,k+1(i,j-1)。

其中，根据差分方程稳定性的冯·牛曼分析法(傅里叶分析方法)，结合欧拉公式，可知显式差分格式成立所要满足的条件为：λ_x(i,j)+λ_y(i,j)＜1/2且λ_x(i,j)和λ_y(i,j)均大于零。

故依据高斯--赛德尔迭代法，在每一个时间步长Δt内，将显式表达式(式(21))改写为最终显式差分格式表达式(式(22))：

$>\left(\begin{array}{c}{P_o}^{n+1,k+1}(i,j)=[1-2{\lambda }_x(i,j)-2{\lambda }_y(i,j)]{P_o}^{n,k}(i,j)\\ +{\lambda }_x(i-1,j)·{P_o}^{n,k+1}(i-1,j)+{\lambda }_x(i+1,j)·{P_o}^{n,k}(i+1,j)\\ +{\lambda }_y(i,j-1)·{P_o}^{n,k+1}(i,j-1)+{\lambda }_y(i,j+1)·{P_o}^{n,k}(i,j+1)\\ -{\alpha }_x(i,j)·[D(i-1,j)-2D(i,j)+D(i+1,j)]·\mathrm{\Delta t}\\ -{\alpha }_y(i,j)·[D(i-1,j)-2D(i,j)+D(i,j+1)]·\mathrm{\Delta t}\\ -{\beta }_x(i,j)·[P_{\mathrm{cow}}(i-1,j)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i+1,j)]·\mathrm{\Delta t}\\ -{\beta }_y(i,j)·[P_{\mathrm{cow}}(i,j-1)-2P_{\mathrm{cow}}(i,j)+P_{\mathrm{cow}}(i,j+1)]·\mathrm{\Delta t}\end{array}\right)---\left(22\right)$>

其中，P_o^n+1,k+1(i,j)为第(i,j)个网格点在第n到第n+1个时间步长内、第k+1次迭代后的孔压，式(22)是本发明实施例最终使用的显式差分格式表达式。

参见图2，对最终显式差分格式表达式进行求解得到孔压，包括：

步骤(1)：“开始”后，获得初始参数，设定初始条件和边界条件，设置空间步长、时间步长和时间步长数；

其中，在本实施例中，初始参数包括：孔压、目的层位高程(深度)、储层物性和流体物性等参数；模型的边界条件和初始条件，设定如下：边界条件：常压边界(静水压力边界)，即需要研究的目的层位边界的油和水流体势均为零，Φ_o＝Φ_w＝0；其中，Φ_o为目的层位边界的油的流体势，Φ_w为目的层位边界的水的流体势；

初始条件：某一关键成藏期原油二次运移初始时刻第(i,j)个网格点在t(0)时的孔压P_o⁰(i,j)、水压P_w⁰(i,j)。

时间步长和空间步长的设定要使2个无量纲的传导系数λ_x(i,j)和λ_y(i,j)满足λ_x(i_,j)+λ_y(i,j)＜1/2且λ_x(i,j)和λ_y(i，j)均大于零；时间步长数根据地层原油成藏期时间设定，用N表示。

步骤(2)：把t(n)时刻的孔压赋给指标，即oldP_oⁿ(i,j)＝P_oⁿ(i,j)，并将时间t(n)增加一个时间步长，即t(n+1)＝t(n)+Δt，进入下一个时间步长；

步骤(3)：计算2个无量纲的传导系数和四个中间参数；

步骤(4)：计算Δt时间后所有网格点第k+1次迭代后的孔压P_o^n+1,k+1(i,j)；

其中，在计算第n+1个时间步长时、第(i,j)个网格点在第k+1次迭代后的孔压P_o^n+1,k+1(i,j)时，将显式表达式(式(21))中第(i-1,j)和第(i,j-1)个网格点的压力值P_o^n,k(i-1,j)和P_o^n,k(i,j-1)改为当前次(第k+1次)迭代已经产生的新值P_o^n,k+1(i-1,j)和P_o^n,k+1(i,j-1)；因为在计算第(i,j)个网格点的孔压时，它前面的网格点(第(i-1,j)和第(i,j-1)个网格点)第k+1次迭代后的孔压P_o^n,k+1(i-1,j)和P_o^n,k+1(i,j-1)已经求出；

步骤(5)：判断所有网格点在该次(第k+1次)迭代后得到的孔压P_o^n+1,k+1(i,j)与指标oldP_oⁿ(i,j)之间的相对误差的最大值是否在设定精度eps之内，即其是否小于eps；

步骤(6)：若否，将其赋值给指标，即oldP_oⁿ(i,j)＝P_o^n+1,k+1(i,j)，同时将迭代次数k加1，即k＝k+1，进入下一次(第(k+2))迭代，然后转到步骤(4)，依次执行步骤(4)和步骤(5)；

步骤(7)：若是，则将其作为t(n+1)时刻的孔压P_oⁿ⁺¹(i,j)，即P_oⁿ⁺¹(i,j)＝P_o^n+1,k+1(i,j)，并判断此时的时间步长数n+1是否大于或等于设定的时间步长数N，即判断n+1≥N是否成立；

步骤(8)：如果否，说明当前时间并未达到成藏期时间，并将t(n+1)时刻的孔压P_oⁿ⁺¹(i,j)覆盖上一个时间步长(t(n)时刻)的孔压，即P_oⁿ(i,j)＝P_oⁿ⁺¹(i,j)，然后转到步骤(2)，依次执行步骤(2)、(3)、(4)、(5)、(6)和步骤(7)；

步骤(9)：如果是，说明当前时间达到成藏期时间，则t(n+1)时刻的孔压P_oⁿ⁺¹(i,j)为所求的孔压，“结束”。

基于孔压计算得到超压驱动下原油的二次运移速率，包括：

基于孔压，结合达西定律和流体势方程，计算得到超压驱动下原油的二次运移速率。具体地，在通过以上步骤求解出第(i,j)个网格点在第n+1个时间步长时的孔压P_oⁿ⁺¹(i,j)后，根据流体势方程(式(3))，该时刻的油势Φ_oⁿ⁺¹(i,j)也被解出。采用达西定律(式(4))，结合中心差商法计算得到超压驱动下原油的二次运移速率(如式(23)、(24)和(25))：

$>{v_{\mathrm{ox}}}^{n+1}(i,j)=-8.64\times {10}^{-8}{\rho }_o\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}\frac{{{\Phi }_o}^{n+1}(i+1,j)-{{\Phi }_o}^{n+1}(i-1,j)}{2\mathrm{\Delta x}}---\left(23\right)$>

$>{v_{\mathrm{oy}}}^{n+1}(i,j)=-8.64\times {10}^{-8}{\rho }_o\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}\frac{{{\Phi }_o}^{n+1}(i,j+1)-{{\Phi }_o}^{n+1}(i,j-1)}{2\mathrm{\Delta y}}---\left(24\right)$>

$>{v_{\mathrm{oxy}}}^{n+1}(i,j)=3.65\times {10}^5\sqrt{{\left[{v_{\mathrm{ox}}}^{n+1}(i,j)\right]}^2+{\left[{v_{\mathrm{oy}}}^{n+1}(i,j)\right]}^2}---\left(25\right)$>

其中，为第(i,j)个网格点在第n+1个时间步长时的油势；P_oⁿ⁺¹(i,j)的单位为MPa；式(23)为第(i,j)个网格点在第n+1个时间步长时在x方向的运移速度，单位为m/d；式(24)为第(i,j)个网格点在第n+1个时间步长时在y方向的运移速度，单位为m/d；式(25)为第(i,j)个网格点在第n+1个时间步长时在x、y两个方向叠加后的运移速率，即为所求的超压驱动原油二次运移的速率，单位为千米/百万年，符号为km/Ma。

以渤海湾盆地沾化凹陷渤南洼陷超压储层为例

为检验本发明实施例提供的模型和方法的实用性，选取我国东部的渤海湾盆地沾化凹陷中的渤南洼陷作为研究区，将其内的超压储层--沙河街组第三段下亚段(Es3x)作为研究层位，使用该模型和方法，求取超压驱动下，Es3x亚段原油二次运移的速率。

1、研究区概况

渤南洼陷位于渤海湾盆地济阳坳陷沾化凹陷东部，是胜利油田沾化凹陷中地层埋藏深度最大的洼陷，面积约为600km²。渤海湾盆地是一个富油气新生代盆地，盆地内普遍发育超压。据已有的研究，渤南洼陷超压出现的深度范围为2300-41000m，超压主要发育在沙河街组沙三段(Es3)和沙四段(Es4)，尤其是沙三段下亚段(Es3x)发育强超压，最大压力系数可达1.8。且Es3段底部有致密的膏岩层封闭，其顶部也被致密的泥岩封隔，使得该段内的原油在垂向上运移几乎被遮挡。其中的原油在其超压的驱动下、克服毛细管力，沿着地层做侧向运移。

对渤南洼陷包裹体均一温度测定和充注时期分析认为，该洼陷有3期油气充注的时期：第一期36～24Ma，发生于沙一段沉积早期；第二期16～10Ma，发生于东营组沉积末～馆陶组沉积早期；第三期5～2Ma，发生于明化镇组沉积末期。但主要的充注期次是第二期和第三期。根据本次实例研究所取的数据和假设条件，在地质空间、时间尺度内，本次模拟计算的运移速率可应用于最晚一期，即Es3x亚段第三期油气充注时的原油运移速率。

2、网格剖分和基础参数的取值

渤南洼陷大地坐标范围为：x(20615000，20653000)、y(4178000，4210000)，取空间步长Δx＝380m、Δy＝320m，可知网格数为100*100＝10000，取时间步长Δt＝30天，为确定两个大于零的无量纲参数λ_x(i,j)和λ_y(i,j)、毛细管力P_cow(i,j)和四个中间参数α_x(i,j)、α_y(i,j)、β_x(i,j)、β_y(i,j)，对所需的相关参数做如下说明：

P_cow＝0.001(2γ_hc-w/r_t)

γ_ow＝38.379(ρ_w-ρ_o)^0.09935

log(r_t)＝0.459+0.5log(k)-0.385log(100φ)

其中：P_cow为毛细管力，单位为MPa；ρ_w为水的密度，单位为kg/m³；ρ_o为烃的密度，单位为kg/m³；γ_ow为油水界面张力，单位为dynes/cm；r_t为孔喉半径，10^-6m；k为储层的绝对渗透率，单位为mD；φ为储层孔隙度，小数。

通过对研究区实际资料进行研究，地层原油平均密度取795kg/m³，地层水密度取1015kg/m³；地层原油和地层水动力粘度分别取μ_o＝1.0mPa·s和μ_w＝0.15mPa·s；g为重力加速度常数，取9.81N/kg。

在考虑束缚水饱和度S_wi和含水饱和度S_w的条件下，根据公式

S_wi＝0.425-0.0862log(k)

$>k_{\mathrm{ro}}=0.567\exp \left(0.218\log \left(k\right)\right){\left[\frac{1-S_w}{1-S_{\mathrm{wi}}}\right]}^2$>

计算得到地层原油运移的相对渗透率k_ro；

由于可知当含油饱和度达到临界值时，才能渗流；且一般认为以含油饱和度等于50％是油层出现的下限条件。因此，这里取S_o＝S_w＝50％，k_rw＝0.20。其中，S_wi为地层孔隙束缚水饱和度，即地层岩石空隙中不能流动的水；k_ro为地层原油的相对渗透率；

结合以上各参数的单位，在本实例中，将无量纲的传导系数λ_x(i,j)和λ_y(i,j)的表达式化为：

$>{\lambda }_x(i,j)=8.64\times {10}^{-2}\frac{(\frac{k_x(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+\frac{k_x(i,j)·k_{\mathrm{rw}}}{{\mu }_w})·\mathrm{\Delta t}}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta x}\right)}^2}$>

$>{\lambda }_y(i,j)=8.64\times {10}^{-2}\frac{(\frac{k_y(i,j)·k_{\mathrm{ro}}}{{\mu }_o}+\frac{k_y(i,j)·k_{\mathrm{rw}}}{{\mu }_w})·\mathrm{\Delta t}}{\phi [S_o{C}_o+(1-S_o)C_w+C_{\phi }]{\left(\mathrm{\Delta y}\right)}^2};$>

通过分析研究区502口井的岩心可知，Es3x亚段孔隙度φ为φ∈(0.08,0.15)，小数，无量纲；绝对渗透率k为：k∈(0.5,12.0)，mD。通过对以上参数的分析，计算λ_x(i,j)和λ_y(i,j)的大小范围，可知其满足λ_x(i,j)+λ_y(i,j)＜1/2且λ_x(i,j)和λ_y(i,j)均大于0的模型收敛的条件。

将渤南洼陷502口井的1000余个地层试孔压力实测点数据、孔渗数据、经地震勘探得到的层位高程数据和相关流体参数代入计算模型，利用高斯--赛德尔迭代法，通过Visual Fortran计算机语言和Visual FortranComposer XE 2011编译器，编辑调用数据和核心迭代程序代码，求解出孔压P_oⁿ⁺¹(i,j)，最后计算出Es3x亚段的超压释放驱动原油二次运移的速率。

3、渤南洼陷Es3x亚段超压驱动原油二次运移速率计算结果

通过自行编写的Visual Fortran程序计算出Es3x亚段的二次运移速率ν_oxyⁿ⁺¹(i,j)后，将数据导入绘图软件(Golden Software Surfer 8.0)绘制Es3x亚段原油在晚期(第三期)油气充注时、在自身超压驱动下，石油二次运移速率的等值线图，如图3(渤南洼陷Es3x亚段超压驱动原油二次运移速率的等值线图，km/Ma)所示,单位km/Ma表示千米/百万年。

由图3可以看出：渤南洼陷Es3x亚段原油在晚期(第三期)油气充注时、超压驱动下在相对强超压区石油二次运移的速率最大可达300km/Ma、大多数范围处在40-80km/Ma。又可知，若在储层中原油运移速率较小的地方存在圈闭，则原油开始在此聚集成藏。因此，在超压的驱动下、克服毛细管压力，Es3x亚段原油第三期成藏时间大约只要1Ma(1百万年)。呈现出：超压驱动快速运移、在运移速率较小且运移方向上存在圈闭处可发生油气聚集成藏。

本发明实施例提供的超压驱动下原油二次运移速率的计算方法，结合超压储层中流体流动特点和物性特点，经地质建模、数学推导和求解，计算得出超压储层内超压释放过程中驱动原油二次运移速率。其中，本发明考虑了超压驱动下油气二次运移的达西流的动力学过程和油-水两相的毛细管阻力和重力作用，从而实现了超压的动态释放和原油二次运移速率的计算，并针对渤海湾盆地沾化凹陷大尺度的超压储层进行了石油运移速率的实际应用计算，从而可以进一步得到超压驱动油气运移的距离和原油分布范围，可为油气地质评价和勘探以及新的油气储量的发现提供科学依据。

最后所应说明的是，以上具体实施方式仅用以说明本发明的技术方案而非限制，尽管参照实例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。