CN117113873B - 一种多相流体地层渗流的数值模拟方法及应用 - Google Patents

一种多相流体地层渗流的数值模拟方法及应用 Download PDF

Info

Publication number
CN117113873B
CN117113873B CN202311027224.7A CN202311027224A CN117113873B CN 117113873 B CN117113873 B CN 117113873B CN 202311027224 A CN202311027224 A CN 202311027224A CN 117113873 B CN117113873 B CN 117113873B
Authority
CN
China
Prior art keywords
fluid
phase
multiphase fluid
multiphase
seepage
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
CN202311027224.7A
Other languages
English (en)
Other versions
CN117113873A (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 CN202311027224.7A priority Critical patent/CN117113873B/zh
Publication of CN117113873A publication Critical patent/CN117113873A/zh
Application granted granted Critical
Publication of CN117113873B publication Critical patent/CN117113873B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/02Agriculture; Fishing; Forestry; Mining
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • 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
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Business, Economics & Management (AREA)
  • Computational Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Fluid Mechanics (AREA)
  • Databases & Information Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Agronomy & Crop Science (AREA)
  • Animal Husbandry (AREA)
  • Marine Sciences & Fisheries (AREA)
  • Mining & Mineral Resources (AREA)
  • Computer Graphics (AREA)
  • Health & Medical Sciences (AREA)
  • Economics (AREA)
  • General Health & Medical Sciences (AREA)
  • Human Resources & Organizations (AREA)
  • Marketing (AREA)
  • Primary Health Care (AREA)
  • Strategic Management (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请提供了一种多相流体地层渗流的数值模拟方法及应用,涉及数值模拟技术领域,该方法包括:S1:获取岩心样品的孔隙度三维数据体和渗透率三维数据体,构建高精度数字岩心模型;S2:通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型;S3:根据高精度数字岩心模型和渗流模拟数学模型建立三维流体力学模型,根据三维流体力学模型求解多相流体的压力场;其中,在建立渗流模拟数学模型中,引入多相流体的压缩系数、混合流体的混合粘度和混合密度。本申请可以改善现有技术中多相流体地层渗流模拟结果与实际渗流情况相距甚远的技术问题。

Description

一种多相流体地层渗流的数值模拟方法及应用
技术领域
本申请涉及数值模拟技术领域,尤其是涉及一种多相流体地层渗流的数值模拟方法及应用。
背景技术
在过去的多相渗流模拟中,经典的VOF(Volume of Fluid)方法已经被广泛应用。该方法基于对流体在空间中的分布进行体积分数的描述,能够有效地模拟多相流体的界面行为。然而,传统的VOF方法存在一些限制,例如对于复杂的多孔介质结构和非均质介质的处理相对困难。
在传统的流体力学模型中,常常假设液体和气体是不可压缩的,即密度保持不变。然而,这种假设在某些情况下存在一定的局限性,因此其流体的渗流模拟结果常与实际渗流情况相距甚远,影响了模拟结果的有效性,难以为油气开发、地质工程等需要考虑渗流情况的项目提供有用的数据支持,尤其是多相流体渗流模拟中。
发明内容
本申请的目的在于提供一种多相流体地层渗流的数值模拟方法及应用,以改善现有技术中多相流体地层渗流模拟结果与实际渗流情况相距甚远的技术问题,可以更为准确的模拟并描述多相流体的相互作用、界面行为和动态演化过程。
本申请的另一目的在于提供一种多相流体地层渗流的数值模拟方法的应用。
为了上述目的,本申请提供以下技术方案:
第一方面,本申请提供一种多相流体地层渗流的数值模拟方法,包括:
S1:获取岩心样品的孔隙度三维数据体和渗透率三维数据体,构建高精度数字岩心模型;
S2:通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型;
S3:根据高精度数字岩心模型和渗流模拟数学模型建立三维流体力学模型,根据所述三维流体力学模型求解所述多相流体的压力场;
其中,在通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型中,引入所述多相流体的压缩系数、混合流体的混合粘度和混合密度。
进一步地,在本申请的一些实施例中,所述渗流模拟数学模型根据多相流体动量守恒方程、多相流体质量守恒方程以及多相流体相方程构建。
进一步地,在本申请的一些实施例中,其中多相流体质量守恒方程为:
dA为微元体的面积,单位为m2;dV为微元体的体积(dV=dxdydz,x、y、z为微元体的三维坐标的三个方向),单位为m3;m为微元体流体的质量,其单位为kg;t为时间,单位为s;qogw为源汇相,单位为kg/s;此时vm为混合流体的速度,其单位为m/s;ρm为混合流体的密度,其单位为kg/m3,其表达式如下:
ρm=Soρo+Sgρg+Swρw
其中,So、Sg和Sw分别为油、气、水三相的饱和度;ρo、ρg、ρw为油、气、水三相流体的密度,单位为kg/m3
进一步地,在本申请的一些实施例中,多相流体动量守恒方程为:
其中,φ为岩石的孔隙度,孔隙度无量纲;ρm为混合流体的密度,单位为kg/m3;k为单相流体的渗透率,单位为m2为压力梯度,单位为Pa/m;pc为三相毛管压力,单位为Pa;μm为混合流体的粘度,表达式如下:
μm=Soμo+Sgμg+Swμw
其中,μo、μg、μw分别为油气水三相的粘度,粘度单位为Pa·s;
进一步地,在本申请的一些实施例中,多相流体相方程为:
其中,qo、qg、qw为油、气、水的源汇相,单位为kg/s。
进一步地,在本申请的一些实施例中,所述三相毛管压力pc,根据多相流体的优势相,分别通过下式获得:
当油相为优势相,则三相毛管压力pc为:
当气相为优势相,则三相毛管压力pc为:
当水相为优势相,则三相毛管压力pc为:
其中pcog、pcow、pcwg分别为油气、油水以及气水直接的毛管压力,其表达式如下:
其中,σog、σow、σwg分别为油气、油水、气水之间的界面张力,单位为N/m;θog、θow、θwg分别为油气、油水、气水之间的润湿接触角,无量纲;r为喉道半径,单位为m。
进一步地,在本申请的一些实施例中,在步骤S3中,
利用快速迭代方法或多重网格方法求解所述多相流体的压力场。进一步地,在本申请的一些实施例中,所述多相流体的压力场求解包括:
求解多相流体质量方程中的第一项:
其中,Co、Cg、Cw分别为油、气、水三相的压缩系数,压缩系数无量纲;t为当前时刻,t+Δt为下一时刻,单位均为s;
求解多相流体质量方程中的第二项:
其中,zc为相邻网格数,该参数无量纲;A为相邻网格的面积,单位为m2,i为网格i,j为与网格i相邻的网格j;Gm为传导率;
联立多相流体质量方程、多相流体质量方程中的第一项的求解公式和第二项的求解公式,得到:
将当前时刻t和下一时刻t+Δt的不同参数代入上式,得到N个在不同时刻的上式,形成上式的矩阵方程At+ΔtXt+Δt=Bt,其中:
其中,At+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵(N为模型的网格数量),Xt+Δt和Bt是长度为N的两个向量,Xt+Δt为下一时刻t+Δt的压力场向量,Bt为与当前时刻t压力场和边界条件有关的向量;
求解每一时刻的矩阵方程,得到所述多相流体的压力场。
进一步地,在本申请的一些实施例中,所述高精度数字岩心模型的建立包括:
利用岩心样品和扫描技术,获取所述岩心样品的孔隙度三维数据体和渗透率三维数据体;
建立所述岩心样品的三维几何模型,根据所述孔隙度三维数据体和渗透率三维数据体赋值所述三维几何模型,得到所述岩心样品的高精度数字岩心模型;
其中,所述孔隙度三维数据体和渗透率三维数据体的获得包括以下步骤:
S11:获取待模拟储层的岩心样品,并利用所述岩心样本获取所述岩心样品的实际孔隙度φ和实际渗透率k;
S12:利用扫描技术获取所述岩心样本在不同断面上的二维图像,获得二维图像数据;
S13:根据多个所述断面上的二维图像数据,获得岩心三维数据体;
S14:根据岩心的实际孔隙度φ与二维图像数据的关系,获得所述岩心样品的孔隙度转化系数;根据所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系,获取所述岩心样品的渗透率转化系数;
S15:在所述岩心三维数据体上应用所述孔隙度转化系数和所述渗透率转化系数,以及实际孔隙度与二维图像数据的关系、实际渗透率与实际孔隙度及二维图像数据的关系,得到岩心孔隙度三维数据体和岩心渗透率三维数据体。
第二方面,本申请还提供了一种多相流体地层渗流的数值模拟方法在油田开发、地下水管理、地质工程和环境保护技术领域中的应用。
本申请提供一种多相流体地层渗流的数值模拟方法及应用,在VOF模型建立的渗流模拟数学模型中引入可压缩性项,修正传统的流体力学方程,使其能够考虑流体在压缩或膨胀过程中的体积变化;例如,在高速气体流动中,考虑到气体的可压缩性可以更准确地模拟激波、压力波和射流的行为;在液体喷射问题中,液体的可压缩性可以影响喷射流的速度和喷雾的形成过程。同时还考虑了多相流体因为优势相的不同而导致流体的密度和粘度变化的问题,引入了混合粘度和混合密度,进而建立更为可靠的CFD模型,可以更准确地捕捉到多相流体或混合流体中不同组分之间的相互作用和影响;该CFD模型可以更精确地预测流体的速度场、压力分布、涡旋结构等流动特性,从而提高模拟结果的准确性和可靠性这样可以更精确地描述流体的运动和相互作用。此外,引入混合密度和混合粘度的概念还可以改善对于流体界面和相变过程的建模。在多相流体中,流体界面的位置和形态对于流动行为和传热过程有着重要影响。通过考虑不同组分的密度和粘度差异,可以更准确地描述流体界面的动态变化和相变现象,如气液界面的蒸发和凝结、液滴的形成和破裂等。
本申请提供的数值模拟方法可以使流体力学方程模型更加准确地描述物质的行为和性质。这对于模拟和分析涉及液体和气体的复杂流动和相互作用具有重要意义,为各种工程和科学领域提供了更可靠的数值模拟工具。同时,该数值模拟方法可以改善传统CFD模型对于多相流体或混合流体的建模和预测能力。这对于涉及复杂流体***的工程和科学研究具有重要意义,提高了对于多相流体行为的理解和优化设计的能力。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图示出的结构获得其他的附图。
图1为本申请提供的基于VOF的多相流体渗流过程模拟方法的流程示意图。
具体实施方式
下面将结合实施例对本申请方案进行清楚、完整地描述,显然,所描述的实施例是本申请一部分实施例,而不是全部的实施例。基于本申请中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本申请保护的范围。
发明人在本申请中提出了一种多相流体地层渗流的数值模拟方法,参阅图1,包括:
S1:获取岩心样品的孔隙度三维数据体和渗透率三维数据体,构建高精度数字岩心模型;
S2:通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型;
S3:根据高精度数字岩心模型和渗流模拟数学模型建立三维流体力学模型,根据所述三维流体力学模型求解所述多相流体的压力场;
其中,在通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型中,引入所述多相流体的压缩系数、混合流体的混合粘度和混合密度。
在本申请中,所述渗流模拟数学模型根据多相流体动量守恒方程、多相流体质量守恒方程以及多相流体相方程构建。
其中,其中多相流体质量守恒方程为:
dA为微元体的面积,单位为m2;dV为微元体的体积(dV=dxdydz,x、y、z为微元体的三维坐标的三个方向),单位为m3;m为微元体流体的质量,其单位为kg;t为时间,单位为s;qogw为源汇相,单位为kg/s;此时vm为混合流体的速度,其单位为m/s;ρm为混合流体的密度,其单位为kg/m3,其表达式如下:
ρm=Soρo+Sgρg+Swρw
其中,So、Sg和Sw分别为油、气、水三相的饱和度;ρo、ρg、ρw为油、气、水三相流体的密度,单位为kg/m3
多相流体动量守恒方程为:
其中,φ为岩石的孔隙度,孔隙度无量纲;ρm为混合流体的密度,单位为kg/m3;k为单相流体的渗透率,单位为m2为压力梯度,单位为Pa/m;pc为三相毛管压力,单位为Pa;μm为混合流体的粘度,表达式如下:
μm=Soμo+Sgμg+Swμw
其中,μo、μg、μw分别为油气水三相的粘度,粘度单位为Pa·s。
通过多相流体动量守恒方程,可以计算得到多相流体的速度,即可得到所有网格点的多相流体的速度,即可得到多相流体的速度场,用于计算后续其他如压力等参数。本申请中,多相流体的速度场在三维流体力学模型中通过多相流体动量守恒方程计算得到,其可以确定流体的流动方向,了解流体在三维空间中的流动路径;此外还可以计算流体流动速率,提供流体在各个点的流动速率信息,了解流体在多孔介质中的运动状态,以及进一步预测其未来的流动行为。通过对速度场的分析,还可以评估流体的动态行为,如流体的湍动情况、流体是否存在旋转等。这有助于我们更好地理解流体的物理性质,以及流体在复杂条件下的行为。速度场信息也是计算其他重要物理量,如动能、动量等的基础。
多相流体相方程为:
其中,qo、qg、qw为油、气、水的源汇相,单位为kg/s,So、Sg和Sw分别为油相、气相、水相的饱和度;ρo、ρg、ρw分别为油相、气相、水相流体的密度,kg/m3;vo、vg、vw分别为油相、气相、水相流体的速度,m/s。
其具体建立过程如下:
建立岩石样品的三维几何模型,其中通过网格法将岩石样品的三维几何模型分为若干个网格单元,每一网格单元即为1个微元体,该微元体的三维坐标即可表示为(x,y,z)。
因此当流体为单相流体时,岩石在渗流过程中,根据牛顿第二定律,有如下表达式:
其中,Fpres、Fgrav、Fvisc、Fcap分别为流体微元受到周围流体的净力、重力、粘滞力以及界面力(毛管压力Pc产生),其中压力作为动力,粘滞力、惯性力为阻力,而重力和界面力既有可能作为动力,也有可能作为阻力,所有力的单位均为N。单相流体没有界面力,此时Fc=0,m为微元体流体的质量,其单位为kg;a为加速度,其单位为m/s2。
以x方向为例,其中流体微元受到周围流体净受力Fpres可以表示为如下公式:
其中p为压力,单位为Pa;dA为微元体的面积,单位为m2;dV为微元体的体积(dV=dxdydz),单位为m3;忽略流体的重力项。
此时粘滞力作为阻力,单位为N;粘滞力(Fvisc)可以表示为如下公式:
μ为单相流体的粘度,单位为Pa·s;kx为x方向的绝对渗透率,单位为m2;vx为x方向的达西渗流速度,单位为m/s;ax为x方向上的加速度,其单位为m/s2
对于x,y,z三个方向,可以用如下的通式来表达:
其中,v为达西渗流速度,单位为m/s;a加速度,其单位为m/s2
其中压力梯度可以用如下公式来表示:
对于公式(5),两边约去dxdydz三项,可以进一步推导为:
其中为达西渗流速度,单位为m/s;而/>为真实渗流速度,单位为m/s;其中两者有如下的速度关系:
其中,真实渗流速度的物质导数/>可以表示为:
引入连续性方程,其表达式如下:
上述方程两边同时乘以速度
联立公式(7)、公式(9)、公式(11)可以得到:
而岩石中,单相流体的动量守恒方程(Ledda P G,Siconolfi L,Viola F,etal.Suppression of von Kármán vortex streets past porous rectangular cylinders[J].Physical Review Fluids,2018,3(10).)如下:
式中,φ为岩石的孔隙度,孔隙度无量纲;ρ为单相流体的密度,单位为kg/m3;μ为单相流体的粘度,单位Pa·s;k为单相流体的渗透率,单位为m2
而对于多相流体而言在岩石中的渗流过程,可以通过微元体控制面净流入该物理量的通量与微元体单位时间内流体质量的增加量相等,方向是流出为正,流入为负,得到单位时间内油、气、水质量增加量的公式,即为多相流体质量守恒方程,如下:
其中,qogw为源汇相,单位为kg/s;此时vm为混合流体的速度,其单位为m/s;ρm为混合流体的密度,其单位为kg/m3,其表达式如下:
ρm=Soρo+Sgρg+Swρw (15)
对于油、气、水三相流体的多相动量守恒方程,此时需要考虑毛管压力Pc的作用;同时,由于多相流体或混合流体的不同组分之间存在着密度和粘度的差异,因此为了更准确地模拟和预测这些复杂流体***的行为,引入混合密度和混合粘度的概念,结合单相流体动量守恒方程(13),多相流体的动量守恒方程如下:
μw为混合流体的粘度,表达式如下:
m=Soμo+Sgμg+Swμw (17)
其中μo、μg、μw分别为油气水三相的粘度,粘度单位为Pa·s。对于毛管压力Pc,如果油相占优,那么此时三相毛管压力可以表示为:
如果气相占优,那么此时三相毛管压力可以表示为:
如果水相占优,那么此时三相毛管压力可以表示为:
pcog、pcow、pcwg分别为油气、油水以及气水直接的毛管压力,其表达式如下:
σog、σow、σwg分别为油气、油水、气水之间的界面张力,单位为N/m;θog、θow、θwg分别为油气、油水、气水之间的润湿接触角,无量纲;r为喉道半径,单位为m。
其中,混合密度是指考虑到不同组分之间的密度差异后的整体平均密度,而混合粘度则是考虑到不同组分之间的粘度差异后的整体平均粘度。
而多相流体这样的多相***,还需要考虑相方程,即连续方程(每一相的质量守恒方程),以求解各组分的体积分数(Lagrée B,Zaleski S,Bondino I.Simulation ofViscous Fingering in Rectangular Porous Media with Lateral Injection and Two-and Three-Phase Flows[J].Transport in Porous Media,2016,113(3).)。相方程的表达式如下:
qo、qg、qw为油、气、水的源汇相,单位为kg/s。
联立上述多相渗流动量守恒方程、多相渗流质量守恒方程、多相渗流相方程即可得到多相渗流数学模型。
然后基于岩心样品的数字岩心模型对多相渗流数学模型进行求解,即可得到多相渗流的岩心样品室内渗流模拟结果,其求解包括压力场求解,其具体求解过程包括以下步骤:
对于式(14)的第一项,
其中,Co、Cg、Cw分别为油、气、水三相的压缩系数,压缩系数无量纲;这一项可以进一步简化为:
其中,t为当前时刻,t+Δt为下一时刻,单位均为s。对于获取公式(14)的第二项,为了得到需要基于公式(20),将其进一步差分后可以得到油、气、水的分相渗流方程,用以计算网格中每一相流体的饱和度:
其中,A为面积,单位为m2;So t、Sg t、Sw t分别为当前时刻油、气、水三相的饱和度,So t +Δt、Sg t+Δt、Sw t+Δt为下一时刻油、气、水三相的饱和度,饱和度均无量纲;ρo t、ρg t、ρw t分别为当前时刻油、气、水三相的密度,ρo t+Δt、ρg t+Δt、ρw t+Δt为下一时刻油、气、水三相的密度,密度单位均为kg/m3;qo t、qg t、qw t为油、气、水三相的源汇相,单位为kg/s;vo t+Δt、vg t+Δt、vw t+Δt分别为油气水三相下一时刻的速度,单位为m/s;|φt+Δt为下一时刻的孔隙度,孔隙度无量纲。
对公式(16)进行差分离散,可以表示为如下:
两边同时乘以
可以求解公式得到
基于公式(26)可以得到公式(27):
所以,传导率Gm可以表示为:
基于公式(27),公式(14)的第二项可以表示为,
其中,zc为相邻网格数,该参数无量纲;A为相邻网格的面积,单位为m2
联立公式(14)、公式(22)、公式(29),可以得到:
移项可以得到:
所有网格均可将当前时刻t和下一时刻t+Δt的不同参数代入公式(31)中,从而可得到由N个与上述公式相似的方程组成的方程组。根据方程组的下标,可将方程组整理为矩阵形式,即At+ΔtXt+Δt=Bt
式中,At+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵(N为模型的网格数量),Xt+Δt和Bt是长度为N的两个向量,Xt+Δt为下一时刻t+Δt的压力场向量,Bt为与当前时刻t压力场和边界条件有关的向量,求解以上矩阵方程即得到下一时刻模型中流体的压力场。
在流体力学(CFD)模型求解过程中,对CFD求解过程进行了改进,其采用了更精确和稳定的离散化方法来计算速度场,例如,使用高阶精度格式或网格自适应技术可以减少数值耗散和扩散,并提高速度场的解析精度。这有助于更准确地捕捉流体流动的细节和涡旋结构。尤其是在处理流体在复杂多孔介质中的传播问题时,高阶精度格式和网格自适应技术对于空间离散化和时间积分步长的选择有极高的自适应性,能够根据流体流动的复杂性和流动特性自动调整,进而动态调整网格和时间步长,我们实现了加速收敛和提高求解稳定性。
同时在压力场求解过程中,为了加快压力场的收敛速度并提高求解的稳定性,引入了快速迭代方法或多重网格方法这样可以减少迭代次数并提高求解效率的高效的压力场求解算法,提高压力场的求解效率和计算量。
其中,本申请中快速迭代方法(Fast Iterative Method)是指一类可以加快迭代收敛速度的方法。这类方法主要利用方程的一些性质或者引入预条件子来改进传统的迭代方法(如雅可比方法或高斯-赛德尔方法),使得每次迭代后,解的改进更为显著,从而加快收敛的速度。快速迭代方法的种类很多,包括共轭梯度法、GMRES方法、BICGSTAB方法等。
多重网格方法:
本申请中多重网格方法(Multigrid Method)是指一类高效解决偏微分方程离散化所得到的大规模线性方程组的方法。其主要是通过在不同的网格尺度上进行迭代,充分利用了不同尺度上的信息,从而达到更快的收敛速度。多重网格方法主要包括两个步骤:平滑步骤和校正步骤。平滑步骤主要是在当前网格尺度上进行局部的迭代,消除高频误差;校正步骤则是将问题转移到更粗的网格尺度上,进一步进行迭代,消除低频误差。这两个步骤的交替执行,使得多重网格方法可以在很少的迭代次数内,就得到很好的收敛效果。
同时,在求解过程中还引入了高精度数字岩心模型以应对具有复杂边界条件或流动现象的问题,可以更好地模拟流体流动的特征,提高模拟结果的准确性。同时,其扩展的多相渗流模型,可以通过VOF体积分率方程追踪流体间的相界面来描述油、气和水相的传输行为,实现界面跟踪和相与相之间的相互作用描述,能够更准确地模拟多孔介质中液气界面的行为。此外,本发明的模拟方法还具备多尺度模拟能力,能够处理不同尺度的渗流问题,提供更全面的模拟和分析。
其中,需要说明的是,岩心样品的高精度数字岩心模型的构建包括以下步骤:
S11:获取待模拟储层的岩心样品,并利用所述岩心样本获取所述岩心样品的实际孔隙度φ和实际渗透率k;
S12:利用扫描技术获取所述岩心样本在不同断面上的二维图像,获得二维图像数据;
S13:根据多个所述断面上的二维图像数据,获得岩心三维数据体;
S14:根据岩心的实际孔隙度φ与二维图像数据的关系,获得所述岩心样品的孔隙度转化系数;根据所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系,获取所述岩心样品的渗透率转化系数;
S15:在所述岩心三维数据体上应用所述孔隙度转化系数和所述渗透率转化系数,以及实际孔隙度与二维图像数据的关系、实际渗透率与实际孔隙度及二维图像数据的关系,得到岩心孔隙度三维数据体和岩心渗透率三维数据体;
S16:建立所述岩心样品的三维几何模型,根据所述岩心孔隙度三维数据体和岩心渗透率三维数据体赋值所述三维几何模型,得到所述岩心样品的高精度数字岩心模型。
需要说明的是,所述岩心样品经过洗油洗盐并且在一定温度下,如80℃下烘干后的岩心样品。其中,所述扫描技术可以为核磁共振扫描技术或者CT扫描技术中的任一种,所述岩心样品在饱和地层水之后进行断面扫描。其岩心样品进行扫描时,其断面位置和断面数可以根据扫描时所采用的仪器的精度进行调整,精度高,则其断面数可以适当提高;精度低,则其断面数可以适当降低。其中获得的二维图像数据包括图像数据,即图像的像素体数据,和断面坐标,即切片位置坐标。获得的二维图像数据,可以保存为TXT文本,以便于进行后续数据处理。
当所述扫描技术为核磁共振扫描技术时,所述实际孔隙度φ与二维图像数据的关系为φ=αMRI×VMRI;其中,VMRI为图像数据,αMRI为转换系数;所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系:k=βMRI×VMRI 2×φ;βMRI为转换系数;其中VMRI为岩石样品每一个扫描区域i的扫描数据VMRIi的平均值。
当所述扫描技术为CT扫描技术时,所述实际孔隙度φ与二维图像数据的关系为φ=αCT×1/VCT;其中,VCT为图像数据,αCT为转换系数;所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系:k=βCT×(1/VCT 2)×φ;βCT为转换系数;其中VCT为岩石样品每一个扫描区域i的扫描数据VCTi的平均值。
其中,所述图像数据,即像素体数据可以为岩石密度或者灰度或者像素数中的任一种。由于图像的灰度或者像素数可以用于表示岩石密度的大小,因此在本申请中直接用灰度或者像素数直接表示相对的岩石密度。
由于扫描得到的数据体难以正好落在便于处理的范围内,因此在所述步骤S3中,在根据多个所述断面上的二维图像数据,获得岩心三维数据体之前,先对所述二维图像数据进行粗化处理和/或插值处理。如,当数据体中的节点数达到上千万甚至上亿时,对所述二维数据图像进行粗化。其中,所述二维图像数据经所述粗化处理后,其网格数量不低于100万且不高于300万;而当所述二维图像数据的节点数低于100万个时,对所述二位图像数据进行插值处理,利用已知邻近像素点的灰度值(或rgb图像中的三色值)来产生未知像素点的灰度值,以便由原始图像再生出具有更高分辨率的图像。
其中,所述实际孔隙度φ和实际渗透率k为任一所述利用扫描技术扫描的岩心样品的切片的实际孔隙度φ和实际渗透率k;或者;所述实际孔隙度φ和实际渗透率k为整个岩心样品的测试得到的实际孔隙度和实际渗透率。
在一些实施例中,所述插值算法包括三线性插值算法或克里金插值算法。
在一些实施例中,其特征在于,在步骤S11中,当所述扫描技术为CT扫描技术时,所述岩心样品通过洗油、洗盐以及烘干处理后并测量实际孔隙度φ和实际渗透率k以及岩石样本的长度、直径;
当所述扫描技术为核磁共振扫描技术时,所述岩心样品通过洗油、洗盐以及烘干处理后,再进行饱和地层水处理,并测量实际孔隙度φ和实际渗透率k以及岩石样本的长度、直径。
综上,本申请引入液体和气体的可压缩性、考虑混合密度和混合粘度,建立了多相渗流数学模型,同时在求解过程中改进了传统的CFD模型求解过程中最常用的投影法,使得CFD模型可以更准确地描述液体和气体的行为和性质,提高模拟结果的准确性和可靠性,为工程设计和科学研究提供了有力的工具。
最后应说明的是:以上各实施例仅用以说明本申请的技术方案,而非对其限制;尽管参照前述各实施例对本申请进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或对其中部分或全部技术特征进行等同替换;而这些修改或替换,并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。

Claims (10)

1.一种多相流体地层渗流的数值模拟方法,其特征在于,包括:
S1:获取岩心样品的孔隙度三维数据体和渗透率三维数据体,构建高精度数字岩心模型;
S2:通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型;
S3:根据高精度数字岩心模型和渗流模拟数学模型建立三维流体力学模型,根据所述三维流体力学模型求解所述多相流体的压力场;
其中,在通过VOF模型追踪多相流体间的相界面,描述多相流体的传输行为,建立渗流模拟数学模型中,引入所述多相流体的压缩系数、混合流体的混合粘度和混合密度。
2.根据权利要求1所述的一种多相流体地层渗流的数值模拟方法,其特征在于,所述渗流模拟数学模型根据多相流体动量守恒方程、多相流体质量守恒方程以及多相流体相方程构建。
3.根据权利要求2所述的一种多相流体地层渗流的数值模拟方法,其特征在于,其中多相流体质量守恒方程为:
其中,A为微元体的面积,m2;V为微元体的体积(dV=dxdydz,x、y、z为微元体的三维坐标 的三个方向),m3;m为微元体流体的质量,kg;t为时间,s;qogw为源汇相,kg/s;为孔隙度;vm 为混合流体的速度,m/s;ρm为混合流体的密度,kg/m3,其表达式如下:
ρm=Soρo+Sgρg+Swρw
其中,So、Sg和Sw分别为油相、气相、水相的饱和度;ρo、ρg、ρw分别为油相、气相、水相流体的密度,kg/m3
4.根据权利要求2所述的一种多相流体地层渗流的数值模拟方法,其特征在于,多相流体动量守恒方程为:
其中,孔隙度无量纲;ρm为混合流体的密度,kg/m3;k为单相流体的渗透率,m2为压力 梯度,Pa/m;pc为三相毛管压力,毛管压力Pa;l为长度,m;▽表示梯度运算;t为时间,s;为 孔隙度;vm为混合流体的速度,m/s;μm为混合流体的粘度,表达式如下:
μm=Soμo+Sgμg+Swμw
其中,μo、μg、μw分别为油气水三相的粘度,粘度单位为Pa·s;So、Sg和Sw分别为油相、气相、水相的饱和度。
5.根据权利要求2所述的一种多相流体地层渗流的数值模拟方法,其特征在于,多相流体相方程为:
其中,qo、qg、qw为油、气、水的源汇相,kg/s;So、Sg和Sw分别为油相、气相、水相的饱和度; ρo、ρg、ρw分别为油相、气相、水相流体的密度,kg/m3;vo、vg、vw分别为油相、气相、水相流体的 速度,m/s;dV为微元体的体积(dV=dxdydz,x、y、z为微元体的三维坐标的三个方向),m3;A 为微元体的面积,m2;t为时间,s;为孔隙度。
6.根据权利要求4所述的一种多相流体地层渗流的数值模拟方法,其特征在于,所述三相毛管压力pc根据多相流体的优势相,分别通过下式获得:
当油相为优势相,则三相毛管压力pc为:
当气相为优势相,则三相毛管压力pc为:
当水相为优势相,则三相毛管压力pc为:
其中:So、Sg和Sw分别为油相、气相、水相的饱和度;▽表示梯度运算;pcog、pcow、pcwg分别为油气、油水以及气水直接的毛管压力,其表达式如下:
其中,σog、σow、σwg分别为油气、油水、气水之间的界面张力,单位为N/m;θog、θow、θwg分别为油气、油水、气水之间的润湿接触角,无量纲;r为喉道半径,单位为m。
7.根据权利要求1所述的一种多相流体地层渗流的数值模拟方法,其特征在于,在步骤S3中,利用快速迭代方法或多重网格方法求解所述多相流体的压力场。
8.根据权利要求1~7任一项所述的一种多相流体地层渗流的数值模拟方法,其特征在于,所述多相流体的压力场求解包括:
求解多相流体质量方程中的第一项:
其中,Co、Cg、Cw分别为油、气、水三相的压缩系数,压缩系数无量纲;t为当前时刻,t+Δt 为下一时刻,单位均为s;So、Sg和Sw分别为油相、气相、水相的饱和度;ρo、ρg、ρw分别为油相、 气相、水相流体的密度,kg/m3;V为体积,m3为孔隙度;
求解多相流体质量方程中的第二项:
其中,zc为相邻网格数,该参数无量纲;A为微元体的面积,单位为m2;每一网格单元即为1个微元体,i为网格i,j为与网格i相邻的网格j;Gm为传导率;ρm为混合流体的密度,kg/m3;vm为混合流体的速度,m/s;p为压力,单位为Pa;
联立多相流体质量方程、多相流体质量方程中的第一项的求解公式和第二项的求解公式,得到:
将当前时刻t和下一时刻t+Δt的不同参数代入上式,得到N个在不同时刻的上式,形成上式的矩阵方程At+ΔtXt+Δt=Bt,其中:
其中,At+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵(N为模型的网格数量),Xt +Δt和Bt是长度为N的两个向量,Xt+Δt为下一时刻t+Δt的压力场向量,Bt为与当前时刻t压力场和边界条件有关的向量;qogw为源汇相,kg/s;
求解每一时刻的矩阵方程,得到所述多相流体的压力场。
9.根据权利要求1所述的一种多相流体地层渗流的数值模拟方法,所述高精度数字岩心模型的建立包括:
利用岩心样品和扫描技术,获取所述岩心样品的孔隙度三维数据体和渗透率三维数据体;
建立所述岩心样品的三维几何模型,根据所述孔隙度三维数据体和渗透率三维数据体赋值所述三维几何模型,得到所述岩心样品的高精度数字岩心模型;
其中,所述孔隙度三维数据体和渗透率三维数据体的获得包括以下步骤:
S11:获取待模拟储层的岩心样品,并利用所述岩心样品获取所述岩心样品的实际孔隙度φ和实际渗透率k;
S12:利用扫描技术获取所述岩心样品在不同断面上的二维图像,获得二维图像数据;
S13:根据多个所述断面上的二维图像数据,获得岩心三维数据体;
S14:根据岩心的实际孔隙度φ与二维图像数据的关系,获得所述岩心样品的孔隙度转化系数;根据所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系,获取所述岩心样品的渗透率转化系数;
S15:在所述岩心三维数据体上应用所述孔隙度转化系数和所述渗透率转化系数,以及实际孔隙度与二维图像数据的关系、实际渗透率与实际孔隙度及二维图像数据的关系,得到岩心孔隙度三维数据体和岩心渗透率三维数据体。
10.一种在油田开发、地下水管理、地质工程和环境保护技术领域中的应用,具体包括权利要求1~9任一项所述一种多相流体地层渗流的数值模拟方法。
CN202311027224.7A 2023-08-15 2023-08-15 一种多相流体地层渗流的数值模拟方法及应用 Active CN117113873B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311027224.7A CN117113873B (zh) 2023-08-15 2023-08-15 一种多相流体地层渗流的数值模拟方法及应用

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311027224.7A CN117113873B (zh) 2023-08-15 2023-08-15 一种多相流体地层渗流的数值模拟方法及应用

Publications (2)

Publication Number Publication Date
CN117113873A CN117113873A (zh) 2023-11-24
CN117113873B true CN117113873B (zh) 2024-04-09

Family

ID=88812081

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311027224.7A Active CN117113873B (zh) 2023-08-15 2023-08-15 一种多相流体地层渗流的数值模拟方法及应用

Country Status (1)

Country Link
CN (1) CN117113873B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103617291A (zh) * 2013-12-14 2014-03-05 中国海洋石油总公司 储层成因单元界面等效表征方法
CN105021506A (zh) * 2015-07-09 2015-11-04 中国石油大学(华东) 一种基于孔隙网络模型的三相相对渗透率的计算方法
CN111624147A (zh) * 2020-04-16 2020-09-04 中国石油天然气股份有限公司 岩心的相对渗透率测定方法及装置
WO2021180189A1 (zh) * 2020-03-13 2021-09-16 重庆科技学院 一种多元热流体热采油藏数值模拟方法
CN114239367A (zh) * 2021-12-31 2022-03-25 西南石油大学 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
CN114386302A (zh) * 2021-12-31 2022-04-22 西南石油大学 一种非定常流固耦合多相渗流模型构建方法
CN115329616A (zh) * 2022-08-08 2022-11-11 中国科学院武汉岩土力学研究所 一种含相变的多孔介质热流固化耦合多相渗流模拟方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180252087A1 (en) * 2015-03-27 2018-09-06 Schlumberger Technology Corporation A method and a system for performing chemical treatment of a near wellbore area

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103617291A (zh) * 2013-12-14 2014-03-05 中国海洋石油总公司 储层成因单元界面等效表征方法
CN105021506A (zh) * 2015-07-09 2015-11-04 中国石油大学(华东) 一种基于孔隙网络模型的三相相对渗透率的计算方法
WO2021180189A1 (zh) * 2020-03-13 2021-09-16 重庆科技学院 一种多元热流体热采油藏数值模拟方法
CN111624147A (zh) * 2020-04-16 2020-09-04 中国石油天然气股份有限公司 岩心的相对渗透率测定方法及装置
CN114239367A (zh) * 2021-12-31 2022-03-25 西南石油大学 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
CN114386302A (zh) * 2021-12-31 2022-04-22 西南石油大学 一种非定常流固耦合多相渗流模型构建方法
CN115329616A (zh) * 2022-08-08 2022-11-11 中国科学院武汉岩土力学研究所 一种含相变的多孔介质热流固化耦合多相渗流模拟方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Three-dimensional Digital Core Structure Characterization Method and Application of Thick Limestone in Baotu Spring Basin of Jinan;Chuanlei Li等;《IEEE》;全文 *
低渗透砂岩油水两相流动压力波动特征;王建忠等;《石油地质与工程 》;全文 *

Also Published As

Publication number Publication date
CN117113873A (zh) 2023-11-24

Similar Documents

Publication Publication Date Title
Zhang Pore-scale modelling of relative permeability of cementitious materials using X-ray computed microtomography images
CN107060746B (zh) 一种复杂裂缝性油藏流动模拟的方法
Laloui et al. Solid–liquid–air coupling in multiphase porous media
Aghaei et al. Direct pore-to-core up-scaling of displacement processes: Dynamic pore network modeling and experimentation
Zhu et al. A pore‐scale numerical model for flow through porous media
Alpak et al. Direct simulation of pore-scale two-phase visco-capillary flow on large digital rock images using a phase-field lattice Boltzmann method on general-purpose graphics processing units
CN114239367B (zh) 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
Kuila et al. A Robust and accurate Riemann solver for a compressible two-phase flow model
Van Marcke et al. An improved pore network model for the computation of the saturated permeability of porous rock
Das et al. A finite volume model for the hydrodynamics of combined free and porous flow in sub-surface regions
Verma et al. Effect of wettability on two-phase quasi-static displacement: Validation of two pore scale modeling approaches
Yang et al. Modelling of fluid–structure interaction with multiphase viscous flows using an immersed-body method
CN114386302B (zh) 一种非定常流固耦合多相渗流模型构建方法
Khayrat et al. A multi-scale network method for two-phase flow in porous media
CN110263434A (zh) 一种基于多尺度混合有限元的流动单元数值模拟方法
Wang et al. Hydro-mechanical analysis of piping erosion based on similarity criterion at micro-level by PFC3D
He et al. Numerical estimation and prediction of stress-dependent permeability tensor for fractured rock masses
Liu et al. Partially implicit finite difference scheme for calculating dynamic pressure in a terrain-following coordinate non-hydrostatic ocean model
Tavakoli et al. Numerical simulation of liquid/gas phase flow during mold filling
CN117113873B (zh) 一种多相流体地层渗流的数值模拟方法及应用
Pepper et al. Application of h‐adaptation for environmental fluid flow and species transport
Khayrullin et al. Methods for Studying Two-Phase Flows in Porous Media: Numerical Simulation and Experiments on Microfluidics Chips
Yan et al. Dynamic effect in capillary pressure–saturation relationship using lattice Boltzmann simulation
Zhao et al. A three-dimensional one-layer particle level set method
Khokhar et al. A numerical analysis of mixing and separating of Newtonian fluids in a channel filled with porous materials

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