CN102339325A - 一种分析离散裂缝性油藏流体流动的方法 - Google Patents

Abstract

本发明提供了一种离散裂缝网络模型数值模拟有限元数值求解方法，属于油气田开发领域。本方法根据离散裂缝网络模型，对缝洞型油藏的基岩和裂缝这两种不同的连续性介质采取不同的离散化方法，得到数学模型；然后对所述数学模型进行数值模拟，得到离散裂缝网络数值模拟；最后分析所述离散裂缝网络数值模拟的有限元单元的特性，对所述数学模型进行求解。本发明提高了数值模拟的计算速度、增加了计算的收敛性，实现了对缝洞型油藏的科学处理，为找准缝洞型油藏剩余油的分布位置、定量确定油藏的储量丰度、科学合理地开发这类油田提供依据，最终达到提高采收率的目的。

Description

一种分析离散裂缝性油藏流体流动的方法

技术领域

本发明属于油气田开发领域，具体涉及一种分析离散裂缝性油藏流体流动的方法。

背景技术

目前由于离散裂缝网络模型在一定程度上比较真实地描述了缝洞型油藏中裂缝的地质分布特征及其对流体流动的影响，所以采用离散裂缝网络模型数值模拟方法。同时因为有限差分方法数学概念直观、表达简单、占用内存小等特点，多在油藏数值模拟的求解方法中使用，但是有限单元法虽计算精度高、网格取向性弱且适合复杂边界问题，但计算量大、收敛性差使得油藏数值模拟中基本不采用这种方法。目前存在的主要难题是：(1)模型不能对裂缝、洞穴和基岩之间的窜流进行精确的评估且不能对大尺度洞穴和控制着油藏流体流动的大裂缝进行模拟；(2)有限差分法对于复杂边界问题、高阶微分方程、高梯度问题求解比较困难，而有限元法计算量大、收敛性差。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种用离散裂缝网络模型和有限单元法对缝洞型油藏油分布进行分析的方法。

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

一种有限元法分析离散裂缝油藏流体流动的方法，所述方法包括以下步骤：

A缝洞型油藏物探数据输入步骤，

B对缝洞型油藏的基岩和裂缝两种不同的连续性介质采取不同的离散化方法，即对裂缝采用线单元离散，对基岩采用三角形单元离散：即对每个待测区域V离散成E个单元，每个单元有N个结点；

C输入结点信息、地层属性数据、相对渗透率和PVT数据：即岩石和流体的地层条件下的物性数据

D分析得到基岩系统和裂缝系统压力的有限元单元特性函数及其饱和度的有限元单元特性函数；

E形成各单元特性矩阵和列阵步骤：形成基岩系统压力函数单元特征矩阵(方程(1-8)中矩阵[H^(e)]和方程(1-9)中列阵

)和裂缝系统压力函数单元特征矩阵(方程(2-6)中矩阵

和方程(2-7)中列阵)，基岩系统饱和度函数单元特征矩阵(方程(3-8)中矩阵[H′^(e)]和方程(3-9)中列阵

)和裂缝系统饱和度函数单元特征矩阵(方程(4-8)中矩阵[H′_f ^(e)]和方程(4-9)中列阵)；

F组装有限个基本单元的单元特征矩阵，形成总特征矩阵，即形成压力整体特征矩阵(方程(5-2)中矩阵

和列阵

)和饱和度整体特征矩阵(方程(5-4)中矩阵和列阵

)；

G根据各物理量求取步骤；

(1)对压力方程求解，即根据初始条件或上一步计算结果采用高斯多重积分方法对压力方程中各项积分的进行计算，形成压力方程的单元特性矩阵和列阵，然后组合压力方程的整体特性矩阵和列阵，最后采取预处理共轭梯度法对压力方程组进行求解；

(2)对饱和度方程求解，即根据这一步求出的压力和初始条件或上一步计算结果对饱和度方程中各项积分进行计算，形成饱和度方程的单元特性矩阵或列阵，然后组合饱和度方程的整体特性矩阵和列阵，最后对饱和度方程组进行求解；

H判断步骤：

(1)计算线性方程组中总装矩阵的条件数，如条件数数量级大于16，则认为不收敛，否则就认为收敛，从而进行下一时间步的计算；

(2)如果计算时间大约设定时间计算结束，否则进人下一时间步的计算；

I输出显示步骤：

输出显示压力场、饱和度场等的图像及动画以及收敛性曲线和生产指数曲线。

所述方法对裂缝进行了单独处理，根据地质上所描述的裂缝发育分布情况，在数值模拟中；在所述数值模拟中把大洞穴看做裂缝来进行处理，油藏模拟假设基岩中的流体和岩石均微可压缩，且压缩系数为常数，裂缝及裂缝中的流体为不可压缩；所述数学模型为：

基岩系统：

${-\mathrm{\φC}}_t\frac{{\∂p}_o}{\∂t}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-(q_{\mathrm{wv}}+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}(q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*}))---(5-5)$

$\mathrm{\φ}\frac{{\∂S}_o}{\∂t}+S_o{\mathrm{\φC}}_{\mathrm{to}}\frac{{\∂p}_o}{\∂t}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}{q}_{\mathrm{ovnf}}^{*}---(5-6)$

S_o+S_w＝1            (5-7)

P_c＝P_o-P_w           (5-8)

裂缝系统：

$\frac{\∂}{\∂l}{\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)}_{\mathrm{nf}}+(q_{\mathrm{wvnf}}^{*}+q_{\mathrm{ovnf}}^{*})=0---(5-9)$

$\frac{\∂}{\∂l}{\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right)}_{\mathrm{nf}}={\mathrm{\φ}}_{\mathrm{nf}}\frac{{\∂S}_{\mathrm{onf}}}{\∂t}-q_{\mathrm{ovnf}}^{*}---(5-10)$

S_onf+S_wnf＝1        (5-11)

基岩-裂缝间的窜流：

$q_{\mathrm{onf}}^{*}={\left({\mathrm{\ξ\λ}}_o(p_o-p)\right)}_{\mathrm{nf}}---(5-12)$

$q_{\mathrm{wnf}}^{*}={\left({\mathrm{\ξ\λ}}_w(p_w-p)\right)}_{\mathrm{nf}}---(5-13)$

其中， $\mathrm{\ξ}=\frac{\mathrm{lgh}}{d}.$

所述步骤D包括以下步骤：

I基岩系统压力有限元单元特性分析

基岩系统压力方程：

${-\mathrm{\φC}}_t\frac{{\∂p}_o}{\∂t}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-q_v---(1-1)$

其中， $q_v=q_{\mathrm{wv}}+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}(q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*})$

对时间项的处理：

${-\mathrm{\φC}}_t\frac{p_o-p_o^{(n-1)}}{\mathrm{\Δt}}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-q_v---(1-2)$

取

移项后得：

$-\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{p}_o+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-\mathrm{\γ}---(1-3)$

根据虚位移原理将式(1-1)化为基岩系统微分方程的等效积分“弱”形式并化简得：

${-\∫\∫}_{\mathrm{\Ω}}\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{p}_o{\mathrm{\δp}}_o\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}\▿\·\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\▿p}_o\right){\mathrm{\δp}}_o\mathrm{d\Ω}={\∫\∫}_{\mathrm{\Ω}}\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right){\mathrm{\δp}}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}\mathrm{\γ\δ}{p}_o\mathrm{d\Ω}---(1-4)$

采用分部积分法得：

根据边界条件且边界处在此暂不考虑毛管力，则

$I={\∫\∫}_{\mathrm{\Ω}}\frac{{\mathrm{\φC}}_t}{2\mathrm{\Δt}}{\left(p_o\right)}^2\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}\frac{({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)}{2}\left({\▿p}_o\right)\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\λ}}_w{\▿p}_c{\▿p}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γp}}_o\mathrm{d\Ω}---(1-6)$

所述步骤D中，用来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I取极值，则：

$\frac{\∂\left(I^{\left(e\right)}\right)}{\∂\left\{{p_o}^{\left(e\right)}\right\}}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{\left[N\right]}^T\left[N\right]\left\{{p_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂x}\right]}^T\right[\frac{\∂N}{\∂x}\left]\right\{{p_o}^{\left(e\right)}\}+({\mathrm{\λ}}_o+{\mathrm{\λ}}_e){\left[\frac{\∂N}{\∂y}\right]}^T[\frac{\∂N}{\∂y}\left]\right\{{p_o}^{\left(e\right)}\left\}\right)\mathrm{d\Ω}$

$-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_w\frac{\∂\left({p_c}^{\left(e\right)}\right)}{\∂x}{\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_w\frac{\∂\left({p_c}^{\left(e\right)}\right)}{\∂y}{\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\γ}{\left[N\right]}^T\mathrm{d\Ω}---(1-7)$

其中：

$\left[{H_1}^{\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{\left[N\right]}^T\right[N\left]\right)\mathrm{d\Ω}$

$\left[{H_2}^{\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂x}\right]}^T\right[\frac{\∂N}{\∂x}]+({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂y}\right]}^T[\frac{\∂N}{\∂y}\left]\right)\mathrm{d\Ω}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{d\Ω}$

$\left[B\right]=\left[\begin{array}{cccc}\frac{{\∂N}_1}{\∂x}& \frac{{\∂N}_2}{\∂x}& ...& \frac{{\∂N}_n}{\∂x}\\ \frac{{\∂N}_1}{\∂y}& \frac{{\∂N}_2}{\∂y}& ...& \frac{{\∂N}_n}{\∂y}\end{array}\right],\left[D\right]=\left[\begin{array}{cc}({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)& \\ & ({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\end{array}\right]$

$\left[H^{\left(e\right)}\right]=\left[H_1^{\left(e\right)}\right]+\left[{H_2}^{\left(e\right)}\right]---(1-8)$

最终得到：

II裂缝系统压力有限元单元特性分析

裂缝系统压力方程：

$\frac{\∂}{\∂l}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)+q_v^{*}=0---(2-1)$

其中， $q_v^{*}=q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*}$

根据虚位移原理将(2-1)化为裂缝系统微分方程的等效积分“弱”形式并化简得：

${\∫}_l\frac{\∂}{\∂l}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)\mathrm{\δpdl}+{\∫}_l{q}_v^{*}\mathrm{\δpdl}=0---(2-2)$

采用分部积分法得：

结合边界条件可得：

$I={\∫}_l\frac{1}{2}({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left(\frac{\∂p}{\∂l}\right)}^2\mathrm{dl}-{\∫}_l{q}_v^{*}\mathrm{pdl}---(2-4)$

根据所述步骤D中，分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I_f取极值，则：

$\frac{\∂\left({I_f}^{\left(e\right)}\right)}{\∂\left\{p^{\left(e\right)}\right\}}={\∫}_{l^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{{\∂N}^{\′}}{\∂l}\right]}^T\right[\frac{{\∂N}^{\′}}{\∂l}\left]\right\{p^e\left\}\right)\mathrm{dl}-{\∫}_{l^{\left(e\right)}}{q}_v^{*}{\left[N\right]}^T\mathrm{dl}---(2-5)$

其中：

$\left[{H_f}^{\left(e\right)}\right]={\∫}_{l^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N^{\′}}{\∂l}\right]}^T\right[\frac{\∂N^{\′}}{\∂l}\left]\right\{p^{\left(e\right)}\left\}\right)\mathrm{dl}={\∫}_l{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{dl}---(2-6)$

$\left[B\right]=\left[\begin{array}{cc}\frac{{\∂N}_1^{\′}}{\∂l}& \frac{{\∂N}_2^{\′}}{\∂l}\end{array}\right],\left[D\right]=({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)$

最终得到：

III基岩系统饱和度有限元单元特性分析

基岩系统饱和度方程：

$\mathrm{\φ}\frac{{\∂S}_o}{\∂t}+S_o\mathrm{\φ}{C}_{\mathrm{to}}\frac{{\∂p}_o}{\∂t}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}^{\′}---(3-1)$

其中： $q_{\mathrm{ov}}^{\′}=q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}{q}_{\mathrm{ovnf}}^{*}$

对时间项处理后得：

$\mathrm{\φ}\frac{S_o-S_o^{(n-1)}}{\mathrm{\Δt}}+S_o{\mathrm{\φC}}_{\mathrm{to}}\frac{p_o-p_o^{(n-1)}}{\mathrm{\Δt}}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}^{\′}---(3-2)$

取

移项后得：

${\mathrm{\ψS}}_o=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+{\mathrm{\γ}}_o^{\′}---(3-3)$

根据虚位移原理将(3-1)化为基岩系统微分方程的等效积分“弱”形式并化简得：

${\∫\∫}_{\mathrm{\Ω}}{\mathrm{\ψS}}_o{\mathrm{\δS}}_o\mathrm{d\Ω}={\∫\∫}_{\mathrm{\Ω}}\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right){\mathrm{\δS}}_o\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γ}}_o^{\′}{\mathrm{\δS}}_o\mathrm{d\Ω}---(3-4)$

采用分部积分法得：

结合边界条件，可得：

$I^{\′}={\∫\∫}_{\mathrm{\Ω}}\frac{1}{2}\mathrm{\ψ}{S_o}^2\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\λ}}_o{\▿p}_o{\▿S}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γ}}_o^{\′}{S}_o\mathrm{d\Ω}---(3-6)$

根据步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I’取极值，则：

$\frac{\∂\left(I^{\′\left(e\right)}\right)}{\∂\left\{{S_o}^{\left(e\right)}\right\}}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\left\{{S_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_o\frac{\∂\left({p_o}^{\left(e\right)}\right)}{\∂x}{\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_o\frac{\∂\left({p_o}^{\left(e\right)}\right)}{\∂y}{\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\mathrm{\γ}}_o^{\′}{\left[N\right]}^T\mathrm{d\Ω}$

$={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\left\{{S_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_o\left(\right[\frac{\∂N}{\∂x}\left]\right\{{p_o}^{\left(e\right)}\left\}\right){\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_o\left(\right[\frac{\∂N}{\∂y}\left]\right\{{p_o}^{\left(e\right)}\left\}\right){\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\mathrm{\γ}}_o^{\′}{\left[N\right]}^T\mathrm{d\Ω}---(3-7)$

其中：

$\left[H^{\′\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\mathrm{d\Ω}---(3-8)$

最终得到：

IV裂缝系统饱和度有限元单元特性分析

裂缝系统饱和度方程：

$\frac{\∂}{\∂l}\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right)=\mathrm{\φ}\frac{{\∂S}_o}{\∂t}-q_{\mathrm{ov}}^{*}---(4-1)$

其中： $q_{\mathrm{ov}}^{*}=q_{\mathrm{ovnf}}^{*}$

对时间项处理得：

取

移项后得：

根据虚位移原理将式(4-1)化为裂缝系统微分方程的等效积分“弱”形式并化简得：

${\∫}_l\frac{\mathrm{\φ}}{\mathrm{\Δt}}{S}_o{\mathrm{\δS}}_o\mathrm{dl}={\∫}_l\frac{\∂}{\∂l}\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right){\mathrm{\δS}}_o\mathrm{dl}+{\∫}_l{\mathrm{\γ}}_{\mathrm{fo}}^{\′}{\mathrm{\δS}}_o\mathrm{dl}---(4-4)$

采用分部积分法得：

结合边界条件，可得：

${\∫}_l\frac{\mathrm{\φ}}{2\mathrm{\Δt}}{\left(S_o\right)}^2\mathrm{dl}=-{\∫}_l{\mathrm{\λ}}_o\frac{\∂p}{\∂l}\frac{\∂S_o}{\∂l}\mathrm{dl}+{\∫}_l{\mathrm{\γ}}_{\mathrm{fo}}^{\′}{S}_o\mathrm{dl}---(4-6)$

根据步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I_f’取极值，则：

$\frac{\∂\left(I_f^{\′\left(e\right)}\right)}{\∂\left({S_o}^{\left(e\right)}\right)}={\∫}_{l^{\left(e\right)}}\left(\frac{\mathrm{\φ}}{2\mathrm{\Δt}}{\left[N^{\′}\right]}^T\right[N^{\′}\left]\right\{{S_o}^{\left(e\right)}\left\}\right)\mathrm{dl}+{\∫}_{l^{\left(e\right)}}\left({\mathrm{\λ}}_o\frac{\∂\left(p^{\left(e\right)}\right)}{\∂l}{\left[N^{\′}\right]}^T\right)\mathrm{dl}-{\∫}_{l^{\left(e\right)}}\left({\mathrm{\γ}}_{\mathrm{fo}}^{\′}{\left[N^{\′}\right]}^T\right)\mathrm{dl}---(4-7)$

其中：

$\left[H_f^{\′\left(e\right)}\right]={\∫}_{l^{\left(e\right)}}\left(\frac{\mathrm{\φ}}{\mathrm{\Δt}}{\left[N^{\′}\right]}^T\right[N^{\′}\left]\right)\mathrm{dl}---(4-8)$

最终得到：

V形成总特征矩阵步骤：

根据整体研究区域插值单元的划分来组装整体特性矩阵，即把有限个基本单元的单元特征矩阵组装在一起，形成总特征矩阵；即

压力求解的整体方程为：

其中：

饱和度求解的整体方程为：

其中：

所述方法的离散裂缝网络数值模拟的过程中采用了动态分配存储空间技术。

与现有技术相比，本发明的有益效果是：提高了数值模拟的计算速度、增加了计算的收敛性，实现了对缝洞型油藏的科学处理，为找准缝洞型油藏剩余油的分布位置、定量确定油藏的储量丰度、科学合理地开发这类油田提供依据，最终达到提高采收率的目的。

附图说明

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

图1是本发明离散介质的离散化方法示意图。

图2是本发明基岩-裂缝窜流计算示意图。

图3是本发明离散裂缝网络数值模拟有限元单元特性分析过程示意图。

图4是本发明离散裂缝网络数值模拟程序流程图。

图5是本发明实施例中裂缝性油藏不采用离散裂缝网络模型的概念模型的有限元网格。

图6是本发明实施例中裂缝性油藏采用了离散裂缝网络模型的概念模型的有限元网格。

图7是本发明实施例中含水率和时间的关系曲线图。

图8是本发明实施例中计算时间步和收敛性关系图。

图9是本发明实施例中生产和模拟的累计产油曲线。

具体实施方式

一种离散裂缝网络模型数值模拟有限元数值求解方法，所述方法包括以下步骤：

(1)根据离散裂缝网络模型，对缝洞型油藏的基岩和裂缝这两种不同的连续性介质采取不同的离散化方法，即对裂缝采用线单元离散，对基岩采用三角形单元离散，得到数学模型(如图1所示)；

(2)对所述数学模型进行数值模拟，得到离散裂缝网络数值模拟；

(3)分析所述离散裂缝网络数值模拟的有限元单元的特性，对所述数学模型进行求解。

所述方法对裂缝进行了显式的处理，根据地质上所描述的裂缝发育分布情况，在数值模拟中给予尽量符合实际的考虑；在所述数值模拟中把大洞穴看做裂缝来进行处理，油藏模拟假设基岩中的流体和岩石均微可压缩，且压缩系数为常数，裂缝及裂缝中的流体为不可压缩；所述数学模型为：

基岩系统：

${-\mathrm{\φC}}_t\frac{{\∂p}_o}{\∂t}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-(q_{\mathrm{wv}}+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}(q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*}))---(0-1)$

$\mathrm{\φ}\frac{{\∂S}_o}{\∂t}+S_o{\mathrm{\φC}}_{\mathrm{to}}\frac{{\∂p}_o}{\∂t}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}{q}_{\mathrm{ovnf}}^{*}---(0-2)$

S_o+S_w＝1     (0-3)

P_c＝P_o-P_w    (0-4)

裂缝系统：

$\frac{\∂}{\∂l}{\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)}_{\mathrm{nf}}+(q_{\mathrm{wvnf}}^{*}+q_{\mathrm{ovnf}}^{*})=0---(0-5)$

$\frac{\∂}{\∂l}{\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right)}_{\mathrm{nf}}={\mathrm{\φ}}_{\mathrm{nf}}\frac{{\∂S}_{\mathrm{onf}}}{\∂t}-q_{\mathrm{ovnf}}^{*}---(0-6)$

S_onf+S_wnf＝1 (0-7)

基岩-裂缝间的窜流：

$q_{\mathrm{onf}}^{*}={\left({\mathrm{\ξ\λ}}_o(p_o-p)\right)}_{\mathrm{nf}}---(0-8)$

$q_{\mathrm{wnf}}^{*}={\left({\mathrm{\ξ\λ}}_w(p_w-p)\right)}_{\mathrm{nf}}---(0-9)$

其中，

图2说明了基岩和裂缝间窜流的计算方法，即每条裂缝的每个端点处P_mnf＝P_fnf。

所述数学模型可以准确描述非均质性，也可以应用到水润湿和混合润湿介质中。

所述步骤(3)包括以下步骤：

(31)基岩系统压力有限元单元特性分析

基岩系统压力方程：

${-\mathrm{\φC}}_t\frac{{\∂p}_o}{\∂t}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-q_v---(1-1)$

其中， $q_v=q_{\mathrm{wv}}+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}(q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*})$

对时间项的处理：

${-\mathrm{\φC}}_t\frac{p_o-p_o^{(n-1)}}{\mathrm{\Δt}}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-q_v---(1-2)$

取

移项后得：

$-\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{p}_o+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-\mathrm{\γ}---(1-3)$

根据虚位移原理将式(1-1)化为基岩系统微分方程的等效积分“弱”形式并化简得：

${-\∫\∫}_{\mathrm{\Ω}}\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{p}_o{\mathrm{\δp}}_o\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}\▿\·\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\▿p}_o\right){\mathrm{\δp}}_o\mathrm{d\Ω}={\∫\∫}_{\mathrm{\Ω}}\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right){\mathrm{\δp}}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}\mathrm{\γ\δ}{p}_o\mathrm{d\Ω}---(1-4)$

采用分部积分法得：

根据边界条件且边界处在此暂不考虑毛管力，则

$I={\∫\∫}_{\mathrm{\Ω}}\frac{{\mathrm{\φC}}_t}{2\mathrm{\Δt}}{\left(p_o\right)}^2\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}\frac{({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)}{2}\left({\▿p}_o\right)\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\λ}}_w{\▿p}_c{\▿p}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γp}}_o\mathrm{d\Ω}---(1-6)$

按照步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I取极值，则：

$\frac{\∂\left(I^{\left(e\right)}\right)}{\∂\left\{{p_o}^{\left(e\right)}\right\}}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{\left[N\right]}^T\left[N\right]\left\{{p_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂x}\right]}^T\right[\frac{\∂N}{\∂x}\left]\right\{{p_o}^{\left(e\right)}\}+({\mathrm{\λ}}_o+{\mathrm{\λ}}_e){\left[\frac{\∂N}{\∂y}\right]}^T[\frac{\∂N}{\∂y}\left]\right\{{p_o}^{\left(e\right)}\left\}\right)\mathrm{d\Ω}$

$-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_w\frac{\∂\left({p_c}^{\left(e\right)}\right)}{\∂x}{\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_w\frac{\∂\left({p_c}^{\left(e\right)}\right)}{\∂y}{\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\γ}{\left[N\right]}^T\mathrm{d\Ω}---(1-7)$

其中：

$\left[{H_1}^{\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{\left[N\right]}^T\right[N\left]\right)\mathrm{d\Ω}$

$\left[{H_2}^{\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂x}\right]}^T\right[\frac{\∂N}{\∂x}]+({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂y}\right]}^T[\frac{\∂N}{\∂y}\left]\right)\mathrm{d\Ω}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{d\Ω}$

$\left[B\right]=\left[\begin{array}{cccc}\frac{{\∂N}_1}{\∂x}& \frac{{\∂N}_2}{\∂x}& ...& \frac{{\∂N}_n}{\∂x}\\ \frac{{\∂N}_1}{\∂y}& \frac{{\∂N}_2}{\∂y}& ...& \frac{{\∂N}_n}{\∂y}\end{array}\right],\left[D\right]=\left[\begin{array}{cc}({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)& \\ & ({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\end{array}\right]$

$\left[H^{\left(e\right)}\right]=\left[H_1^{\left(e\right)}\right]+\left[{H_2}^{\left(e\right)}\right]---(1-8)$

最终得到：

(32)裂缝系统压力有限元单元特性分析

裂缝系统压力方程：

$\frac{\∂}{\∂l}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)+q_v^{*}=0---(2-1)$

其中， $q_v^{*}=q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*}$

根据虚位移原理将(2-1)化为裂缝系统微分方程的等效积分“弱”形式并化简得：

${\∫}_l\frac{\∂}{\∂l}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)\mathrm{\δpdl}+{\∫}_l{q}_v^{*}\mathrm{\δpdl}=0---(2-2)$

采用分部积分法得：

结合边界条件可得：

$I={\∫}_l\frac{1}{2}({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left(\frac{\∂p}{\∂l}\right)}^2\mathrm{dl}-{\∫}_l{q}_v^{*}\mathrm{pdl}---(2-4)$

按照步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I_f取极值，则：

$\frac{\∂\left({I_f}^{\left(e\right)}\right)}{\∂\left\{p^{\left(e\right)}\right\}}={\∫}_{l^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{{\∂N}^{\′}}{\∂l}\right]}^T\right[\frac{{\∂N}^{\′}}{\∂l}\left]\right\{p^e\left\}\right)\mathrm{dl}-{\∫}_{l^{\left(e\right)}}{q}_v^{*}{\left[N\right]}^T\mathrm{dl}---(2-5)$

其中：

$\left[{H_f}^{\left(e\right)}\right]={\∫}_{l^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N^{\′}}{\∂l}\right]}^T\right[\frac{\∂N^{\′}}{\∂l}\left]\right\{p^{\left(e\right)}\left\}\right)\mathrm{dl}={\∫}_l{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{dl}---(2-6)$

$\left[B\right]=\left[\begin{array}{cc}\frac{{\∂N}_1^{\′}}{\∂l}& \frac{{\∂N}_2^{\′}}{\∂l}\end{array}\right],\left[D\right]=({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)$

最终得到：

(33)基岩系统饱和度有限元单元特性分析

基岩系统饱和度方程：

$\mathrm{\φ}\frac{{\∂S}_o}{\∂t}+S_o\mathrm{\φ}{C}_{\mathrm{to}}\frac{{\∂p}_o}{\∂t}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}^{\′}---(3-1)$

其中： $q_{\mathrm{ov}}^{\′}=q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}{q}_{\mathrm{ovnf}}^{*}$

对时间项处理后得：

$\mathrm{\φ}\frac{S_o-S_o^{(n-1)}}{\mathrm{\Δt}}+S_o{\mathrm{\φC}}_{\mathrm{to}}\frac{p_o-p_o^{(n-1)}}{\mathrm{\Δt}}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}^{\′}---(3-2)$

取

移项后得：

${\mathrm{\ψS}}_o=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+{\mathrm{\γ}}_o^{\′}---(3-3)$

根据虚位移原理将(3-1)化为基岩系统微分方程的等效积分“弱”形式并化简得：

${\∫\∫}_{\mathrm{\Ω}}{\mathrm{\ψS}}_o{\mathrm{\δS}}_o\mathrm{d\Ω}={\∫\∫}_{\mathrm{\Ω}}\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right){\mathrm{\δS}}_o\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γ}}_o^{\′}{\mathrm{\δS}}_o\mathrm{d\Ω}---(3-4)$

采用分部积分法得：

结合边界条件，可得：

$I^{\′}={\∫\∫}_{\mathrm{\Ω}}\frac{1}{2}\mathrm{\ψ}{S_o}^2\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\λ}}_o{\▿p}_o{\▿S}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γ}}_o^{\′}{S}_o\mathrm{d\Ω}---(3-6)$

-按照步骤D中所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I’取极值，则：

$\frac{\∂\left(I^{\′\left(e\right)}\right)}{\∂\left\{{S_o}^{\left(e\right)}\right\}}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\left\{{S_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_o\frac{\∂\left({p_o}^{\left(e\right)}\right)}{\∂x}{\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_o\frac{\∂\left({p_o}^{\left(e\right)}\right)}{\∂y}{\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\mathrm{\γ}}_o^{\′}{\left[N\right]}^T\mathrm{d\Ω}$

$={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\left\{{S_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_o\left(\right[\frac{\∂N}{\∂x}\left]\right\{{p_o}^{\left(e\right)}\left\}\right){\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_o\left(\right[\frac{\∂N}{\∂y}\left]\right\{{p_o}^{\left(e\right)}\left\}\right){\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\mathrm{\γ}}_o^{\′}{\left[N\right]}^T\mathrm{d\Ω}---(3-7)$

其中：

$\left[H^{\′\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\mathrm{d\Ω}---(3-8)$

最终得到：

(34)裂缝系统饱和度有限元单元特性分析

裂缝系统饱和度方程：

$\frac{\∂}{\∂l}\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right)=\mathrm{\φ}\frac{{\∂S}_o}{\∂t}-q_{\mathrm{ov}}^{*}---(4-1)$

其中： $q_{\mathrm{ov}}^{*}=q_{\mathrm{ovnf}}^{*}$

对时间项处理得：

取

移项后得：

根据虚位移原理将式(4-1)化为裂缝系统微分方程的等效积分“弱”形式并化简得：

${\∫}_l\frac{\mathrm{\φ}}{\mathrm{\Δt}}{S}_o{\mathrm{\δS}}_o\mathrm{dl}={\∫}_l\frac{\∂}{\∂l}\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right){\mathrm{\δS}}_o\mathrm{dl}+{\∫}_l{\mathrm{\γ}}_{\mathrm{fo}}^{\′}{\mathrm{\δS}}_o\mathrm{dl}---(4-4)$

采用分部积分法得：

结合边界条件，可得：

${\∫}_l\frac{\mathrm{\φ}}{2\mathrm{\Δt}}{\left(S_o\right)}^2\mathrm{dl}=-{\∫}_l{\mathrm{\λ}}_o\frac{\∂p}{\∂l}\frac{\∂S_o}{\∂l}\mathrm{dl}+{\∫}_l{\mathrm{\γ}}_{\mathrm{fo}}^{\′}{S}_o\mathrm{dl}---(4-6)$

按照步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I_f’取极值，则：

$\frac{\∂\left(I_f^{\′\left(e\right)}\right)}{\∂\left({S_o}^{\left(e\right)}\right)}={\∫}_{l^{\left(e\right)}}\left(\frac{\mathrm{\φ}}{2\mathrm{\Δt}}{\left[N^{\′}\right]}^T\right[N^{\′}\left]\right\{{S_o}^{\left(e\right)}\left\}\right)\mathrm{dl}+{\∫}_{l^{\left(e\right)}}\left({\mathrm{\λ}}_o\frac{\∂\left(p^{\left(e\right)}\right)}{\∂l}{\left[N^{\′}\right]}^T\right)\mathrm{dl}-{\∫}_{l^{\left(e\right)}}\left({\mathrm{\γ}}_{\mathrm{fo}}^{\′}{\left[N^{\′}\right]}^T\right)\mathrm{dl}---(4-7)$

其中：

$\left[H_f^{\′\left(e\right)}\right]={\∫}_{l^{\left(e\right)}}\left(\frac{\mathrm{\φ}}{\mathrm{\Δt}}{\left[N^{\′}\right]}^T\right[N^{\′}\left]\right)\mathrm{dl}---(4-8)$

最终得到：

(35)根据整体研究区域插值单元的划分来组装整体特性矩阵，即把有限个基本单元的单元特征矩阵组装在一起，形成总特征矩阵。

有限元是用有限个基本单元来逼近整体模型，步骤(31)至(34)对每个插值单元形成单元特性矩阵和列阵，这样就可以根据整体研究区域插值单元的划分来组装整体特性矩阵，即把有限个基本单元的单元特征矩阵组装在一起，形成总特征矩阵。即

压力求解的整体方程为：

其中：

饱和度求解的整体方程为：

其中：

上述离散裂缝网络数值模拟有限元单元特性分析过程如图3所示，图中对裂缝性油藏数值模拟有限元方法采用示意图的形式进行了说明。图4是本发明离散裂缝网络数值模拟程序流程图，图中说明了利用该方法进行数值模拟的流程，其中结点信息、地层属性数据、相对渗透率和PVT数据以文件的形式存储，易于管理和修改。由于油藏数值模拟的数据量较大，为了合理安排和节省存储空间，在程序编制中采用了动态分配存储空间技术。

1)模型验证

为了验证方法的正确性，根据图4所示的流程编制了该方法的数值模拟程序，并进行了模拟计算。裂缝的发育往往和断层及构造的发育密切相关，主要产生在断层附近和背斜顶部，并且研究表明在深于3000m的地层中主要发育高角度的裂缝，所以实施例中假设裂缝均为垂直裂缝，因此在油藏中就存在裂缝及其以外的区域，分别称之为裂缝和基岩。在数学表达中，对裂缝采用一维实体进行表示，而对于裂缝以外的基岩区域则仍采用二维平面进行表示。本实施例采用的概念模型中设置一50m×50m二维油藏，并布置一口注水井和一口生产井，左下角注水，右上角采出，注入量为150m。图5所示的是实施例中不采用离散裂缝网络模型时离散的有限元网格，图6所示的是实施例中采用图1所示的方法进行离散后得到的裂缝性油藏概念模型的有限元网格；从图5到图6可以看出采用离散裂缝网络模型后剖分的有限元网格数大大减小，从而使计算量也大大减小；图7给出了实施例中含水率和时间的关系曲线图。

本发明准确的进行了一简单概念模型的模拟，其结果是注入水进入油藏后先是在基岩中推挤油向前缓慢推进，当注入水遇到裂缝后，便沿着裂缝方向向前推进，其水驱油速度大大加快；其主要原因是由于基岩的渗透率很小，裂缝的渗透率远远大于基岩的渗透率，油藏注入水后水挤压着油必然主要从渗透率高的裂缝中流到生产井，最后得到的图7所示含水率曲线也说明生产井在投产初期含水率很小，且增加很小，当生产一段时间后，含水迅速增加，其也是由于裂缝其主导作用，从而理论分析和实际结果的一致性说明了该方法的正确性。

图8所示为时间步和迭代次数的关系图，由于每一时间步都需要经过很多次迭代才能够得到精确值，也就是说每一时间步的计算结果只有达到给定精度迭代才进入下一时间步的计算，直至计算结束。迭代次数越多，收敛性越差，从图可以看出采用传统有限元方法进行计算的迭代次数明显比采用本专利所述的方法少得多，收敛性也就好的多，经对比，收敛性比普通方法高5倍左右，迭代次数的减少导致计算速度也有了幅度的提高。

2)裂缝性油藏的模拟

某一缝洞型油藏，溶洞近似处理为裂缝，其渗透率为1D，基岩系统(包含溶蚀孔洞)渗透率为0.1mD，基岩系统孔隙度为0.05，裂缝系统的孔隙度为0.001；油藏温度124.1摄氏度，油藏原始压力59MPa，油藏饱和压力20.2MPa，地饱压差为38.55MPa；地层原油粘度21.703mPa咣地面原油密度0.96g/cm³，原油体积系数1.46。地层水密度1.147g/cm³，原始气油比60m³/m³；油水界面-5610m，原始石油地质储量3069？0⁴t。

通过采用本专利所述的方法进行缝洞型油藏的数值模拟，模拟结果如图9所示，该图为某油田模拟后的累计产量曲线与实际累积产量曲线，计算结果表明该方法不仅比较准确的描述了缝洞型油藏中流体的流动，还揭示了缝洞型油藏裂缝为流体流动通道等规律，为下一步生产预测奠定了基础。

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

Claims (4)

1.一种分析离散裂缝性油藏流体流动的方法，其特征在于，所述方法包括以下步骤：

A缝洞型油藏物探数据输入步骤，

B对缝洞型油藏的基岩和裂缝两种不同的连续性介质采取不同的离散化方法，即对裂缝采用线单元离散，对基岩采用三角形单元离散：即对每个待测区域V离散成E个单元，每个单元有N个结点；

C输入结点信息、地层属，性数据、相对渗透率和PVT数据；

D分析得到基岩系统和裂缝系统压力的有限元单元特性函数及其饱和度的有限元单元特性函数；

E形成各单元特性矩阵和列阵步骤：形成基岩系统压力函数单元特征矩阵：矩阵[H^(e)]和列阵

以及裂缝系统压力函数单元特征矩阵：矩阵

和列阵

)，基岩系统饱和度函数单元特征矩阵：矩阵[H′^(e)]和列阵

和裂缝系统饱和度函数单元特征矩阵：矩阵[H′_f ^(e)]和列阵)；

F组装有限个基本单元的单元特征矩阵，形成总特征矩阵，即形成压力整体特征矩阵：矩阵和列阵

和饱和度整体特征矩阵：矩阵和列阵

G根据各物理量求取步骤；

(1)对压力方程求解，即根据初始条件或上一步计算结果采用高斯多重积分方法对压力方程中各项积分的进行计算，形成压力方程的单元特性矩阵和列阵，然后组合压力方程的整体特性矩阵和列阵，最后采取预处理共轭梯度法对压力方程组进行求解；

(2)对饱和度方程求解，即根据这一步求出的压力和初始条件或上一步计算结果对饱和度方程中各项积分进行计算，形成饱和度方程的单元特性矩阵或列阵，然后组合饱和度方程的整体特性矩阵和列阵，最后对饱和度方程组进行求解；

H判断步骤：

(1)计算线性方程组中总装矩阵的条件数，如条件数数量级大于16，则认为不收敛，否则就认为收敛，从而进行下一时间步的计算；

(2)如果计算时间大约设定时间计算结束，否则进人下一时间步的计算；

I输出显示步骤：

输出显示压力场、饱和度场等的图像及动画以及收敛性曲线和生产指数曲线。

2.根据权利要求1所述的方法，其特征在于，

所述方法对裂缝进行了单独处理，根据地质上所描述的裂缝发育分布情况，在数值模拟中；在所述数值模拟中把大洞穴作为裂缝来进行处理，油藏模拟假设基岩中的流体和岩石均微可压缩，且压缩系数为常数，裂缝及裂缝中的流体为不可压缩；所述数学模型为：

基岩系统：

${-\mathrm{\φC}}_t\frac{{\∂p}_o}{\∂t}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-(q_{\mathrm{wv}}+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}(q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*}))---(0-1)$

$\mathrm{\φ}\frac{{\∂S}_o}{\∂t}+S_o{\mathrm{\φC}}_{\mathrm{to}}\frac{{\∂p}_o}{\∂t}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}{q}_{\mathrm{ovnf}}^{*}---(0-2)$

S_o+S_w＝1           (0-3)

P_c＝P_o-P_w          (0-4)

裂缝系统：

$\frac{\∂}{\∂l}{\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)}_{\mathrm{nf}}+(q_{\mathrm{wvnf}}^{*}+q_{\mathrm{ovnf}}^{*})=0---(0-5)$

$\frac{\∂}{\∂l}{\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right)}_{\mathrm{nf}}={\mathrm{\φ}}_{\mathrm{nf}}\frac{{\∂S}_{\mathrm{onf}}}{\∂t}-q_{\mathrm{ovnf}}^{*}---(0-6)$

S_onf+S_wnf＝1       (0-7)

基岩-裂缝间的窜流：

$q_{\mathrm{onf}}^{*}={\left({\mathrm{\ξ\λ}}_o(p_o-p)\right)}_{\mathrm{nf}}---(0-8)$

$q_{\mathrm{wnf}}^{*}={\left({\mathrm{\ξ\λ}}_w(p_w-p)\right)}_{\mathrm{nf}}---(0-9)$

其中， $\mathrm{\ξ}=\frac{\mathrm{lgh}}{d}.$

3.根据权利要求1或2所述的方法，其特征在于，所述步骤D包括以下步骤：

I基岩系统压力有限元单元特性分析步骤：

基岩系统压力方程：

${-\mathrm{\φC}}_t\frac{{\∂p}_o}{\∂t}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-q_v---(1-1)$

其中， $q_v=q_{\mathrm{wv}}+q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}(q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*})$

对时间项的处理：

${-\mathrm{\φC}}_t\frac{p_o-p_o^{(n-1)}}{\mathrm{\Δt}}+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-q_v---(1-2)$

取

移项后得：

$-\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{p}_o+\▿\·\left(({\mathrm{\λ}}_w+{\mathrm{\λ}}_o){\▿p}_o\right)=\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right)-\mathrm{\γ}---(1-3)$

根据虚位移原理将式(1-1)化为基岩系统微分方程的等效积分“弱”形式并化简得：

${-\∫\∫}_{\mathrm{\Ω}}\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{p}_o{\mathrm{\δp}}_o\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}\▿\·\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\▿p}_o\right){\mathrm{\δp}}_o\mathrm{d\Ω}={\∫\∫}_{\mathrm{\Ω}}\▿\·\left({\mathrm{\λ}}_w{\▿p}_c\right){\mathrm{\δp}}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}\mathrm{\γ\δ}{p}_o\mathrm{d\Ω}---(1-4)$

采用分部积分法得：

根据边界条件且边界处在此暂不考虑毛管力，则

$I={\∫\∫}_{\mathrm{\Ω}}\frac{{\mathrm{\φC}}_t}{2\mathrm{\Δt}}{\left(p_o\right)}^2\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}\frac{({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)}{2}\left({\▿p}_o\right)\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\λ}}_w{\▿p}_c{\▿p}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γp}}_o\mathrm{d\Ω}---(1-6)$

所述步骤D中，用来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I取极值，则：

$\frac{\∂\left(I^{\left(e\right)}\right)}{\∂\left\{{p_o}^{\left(e\right)}\right\}}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{\left[N\right]}^T\left[N\right]\left\{{p_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂x}\right]}^T\right[\frac{\∂N}{\∂x}\left]\right\{{p_o}^{\left(e\right)}\}+({\mathrm{\λ}}_o+{\mathrm{\λ}}_e){\left[\frac{\∂N}{\∂y}\right]}^T[\frac{\∂N}{\∂y}\left]\right\{{p_o}^{\left(e\right)}\left\}\right)\mathrm{d\Ω}$

$-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_w\frac{\∂\left({p_c}^{\left(e\right)}\right)}{\∂x}{\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_w\frac{\∂\left({p_c}^{\left(e\right)}\right)}{\∂y}{\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\γ}{\left[N\right]}^T\mathrm{d\Ω}---(1-7)$

其中：

$\left[{H_1}^{\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(\frac{{\mathrm{\φC}}_t}{\mathrm{\Δt}}{\left[N\right]}^T\right[N\left]\right)\mathrm{d\Ω}$

$\left[{H_2}^{\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂x}\right]}^T\right[\frac{\∂N}{\∂x}]+({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N}{\∂y}\right]}^T[\frac{\∂N}{\∂y}\left]\right)\mathrm{d\Ω}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{d\Ω}$

$\left[B\right]=\left[\begin{array}{cccc}\frac{{\∂N}_1}{\∂x}& \frac{{\∂N}_2}{\∂x}& ...& \frac{{\∂N}_n}{\∂x}\\ \frac{{\∂N}_1}{\∂y}& \frac{{\∂N}_2}{\∂y}& ...& \frac{{\∂N}_n}{\∂y}\end{array}\right],\left[D\right]=\left[\begin{array}{cc}({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)& \\ & ({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\end{array}\right]$

$\left[H^{\left(e\right)}\right]=\left[H_1^{\left(e\right)}\right]+\left[{H_2}^{\left(e\right)}\right]---(1-8)$

最终得到：

II裂缝系统压力有限元单元特性分析步骤：

裂缝系统压力方程：

$\frac{\∂}{\∂l}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)+q_v^{*}=0---(2-1)$

其中， $q_v^{*}=q_{\mathrm{ovnf}}^{*}+q_{\mathrm{wvnf}}^{*}$

根据虚位移原理将(2-1)化为裂缝系统微分方程的等效积分“弱”形式并化简得：

${\∫}_l\frac{\∂}{\∂l}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\frac{\∂p}{\∂l}\right)\mathrm{\δpdl}+{\∫}_l{q}_v^{*}\mathrm{\δpdl}=0---(2-2)$

采用分部积分法得：

结合边界条件可得：

$I={\∫}_l\frac{1}{2}({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left(\frac{\∂p}{\∂l}\right)}^2\mathrm{dl}-{\∫}_l{q}_v^{*}\mathrm{pdl}---(2-4)$

根据所述步骤D中，分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I_f取极值，则：

$\frac{\∂\left({I_f}^{\left(e\right)}\right)}{\∂\left\{p^{\left(e\right)}\right\}}={\∫}_{l^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{{\∂N}^{\′}}{\∂l}\right]}^T\right[\frac{{\∂N}^{\′}}{\∂l}\left]\right\{p^e\left\}\right)\mathrm{dl}-{\∫}_{l^{\left(e\right)}}{q}_v^{*}{\left[N\right]}^T\mathrm{dl}---(2-5)$

其中：

$\left[{H_f}^{\left(e\right)}\right]={\∫}_{l^{\left(e\right)}}\left(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w){\left[\frac{\∂N^{\′}}{\∂l}\right]}^T\right[\frac{\∂N^{\′}}{\∂l}\left]\right\{p^{\left(e\right)}\left\}\right)\mathrm{dl}={\∫}_l{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{dl}---(2-6)$

$\left[B\right]=\left[\begin{array}{cc}\frac{{\∂N}_1^{\′}}{\∂l}& \frac{{\∂N}_2^{\′}}{\∂l}\end{array}\right],\left[D\right]=({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)$

最终得到：

III基岩系统饱和度有限元单元特性分析步骤：

基岩系统饱和度方程：

$\mathrm{\φ}\frac{{\∂S}_o}{\∂t}+S_o\mathrm{\φ}{C}_{\mathrm{to}}\frac{{\∂p}_o}{\∂t}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}^{\′}---(3-1)$

其中： $q_{\mathrm{ov}}^{\′}=q_{\mathrm{ov}}-{\mathrm{\δ}}_{\mathrm{nf}}{q}_{\mathrm{ovnf}}^{*}$

对时间项处理后得：

$\mathrm{\φ}\frac{S_o-S_o^{(n-1)}}{\mathrm{\Δt}}+S_o{\mathrm{\φC}}_{\mathrm{to}}\frac{p_o-p_o^{(n-1)}}{\mathrm{\Δt}}=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+q_{\mathrm{ov}}^{\′}---(3-2)$

取移项后得：

${\mathrm{\ψS}}_o=\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right)+{\mathrm{\γ}}_o^{\′}---(3-3)$

根据虚位移原理将(3-1)化为基岩系统微分方程的等效积分“弱”形式并化简得：

${\∫\∫}_{\mathrm{\Ω}}{\mathrm{\ψS}}_o{\mathrm{\δS}}_o\mathrm{d\Ω}={\∫\∫}_{\mathrm{\Ω}}\▿\·\left({\mathrm{\λ}}_o{\▿p}_o\right){\mathrm{\δS}}_o\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γ}}_o^{\′}{\mathrm{\δS}}_o\mathrm{d\Ω}---(3-4)$

采用分部积分法得：

结合边界条件，可得：

$I^{\′}={\∫\∫}_{\mathrm{\Ω}}\frac{1}{2}\mathrm{\ψ}{S_o}^2\mathrm{d\Ω}+{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\λ}}_o{\▿p}_o{\▿S}_o\mathrm{d\Ω}-{\∫\∫}_{\mathrm{\Ω}}{\mathrm{\γ}}_o^{\′}{S}_o\mathrm{d\Ω}---(3-6)$

根据步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I’取极值，则：

$\frac{\∂\left(I^{\′\left(e\right)}\right)}{\∂\left\{{S_o}^{\left(e\right)}\right\}}={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\left\{{S_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_o\frac{\∂\left({p_o}^{\left(e\right)}\right)}{\∂x}{\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_o\frac{\∂\left({p_o}^{\left(e\right)}\right)}{\∂y}{\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\mathrm{\γ}}_o^{\′}{\left[N\right]}^T\mathrm{d\Ω}$

$={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\left\{{S_o}^{\left(e\right)}\right\}\mathrm{d\Ω}+{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}({\mathrm{\λ}}_o\left(\right[\frac{\∂N}{\∂x}\left]\right\{{p_o}^{\left(e\right)}\left\}\right){\left[\frac{\∂N}{\∂x}\right]}^T+{\mathrm{\λ}}_o\left(\right[\frac{\∂N}{\∂y}\left]\right\{{p_o}^{\left(e\right)}\left\}\right){\left[\frac{\∂N}{\∂y}\right]}^T)\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}{\mathrm{\γ}}_o^{\′}{\left[N\right]}^T\mathrm{d\Ω}---(3-7)$

其中：

$\left[H^{\′\left(e\right)}\right]={\∫\∫}_{{\mathrm{\Ω}}^{\left(e\right)}}\mathrm{\ψ}{\left[N\right]}^T\left[N\right]\mathrm{d\Ω}---(3-8)$

最终得到：

IV  裂缝系统饱和度有限元单元特性分析步骤：

裂缝系统饱和度方程：

$\frac{\∂}{\∂l}\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right)=\mathrm{\φ}\frac{{\∂S}_o}{\∂t}-q_{\mathrm{ov}}^{*}---(4-1)$

其中： $q_{\mathrm{ov}}^{*}=q_{\mathrm{ovnf}}^{*}$

对时间项处理得：

取

移项后得：

根据虚位移原理将式(4-1)化为裂缝系统微分方程的等效积分“弱”形式并化简得：

${\∫}_l\frac{\mathrm{\φ}}{\mathrm{\Δt}}{S}_o{\mathrm{\δS}}_o\mathrm{dl}={\∫}_l\frac{\∂}{\∂l}\left({\mathrm{\λ}}_o\frac{\∂p}{\∂l}\right){\mathrm{\δS}}_o\mathrm{dl}+{\∫}_l{\mathrm{\γ}}_{\mathrm{fo}}^{\′}{\mathrm{\δS}}_o\mathrm{dl}---(4-4)$

采用分部积分法得：

结合边界条件，可得：

${\∫}_l\frac{\mathrm{\φ}}{2\mathrm{\Δt}}{\left(S_o\right)}^2\mathrm{dl}=-{\∫}_l{\mathrm{\λ}}_o\frac{\∂p}{\∂l}\frac{\∂S_o}{\∂l}\mathrm{dl}+{\∫}_l{\mathrm{\γ}}_{\mathrm{fo}}^{\′}{S}_o\mathrm{dl}---(4-6)$

根据步骤D中，所述来分析基岩系统压力的有限元特性函数，对其求解区域进行离散，并对泛函I，’取极值，则：

$\frac{\∂\left(I_f^{\′\left(e\right)}\right)}{\∂\left({S_o}^{\left(e\right)}\right)}={\∫}_{l^{\left(e\right)}}\left(\frac{\mathrm{\φ}}{2\mathrm{\Δt}}{\left[N^{\′}\right]}^T\right[N^{\′}\left]\right\{{S_o}^{\left(e\right)}\left\}\right)\mathrm{dl}+{\∫}_{l^{\left(e\right)}}\left({\mathrm{\λ}}_o\frac{\∂\left(p^{\left(e\right)}\right)}{\∂l}{\left[N^{\′}\right]}^T\right)\mathrm{dl}-{\∫}_{l^{\left(e\right)}}\left({\mathrm{\γ}}_{\mathrm{fo}}^{\′}{\left[N^{\′}\right]}^T\right)\mathrm{dl}---(4-7)$

其中：

$\left[H_f^{\′\left(e\right)}\right]={\∫}_{l^{\left(e\right)}}\left(\frac{\mathrm{\φ}}{\mathrm{\Δt}}{\left[N^{\′}\right]}^T\right[N^{\′}\left]\right)\mathrm{dl}---(4-8)$

最终得到：

V形成总特征特征矩阵步骤：

根据整体研究区域插值单元的划分来组装整体特性矩阵，即把有限个基本单元的单元特征矩阵组装在一起，形成总特征矩阵；即

压力求解的整体方程为：

其中：

饱和度求解的整体方程为：

其中：

4.根据权利要求1所述的方法，其特征在于，所述方法的离散裂缝网络数值模拟的过程中采用了动态分配存储空间技术。

Images