CN102339326A - 一种分析模拟缝洞型油藏流体流动的方法 - Google Patents

Abstract

本发明提供了一种模拟缝洞型油藏流体流动的方法，属于油藏数值模拟领域。本方法提出了基于Navier-Stokes方程与Darcy方程耦合的油藏数值模拟方法，主要包括：(1)建立复杂介质Navier-Stokes与Darcy耦合数学模型；(2)Navier-Stokes与Darcy耦合模型的数值计算方法；(3)将岩体与裂缝分开造型，各自造型不存在困难，对于裂缝的描述也清晰容易，裂缝数量及空间产状不受限制，特别是不要求裂缝上的节点与基岩节点完全重合，使得三维工作容易很多。并编写了数值模拟程序来实现本方法，本方法发展了缝洞型油藏数值模拟理论和方法，科学地实现了对缝洞型油藏的模拟。

Description

一种分析模拟缝洞型油藏流体流动的方法

技术领域

本发明属于油气田开发领域，尤其针对缝洞型油藏，分析其时空下地层压力及油水分布的技术内容，具体涉及一种采用耦合型油藏数值模拟方法。

背景技术

目前国内外有关缝洞型油藏数值模拟研究才刚刚开始。美国Stanford大学针对裂缝型油藏开展了离散裂缝网络(DFN)模型研究，初步解决了小型计算能力的离散裂缝型油藏的数值模拟问题，但针对缝洞型油藏中的大型洞穴(长度、宽度、高度均大于0.5米)及大型裂缝(也叫特大缝，裂缝开度大于10毫米，延伸长度大于500米)没有考虑；美国A&M大学利用Stokes-Brinkman方程研究缝洞型油藏的网格粗化问题，开展了二维单相流数值模拟研究，三维及两相流的方法没有考虑、计算方法没有研究。有关缝洞型油藏动力学的理论研究，包括一般性数学模型和数值模拟方法技术还未研发和公开。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种耦合型油藏数值模拟方法，具体的是一种分析判断缝洞型油藏其时空下地层压力及油水分布的方法，实现对缝洞型油藏的科学开发，找准缝洞型油藏剩余油的分布位置，为科学合理地开发这类油田提供依据，达到提高采收率的目的。

本发明是通过以下技术方案实现的：

一种模拟缝洞型油藏流体流动的方法，所述分析方法包括如下过程，

A确定缝洞型油藏的基础参数，包括：测试油藏目标层的构造形态、缝洞体空间发育形态，充填溶洞内物理性质参数；和油藏内油、水流体性质参数；以及油藏内岩石压缩系数，生产井、注水井的钻井位置及完井方式，油藏的边界条件参数和地层初始参数；

B创建缝洞型油藏耦合数学模型，根据缝洞型油藏油水的流动特征分为三个区域：大型洞穴区、多孔介质区域和裂缝带；所述大型洞穴区包括大型洞穴和大裂缝；所述多孔介质区域包括基岩、溶孔和微裂缝；所述裂缝包括岩石受外力或内应力作用时丧失结合力产生破裂但没有产生位移的区域，主要为网状裂缝带；

根据上述划分的不同地质区域，建立复杂介质Navier-Stokes方程和Darcy方程耦合数学模型：用Navier-Stokes方程表示流体在洞穴区中流动的运动方程；用Darcy定律表示流体在多孔介质区中流动的运动方程；结合质量守恒定律、动量守恒定律和状态方程建立数学模型；根据界面处法向速度连续和法向方向力的平衡条件建立界面之间的耦合条件，再根据界面处的耦合条件确定界面处的耦合模型；

C初始化并设置参数步骤：用于依次设置非线性分析控制参数，初始化计算数据和时间参数，对时间参数更新和初始化迭代初值；

D依据所述步骤B的模型，求解洞穴区油水两相流动方程和多孔介质区油水两相渗流方程；其中，洞穴区油水两相流动方程采用算子分裂法CBS；多孔介质区域油水两相渗流方程采用解耦全隐式算法；对裂缝采用锚杆模型算法；

E判断结果步骤：

(1)根据界面条件判断迭代是否收敛，收敛则存储计算结果；不收敛则重复步骤D；

(2)判断计算时间是否结束，是则进入下面判断过程，否则进入步骤C；

(3)判断分析结果；

F根据缝洞型油藏地质模型，建立模拟前需要的三维定量数据体，导入油田生产井产液量、产油量及测压数据，进行流体流动预测与分析，显示待测缝洞型油藏区域油水分布结果。

所述步骤B中包括，

流体在洞穴区中流动的运动方程，如下：

油水两相在溶洞内流动过程中，设油和水为不可压缩流体，流体的密度不随压力发生变化，即ρ＝常数；

连续性方程：

$\▿\·\left(S_w{u}_w\right)=0---\left(1.1\right)$

$\▿\·\left(S_o{u}_o\right)=0---(1.2)$

运动方程：

$\frac{\∂}{\∂t}\left(S_w{\mathrm{\ρ}}_w{u}_w\right)=S_w{\mathrm{\ρ}}_{w}g-{\▿p}_w+{\mathrm{\μ}}_w{\▿}^2{u}_w+S_w{F}_w---\left(1.3\right)$

$\frac{\∂}{\∂t}\left(S_o{\mathrm{\ρ}}_o{u}_o\right)=S_o{\mathrm{\ρ}}_{o}g-{\▿p}_o+{\mathrm{\μ}}_o{\▿}^2{u}_o+S_o{F}_o---(1.4)$

忽略掉重力： $\frac{\∂}{\∂t}\left(S_w{\mathrm{\ρ}}_w{u}_w\right)=-{\▿p}_w+{\mathrm{\μ}}_w{\▿}^2{u}_w+S_w{F}_w---\left(1.3a\right)$

$\frac{\∂}{\∂t}\left(S_o{\mathrm{\ρ}}_o{u}_o\right)=-{\▿p}_o+{\mathrm{\μ}}_o{\▿}^2{u}_o+S_o{F}_o---(1.4a)$

油相和水相的体积分数有如下关系：

                          S_w+S_o＝1    (1.5)

在油水接触面上，压力相等：P_w＝P_o      (1.6)

式中：S_o为油占据的体积百分数；S_w为水占据的体积百分数；u_k为k相流速；方程中u_o、u_w、p_o、p_w、S_o和S_w是未知数。

所述步骤B中，流体在多孔介质区域中流动的运动方程如下：

设油和水为不可压缩流体，那么流体密度是常量，不用考虑流体的状态方程；

流体运动方程：

$u_{\mathrm{pi}}=-\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}({\▿p}_i-{\mathrm{\ρ}}_{i}F),i=o,w---\left(1.7\right)$

若忽略单位质量上的力，其形式为：

$u_{\mathrm{pi}}=-\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}{\▿p}_i,i=o,w---(1.8)$

根据质量守恒方程可推得连续性方程为：

$\▿\·u_{\mathrm{pi}}=0,i=o,w---\left(1.9\right)$

联立运动方程、连续性方程，得到不可压缩油水两相流体的渗流控制方程，即：

$\▿\·\left(\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}{\▿p}_i\right)=\frac{\∂}{\∂t}\left({\mathrm{\φS}}_i\right),i=o,w---\left(1.10\right)$

p_c＝p_o-p_w       (1.11)

S_w+S_o＝1        (1.12)

S_i是各相的体积百分数，K_ri是各相的相对渗透率，p_c是毛管力。

此方程组的使用条件是：(1)彼此不互溶和不起化学反应的油水两相同时流动；(2)岩石和流体均不可压缩，服从线性渗流规律。

所述步骤B中，界面之间的耦合条件为：

法向速度连续性条件：

在界面上满足力的平衡条件：

$({\mathrm{\μ}}_o\frac{\∂u}{\∂n}-P\·n){|}_{{\mathrm{\Σ}}_f}=P_{p,o}\·n{|}_{{\mathrm{\Σ}}_P}$

$({\mathrm{\μ}}_w\frac{\∂u}{\∂n}-P\·n){|}_{{\mathrm{\Σ}}_f}=P_{p,w}\·n{|}_{{\mathrm{\Σ}}_P}$

其中，u为法向流速；P为法向压力；n为法向方向；μ_o为油的粘度；μ_w为水的粘度。

所述步骤D中，算子分裂法CBS算法如下：

分离成对流项和扩散项(0＜θ₁＜θ₂≤1)：

对流项： $\frac{u^{{n+\mathrm{\θ}}_1}-u^n}{\mathrm{dt}}-{\mathrm{\μ\Δu}}^{n+{\mathrm{\θ}}_1}+{\▿p}^{n+{\mathrm{\θ}}_1}=0$

扩散项： $\frac{u^{{n+\mathrm{\θ}}_2}-u^{n+{\mathrm{\θ}}_1}}{\mathrm{dt}}-u^{n+{\mathrm{\θ}}_2}+{\▿u}^{n+{\mathrm{\θ}}_2}=0$

将对流项方程写成时间与空间的形式，为

其中φ为通用符号包括u和p，对φ做Taylor展开，得到：

$\mathrm{\φ}(t+\mathrm{dt},x)=\mathrm{\φ}(t,x)+\mathrm{dt}\frac{\∂\mathrm{\φ}(t,x)}{\∂t}+\frac{{\mathrm{dt}}^2}{2}\frac{{\∂}^2\mathrm{\φ}(t,x)}{{\∂t}^2}+O\left({\mathrm{dt}}^3\right)$

由

上式中

均用

替换，于是，上式改写成

$\mathrm{\φ}(t+\mathrm{dt},x)=\mathrm{\φ}(t,x)-\mathrm{dt}(u\·\▿\mathrm{\φ}(t,x))-\frac{{\mathrm{dt}}^2}{2}\frac{\∂}{\∂t}(u\·\▿\mathrm{\φ}(t,x))$

其中u为上一时间步的速度；

所述解耦全隐式算法是把压力和饱和度分开来顺序求解，每次只解一个未知变量，总体交替迭代，直至收敛。

所述步骤D中，锚杆模型算法的具体算法如下：

(1)基本方程

基岩系统渗流方程：

$\▿\left({\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_m\right)+q_m={\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_m}{\∂t}---\left(2.1\right)$

裂缝系统渗流方程：

$\▿\left({\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_f\right)+q_f={\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_f}{\∂t}---(2.2)$

式中，下标m表示与基岩系统有关的参数与变量，下标f表示与裂缝系统有关的参数与变量；

(2)有限元弱形式

裂缝系统的弱形式如下：

${\∫}_F{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_f}{\∂t}{\mathrm{\δP}}_f\mathrm{dF}+{\∫}_F{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_f\▿\mathrm{\δ}{P}_f\mathrm{dF}-{\∫}_F\tilde{k}(P_m-P_f){\mathrm{\δP}}_f\mathrm{dF}$

$-{\∫}_{{\mathrm{\Γ}}_f}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_m-P_f){\mathrm{\δP}}_{f}d{\mathrm{\Γ}}_f={\∫}_F{q}_f{\mathrm{\δP}}_f\mathrm{dF}---\left(2.3\right)$

将裂缝边界的流入流出量添加到基岩系统的方程中，得到最终基岩系统的弱形式如下：

${\∫}_{\mathrm{\Ω}}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_m}{\∂t}{\mathrm{\δP}}_m\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_m\▿\mathrm{\δ}{P}_m\mathrm{d\Ω}-{\∫}_F\tilde{k}(P_m-P_f){\mathrm{\δP}}_m\mathrm{dF}$

${-\∫}_{{\mathrm{\Γ}}_m}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_m}{\∂n}{\mathrm{\δP}}_{m}d{\mathrm{\Γ}}_m+{\∫}_{{\mathrm{\Γ}}_f}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_m-P_f){\mathrm{\δP}}_{m}d{\mathrm{\Γ}}_f={\∫}_{\mathrm{\Ω}}{q}_m{\mathrm{\δP}}_m\mathrm{d\Ω}---\left(2.4\right)$

(3)有限元解及基岩与裂缝交换项的处理

令mh表示基岩区域Ω上的网格剖分，fh表示裂缝区域上F的网格剖分，由于基岩内为二维流动，裂缝内为一维流动，此时mh表示一个二维的网格剖分，fh表示一个一维的网格剖分；

令P_mh表示基岩系统的有限元解，P_fh表示裂缝系统的有限元解，Ω_h为剖分后的基岩区域，F_h为剖分后的裂缝区域，Γ_mh为剖分后基岩区域的边界，Γ_fh为剖分后的裂缝区域的边界，寻找P_mh及P_fh使得其满足方程：

$\begin{array}{c}{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_{\mathrm{mh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_{\mathrm{mh}}\▿\mathrm{\δ}{P}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{mh}}{\mathrm{dF}}_h\\ -{\∫}_{{\mathrm{\Γ}}_{\mathrm{mh}}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_{\mathrm{mh}}}{\∂n}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Γ}}_{\mathrm{mh}}-{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Γ}}_{\mathrm{fh}}={\∫}_{{\mathrm{\Ω}}_h}{q}_m{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h\end{array}---\left(2.5\right)$

$\begin{array}{c}{\∫}_{F_h}{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_{\mathrm{fh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h+{\∫}_{F_h}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_{\mathrm{fh}}{\▿\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h-{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h\\ +{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{fh}}d{\mathrm{\Γ}}_{\mathrm{fh}}={\∫}_{F_h}{q}_f{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h\end{array}---\left(2.6\right)$

合并方程(2.5)、(2.6)，整理后得到最终求解的弱形式：

${\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_{\mathrm{mh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_{\mathrm{mh}}{\▿\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h-{\∫}_{{\mathrm{\Γ}}_{\mathrm{mh}}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_{\mathrm{mh}}}{\∂n}{\mathrm{\δP}}_{\mathrm{mh}}{\mathrm{d\Γ}}_{\mathrm{mh}}$

$+{\∫}_{F_h}{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_{\mathrm{fh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h+{\∫}_{F_h}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_{\mathrm{fh}}{\▿\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h---\left(2.7\right)$

$+{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}})({\mathrm{\δP}}_{\mathrm{mh}}-{\mathrm{\δP}}_{\mathrm{fh}}){\mathrm{dF}}_h-{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}})({\mathrm{\δP}}_{\mathrm{mh}}-{\mathrm{\δP}}_{\mathrm{fh}})d{\mathrm{\Γ}}_{\mathrm{fh}}$

$={\∫}_{{\mathrm{\Ω}}_h}{q}_m{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{F_h}{q}_f{\mathrm{\δP}}_{\mathrm{fh}}d{F}_h$

写成矩阵的形式，即求解如下线性代数方程组：

$\left(\begin{array}{cc}{A}_m& B_{\mathrm{mf}}\\ B_{\mathrm{mf}}^{\mathrm{\τ}}& A_f\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{mh}}\\ P_{\mathrm{fh}}\end{array}\right)=\left(\begin{array}{c}{F}_m\\ F_f\end{array}\right)---\left(2.8\right)$

后得到基岩系统的有限元解P_mh与裂缝系统的有限元解P_fh；

(4)网格的处理

对基岩网格和裂缝网格独立进行剖分，具体方法如下：

1)、将不考虑裂缝存在的基岩单独造型划分网格，按照边界条件和无裂缝的物理参数计算基岩渗透矩阵；

2)、裂缝单独造型和划分网格，根据裂缝网格在基岩网格中的位置，将基岩变量在裂缝网格点上的值，用基岩网格上的节点的线性插值描达。

在所述步骤A中，探测得到该油藏区域的各项参数过程分别包括，

测试油藏的目标层的构造形态、缝洞体空间发育形态，充填溶洞内的物理性质参数：孔隙度、渗透率、孔喉半径；

油藏内油、水流体性质参数：不同压力下的密度、粘度、重度；

油藏内的岩石压缩系数，生产井、注水井的钻井位置及完井方式，油藏的边界条件参数指油藏的内外边界的性质，存在恒压力边界、等流量边界等；

地层初始参数：包括初始压力、初始油水分布。

所述步骤C中，初始化并设置的参数包括，

设置非线性分析控制参数：设置计算精度、迭代求解次数初始化计算数据：初始化地质模型的三维压力及饱和度数据和时间参数：即初始化设定最初的模拟时间。时间参数更新及初始化赋值：指计算完某一时间步后，把这一时间步计算的压力与饱和度数据赋值到下一个时间步，作为下一时间步的初始值。

所述步骤E，判断结果步骤中

(1)根据界面条件即：油水两种流体在界面两边的压力变化及速度变化相等条件判断迭代是否收敛，收敛则存储计算结果；不收敛则重复步骤D；

(2)判断分析结果：判断油水两种流体在界面两边的压力变化及速度变化是否相等，如相等则结果满意。

本发明是一种耦合型油藏数值模拟方法，是分析模拟缝洞型油藏流体流动的方法，所述方法根据油水的流动特征将缝洞型油藏分为三个区域：大型洞穴区(以下简称洞穴区)、多孔介质区域和裂缝；

所述多孔介质区域包括孔洞：长度、宽度、高度均小于0.5米，和微小裂缝，裂缝开度小于0.1毫米，延伸长度小于0.5米；所述裂缝是指岩石受外力或内应力作用时丧失结合力产生破裂但没有产生位移的区域，包括网状裂缝，裂缝开度在0.1毫米和10毫米之间，延伸长度小于0.5米和500米之间，和大型裂缝，也叫特大缝，裂缝开度大于10毫米，延伸长度大于500米。结合缝洞体的形成机制分析，大型洞穴与大裂缝为同一组合体特征，流体在其空间内流动特征密不可分，为此计算流体流动时，大型裂缝融合在洞穴区进行统一研究计算。

所述方法包括以下步骤：

(1)建立复杂介质Navier-Stokes方程和Darcy方程耦合数学模型；

(2)计算所述耦合数学模型的数值；

所述步骤(1)包括：

用Navier-Stokes方程表示流体在洞穴区中流动的运动方程；

用Darcy定律表示流体在多孔介质区域中流动的运动方程；

结合质量守恒定律、动量守恒定律和状态方程建立数学模型；

根据界面处法向速度连续和法向方向力的平衡条件建立界面之间的耦合条件，再根据界面处的耦合条件确定界面处的耦合模型；

所述步骤(2)包括：

洞穴区采用算子分裂法CBS；

多孔介质区域采用解耦全隐式算法；

裂缝采用“锚杆模型”算法。

所述方法采用计算机流程实现如下：

(1)程序启动；

(2)前处理建立几何模型、网格模型；

(3)设置岩体材料模型和材料参数；

(4)设定边值条件，给出地层初值信息；

(5)设定非线性分析控制参数；

(6)初始化计算数据和时间参数；

(7)时间参数更新；

(8)初始化迭代初值；

(9)求解洞穴区油水两相流动方程；

(10)求解多孔介质区油水两相渗流方程；

(11)由界面条件判断迭代是否收敛，如果是则转入(12)，否则更新迭代值，然后转入(9)；

(12)存储当前步时间计算结果；

(13)判断是否达到计算结束时间，如果达到则转入(14)，否则转入(7)；

(14)判断计算分析结果是否满意，如果满意则转入(15)，否则转入(3)；

(15)后处理显示结果；

(16)关闭程序。

本发明与现有技术相比，有益效果是：(1)针对缝洞型油藏首次提出了基于Navier-Stokes方程和Darcy方程的耦合型油藏数值模拟方法，发展了油藏数值模拟理论和方法；(2)将岩体与裂缝分开造型，各自造型不存在困难，对于裂缝的描述也清晰容易，裂缝数量及空间产状不受限制，特别是不要求裂缝上的节点与基岩节点完全重合，使得三维工作容易很多，科学地实现了对缝洞型油藏的模拟；(3)实现了一种模拟缝洞型油藏流体流动的方法。

附图说明

下面结合附图对本发明作进一步详细描述：

图1是本发明中缝洞型油藏的示意图。图中说明了缝洞型油藏存在大型洞穴区

多孔介质区

与裂缝图中箭头表示流体向井筒流动的方向。

图2是本发明中的计算流程图。

图3是本发明中裂缝处理方法示意图。

图4是本发明实施例的模型示意图。

图5是本发明实施例的洞穴区压力图。

图6是本发明实施例的多孔介质区压力图。

图7是本发明实施例的压力耦合结果图。

图8是本发明实施例的整体区域速度矢量图。

图9是本发明实施例的整体区域速度分布图。

图10是发明实施例的油藏模型饱和度分布图和时间步与迭代步长的关系图。

图11是实施例油藏模型饱和度场分布图。

具体实施方式

(1)建立复杂介质Navier-Stokes方程和Darcy方程耦合数学模型

对缝洞型碳酸盐岩油藏，流体存在两种流动形式：流体在基岩、溶孔和微裂缝中的流动表现为渗流特征；在大裂缝和溶洞中表现为自由流动特征。前者用Darcy定律表示其运动方程，后者用Navier-Stokes(简称N-S)方程表示其运动方程。结合质量守恒定律、动量守恒定律和状态方程建立数学模型。根据界面处法向速度连续和法向方向力的平衡条件建立界面之间的耦合条件，再根据界面处的耦合条件确定界面处的耦合模型。

A，洞穴区不可压缩流体(油和水)两相流动方程

如图1所示，Ω_v为洞穴区。油水两相在溶洞中的流动过程中，设油和水为不可压缩流体，流体的密度不随压力发生变化，即ρ＝常数。

1、连续性方程：

$\▿\·\left(S_w{u}_w\right)=0---\left(1.1\right)$

$\▿\·\left(S_o{u}_o\right)=0---(1.2)$

2、运动方程：

$\frac{\∂}{\∂t}\left(S_w{\mathrm{\ρ}}_w{u}_w\right)=S_w{\mathrm{\ρ}}_{w}g-{\▿p}_w+{\mathrm{\μ}}_w{\▿}^2{u}_w+S_w{F}_w---\left(1.3\right)$

$\frac{\∂}{\∂t}\left(S_o{\mathrm{\ρ}}_o{u}_o\right)=S_o{\mathrm{\ρ}}_{o}g-{\▿p}_o+{\mathrm{\μ}}_o{\▿}^2{u}_o+S_o{F}_o---(1.4)$

忽略掉重力： $\frac{\∂}{\∂t}\left(S_w{\mathrm{\ρ}}_w{u}_w\right)=-{\▿p}_w+{\mathrm{\μ}}_w{\▿}^2{u}_w+S_w{F}_w---\left(1.3a\right)$

$\frac{\∂}{\∂t}\left(S_o{\mathrm{\ρ}}_o{u}_o\right)=-{\▿p}_o+{\mathrm{\μ}}_o{\▿}^2{u}_o+S_o{F}_o---(1.4a)$

3，油相和水相的体积分数有如下关系：

S_w+S_o＝1                                (1.5)

4，在油水接触面上，压力相等：P_w＝P_o     (1.6)

式中：S_o为油占据的体积百分数；S_w为水占据的体积百分数；u_k为k相流速

方程中u_o、u_w、p_o、p_w、S_o和S_w是未知数。六个方程，六个未知数，因此方程是封闭的。

B，多孔介质中不可压缩流体油水两相的渗流方程

如图1所示，Ω_φ为多孔介质区域。设油和水为不可压缩流体，那么流体密度是常量，不用考虑流体的状态方程。

流体运动方程仍然为达西(Darcy)定律：

$u_{\mathrm{pi}}=-\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}({\▿p}_i-{\mathrm{\ρ}}_{i}F),i=o,w---\left(1.7\right)$

若忽略单位质量上的力(其中包括重力等)，其形式为：

$u_{\mathrm{pi}}=-\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}{\▿p}_i,i=o,w---(1.8)$

根据质量守恒方程可推得连续性方程为：

$\▿\·u_{\mathrm{pi}}=0,i=o,w---\left(1.9\right)$

联立运动方程、连续性方程，得到不可压缩油水两相流体的渗流控制方程，即：

$\▿\·\left(\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}{\▿p}_i\right)=\frac{\∂}{\∂t}\left({\mathrm{\φS}}_i\right),i=o,w---\left(1.10\right)$

p_c＝p_o-p_w    (1.11)

S_w+S_o＝1     (1.12)

S_i是各相的体积百分数，K_ri是各相的相对渗透率，p_c是毛管力。

此方程组的使用条件是：(1)彼此不互溶和不起化学反应的油水两相同时流动；(2)岩石和流体均不可压缩，服从线性渗流规律。

方程未知量的个数与方程数相等，因此方程组是封闭的。

C，界面耦合条件

如图1所示，Ω_f为裂缝。

法向速度连续性条件：

在界面上满足力的平衡条件：

$({\mathrm{\μ}}_o\frac{\∂u}{\∂n}-P\·n){|}_{{\mathrm{\Σ}}_f}=P_{p,o}\·n{|}_{{\mathrm{\Σ}}_P}$

$({\mathrm{\μ}}_w\frac{\∂u}{\∂n}-P\·n){|}_{{\mathrm{\Σ}}_f}=P_{p,w}\·n{|}_{{\mathrm{\Σ}}_P}$

(2)N-S和Darcy耦合模型数值计算方法

A，洞穴区采用算子分裂法

采用算子分裂法将方程组分离成对流项和扩散项(0＜θ₁＜θ₂≤1)：

对流项： $\frac{u^{{n+\mathrm{\θ}}_1}-u^n}{\mathrm{dt}}-{\mathrm{\μ\Δu}}^{n+{\mathrm{\θ}}_1}+{\▿p}^{n+{\mathrm{\θ}}_1}=0$

扩散项： $\frac{u^{{n+\mathrm{\θ}}_2}-u^{n+{\mathrm{\θ}}_1}}{\mathrm{dt}}-u^{n+{\mathrm{\θ}}_2}+{\▿u}^{n+{\mathrm{\θ}}_2}=0$

将对流方程(写成一般对流方程形式

以便通用，φ为通用符号)

对φ做Taylor展开，得到：

$\mathrm{\φ}(t+\mathrm{dt},x)=\mathrm{\φ}(t,x)+\mathrm{dt}\frac{\∂\mathrm{\φ}(t,x)}{\∂t}+\frac{{\mathrm{dt}}^2}{2}\frac{{\∂}^2\mathrm{\φ}(t,x)}{{\∂t}^2}+O\left({\mathrm{dt}}^3\right)$

由

上式中

均用

替换。于是，上式可改写成

$\mathrm{\φ}(t+\mathrm{dt},x)=\mathrm{\φ}(t,x)-\mathrm{dt}(u\·\▿\mathrm{\φ}(t,x))-\frac{{\mathrm{dt}}^2}{2}\frac{\∂}{\∂t}(u\·\▿\mathrm{\φ}(t,x))$

其中u为上一时间步的速度。

算子分裂法是求解流体问题的一种经典的数值模拟方法，由于是显式求解，能大大节约计算时间。

B，多孔介质区采用解耦全隐式算法

解耦全隐式算法是把压力和饱和度分开来顺序求解，每次只解一个未知变量，总体交替迭代，直至收敛。该解耦的全隐式算法比普通不解耦的全隐式算法计算速度快7-8倍，既保证了收敛性和精度，又能节省大量的存储量。。

C.裂缝处理的数值算法

本发明采用“锚杆模型”来处理裂缝。“锚杆模型”是一种适用于渗透矩阵加强型的组合网格模型，专门适用于解决裂缝的模拟问题，如图3所示，图中左边是基岩网格，中间是裂缝网格，右边是锚杆模型原理示意。采用“锚杆模型”，岩体与裂缝分开进行造型，裂缝的数量及空间产状不受限制，特别是不要求裂缝上的节点与基岩节点完全重合，使得三维工作容易很多。经过实际应用已经证明该方法能够较好地模拟裂缝。

锚杆模型具体算法如下：

(1)基本方程

为描述简便，以下所有的方程中，我们用下标m表示与基岩系统有关的参数与变量，用下标f表示与裂缝系统有关的参数与变量。

基岩系统渗流方程：

$\▿\left({\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_m\right)+q_m={\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_m}{\∂t}---\left(2.1\right)$

裂缝系统渗流方程：

$\▿\left({\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_f\right)+q_f={\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_f}{\∂t}---(2.2)$

(2)有限元弱形式

在进行有限元计算前我们需要得到方程的弱形式，

裂缝系统的弱形式如下：

${\∫}_F{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_f}{\∂t}{\mathrm{\δP}}_f\mathrm{dF}+{\∫}_F{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_f\▿\mathrm{\δ}{P}_f\mathrm{dF}-{\∫}_F\tilde{k}(P_m-P_f){\mathrm{\δP}}_f\mathrm{dF}$

$-{\∫}_{{\mathrm{\Γ}}_f}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_m-P_f){\mathrm{\δP}}_{f}d{\mathrm{\Γ}}_f={\∫}_F{q}_f{\mathrm{\δP}}_f\mathrm{dF}---\left(2.3\right)$

为了保证基岩系统与裂缝系统的流量守恒，我们还需将裂缝边界的流入(出)量添加到基岩系统的方程中，即最终基岩系统的弱形式如下：

${\∫}_{\mathrm{\Ω}}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_m}{\∂t}{\mathrm{\δP}}_m\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_m\▿\mathrm{\δ}{P}_m\mathrm{d\Ω}-{\∫}_F\tilde{k}(P_m-P_f){\mathrm{\δP}}_m\mathrm{dF}$

${-\∫}_{{\mathrm{\Γ}}_m}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_m}{\∂n}{\mathrm{\δP}}_{m}d{\mathrm{\Γ}}_m+{\∫}_{{\mathrm{\Γ}}_f}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_m-P_f){\mathrm{\δP}}_{m}d{\mathrm{\Γ}}_f={\∫}_{\mathrm{\Ω}}{q}_m{\mathrm{\δP}}_m\mathrm{d\Ω}---\left(2.4\right)$

(3)有限元解及基岩与裂缝交换项的处理

令mh表示基岩区域Ω上的网格剖分，fh表示裂缝区域上F的网格剖分，比如二维问题，基岩内为二维流动，裂缝内为一维流动，这时mh表示一个二维的网格剖分，fh表示一个一维的网格剖分。

令P_mh表示基岩系统的有限元解，P_fh表示裂缝系统的有限元解，Ω_h为剖分后的基岩区域，F_h为剖分后的裂缝区域，Γ_mh为剖分后基岩区域的边界，Γ_fh为剖分后的裂缝区域的边界，我们需要寻找P_mh及P_fh使得其满足方程：

$\begin{array}{c}{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_{\mathrm{mh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_{\mathrm{mh}}\▿\mathrm{\δ}{P}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{mh}}{\mathrm{dF}}_h\\ -{\∫}_{{\mathrm{\Γ}}_{\mathrm{mh}}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_{\mathrm{mh}}}{\∂n}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Γ}}_{\mathrm{mh}}-{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Γ}}_{\mathrm{fh}}={\∫}_{{\mathrm{\Ω}}_h}{q}_m{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h\end{array}---\left(2.5\right)$

$\begin{array}{c}{\∫}_{F_h}{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_{\mathrm{fh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h+{\∫}_{F_h}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_{\mathrm{fh}}{\▿\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h-{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h\\ +{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{fh}}d{\mathrm{\Γ}}_{\mathrm{fh}}={\∫}_{F_h}{q}_f{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h\end{array}---\left(2.6\right)$

合并方程(2.5)、(2.6)，整理后得到最终求解的弱形式：

${\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_{\mathrm{mh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_{\mathrm{mh}}{\▿\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h-{\∫}_{{\mathrm{\Γ}}_{\mathrm{mh}}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_{\mathrm{mh}}}{\∂n}{\mathrm{\δP}}_{\mathrm{mh}}{\mathrm{d\Γ}}_{\mathrm{mh}}$

$+{\∫}_{F_h}{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_{\mathrm{fh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h+{\∫}_{F_h}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_{\mathrm{fh}}{\▿\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h---\left(2.7\right)$

$+{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}})({\mathrm{\δP}}_{\mathrm{mh}}-{\mathrm{\δP}}_{\mathrm{fh}}){\mathrm{dF}}_h-{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}})({\mathrm{\δP}}_{\mathrm{mh}}-{\mathrm{\δP}}_{\mathrm{fh}})d{\mathrm{\Γ}}_{\mathrm{fh}}$

$={\∫}_{{\mathrm{\Ω}}_h}{q}_m{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{F_h}{q}_f{\mathrm{\δP}}_{\mathrm{fh}}d{F}_h$

写成矩阵的形式，即求解如下线性代数方程组：

$\left(\begin{array}{cc}{A}_m& B_{\mathrm{mf}}\\ B_{\mathrm{mf}}^{\mathrm{\τ}}& A_f\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{mh}}\\ P_{\mathrm{fh}}\end{array}\right)=\left(\begin{array}{c}{F}_m\\ F_f\end{array}\right)---\left(2.8\right)$

后得到基岩系统的有限元解P_mh与裂缝系统的有限元解P_ft。

(4)网格的处理

在基岩中存在裂缝时，构造网格普遍比较困难，而且构造出来的网格质量也难以保证，所以基岩和裂缝一同构造网格的做法在应用中常常受到限制。本方法中，对基岩网格和裂缝网格独立进行剖分，适用于解决裂缝的模拟问题。

具体做法如下：

1)、将不考虑裂缝存在的基岩单独造型划分网格如图3，按照边界条件和无裂缝的物理参数计算基岩渗透矩阵。

2)、裂缝单独造型和划分网格，根据裂缝网格在基岩网格中的位置，将基岩变量在裂缝网格点上的值，用基岩网格上的节点的线性插值描达，如基岩变量在裂缝节点①上的自由度用基岩节点5、7、12、10四个节点的自由度来组合书写，基岩变量在裂缝节点②的自由度用基岩节点10、12、17、15四个节点的自由度来组合书写，基岩变量在裂缝节点③的自由度用基岩节点12、13、18、17四个节点的自由度来组合书写，基岩变量在裂缝节点④的自由度用基岩节点18、19、24、23四个节点的自由度来组合书写。

上述做法的优势是将岩体与裂缝分开造型，各自造型不存在困难，对于裂缝的描述也清晰容易，裂缝数量及空间产状不受限制，特别是不要求裂缝上的节点与基岩节点完全重合，使得三维工作容易很多。

N-S和Darcy耦合模型数值计算方法的计算流程如图2所示，包括以下步骤：

(1)程序启动；

(2)前处理建立几何模型、网格模型；

(3)设置岩体材料模型和材料参数；

(4)设定边值条件，给出地层初值信息；

(5)设定非线性分析控制参数；

(6)初始化计算数据和时间参数；

(7)时间参数更新；

(8)初始化迭代初值；

(9)求解洞穴区油水两相流动方程；

(10)求解多孔介质区油水两相渗流方程；

(11)由界面条件判断迭代是否收敛，如果是则转入(12)，否则更新迭代值，然后转入(9)；

(12)存储当前步时间计算结果；

(13)判断是否达到计算结束时间，如果达到则转入(14)，否则转入(7)；

(14)判断计算分析结果是否满意，如果满意则转入(15)，否则转入(3)；

(15)后处理显示结果；

(16)关闭程序。

为了验证方法的正确性，根据图2所示的流程编制了该方法的数值模拟程序，并进行了模拟计算。模型为1m×1m的正方形，中心是0.01m×0.01m的洞，其余部分为多孔介质区。模型左边界压力为10Mpa，右边界压力为1Mpa，流体(油和水)从左向右流动，如图4所示。模拟结果如图5-图9所示：图5给出了自由流动区的压力模拟结果；图6给出了多孔介质区的压力模拟结果；图7给出了压力耦合的模拟结果；图8给出了整体区域速度矢量的模拟结果；图9给出了整体区域速度分布的模拟结果。

建立油藏地质模型如图10所示，模型主体尺寸：3X0.3m；出口宽度：0.1m；其假设容器左右两侧、底部为入口，左侧上半部分(0.15m)为油相入口，下半部分(0.15m)为水相入口，右侧与左侧相同；底部入口为水相入口。油、水流入速度均为0.017m/s；油的密度为960kg/m^3，粘度为0.02kg/m.s；水的密度为1140kg/m^3，粘度为0.001kg/m.s。井打在洞穴上，后期油水同出。

建立了如图11所示的模型，该模型为两个裂缝连通的溶洞，底部都有水注入，观测溶洞内原油流体情况。其假设两个溶洞内没有填充物，左边溶洞顶部有一采油井、底水沿底部裂缝窜进。初始时两溶洞充满油。模拟结果说明沟通位置决定了洞内剩余油的分布

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (9)

1.一种分析模拟缝洞型油藏流体流动的方法，其特征在于，所述分析方法包括如下过程，

A确定缝洞型油藏的基础参数，包括：测试油藏的目标层的构造形态、缝洞体空间发育形态，充填溶洞内的物理性质参数；和油藏内油、水流体性质参数；以及油藏内的岩石压缩系数，生产井、注水井的钻井位置及完井方式，油藏的边界条件参数和地层初始参数；

B创建缝洞型油藏耦合数学模型，根据缝洞型油藏油水的流动特征分为三个区域：大型洞穴区、多孔介质区域和裂缝带；所述大型洞穴区包括大型洞穴和大裂缝；所述多孔介质区域包括基岩、溶孔和微裂缝；所述裂缝包括岩石受外力或内应力作用时丧失结合力产生破裂但没有产生位移的区域，主要为网状裂缝带；

根据上述划分的不同地质区域，建立复杂介质Navier-Stokes方程和Darcy方程耦合数学模型：用Navier-Stokes方程表示流体在洞穴区中流动的运动方程；用Darcy定律表示流体在多孔介质区域中流动的运动方程；结合质量守恒定律、动量守恒定律和状态方程建立数学模型；根据界面处法向速度连续和法向方向力的平衡条件建立界面之间的耦合条件，再根据界面处的耦合条件确定界面处的耦合模型；

C初始化并设置参数步骤：用于依次设置非线性分析控制参数，初始化计算数据和时间参数，对时间参数更新和初始化迭代初值；

D依据所述步骤B的模型，求解洞穴区油水两相流动方程和多孔介质区油水两相渗流方程；其中，洞穴区油水两相流动方程采用算子分裂法CBS；多孔介质区域油水两相渗流方程采用解耦全隐式算法；对裂缝采用锚杆模型算法；

E判断结果步骤：

(1)根据界面条件判断迭代是否收敛，收敛则存储计算结果；不收敛则重复步骤D；

(2)判断计算时间是否结束，是则进入下面判断过程，否则进入步骤C；

(3)判断分析结果；

F根据缝洞型油藏地质模型，建立模拟前需要的三维定量数据体，导入油田生产井产液量、产油量及测压数据，进行流体流动预测与分析，显示待测缝洞型油藏区域油水分布结果。

2.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤B中包括，

流体在洞穴区中流动的运动方程，如下：

油水两相在溶洞中流动过程中，设油和水为不可压缩流体，流体的密度不随压力发送变化，即ρ＝常数；

连续性方程：

$\▿\·\left(S_w{u}_w\right)=0---\left(1.1\right)$

$\▿\·\left(S_o{u}_o\right)=0---(1.2)$

运动方程：

$\frac{\∂}{\∂t}\left(S_w{\mathrm{\ρ}}_w{u}_w\right)=S_w{\mathrm{\ρ}}_{w}g-{\▿p}_w+{\mathrm{\μ}}_w{\▿}^2{u}_w+S_w{F}_w---\left(1.3\right)$

$\frac{\∂}{\∂t}\left(S_o{\mathrm{\ρ}}_o{u}_o\right)=S_o{\mathrm{\ρ}}_{o}g-{\▿p}_o+{\mathrm{\μ}}_o{\▿}^2{u}_o+S_o{F}_o---(1.4)$

忽略掉重力： $\frac{\∂}{\∂t}\left(S_w{\mathrm{\ρ}}_w{u}_w\right)=-{\▿p}_w+{\mathrm{\μ}}_w{\▿}^2{u}_w+S_w{F}_w---\left(1.3a\right)$

$\frac{\∂}{\∂t}\left(S_o{\mathrm{\ρ}}_o{u}_o\right)=-{\▿p}_o+{\mathrm{\μ}}_o{\▿}^2{u}_o+S_o{F}_o---(1.4a)$

油相和水相的体积分数有如下关系：

                          S_w+S_o＝1                      (1.5)

在油水接触面上，压力相等：P_w＝P_o                        (1.6)

式中：S_o为油占据的体积百分数；S_w为水占据的体积百分数；u_k为k相流速；方程中u_o、u_w、p_o、p_w、S_o和S_w是未知数。

3.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤B中，流体在多孔介质区域中流动的运动方程如下：

设油和水为不可压缩流体，那么流体密度是常量，不用考虑流体的状态方程；

流体运动方程：

$u_{\mathrm{pi}}=-\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}({\▿p}_i-{\mathrm{\ρ}}_{i}F),i=o,w---\left(1.7\right)$

若忽略单位质量上的力，其形式为：

$u_{\mathrm{pi}}=-\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}{\▿p}_i,i=o,w---(1.8)$

根据质量守恒方程可推得连续性方程为：

$\▿\·u_{\mathrm{pi}}=0,i=o,w---\left(1.9\right)$

联立运动方程、连续性方程，得到不可压缩油水两相流体的渗流控制方程，即：

$\▿\·\left(\frac{{\mathrm{KK}}_{\mathrm{ri}}}{{\mathrm{\μ}}_i}{\▿p}_i\right)=\frac{\∂}{\∂t}\left({\mathrm{\φS}}_i\right),i=o,w---\left(1.10\right)$

p_c＝p_o-p_w    (1.11)

S_w+S_o＝1     (1.12)

S_i是各相的体积百分数，K_ri是各相的相对渗透率，p_c是毛管力。

此方程组的使用条件是：(1)彼此不互溶和不起化学反应的油水两相同时流动；(2)岩石和流体均不可压缩，服从线性渗流规律。

4.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤B中，界面之间的耦合条件为：

法向速度连续性条件：

在界面上满足力的平衡条件：

$({\mathrm{\μ}}_o\frac{\∂u}{\∂n}-P\·n){|}_{{\mathrm{\Σ}}_f}=P_{p,o}\·n{|}_{{\mathrm{\Σ}}_P}$

$({\mathrm{\μ}}_w\frac{\∂u}{\∂n}-P\·n){|}_{{\mathrm{\Σ}}_f}=P_{p,w}\·n{|}_{{\mathrm{\Σ}}_P}$

其中，u为法向流速；P为法向压力；n为法向方向；μ_o为油的粘度；μ_w为水的粘度。

5.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤D中，算子分裂法CBS算法如下：

分离成对流项和扩散项(0＜θ₁＜θ₂≤1)：

对流项： $\frac{u^{{n+\mathrm{\θ}}_1}-u^n}{\mathrm{dt}}-{\mathrm{\μ\Δu}}^{n+{\mathrm{\θ}}_1}+{\▿p}^{n+{\mathrm{\θ}}_1}=0$

扩散项： $\frac{u^{{n+\mathrm{\θ}}_2}-u^{n+{\mathrm{\θ}}_1}}{\mathrm{dt}}-u^{n+{\mathrm{\θ}}_2}+{\▿u}^{n+{\mathrm{\θ}}_2}=0$

将对流项方程写成时间与空间的形式，为

其中φ为通用符号包括u和p，对φ做Taylor展开，得到：

$\mathrm{\φ}(t+\mathrm{dt},x)=\mathrm{\φ}(t,x)+\mathrm{dt}\frac{\∂\mathrm{\φ}(t,x)}{\∂t}+\frac{{\mathrm{dt}}^2}{2}\frac{{\∂}^2\mathrm{\φ}(t,x)}{{\∂t}^2}+O\left({\mathrm{dt}}^3\right)$

由

上式中均用

替换，于是，上式改写成

$\mathrm{\φ}(t+\mathrm{dt},x)=\mathrm{\φ}(t,x)-\mathrm{dt}(u\·\▿\mathrm{\φ}(t,x))-\frac{{\mathrm{dt}}^2}{2}\frac{\∂}{\∂t}(u\·\▿\mathrm{\φ}(t,x))$

其中u为上一时间步的速度；

所述解耦全隐式算法是把压力和饱和度分开来顺序求解，每次只解一个未知变量，总体交替迭代，直至收敛。

6.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，所述步骤D中，锚杆模型算法的具体算法如下：

(1)基本方程

基岩系统渗流方程：

$\▿\left({\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_m\right)+q_m={\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_m}{\∂t}---\left(2.1\right)$

裂缝系统渗流方程：

$\▿\left({\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_f\right)+q_f={\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_f}{\∂t}---(2.2)$

式中，下标m表示与基岩系统有关的参数与变量，下标f表示与裂缝系统有关的参数与变量；

(2)有限元弱形式

裂缝系统的弱形式如下：

${\∫}_F{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_f}{\∂t}{\mathrm{\δP}}_f\mathrm{dF}+{\∫}_F{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_f\▿\mathrm{\δ}{P}_f\mathrm{dF}-{\∫}_F\tilde{k}(P_m-P_f){\mathrm{\δP}}_f\mathrm{dF}$

$-{\∫}_{{\mathrm{\Γ}}_f}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_m-P_f){\mathrm{\δP}}_{f}d{\mathrm{\Γ}}_f={\∫}_F{q}_f{\mathrm{\δP}}_f\mathrm{dF}---\left(2.3\right)$

将裂缝边界的流入流出量添加到基岩系统的方程中，得到最终基岩系统的弱形式如下：

${\∫}_{\mathrm{\Ω}}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_m}{\∂t}{\mathrm{\δP}}_m\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_m\▿\mathrm{\δ}{P}_m\mathrm{d\Ω}-{\∫}_F\tilde{k}(P_m-P_f){\mathrm{\δP}}_m\mathrm{dF}$

${-\∫}_{{\mathrm{\Γ}}_m}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_m}{\∂n}{\mathrm{\δP}}_{m}d{\mathrm{\Γ}}_m+{\∫}_{{\mathrm{\Γ}}_f}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_m-P_f){\mathrm{\δP}}_{m}d{\mathrm{\Γ}}_f={\∫}_{\mathrm{\Ω}}{q}_m{\mathrm{\δP}}_m\mathrm{d\Ω}---\left(2.4\right)$

(3)有限元解及基岩与裂缝交换项的处理

令mh表示基岩区域Ω上的网格剖分，fh表示裂缝区域上F的网格剖分，由于基岩内为二维流动，裂缝内为一维流动，此时mh表示一个二维的网格剖分，fh表示一个一维的网格剖分；

令P_mh表示基岩系统的有限元解，P_fh表示裂缝系统的有限元解，Ω_h为剖分后的基岩区域，F_h为剖分后的裂缝区域，Γ_mh为剖分后基岩区域的边界，Γ_fh为剖分后的裂缝区域的边界，寻找P_mh及P_fh使得其满足方程：

$\begin{array}{c}{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_{\mathrm{mh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_{\mathrm{mh}}\▿\mathrm{\δ}{P}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{mh}}{\mathrm{dF}}_h\\ -{\∫}_{{\mathrm{\Γ}}_{\mathrm{mh}}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_{\mathrm{mh}}}{\∂n}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Γ}}_{\mathrm{mh}}-{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Γ}}_{\mathrm{fh}}={\∫}_{{\mathrm{\Ω}}_h}{q}_m{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h\end{array}---\left(2.5\right)$

$\begin{array}{c}{\∫}_{F_h}{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_{\mathrm{fh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h+{\∫}_{F_h}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_{\mathrm{fh}}{\▿\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h-{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h\\ +{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}}){\mathrm{\δP}}_{\mathrm{fh}}d{\mathrm{\Γ}}_{\mathrm{fh}}={\∫}_{F_h}{q}_f{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h\end{array}---\left(2.6\right)$

合并方程(2.5)、(2.6)，整理后得到最终求解的弱形式：

${\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m{c}_m{\mathrm{\φ}}_m\frac{{\∂P}_{\mathrm{mh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{{\mathrm{\Ω}}_h}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}{\▿P}_{\mathrm{mh}}{\▿\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h-{\∫}_{{\mathrm{\Γ}}_{\mathrm{mh}}}{\mathrm{\ρ}}_m\frac{k_m}{{\mathrm{\μ}}_m}\frac{{\∂P}_{\mathrm{mh}}}{\∂n}{\mathrm{\δP}}_{\mathrm{mh}}{\mathrm{d\Γ}}_{\mathrm{mh}}$

$+{\∫}_{F_h}{\mathrm{\ρ}}_f{c}_f{\mathrm{\φ}}_f\frac{{\∂P}_{\mathrm{fh}}}{\∂t}{\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h+{\∫}_{F_h}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}{\▿P}_{\mathrm{fh}}{\▿\mathrm{\δP}}_{\mathrm{fh}}{\mathrm{dF}}_h---\left(2.7\right)$

$+{\∫}_{F_h}\tilde{k}(P_{\mathrm{mh}}-P_{\mathrm{fh}})({\mathrm{\δP}}_{\mathrm{mh}}-{\mathrm{\δP}}_{\mathrm{fh}}){\mathrm{dF}}_h-{\∫}_{{\mathrm{\Γ}}_{\mathrm{fh}}}{\mathrm{\ρ}}_f\frac{k_f}{{\mathrm{\μ}}_f}(P_{\mathrm{mh}}-P_{\mathrm{fh}})({\mathrm{\δP}}_{\mathrm{mh}}-{\mathrm{\δP}}_{\mathrm{fh}})d{\mathrm{\Γ}}_{\mathrm{fh}}$

$={\∫}_{{\mathrm{\Ω}}_h}{q}_m{\mathrm{\δP}}_{\mathrm{mh}}d{\mathrm{\Ω}}_h+{\∫}_{F_h}{q}_f{\mathrm{\δP}}_{\mathrm{fh}}d{F}_h$

写成矩阵的形式，即求解如下线性代数方程组：

$\left(\begin{array}{cc}{A}_m& B_{\mathrm{mf}}\\ B_{\mathrm{mf}}^{\mathrm{\τ}}& A_f\end{array}\right)\left(\begin{array}{c}{P}_{\mathrm{mh}}\\ P_{\mathrm{fh}}\end{array}\right)=\left(\begin{array}{c}{F}_m\\ F_f\end{array}\right)---\left(2.8\right)$

后得到基岩系统的有限元解P_mh与裂缝系统的有限元解P_fh；

(4)网格的处理

对基岩网格和裂缝网格独立进行剖分，具体方法如下：

1)、将不考虑裂缝存在的基岩单独造型划分网格，按照边界条件和无裂缝的物理参数计算基岩渗透矩阵；

2)、裂缝单独造型和划分网格，根据裂缝网格在基岩网格中的位置，将基岩变量在裂缝网格点上的值，用基岩网格上节点的线性插值描述。

7.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤A中，探测得到该油藏区域的各项参数过程分别包括，

测试油藏的目标层的构造形态、缝洞体空间发育形态，充填溶洞内的物理性质参数：孔隙度、渗透率、孔喉半径；

油藏内油、水流体性质参数：不同压力下的密度、粘度、重度；

油藏内岩石压缩系数，生产井、注水井的钻井位置及完井方式，油藏的边界条件参数指油藏的内外边界的性质，存在恒压力边界、等流量边界等；

地层初始参数：包括初始压力、初始油水分布。

8.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤C中，初始化并设置的参数包括，

设置非线性分析控制参数：设置计算精度、迭代求解次数初始化计算数据：初始化地质模型的三维压力及饱和度数据和时间参数：即初始化设定最初的模拟时间。

9.根据权利要求1所述的一种模拟缝洞型油藏流体流动的方法，其特征在于，

所述步骤E，判断结果步骤中

(1)根据界面条件即：油水两种流体在界面两边的压力变化及速度变化相等条件判断迭代是否收敛，收敛则存储计算结果；不收敛则重复步骤D；

(2)判断分析结果：判断油水两种流体在界面两边的压力变化及速度变化是否相等，相等则结果满意。

Images