CN112163379B - 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法 - Google Patents

一种基于孔隙网络模型的非稳态两相相对渗透率计算方法 Download PDF

Info

Publication number
CN112163379B
CN112163379B CN202010992354.4A CN202010992354A CN112163379B CN 112163379 B CN112163379 B CN 112163379B CN 202010992354 A CN202010992354 A CN 202010992354A CN 112163379 B CN112163379 B CN 112163379B
Authority
CN
China
Prior art keywords
phase
relative permeability
unsteady
node
network model
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010992354.4A
Other languages
English (en)
Other versions
CN112163379A (zh
Inventor
唐雁冰
杨鑫
李闽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Petroleum University
Original Assignee
Southwest Petroleum University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Petroleum University filed Critical Southwest Petroleum University
Priority to CN202010992354.4A priority Critical patent/CN112163379B/zh
Publication of CN112163379A publication Critical patent/CN112163379A/zh
Application granted granted Critical
Publication of CN112163379B publication Critical patent/CN112163379B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于孔隙网络模型的非稳态两相相对渗透率计算方法,包括以下步骤:在孔隙网络模型的基础上,通过的非稳态渗流方程计算流体流量及压力;根据计算得到的流体流量及压力,结合相对渗透率计算公式,计算得到所述非稳态两相相对渗透率;根据驱替过程中不同含水饱和度下的两相相对渗透率数据,绘制两相相对渗透率曲线。本发明能够计算得到更符合储层物性特征的相对渗透率,并以此绘制两相相对渗透率曲线,为油气藏开发提供准确的技术指导。

Description

一种基于孔隙网络模型的非稳态两相相对渗透率计算方法
技术领域
本发明涉及油气田开发技术领域,特别涉及一种基于孔隙网络模型的非稳态两相相对渗透率计算方法。
背景技术
油气资源是现全世界使用最多与最重要的能源之一,如何高效合理开发油气资源是每一个油藏工程师面对的难题。两相(油水两相、气水两相)相对渗透率曲线是油气田开发工程中非常重要的基础数据,获取该类数据的方法常常有实验法与数值模拟法,每种方法都各有优缺点。实验方法包括常用的稳态法与非稳态法,非稳态法实验枯燥且难度大,耗时长,非稳态法(常用的JBN实验方法)虽然耗时短,但数据常常难以处理;常用的以黑油模型为基础数值模拟法往往较少考虑储层岩石的物性特征。基于孔隙网络模型的数值模拟方法一般假设岩石内部孔喉具有一定的形状,且可以通过岩心分析实验准确获取这类资料,模型复用性高,相比实验,可变条件(流速、压力等)宽泛,且稳定的孔隙网络模拟方法可无限重复使用,具有强大的经济效益。
目前,国内外许多研究人员通常采用连续介质理论研究多孔介质的多相渗流,但由于粘滞力、毛管力在孔隙尺度上具有非连续性,且考虑裂缝介质后其渗流规律如何变化尚不明确;实际流体在储层内部通常具有可压缩性(尤其是气体),这与流体能瞬时从入口传到出口的稳态渗流理论相悖;实际生产情况中,多相流体的渗流规律极其复杂,常规的渗流理论不能够准确的指导油气藏的开发,预测其生产动态,而目前对非稳态两相渗流模拟研究方法还存在一定的局限性,常规的商业数值模拟类软件能大致模拟两相渗流过程,但是绘制的相对渗透率曲线往往与实际产生的准确性有所偏颇,这极大地限制了油气资源的开采。
发明内容
针对上述问题,本发明旨在提供一种基于孔隙网络模型的非稳态两相相对渗透率计算方法。
本发明的技术方案如下:
一种基于孔隙网络模型的非稳态两相相对渗透率计算方法,包括以下步骤:
通过的非稳态渗流方程计算流体流量及压力,所述非稳态渗流方程为:
Figure GDA0003380905100000011
式中:
Figure GDA0003380905100000012
为哈密尔顿算子;G为传导率;p为压力,MPa;φ0为初始压力下的孔隙度,无量纲;Ct为综合压缩系数,无量纲;t为时间,s;
根据计算得到的流体流量及压力,结合相对渗透率计算公式,计算得到所述非稳态两相相对渗透率。
作为优选,若考虑汇源项,则所述非稳态渗流方程为:
Figure GDA0003380905100000021
式中:q为汇或源的流量,m3/s。
作为优选,若考虑网络管束,则非稳态渗流方程为:
Figure GDA0003380905100000022
式中:Δ'为孔隙网络模型的所有方向,当所述孔隙网络模型为SC网络模型时为6个方向,当所述孔隙网络模型为BCC网络模型时为8个方向,当所述孔隙网络模型为FCC网络模型时为12个方向;Δ”p为每个方向的压差,MPa;Vb为网格体积,cm3;Δt为时间差,s。
作为优选,计算油水两相的相对渗透率时,两相流中水相和油相的传导率分别为:
Figure GDA0003380905100000023
Figure GDA0003380905100000024
式中:gw、go分别为水相传导率和油相传导率;rij为节点i和节点j之间的孔道半径,cm;μw、μo分别为水相粘度、油相粘度,Pa·s;lij为节点i和节点j之间的孔道长度,cm;Bo为油相体积系数,无量纲;
计算气水两相的相对渗透率时,两相流中水相的传导率通过式(4)进行计算,两相流中气体的传导率为:
Figure GDA0003380905100000025
Figure GDA0003380905100000026
式中:gg为气相传导率;Bg为气相体积系数,无量纲;μg为气相粘度,Pa·s;psc为地面大气压,MPa;Zsc为地面气体偏差因子,无量纲;Tsc为地面温度,℃;Z为地下气体偏差因子,无量纲;T表示地下温度,℃;<p>表示地下气体压力,MPa;<p>=(pi+pj)/2,pi和pj为管束两端节点i和节点j的压力,MPa。
作为优选,所述相对渗透率计算公式为:
Figure GDA0003380905100000031
Figure GDA0003380905100000032
Figure GDA0003380905100000033
式中:kro、krg、krw1、krw2分别为油相、气相、油水两相中的水相、气水两相中的水相的相对渗透率,%;qo(t)、qg(t)、qw1(t)、qw2(t)分别为t时刻时油相、气相、油水两相中的水相、气水两相中的水相的流量,m3/s;α1、α2均为与岩心末端及驱替时间末端相关的经验常数;kro_min、krg_min分别为油相、气相的最小渗透率,mD;qomax、qgmax分别为整个流动过程中的油相、气相的最大流量,m3/s。
作为优选,本发明所述的非稳态两相相对渗透率计算方法还包括根据驱替过程中不同含水饱和度下的两相相对渗透率数据,绘制两相相对渗透率曲线的步骤。
作为优选,绘制两相相对渗透率曲线的具体步骤如下:
首先,根据所述非稳态渗流方程求得岩心末端部分两相初始状态的流量及压力,根据所述初始状态的流量及压力,结合相对渗透率计算公式,计算得到两相初始状态的相对渗透率;
然后,进行两相驱替,根据两相流动时的流动控制方程及最小时间步长,获得驱替过程中每个时间步长的含水饱和度及传导率,将所述传导率代入所述非稳态渗流方程,迭代计算获得整个驱替过程中的流量及压力变化数据,以此计算得到驱替过程中的两相相对渗透率数据;
最后,以含水饱和度为横坐标,以相对渗透率为纵坐标,绘制得到所述两相相对渗透率曲线。
作为优选,所述流动控制方程为:
Figure GDA0003380905100000034
μeff=μwlijxijo(lij-lijxij),μeff=Bgμglijxijw(lij-lijxij) (12)
Figure GDA0003380905100000035
式中:qij为两相的流量,cm3/s;gij为两相的传导率;pi和pj为管束两端节点i和节点j的压力,MPa;pcij为圆管内毛管压力,MPa;μeff为单个管道中两相有效粘度,Pa·s;μw、μo、μg分别为水相粘度、油相粘度、气相粘度,Pa·s;lij为节点i和节点j之间的孔道长度,cm;xij为与凹液面位置有关的无量纲数,0≤xij≤1,即凹液面所在位置横坐标除以孔隙吼道的长度;Bg为气相体积系数,无量纲;γ为界面张力,mN/m;θij为两相界面中润湿流体的润湿角,°;rij为节点i和节点j之间的孔道半径,cm。
作为优选,所述最小时间步长的计算方法如下:
对于任一与孔隙相连通的喉道i,假定其中被驱替相部分完全被驱替相流体充满所需时间为ti,则有:
Figure GDA0003380905100000041
式中:Vi,nw为被驱替相流体体积,m3;qi,w为被驱替相流体流量,m3/s;m为总迭代次数;
则驱替相流体充填所有孔隙网络的节点单元所需的最小时间步长tmin可表示为:
Figure GDA0003380905100000042
下一次的时间步长更新,将已经整管占据的管束将其设置为1,存在两相界面的管束设置为0,其时间t计算为:
Figure GDA0003380905100000043
式中:lij为节点i和节点j之间的孔道长度,cm;xij为与凹液面位置有关的无量纲数,0≤xij≤1,即凹液面所在位置横坐标除以孔隙吼道的长度;v为有效粘度下的线速度;
然后再统计时间步长更新这一过程中的最小时间步长,作为下一次最小时间步长。
作为优选,所述含水饱和度的计算方法如下:
Figure GDA0003380905100000044
Vw=∑1Vi (18)
式中:Sw为含水饱和度,%;Vw为孔隙在时间步长中被水相全部占据的管束体积之和,cm3;Vp为总孔隙体积,其值为每个连通管束的体积和,cm3;Vi为第i个被水相全部占据的管束体积,cm3
与现有技术相比,本发明具有如下优点:
本发明能够通过结合两相非稳态数理方程与孔隙网络模型,计算得到更符合储层物性特征的相对渗透率,并以此绘制两相相对渗透率曲线,为油气藏开发提供准确的技术指导。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为孔隙网络的二维局部节点示意图;
图2为岩心末端部分流量示意图;
图3为水驱油时圆管内油水界面示意图;
图4为水驱气时圆管内气水界面示意图;
图5为本发明实施例两相相对渗透率曲线结果示意图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。需要说明的是,在不冲突的情况下,本申请中的实施例及实施例中的技术特征可以相互结合。除非另外定义,本发明公开使用的技术术语或者科学术语应当为本公开所属领域内具有一般技能的人士所理解的通常意义。本发明公开使用的“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同,而不排除其他元件或者物件。
一种基于孔隙网络模型的非稳态两相相对渗透率计算方法,包括以下步骤:
S1:通过的非稳态渗流方程计算流体流量及压力,所述非稳态渗流方程为:
Figure GDA0003380905100000051
式中:
Figure GDA0003380905100000052
为哈密尔顿算子;G为传导率;p为压力,MPa;φ0为初始压力下的孔隙度,无量纲;Ct为综合压缩系数,无量纲;t为时间,s;
若考虑汇源项,则所述非稳态渗流方程为:
Figure GDA0003380905100000053
式中:q为汇或源的流量,m3/s。
若考虑网络管束,则利用非稳态渗流积分方程计算流体流量及压力,所述非稳态渗流积分方程为:
Figure GDA0003380905100000054
式中:Δ'为孔隙网络模型的所有方向,当所述孔隙网络模型为SC网络模型时为6个方向,当所述孔隙网络模型为BCC网络模型时为8个方向,当所述孔隙网络模型为FCC网络模型时为12个方向;Δ”p为每个方向的压差,MPa;Vb为网格体积,cm3;Δt为时间差,s。
公式(1)-(3)中的传导率计算方法如下:
计算油水两相的相对渗透率时,两相流中水相和油相的传导率分别为:
Figure GDA0003380905100000061
Figure GDA0003380905100000062
式中:gw、go分别为水相传导率和油相传导率;rij为节点i和节点j之间的孔道半径,cm;μw、μo分别为水相粘度、油相粘度,Pa·s;lij为节点i和节点j之间的孔道长度,cm;Bo为油相体积系数,无量纲;
计算气水两相的相对渗透率时,两相流中水相的传导率通过式(4)进行计算,两相流中气体的传导率为:
Figure GDA0003380905100000063
Figure GDA0003380905100000064
式中:gg为气相传导率;Bg为气相体积系数,无量纲;μg为气相粘度,Pa·s;psc为地面大气压,MPa;Zsc为地面气体偏差因子,无量纲;Tsc为地面温度,℃;Z为地下气体偏差因子,无量纲;T表示地下温度,℃;<p>表示地下气体压力,MPa;<p>=(pi+pj)/2,pi和pj为管束两端节点i和节点j的压力,MPa。
下面以公式(3)为例,介绍所述非稳态渗流方程的求解方法:
对公式(3)右侧时间进行差分:
Figure GDA0003380905100000065
式中:n为当前时刻状态下的各参数值;Δt为时间步长,s;i、j、k分别表示三维空间的三个方向;pn+1、pn分别为下一时刻和当前时刻的压力值,MPa。
基于质量守恒定律,以图1所示的二维模型对公式(3)左边进行求解,每个节点流入流出流量体积之和为0,每个节点的流量如下:
节点1:
Figure GDA0003380905100000071
节点2:
Figure GDA0003380905100000072
节点3:
Figure GDA0003380905100000073
节点4:
Figure GDA0003380905100000074
节点5:
Figure GDA0003380905100000075
节点6:
Figure GDA0003380905100000076
公式(23)-(28)中下角标代表节点或喉道,例如p1代表节点1处的压力,g10代表节点0和节点1之间的喉道中流体的传导率。
以节点2为例,考虑隐式时间推进过程,变形得到:
Figure GDA0003380905100000077
以隐式方法处理所有节点得到:
Figure GDA0003380905100000081
将上述各节点的守恒方程写成矩阵的形式可以得到:
Figure GDA0003380905100000082
得到如[A]n+1[P]n+1=[B]n形式的矩阵方程,其中:
Figure GDA0003380905100000083
Figure GDA0003380905100000084
Figure GDA0003380905100000085
采用梯度下降法求解矩阵f(P)=AP-B:
f(P)=AP-B (32)
f′(P)=A (33)
Figure GDA0003380905100000091
其中,Pi+1为当前时刻压力场迭代求解过程中下一步压力场的迭代值(由初始值迭代而来),当迭代过程中残差值f(P)误差小于1E-7(10的负7次方)时,认为收敛,此时得到的压力场即为当前时刻压力场。
通过计算所述矩阵即可得到这一时刻所有节点的流量和压力,这些流量和压力会形成流量场和压力场(矩阵形式)。
S2:根据计算得到的流体流量及压力,结合相对渗透率计算公式,计算得到所述非稳态两相相对渗透率。所述相对渗透率计算公式为:
Figure GDA0003380905100000092
Figure GDA0003380905100000093
Figure GDA0003380905100000094
式中:kro、krg、krw1、krw2分别为油相、气相、油水两相中的水相、气水两相中的水相的相对渗透率,%;qo(t)、qg(t)、qw1(t)、qw2(t)分别为t时刻时油相、气相、油水两相中的水相、气水两相中的水相的流量,m3/s;α1、α2均为与岩心末端及驱替时间末端相关的经验常数;kro_min、krg_min分别为油相、气相的最小渗透率,mD;qomax、qgmax分别为整个流动过程中的油相、气相的最大流量,m3/s。
上述的相对渗透率计算公式适用于各种渗透率的岩石,且能获取精确的相对渗透率。公式推导过程如下:
对于常规达西流动,两相相对渗透率可以表示为:
Figure GDA0003380905100000095
式中:krw、kro分别为水相相对渗透率和油相相对渗透率,%;qw、qs、qo分别为两相流中水相流量、相同条件下(压差,流量,界面移动情况均相等)流过岩石内部单相流量、两相流中油相流量,cm3/s;ΔPs、ΔPw、ΔPo为岩心两端稳定压差,MPa;μw、μs、μo分别为水相粘度、单相流体粘度、油相粘度,Pa·s。
根据式(35),只要获取单向流动和两相流动时对应状态下各相流体的压差和流量,则可以计算相对渗透率,但在这个过程中,很难保证同模型中单相流和两相流动时的情况一样(主要是压差),因此较难获取相同条件下的单相渗透率,故该公式计算误差较大。
在中低渗岩石中,如图2所示,任意时刻qin不等于qout,不满足达西渗流无法直接用进出口流量计算相渗,取模型或岩石的末端部分,即近似稳态渗流状态下的流动区域内,得到qout=qout',计算该区域内含水饱和度,以及该区域内的压力、粘度、分流量等数据进行相渗计算。
以水驱油为例,在驱替过程中,岩心末端部分流体流动状态分为三个阶段:第一阶段,单相油流;第二阶段,油水两相流;第三阶段,单相水流,整个计算过程忽略毛管力影响,饱和度计量也只取所述的岩心末端部分。
在第一阶段与第三阶段,qs等于不同时刻t时的油相流量qo与水相流量qw之和(此时管束内流体视为单相流),油相相对渗透率计算公式为:
Figure GDA0003380905100000101
此时:
Figure GDA0003380905100000102
在两相流阶段,可以近似认为油相流量与水相流量之和占总流量的100%,油相流量占比40%~60%,水相流量占比40%~60%,此时单管内两相界面中水相入口处至水驱前缘总长lw占管束总长度40%~60%,lo也占管束总长度40%~60%,由泊肃叶定律,两相流动时,每个存在两相界面的毛管中,单一相流量为:
Figure GDA0003380905100000103
因此:
Figure GDA0003380905100000104
所以油相相对渗透率为:
Figure GDA0003380905100000105
又因为:
Figure GDA0003380905100000111
所以:
Figure GDA0003380905100000112
同理可得:
Figure GDA0003380905100000113
考虑到所取位置越靠近岩心末端,流动越接近达西流,令:
Figure GDA0003380905100000116
得到:
Figure GDA0003380905100000114
油驱水或水驱气或气驱水等理论相同,综上得到式(8)-式(10)所示的能够适用于各种储层且计算结果精确的相对渗透率计算公式。
在另一个具体的实施例中,所述非稳态两相相对渗透率计算方法还包括步骤S3:根据驱替过程中不同含水饱和度下的两相相对渗透率数据,绘制两相相对渗透率曲线,具体步骤如下:
S31:根据所述非稳态渗流方程求得两相初始状态的流量及压力,根据所述初始状态的流量及压力,结合相对渗透率计算公式,计算得到两相初始状态的相对渗透率。
S32:进行两相驱替,根据两相流动时的流动控制方程及最小时间步长,获得驱替过程中每个时间步长(所述时间步长大于等于所述最小时间步长)的含水饱和度及传导率,将所述传导率代入所述非稳态渗流方程,迭代计算获得整个驱替过程中的流量及压力变化数据,以此计算得到驱替过程中的两相相对渗透率数据。
两相驱替过程中,两相界面如图3和图4所示,所述流动控制方程为:
Figure GDA0003380905100000115
μeff=μwlijxijo(lij-lijxij),μeff=Bgμglijxijw(lij-lijxij) (12)
Figure GDA0003380905100000121
式中:qij为两相的流量,cm3/s;gij为两相的传导率;pi和pj为管束两端节点i和节点j的压力,MPa;pcij为圆管内毛管压力,MPa;μeff为单个管道中两相有效粘度,Pa·s;μw、μo、μg分别为水相粘度、油相粘度、气相粘度,Pa·s;lij为节点i和节点j之间的孔道长度,cm;xij为与凹液面位置有关的无量纲数,0≤xij≤1,即凹液面所在位置横坐标除以孔隙吼道的长度;Bg为气相体积系数,无量纲;γ为界面张力,mN/m;θij为两相界面中润湿流体的润湿角,°;rij为节点i和节点j之间的孔道半径,cm。
所述最小时间步长的计算方法如下:
对于任一与孔隙相连通的喉道i,假定其中被驱替相部分完全被驱替相流体充满所需时间为ti,则有:
Figure GDA0003380905100000122
式中:Vi,nw为被驱替相流体体积,m3;qi,w为被驱替相流体流量,m3/s;m为总迭代次数;
则驱替相流体充填所有孔隙网络的节点单元所需的最小时间步长tmin可表示为:
Figure GDA0003380905100000123
下一次的时间步长更新,将已经整管占据的管束将其设置为1,存在两相界面的管束设置为0,其时间t计算为:
Figure GDA0003380905100000124
式中:lij为节点i和节点j之间的孔道长度,cm;xij为与凹液面位置有关的无量纲数,0≤xij≤1,即凹液面所在位置横坐标除以孔隙吼道的长度;v为有效粘度下的线速度;
然后再统计时间步长更新这一过程中的最小时间步长,作为下一次最小时间步长。
所述含水饱和度的计算方法如下:
Figure GDA0003380905100000125
Vw=∑1Vi (18)
式中:Sw为含水饱和度,%;Vw为孔隙在时间步长中被水相全部占据的管束体积之和,cm3;Vp为总孔隙体积,其值为每个连通管束的体积和,cm3;Vi为第i个被水相全部占据的管束体积,cm3
S33:以含水饱和度为横坐标,以相对渗透率为纵坐标,绘制得到所述两相相对渗透率曲线。
在一个具体的实施例中,进行气驱水、水驱油、水驱气三种驱替实验,采用本发明所述的非稳态两相相对渗透率计算方法计算驱替过程中的相对渗透率,并绘制两相相对渗透率曲线,结果如图5所示。
以上所述,仅是本发明的较佳实施例而已,并非对本发明作任何形式上的限制,虽然本发明已以较佳实施例揭露如上,然而并非用以限定本发明,任何熟悉本专业的技术人员,在不脱离本发明技术方案范围内,当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例,但凡是未脱离本发明技术方案的内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰,均仍属于本发明技术方案的范围内。

Claims (9)

1.一种基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,包括以下步骤:
通过的非稳态渗流方程计算流体流量及压力,所述非稳态渗流方程为:
Figure FDA0003380905090000011
式中:
Figure FDA0003380905090000012
为哈密尔顿算子;G为传导率;p为压力,MPa;φ0为初始压力下的孔隙度,无量纲;Ct为综合压缩系数,无量纲;t为时间,s;
根据计算得到的流体流量及压力,结合相对渗透率计算公式,计算得到所述非稳态两相相对渗透率;所述相对渗透率计算公式为:
Figure FDA0003380905090000013
Figure FDA0003380905090000014
Figure FDA0003380905090000015
式中:kro、krg、krw1、krw2分别为油相、气相、油水两相中的水相、气水两相中的水相的相对渗透率,%;qo(t)、qg(t)、qw1(t)、qw2(t)分别为t时刻时油相、气相、油水两相中的水相、气水两相中的水相的流量,m3/s;α1、α2均为与岩心末端及驱替时间末端相关的经验常数;kro_min、krg_min分别为油相、气相的最小渗透率,mD;qomax、qgmax分别为整个流动过程中的油相、气相的最大流量,m3/s。
2.根据权利要求1所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,若考虑汇源项,则所述非稳态渗流方程为:
Figure FDA0003380905090000016
式中:q为汇或源的流量,m3/s。
3.根据权利要求2所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,若考虑网络管束,则所述非稳态渗流方程为:
Figure FDA0003380905090000017
式中:Δ'为孔隙网络模型的所有方向,当所述孔隙网络模型为SC网络模型时为6个方向,当所述孔隙网络模型为BCC网络模型时为8个方向,当所述孔隙网络模型为FCC网络模型时为12个方向;Δ”p为每个方向的压差,MPa;Vb为网格体积,cm3;Δt为时间差,s。
4.根据权利要求1-3中任意一项所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,计算油水两相的相对渗透率时,两相流中水相和油相的传导率分别为:
Figure FDA0003380905090000021
Figure FDA0003380905090000022
式中:gw、go分别为水相传导率和油相传导率;rij为节点i和节点j之间的孔道半径,cm;μw、μo分别为水相粘度、油相粘度,Pa·s;lij为节点i和节点j之间的孔道长度,cm;Bo为油相体积系数,无量纲;
计算气水两相的相对渗透率时,两相流中水相的传导率通过式(4)进行计算,两相流中气相的传导率为:
Figure FDA0003380905090000023
Figure FDA0003380905090000024
式中:gg为气相传导率;Bg为气相体积系数,无量纲;μg为气相粘度,Pa·s;psc为地面大气压,MPa;Zsc为地面气体偏差因子,无量纲;Tsc为地面温度,℃;Z为地下气体偏差因子,无量纲;T表示地下温度,℃;<p>表示地下气体压力,MPa;<p>=(pi+pj)/2,pi和pj为管束两端节点i和节点j的压力,MPa。
5.根据权利要求1-3中任意一项所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,还包括根据驱替过程中不同含水饱和度下的两相相对渗透率数据,绘制两相相对渗透率曲线的步骤。
6.根据权利要求5所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,绘制两相相对渗透率曲线的具体步骤如下:
选取岩心末端部分,根据所述非稳态渗流方程求得该区域内两相初始状态的流量及压力,根据所述初始状态的流量及压力,结合相对渗透率计算公式,计算得到两相初始状态的相对渗透率;
然后,进行两相驱替,根据两相流动时的流动控制方程及最小时间步长,获得驱替过程中每个时间步长的含水饱和度及传导率,将所述传导率代入所述非稳态渗流方程,迭代计算获得整个驱替过程中的流量及压力变化数据,以此计算得到驱替过程中的两相相对渗透率数据;
最后,以含水饱和度为横坐标,以相对渗透率为纵坐标,绘制得到所述两相相对渗透率曲线。
7.根据权利要求6所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,所述流动控制方程为:
Figure FDA0003380905090000031
μeff=μwlijxijo(lij-lijxij),μeff=Bgμglijxijw(lij-lijxij) (12)
Figure FDA0003380905090000032
式中:qij为两相的流量,cm3/s;gij为两相的传导率;pi和pj为管束两端节点i和节点j的压力,MPa;pcij为圆管内毛管压力,MPa;μeff为单个管道中两相有效粘度,Pa·s;μw、μo、μg分别为水相粘度、油相粘度、气相粘度,Pa·s;lij为节点i和节点j之间的孔道长度,cm;xij为与凹液面位置有关的无量纲数,0≤xij≤1,即凹液面所在位置横坐标除以孔隙吼道的长度;Bg为气相体积系数,无量纲;γ为界面张力,mN/m;θij为两相界面中润湿流体的润湿角,°;rij为节点i和节点j之间的孔道半径,cm。
8.根据权利要求6所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,所述最小时间步长的计算方法如下:
对于任一与孔隙相连通的喉道i,假定其中被驱替相部分完全被驱替相流体充满所需时间为ti,则有:
Figure FDA0003380905090000033
式中:Vi,nw为被驱替相流体体积,m3;qi,w为被驱替相流体流量,m3/s;m为总迭代次数;
则驱替相流体充填所有孔隙网络的节点单元所需的最小时间步长tmin可表示为:
Figure FDA0003380905090000034
下一次的时间步长更新,将已经整管占据的管束将其设置为1,存在两相界面的管束设置为0,其时间t计算为:
Figure FDA0003380905090000041
式中:lij为节点i和节点j之间的孔道长度,cm;xij为与凹液面位置有关的无量纲数,0≤xij≤1,即凹液面所在位置横坐标除以孔隙吼道的长度;v为有效粘度下的线速度;
然后再统计时间步长更新这一过程中的最小时间步长,作为下一次最小时间步长。
9.根据权利要求6所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法,其特征在于,所述含水饱和度的计算方法如下:
Figure FDA0003380905090000042
Vw=∑1Vi (18)
式中:Sw为含水饱和度,%;Vw为孔隙在时间步长中被水相全部占据的管束体积之和,cm3;Vp为总孔隙体积,其值为每个连通管束的体积和,cm3;Vi为第i个被水相全部占据的管束体积,cm3
CN202010992354.4A 2020-09-21 2020-09-21 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法 Active CN112163379B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010992354.4A CN112163379B (zh) 2020-09-21 2020-09-21 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010992354.4A CN112163379B (zh) 2020-09-21 2020-09-21 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法

Publications (2)

Publication Number Publication Date
CN112163379A CN112163379A (zh) 2021-01-01
CN112163379B true CN112163379B (zh) 2022-02-15

Family

ID=73862799

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010992354.4A Active CN112163379B (zh) 2020-09-21 2020-09-21 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法

Country Status (1)

Country Link
CN (1) CN112163379B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114021502B (zh) * 2021-11-10 2023-09-05 长江大学 基于迎风gfdm的多孔介质油水两相流计算方法
CN113868930B (zh) * 2021-11-10 2023-09-01 长江大学 基于广义有限差分方法的各向异性储层渗流模拟方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106547938A (zh) * 2015-11-09 2017-03-29 中国地质大学(北京) 裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法
WO2018022956A1 (en) * 2016-07-28 2018-02-01 Board Of Regents, The University Of Texas System Systems and methods for determining unsteady-state two-phase relative permeability
CN109670220A (zh) * 2018-12-05 2019-04-23 西南石油大学 一种基于非结构网格的水平井气水两相数值模拟方法
CN110298048A (zh) * 2018-03-22 2019-10-01 中国石油化工股份有限公司 一种考虑界面相的超临界co2-凝析气多相渗流模拟方法
CN110308495A (zh) * 2018-12-19 2019-10-08 中国石油大学(北京) 地下储层单元线流动数据处理方法及装置
CN110334365A (zh) * 2019-02-27 2019-10-15 中国石油大学(北京) 一种非均质压裂后储层流动数值模拟方法及系统
CN110598167A (zh) * 2019-10-11 2019-12-20 中国石油化工股份有限公司 低渗透油藏油水相对渗透率实验数据的处理方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104533370B (zh) * 2014-11-06 2017-03-15 中国石油大学(北京) 压裂水平井油藏、裂缝、井筒全耦合模拟方法
CN107346518B (zh) * 2016-05-06 2020-10-13 中国石油化工股份有限公司 致密低渗透油藏油水两相流最大渗流阻力梯度的获取方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106547938A (zh) * 2015-11-09 2017-03-29 中国地质大学(北京) 裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法
WO2018022956A1 (en) * 2016-07-28 2018-02-01 Board Of Regents, The University Of Texas System Systems and methods for determining unsteady-state two-phase relative permeability
CN110298048A (zh) * 2018-03-22 2019-10-01 中国石油化工股份有限公司 一种考虑界面相的超临界co2-凝析气多相渗流模拟方法
CN109670220A (zh) * 2018-12-05 2019-04-23 西南石油大学 一种基于非结构网格的水平井气水两相数值模拟方法
CN110308495A (zh) * 2018-12-19 2019-10-08 中国石油大学(北京) 地下储层单元线流动数据处理方法及装置
CN110334365A (zh) * 2019-02-27 2019-10-15 中国石油大学(北京) 一种非均质压裂后储层流动数值模拟方法及系统
CN110598167A (zh) * 2019-10-11 2019-12-20 中国石油化工股份有限公司 低渗透油藏油水相对渗透率实验数据的处理方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Experimental Study on Gas Production from Methane Hydrate-Bearing Sand by Hot-Water Cyclic Injection;Yang, Xin 等;《ENERGY & FUELS》;20101130;第24卷(第11期);第5912-5920页 *
Model test for water inrush caused by karst caves filled with confined water in tunnels;Yang, Weimin 等;《ARABIAN JOURNAL OF GEOSCIENCES》;20191231;第12卷(第24期);第1-11页 *
低渗煤层高压注氮驱替强化抽采技术及应用研究;杨鑫 等;《中国博士学位论文全文数据库 工程科技Ⅰ辑》;20200915(第9期);B021-13 *
低渗透油气藏压裂水平井分形渗流理论研究;谭晓华;《中国博士学位论文全文数据库 工程科技Ⅰ辑》;20170315(第3期);B019-22 *
低渗透油藏非线性渗流基础理论与数值模拟研究;杨仁锋;《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅰ辑》;20110415(第4期);B019-50 *
孔洞型碳酸盐岩储集层中洞对电阻率的影响;何家欢 等;《石油勘探与开发》;20200630;第47卷(第3期);第490-498页 *

Also Published As

Publication number Publication date
CN112163379A (zh) 2021-01-01

Similar Documents

Publication Publication Date Title
CN110334431B (zh) 一种低渗透致密气藏单井控制储量计算及剩余气分析方法
CN112163379B (zh) 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法
Hornung et al. Diffusion models for fractured media
CN104847314B (zh) 高温高压油气直井单相流射孔完井参数优化方法
Gharibshahi et al. Application of CFD to evaluate the pore morphology effect on nanofluid flooding for enhanced oil recovery
CN107045671A (zh) 产水气井积液风险预测方法
Liu et al. Critical parameters of the Jamin effect in a capillary tube with a contracted cross section
CN105021506A (zh) 一种基于孔隙网络模型的三相相对渗透率的计算方法
Sharma et al. Experiments and analysis of multiscale viscous fingering during forced imbibition
CN111425191B (zh) 一种致密储层气-水相对渗透率计算方法及装置
CN113836695B (zh) 一种基于无网格连接元的油藏数值模拟方法
CN104989385B (zh) 基于表皮系数计算的高温高压油气直井射孔参数优化方法
Knudsen et al. Relation between pressure and fractional flow in two-phase flow in porous media
CN113836767A (zh) 一种页岩油储层压后关井时间的优化方法
CN114580100A (zh) 压裂水平井全井筒压力计算方法、设备和计算机可读储存介质
CN110580656B (zh) 一种水平井筒地下极限携液流量预测方法
Lee et al. Numerical modeling of groundwater flow into a radial collector well with horizontal arms
CN109632579B (zh) 一种页岩黏土矿物强制自吸量预测方法
CN115711842A (zh) 考虑水的瞬态渗流的纳米多孔介质页岩渗透率计算方法
CN113255247B (zh) 一种高含水期油藏多尺度数值模拟方法
CN109933951B (zh) 致密油藏体积压裂水平井多尺度、多机理耦合渗流模型的建立方法
CN107133373B (zh) 一种页岩气藏、井筒及地面管网的耦合模拟方法
CN115110945B (zh) 一种低渗透油藏相对渗透率预测方法
CN116050297A (zh) 孔隙网络中多级渗流路径与渗透率等效的提取方法
CN105160056B (zh) 高温高压油气直井两相流射孔完井参数与产能优化方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant