CN114386302B - 一种非定常流固耦合多相渗流模型构建方法 - Google Patents
一种非定常流固耦合多相渗流模型构建方法 Download PDFInfo
- Publication number
- CN114386302B CN114386302B CN202111662690.3A CN202111662690A CN114386302B CN 114386302 B CN114386302 B CN 114386302B CN 202111662690 A CN202111662690 A CN 202111662690A CN 114386302 B CN114386302 B CN 114386302B
- Authority
- CN
- China
- Prior art keywords
- fluid
- seepage
- pore
- core
- 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
Links
- 239000007787 solid Substances 0.000 title claims abstract description 82
- 238000010276 construction Methods 0.000 title claims abstract description 12
- 239000011148 porous material Substances 0.000 claims abstract description 185
- 239000012530 fluid Substances 0.000 claims abstract description 147
- 238000000034 method Methods 0.000 claims abstract description 112
- 238000010168 coupling process Methods 0.000 claims abstract description 70
- 238000005859 coupling reaction Methods 0.000 claims abstract description 70
- 230000008878 coupling Effects 0.000 claims abstract description 66
- 230000008569 process Effects 0.000 claims abstract description 59
- 238000013178 mathematical model Methods 0.000 claims abstract description 55
- 238000004088 simulation Methods 0.000 claims abstract description 55
- 239000007788 liquid Substances 0.000 claims abstract description 31
- 238000006073 displacement reaction Methods 0.000 claims abstract description 20
- 239000011435 rock Substances 0.000 claims description 43
- 230000035699 permeability Effects 0.000 claims description 25
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 18
- 239000011159 matrix material Substances 0.000 claims description 16
- 238000004422 calculation algorithm Methods 0.000 claims description 14
- 230000006835 compression Effects 0.000 claims description 13
- 238000007906 compression Methods 0.000 claims description 13
- 238000006243 chemical reaction Methods 0.000 claims description 12
- 239000013598 vector Substances 0.000 claims description 9
- 238000002347 injection Methods 0.000 claims description 6
- 239000007924 injection Substances 0.000 claims description 6
- 239000000243 solution Substances 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 4
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 239000010779 crude oil Substances 0.000 claims description 3
- 238000001035 drying Methods 0.000 claims description 3
- 239000008398 formation water Substances 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 abstract description 9
- 238000002474 experimental method Methods 0.000 abstract description 8
- 238000012360 testing method Methods 0.000 abstract description 3
- 239000007789 gas Substances 0.000 description 58
- 238000002595 magnetic resonance imaging Methods 0.000 description 32
- 238000005481 NMR spectroscopy Methods 0.000 description 15
- 238000004364 calculation method Methods 0.000 description 14
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 14
- 239000003921 oil Substances 0.000 description 12
- 239000003345 natural gas Substances 0.000 description 7
- 238000011160 research Methods 0.000 description 7
- 238000011161 development Methods 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 6
- 238000005259 measurement Methods 0.000 description 5
- 238000013421 nuclear magnetic resonance imaging Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000009792 diffusion process Methods 0.000 description 4
- 239000003209 petroleum derivative Substances 0.000 description 4
- 238000011084 recovery Methods 0.000 description 4
- 238000009795 derivation Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000002591 computed tomography Methods 0.000 description 2
- 230000001808 coupling effect Effects 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 238000001208 nuclear magnetic resonance pulse sequence Methods 0.000 description 2
- 238000005325 percolation Methods 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 241001270131 Agaricus moelleri Species 0.000 description 1
- 102000004190 Enzymes Human genes 0.000 description 1
- 108090000790 Enzymes Proteins 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013079 data visualisation Methods 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 238000005206 flow analysis Methods 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000003079 shale oil Substances 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000009736 wetting Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种非定常流固耦合多相渗流模型构建方法,包括根据岩心尺度孔隙网络结合非定常流动模型分析流体在孔喉通道中的非定常流动,获得圆形管束中流体轴向速度分布,分别构建非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型;利用非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型,结合孔隙网络模型中混合流体的特性参数构建非混相驱替过程下的流固耦合多相渗流数学模型,将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合,进行岩心流固耦合渗流模拟。本发明可模拟再现室内岩心流动实验过程,实现室内岩心的数字化单相和多相流固耦合渗流模拟分析与测试。
Description
技术领域
本发明涉及油气田开发领域,尤其涉及一种非定常流固耦合多相渗流模型构建方法,适用于常规砂岩油气藏,以及致密油气、页岩油气、天然气水合物等非常规油气藏,也适用于CO2地质封存技术。
背景技术
石油和天然气是保持国民经济高速发展的重要能源之一,如何合理开采石油和天然气并提高其采收率一直是油气田发发过程中的重要问题。实际地下储层岩石内部孔喉结构复杂,很难通过实验手段厘清流体在其中的渗流规律。许多研究人员利用多孔介质模型来模拟不同类型流体在岩石内部流动,从而找到有利于提高石油采收率的方法。数字岩心技术作为多孔介质模型的一个分支,可以用于石油与天然气行业的地质、地震、测井和开发以及提高采收率等各个领域。数字岩心能有效地保留岩心微观物理特征,能确保岩心可以无限次地被使用,是岩石物理实验数值模拟的重要平台,能定量研究岩心内部各种微观因素(如孔隙连通性、润湿性等)对储层渗流过程的影响,能够计算传统物理实验无法直接测量的物理性质,如油气水三相相对渗透率。鉴于其广泛的应用性,开展数字岩心的研究对提高石油和天然气采收率具有重要意义。
近年来,随核磁共振成像技术(MRI)的快速进步和计算机GPU芯片算力的显著提升,为构建与室内实际岩石样品相对应的三维岩心数字化孔隙网络模型,完善孔隙网络模型多相渗流理论、数学模型和数值模拟方法,以及结合岩心数字化孔隙网络模型和孔隙网络模型多相渗流数值模拟方法开展岩心渗流模拟研究奠定了重要的物质基础。同时,也催生了岩心数字化孔隙网络模型和孔隙网络模型多相渗流数值模拟方法这一技术方法。孔隙网络模拟过程中,认为单相或两相流体在任意节点中的流动均遵守Kirchoff定律,即流入和流出节点的流量之和,认为流体压力可以瞬间从入口传递到出口,但是实际的流体通常具有可压缩性,压力波无法瞬间传播,从而无法准确还原压力波对流体流动的影响。因此,需要在传统孔隙网络模型中引入不稳定渗流理论,将拉普拉斯方程转换为非稳态渗流-压力扩散方程,模拟流体流动过程,对石油和天然气相关行业的发展具有重要意义。
发明内容
本发明的目的在于克服现有技术的不足,提供一种非定常流固耦合多相渗流模型构建方法,在构建的岩心数字化孔隙网络模型基础上,为了准确描述多相流体流动过程以及流体-岩石骨架相互作用关系,通过流体力学雷诺输运方程、有限体积法和流固耦合渗流机理,推导出适用于岩心数字化孔隙网络模型的非定常流固耦合多相渗流数学模型,并提出相应的隐式数值模拟方法,形成一套基于核磁共振MRI/T2的岩心数字化建模和非定常流固耦合渗流模拟方法。
本发明的目的是通过以下技术方案来实现的:
一种非定常流固耦合多相渗流模型构建方法,包括:
S1,根据岩心尺度孔隙网络结合非定常流动模型分析流体在孔喉通道中的非定常流动,获得圆形管束中流体轴向速度vr分布,
S2,根据单相液体渗流过程满足的假设条件,采用有限体积法和雷诺运输方程来构建非定常单相液体流固耦合渗流数学模型;
S3,根据单相气体渗流过程中的气体的密度ρg和气体压缩系数Cg,结合雷诺运输方程构建同时满足低压和高压条件的非定常单相气体流固耦合渗流数学模型;
S4,利用非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型,结合孔隙网络模型中混合流体的特性参数构建非混相驱替过程下的流固耦合多相渗流数学模型;
S5,将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合,进行岩心流固耦合渗流模拟。
具体的,非定常单相液体流固耦合渗流数学模型如下列表达式所示:
式中,Ct=Cρ+Cp,Ct为综合压缩系数,Cp为孔隙压缩系数,Pa-1;Cρ为液体压缩系数,Pa-1;Δt为时间步长;qij为管束内流体的体积流量,gij为管束内液体的水力传导率;<p>ij为相邻节点i和j之间的平均压力;Δpij为相邻节点i和j之间管束的流体压力差;Vpi0为节点i在初始时刻的孔隙体积,R0ij为初始时刻节点i和j之间孔喉管束的半径;lij为相邻节点i和节点j之间的孔喉长度;n为与控制体中心节点i相连通的节点数;μ为流体粘度。
具体的,非定常单相气体流固耦合渗流数学模型如下列表达式所示:
式中,ggij为管束中气体流动的水力传导率;μg为气体黏度;p/Zμg为非线性项;p为控制体中心节点处的孔隙流体压力,Rg为气体常数。
具体的,流固耦合多相渗流数学模型如下列表达式所示:
Ct=Cp+CISI+CDSD
式中,gij为管束内流体的水力传导率,pcij为节点i和j之间孔喉通道的毛管力;μeff为节点i和j之间混合流体的有效粘度;CI和CD分别为注入流体和被驱替流体的压缩系数;SI和SD分别为以节点i为中心的控制体内注入流体和被驱替流体的饱和度。
具体的,非定常流固耦合多相渗流数值模拟方法具体包括以下步骤:
S201,首先,对流固耦合多相渗流数学模型中的非线性毛管压力项进行线性化处理:
从而获得线性化后的流固耦合多相渗流数学模型:
S202,采用隐式数值模拟方法,引入源汇项Qi,对线性化后的流固耦合多相渗流数学模型进行离散处理,得到:
式中,上标t表示当前时刻下的流动状态,上标t+Δt表示下一时刻的流动状态;
S203,对步骤7023公式中不同时刻的流动状态进行分离与合并,得到:
S204,遍历孔隙网络模型中的所有网络节点,孔隙网络模型中所有控制体均可将当前时刻和下一时刻的流动状态带入步骤S703的公式中,整理后形成如下矩阵:
[A]t+Δt[X]t+Δt=[B]t,
式中,[A]t+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵,N为孔隙网络模型的节点数,[X]t+Δt和[B]t是长度为N的两个向量,[X]t+Δt为下一时刻的压力场向量,[B]t为与上一时刻压力场和边界条件有关的向量,利用GPU代数多重网格广义最小残差算法求解以上矩阵获得当前时刻孔隙网络模型中的流体压力场分布;
S205,非定常流固耦合多相渗流数值模拟过程中,采用固定时间步长或变时间步长的方式使该步长内流体界面发生位移,每一时间步下,计算出界面移动后新的流体界面位置,并更新整个孔隙网络模型中所有孔喉通道的水力传导率和模型中各相流体的饱和度,再进行流体压力场分布求解,直至整个孔隙网络空间被侵入流体占满或达到某一饱和度数值时停止。
具体的,步骤S5还包括岩心数字化孔隙网络模型的构建过程,包括:
S301,选取岩心样本烘干后测量岩心孔喉长度、孔隙度和渗透率,并对岩心样本进行抽真空饱和度模拟地层水或模拟地层原油;
S302,对岩心样本进行核磁共振MRI/T2扫描测量,获取岩心在不同断面上的二维图像;
S303,对二维图像进行插值,获取岩心MRI/T2三维数据体;
S304,建立三维无序孔隙网络模型,将MRI/T2三维数据体赋值到三维无序孔隙网络模型的节点中;
S305,通过转换系数α计算三维无序孔隙网络模型中各相邻节点的孔喉半径,模拟计算岩心数字化孔隙网络模型的渗透率;调整转换系数α大小,确保孔隙网络渗透率与岩心实测渗透率相近,建立岩心数字化孔隙网络模型。
本发明在构建的岩心数字化孔隙网络模型基础上,为了准确描述多相流体流动过程以及流体-岩石骨架相互作用关系,方法的详细实现流程如下:
①流体在孔喉通道中的非定常流动
孔隙尺度的非定常流动模型是基于假设压力下降是由惯性和流体粘度引起的。假定流动为层流状态,由Navier-Stokes方程简化得到圆形管束内流体的轴向速度vr为:
根据Nguyen和Choi(Nguyen,Q.H.,&Choi,S.B.,A new approach for ananalytical solution of unsteady laminar flow in dispensingprocesses.Proceedings of the Institution of Mechanical Engineers,Part C:Journal of Mechanical Engineering Science,2010,224(6),1231-1243.)的研究,圆形管束中流体轴向速度vr分布为:
式中,R为圆形管束半径,Δp为管束两端压差,l为管束长度,μ流体粘度,t为流动时间,ρ为流体密度。采用式2,分别计算管束半径R=1微米和10微米时,牛顿流体(水,密度ρ=1000kg/m3,粘度μ=1mPa·s)流动时的速度分布可以发现,当时间t分别达到7e-7秒和5e-5秒时,流体的速度剖面达到稳定流动时的抛物线分布。根据式2可以发现,微观孔喉通道中流体流动更快地达到稳定状态,表明泊肃叶方程仍然可用于推导流体在多孔介质中的非定常流动。(Nguyen和Choi没有采用式2分析微纳米级孔喉通道的流体流动)
②非定常单相液体流固耦合渗流数学模型
采用有限体积法进行推导,首先明确孔隙网络中控制体是指以孔隙网络任意节点为中心,与其相连通的所有管束的二分之一处为边界所包含的全部孔隙空间,即孔隙网络的总节点数为N,则控制体就有N个。单相液体渗流过程中满足以下假设:1)流体不能穿过孔隙网络的管束边界到达管束之外,2)流体是连续相,且为牛顿流体,3)忽略流体的质量和能量之间的转换,4)液体和岩石骨架具有一定的可压缩性。实际渗流过程中,由于岩石和流体具有微可压缩性,使得流动过程压力的传播过程不稳定,产生压力波传播。以下是非定常单相液体流固耦合渗流数学模型的推导过程。
流体密度和岩石孔隙空间体积随压力变化满足(忽略孔喉长度随压力变化而变化):
式中,ρ为流体密度,kg/m3;p为孔隙流体压力,Pa;R为管道半径;Vp为管道体积,m3;Cp为孔隙压缩系数,Pa-1;Cρ为液体压缩系数,Pa-1。上述式子分离变量后积分得到:
同理,孔喉半径R满足:
式中:Vp0,R0,ρ0分别为初始时刻(p0压力下)的孔隙体积、孔喉半径和液体密度。将式4按麦克劳林级数展开并取前两项可得:
V≈V0[1+Cp(p-p0)] (5c)
ρ≈ρ0[1+Cρ(p-p0)] (5c)
孔隙网络中,任意控制体内流体流动的质量守恒定律采用雷诺输运方程进行表示:
式中,V表示与任意节点为中心的控制体的体积,ρvdA表示任意时刻流体流入或流出控制体的质量流量。雷诺输运方程描述了在没有汇或源的情况下,通过控制面上净输出的该物理量的通量与控制体内流体质量的当地时间变化率。任意控制体内,
式中,根据式4a和4b,
ρVp≈ρ0Vp0+ρ0Vp0(Cρ+Cp)(p-p0)+ρ0Vp0CρCp(p-p0)2 (7b)
由于Cp和Cρ均为很小的数,忽略式7b中的高阶项可得,
式中,lij为相邻节点i和节点j之间的孔喉长度。
对于控制体中任意方向孔喉通道内流体的流入和流出的净流量表示为:
式中,Aij为相邻节点i和节点j之间孔喉的横截面积。管道变形和压力波动会破坏稳态层流速度分布,但根据式2的分析结果可以发现,微观孔喉通道中流体流动的速度剖面会快速恢复至稳态抛物线速度分布。因此,任意孔喉通道的流量满足:
由式9a和9b可得,
式中,<p>ij为相邻节点i和j之间的平均压力,Δpij为相邻节点i和j之间管束的流体压力差(压差)。根据式6、8和10,可得到孔隙网络模型的非定常单相液体流动数学模型:
式中,Ct=Cρ+Cp(Ct为综合压缩系数),方程两端的流体密度被约去,Δt为时间步长。式11a中有单管流动方程qij=gijΔpij(式中qij为管束内流体的体积流量,gij为管束内液体的水力传导率,)。式11是扩散方程,该方程描述了孔隙网络模型中流体流动过程中的压力扩散(传播)。当Ct=Cρ+Cp=0时,式11退化为基尔霍夫定律和拉普拉斯方程(稳态单相渗流数学模型)。
③非定常单相气体流固耦合渗流数学模型
对于真实气体,气体的密度ρg和气体压缩系数Cg计算式如下:
式中,M为气体分子摩尔质量,Z为气体偏差因子,p为控制体中心节点的压力,Rg为气体常数,T为开氏温度。气体在孔隙网络模型中的流动过程仍然满足雷诺输运方程(式6)。此时,式6的左边第一项可写成如下关系:
由式4、5和12,
其中,
同样地,相较于Cp,CpCρ是很小的数,忽略该项后,根据式14可得,
对于气体流动,雷诺输运方程(式6)左边第二项满足如下关系:
此时,气体流动穿过控制面的净流量方程如下:
联立式6、式15和式17,得到适用于孔隙网络的气体单相非定常流固耦合渗流数学模型通式:
对式18中非线性项p/Zμg的处理如下:低压条件(气体压力小于10MPa)下,假定气体为理想气体,此时气体偏差因子Z为1;气体黏度μg受压力影响较小,可近似认为μg为常数;采用数值解法时,式18中的压力p可将上一时间步的压力值直接带入计算下一时间步的压力场,此时式18两边的压力可约去。高压条件(气体压力大于10MPa)下,p/Zμg近似为常数。当满足以上两种条件时,时18可简化为:
式中,ggij为管束中气体流动的水力传导率。此时,气体压缩系数Cg远大于岩石的孔隙压缩系数Cp,因此上式右端项忽略了孔隙压缩系数。对比式11和式19可以发现,非定常单相液体流固耦合渗流(式11)和非定常单相气体流固耦合渗流(式18、19)在不同的约束条件下,具有相同的数学表达式。
④非定常流固耦合多相渗流数学模型
这里的多相是指,油-水两相,气-水两相和油-气两相。以下以油水两相驱替过程为例进行介绍。水驱油过程满足如下假设:1)任意孔喉中不同流体之间最多只有一个流体界面,2)流动过程中不发生混相,3)孔喉中发生活塞式驱替,4)忽略重力的影响。
驱替过程中,任意两相邻节点i和j之间的孔喉通道内,混合流体的密度计算式为:
任意以节点i为中心的控制体内,混合流体的密度计算式为:
式中,ρI为注入流体的密度,ρD为被驱替流体的密度,VijI和VijD分别为节点i和j之间孔喉通道内注入流体和被驱替流体的体积,ViI和ViD为以节点i为中心的控制体内注入流体和被驱替流体的体积。为简化计算过程,这里假定以节点i为中心的控制体内混合流体的密度与节点i和j之间孔喉通道内混合流体的密度相近,则定义控制体内混合流体的有效密度ρeff为:
ρij≈ρi=ρeff=ρISI+ρDSD (21)
式中,SI和SD分别为以节点i为中心的控制体内注入流体和被驱替流体的饱和度。当控制体内充满注入流体时,则SI=1,控制体内流体密度为ρeff=ρI;若充满被驱替流体时,则控制体内流体密度为ρeff=ρD。孔隙网络模型中,任意控制体内混合流体的压缩性满足如下关系:
Ceffρ=CISI+CDSD (22c)
式中,Ceffρ为混合流体的有效压缩系数,ρeff0为初始条件(压力p0)下混合流体的密度,CI和CD分别为注入流体和被驱替流体的压缩系数。非混相驱替过程中,任意时刻多相流体在孔隙网络模型中的流动过程均服从雷诺输运方程。与非定常单相流体(气体和液体)流固耦合渗流数学模型推导过程中采用的方法类似,可以得到非混相驱替过程下的流固耦合多相渗流数学模型:
Ct=Cp+CISI+CDSD (23c)
式中,gij为管束内流体的水力传导率,pcij为节点i和j之间孔喉通道的毛管力(当该孔喉通道中存在非混相流体界面时,pcij=2γcosθ/Rij,γ为非混相流体之间的界面张力,θ为润湿接触角;当孔喉通道为单相流体时,pcij=0),μeff为节点i和j之间混合流体的有效粘度(μeff=μIXij+μD(1-Xij),Xij为非混相流体界面在孔喉通道内的相对位置(0≤Xij≤1),μI为注入流体的粘度,μD为被驱替流体的粘度)。当控制体内只有注入流体时,μeff=μI;当控制体内只有被驱替流体时,μeff=μD。当某一个存在两相流体界面的孔喉通道中Δpij<pcij,对应孔喉的水力传导率gij=0,此时该流体界面“锁死”。当孔隙网络所有控制体内均只有一相流体(气、油或水)时,式23退化为式11或式19。
式23适用于油水、气水和油气两相渗流。若多渗流过程中某一相流体为气体时,气体的压缩系数Cg应采用式13进行计算;若假定气体为理想气体,则气体压缩系数Cg=1/p。
⑤非定常流固耦合多相渗流数值模拟方法
流固耦合多相渗流数学模型(式23)可以看作是非定常单相渗流数学模型(式11和19)的进一步推广。这里以式23为例介绍数值模拟方法。首先,对式23中非线性毛管压力项进行线性化处理:
从而将式23a改写为:
式25与式23完全等价,但数学模型的形式与非定常单相渗流数学模型(式11和19)更为接近。采用隐式数值模拟方法,引入源汇项Qi(表示某一节点内可能有外部流体注入,或有流体从孔隙网络模型中流出;否则Qi=0),对式25进行离散处理:
式中,上标t表示当前时刻下的流动状态(流体水力传导率和压力),上标t+Δt表示下一时刻的流动状态。对式26不同时刻的流动状态进行分离与合并(将当前时刻的流动状态量放到方程右边,下一时刻的流动状态量放到左边),得到:
遍历孔隙网络模型中的所有网络节点,孔隙网络模型中所有控制体均可将当前时刻和下一时刻的流动状态带入式26b中,从而可得到由N个与式26b相似的方程组成的方程组。根据方程组的下标,可将方程组整理为矩阵形式[A]t+Δt[X]t+Δt=[B]t,式中,[A]t+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵(N为孔隙网络模型的节点数),amn为其中的元素;[X]t+Δt和[B]t是长度为N的两个向量,[X]t+Δt为下一时刻的压力场向量,[B]t为与上一时刻压力场和边界条件有关的向量,bi为其中的元素,且:
其中,
求解以上矩阵方程即得到当前时刻孔隙网络模型中的流体压力场。
由于流动过程中两相界面随时间推进,必须选定一个时间步长使该步长内流体界面均发生了适量的位移Δx=Δt(qij/Aij)(Aij为孔喉通道的横截面积),这里的“适量”是指必须在保证精度的前提下尽量减少运算次数,这里可采用固定时间步长。固定时间步长是指将渗流瞬态划分为长度相等的时间间隔Δtk(k=1,2,…),当Δtk足够小时,相邻时间步长的压力场变化可以认为是稳定、线性的,此时可以使用与单相流部分相同的方法对压力场进行求解。
模拟过程也可以采用变时间步长的方法以加快数值模拟计算速度。驱替过程一直进行直到。每一时间步下,计算出界面移动后新的流体界面位置,并更新整个网络模整个孔隙网络空间被侵入流体占满或达到某一饱和度数值时停止型中所有孔喉通道的水力传导率和模型中各相流体的饱和度,再进行压力场求解。模拟过程中,流体界面发生移动后需要更新水力传导率,甚至发生流体界面“锁死”的情况,因此求解过程中会产生一个严重病态的矩阵(式27),大大增加计算难度。模拟过程中使用AMGCL算法库中的GPU代数多重网格广义最小残差算法(AMG-GMRES)求解该矩阵。整个模拟过程保持注入流量Q不变。
本发明还包括基于核磁共振MRI/T2数据建立岩心数字化孔隙网络模型,具体流程如下:
①选取岩心样本烘干后,测量岩心样本的长度、直径、孔隙度和渗透率;对该岩心样本进行抽真空饱和度模拟地层水或者模拟地层原油,进行核磁共振扫描获取对应岩心样本的核磁共振成像MRI数据和核磁共振T2谱数据,扫描过程可借鉴核磁共振成像学技术实现(俎栋林,《核磁共振成像学》,2004,高等教育出版社)。对岩心不同位置进行MRI测量,获取对应位置的岩心MRI二维图像。根据仪器扫描精度,在仪器可处理范围内选择合适切片位置和断面数,获取岩心不同切片位置的核磁共振MRI二维图像。将切片位置坐标与像素数据保存于TXT文本中。
②获取岩心的核磁共振MRI/T2三维张量数据体。根据步骤①获取的二维图像TXT文本数据和岩心孔喉长度的大致范围对核磁共振成像数据体进行插值,使数据体规模满足反应岩石微观孔喉特征的要求。插值算法可采用三线性插值、克里金插值等算法。插值的基本参数由岩心的实际长度和直径,以及二维图像切片的空间位置确定,插值得到的MRI/T2三维张量数据体。
③构建空间无序结构孔隙网络模型。根据岩心尺度设置模型大小、配位数、平均孔喉长度等参数。首先采用C++语言和矩阵计算库Eigen构造规则立方体网络结构,生成一个X×Y×Z的三维规则立方体网络(X、Y和Z的值由步骤②所得的MRI/T2三维张量数据体规模确定),设置网络模型的总节点数为(X-1×(Y-1)×(Z-1),每个节点代表一个孔隙,节点与节点之间由喉道(圆形管道)相连,其余部分为岩石骨架。建立的网络模型中每个代表孔隙的节点周围都有六个喉道相连,喉道长度为岩石平均孔喉长度<l>。模型的x,y,z方向边长为分别Lx=(X-1)<l>,Ly=(Y-1)<l>,Lz=(Z-1)<l>。模型中所有网格节点之间均通过圆管实现全连接(此时任意节点配位数z=6),设置孔隙与喉道的半径之比为1。设置Ly和Lz为实际岩心的直径,并将每一层的yoz平面中,距离中心点大于Ly的点全部移除,使得模型成为与真实岩心形状一致的柱塞状模型。将核磁共振MRI/T2三维张量数据体中的数值赋予三维规则立方体网络的节点上,2个节点之间连线的值即为孔喉半径R。网络中所有孔喉半径R均可取相邻2节点上MRI/T2值的平均值。对网络各个节点坐标在球面内空间内进行随机移动,生成无序网络空间结构并产生孔喉长度随机变化。随机从网络结构中移除部分连接,可得到不同连通性(配位数)特征的孔隙网络模型。
④岩心数字化孔隙网络模型。岩心样本的核磁共振MRI/T2数据体反应的是岩石内部孔隙空间的相对大小值,而不直接反应岩石的孔喉半径大小。因此,这里可采用试错法估算岩石的孔喉半径Ri与对应MRI/T2数据幅值Ai之间的转换系数α(Ri=αAi),构建岩心数字化孔隙网络模型:假定转换系数的初始值,根据步骤③所述方法并通过转换系数计算孔喉半径Ri得到岩心数字化孔隙网络模型,采用单相稳定渗流孔隙网络模拟算法计算构建的岩心数字化孔隙网络模型的渗透率,检验模拟得到的渗透率与岩心测量的渗透率是否一致;若不一致则调整转换系数,重新计算孔隙网络孔喉半径Ri重新建立岩心数字化孔隙网络模型并计算其渗透率值,直到孔隙网络模型的渗透率与真实岩心的渗透率测量值基本一致,则得到与该实际岩心样本对应的岩心数字化孔隙网络模型。
对于直径2.5cm和长度为5cm的圆柱形岩心样本,对应的岩心数字化孔隙网络模型的节点数将大于100万。此时,常规的CPU稀疏矩阵方程求解算法处理此类问题已较为困难,应采用GPU算法进行计算和求解。
本发明将岩心数字化孔隙网络模型与流固耦合渗流数值模拟方法相结合,即可模拟再现室内岩心流动实验过程,实现室内岩心的数字化单相和多相(油-水,气-水,油-气)流固耦合渗流模拟分析与测试。本发明提出的流固耦合渗流数值模拟方法也可以与基于CT扫描的数字岩心/孔隙网络模型等相关方法相结合,进行多相渗流模拟和分析。
本发明提出的模拟技术通过隐式方法计算孔隙网络模型当前时刻的压力场分布,然后显式计算两相流体界面移动。本发明提出的岩心数字化孔隙网络模型的构建和多相渗流模拟方法可通过GPU加速计算技术提升模型尺度、精度和计算效率,实现基于同一渗流数学模型,从孔隙网络、室内岩心和物理模型实验、到井筒-单井油藏模型尺度的多尺度跨越。
本发明的有益效果:
1.分析了单个微观圆形通道中流体的非定常流动特征,证明泊肃叶公式和抛物线形态的速度剖面仍然适用于孔隙尺度非定常渗流过程;
2.单相液体渗流模型中,同时考虑了液体和岩石骨架的微可压缩性特征及其耦合作用关系,得到单相液体流固耦合渗流数学模型;单相气体渗流模型中,同时考虑岩石骨架的微可压缩性和气体较强的可压缩性,根据气体压缩和固体岩石骨架形变的耦合作用关系,得到了单相气体流固耦合渗流数学模型。进一步完善了孔隙尺度流固耦合渗流理论和方法;
3.孔隙网络模型中,单相气体非定常渗流与单相液体非定常渗流具有相似的数学表达形式,因此可进一步推导多相(油-水,气-水,油-气)渗流数学模型;在孔隙网络多相渗流数学模型中引入了流固耦合(考虑了液体和岩石骨架的微可压缩性,以及气体的可压缩性),得到非定常流固耦合多相渗流数学模型,并导出与数学模型对应的隐式数值模拟方法,完善了孔隙尺度流固耦合多相渗流理论和数值模拟方法;
4.将岩心的数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法相结合,可模拟再现室内岩心流动(驱替、自吸等)真实实验过程和结果,且数值模拟能跟岩心流动实验的磁磁共振MRI图像分析结果进行直接对比和分析,确保模拟结果的真实可靠;
5.通过GPU算法可大幅度提高模型的构建和多相渗流模拟速度。因此,可在同一岩石样品上开展不同驱替速度(或压力)下的渗流模拟研究,相比室内岩心实验具有周期短、速度快、操作方便、结果准确等优势,是室内岩心流动实验的重要补充;
6.部分室内岩心流动实验分析需要在高温高压下进行,操作较为复杂;采用本发明提出的建模和模拟方法只需根据温度和压力调整岩石骨架和流体相关参数,即可实现高温高压下的多相渗流模拟研究和分析,能够方便快速地获得岩心流动分析结果;
7.本发明提出的非定常流固耦合渗流数值模拟方法同样能与基于CT扫描的数字岩心/孔隙网络模型等相关方法相结合,进行多相渗流模拟和分析;
8.本发明引入不稳定渗流理论,将拉普拉斯方程转换为非稳态渗流-压力扩散方程,模拟流体流动过程,准确还原压力波对流体流动的影响。
附图说明
图1是本发明的技术流程图;
图2是核磁共振MRI扫描流程图;
图3是岩心不同位置的MRI二维图像示例图;
图4是基于核磁共振MRI/T2的三维张量数据可视化图像示意图;
图5是空间无序网络结构的构建方法示意图;
图6是岩心数字化孔隙网络模型示意图;
图7是数字化岩心模型中的两相驱替过程示例图一;
图8是数字化岩心模型中的两相驱替过程示例图二;
图9是数字化岩心模型中的两相驱替过程示例图三;
图10是数字化岩心模型中的两相驱替过程示例图四;
图11是数字化岩心模型中的两相驱替过程示例图五。
具体实施方式
为了对本发明的技术特征、目的和有益效果有更加清楚的理解,现对本发明的技术方案精选以下详细说明。显然,所描述的实施案例是本发明一部分实施例,而不是全部实施例,不能理解为对本发明可实施范围的限定。基于本发明的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的其他所有实施例,都属于本发明的保护范围。
实施例一:
本实施例中,如图1所示,一种非定常流固耦合多相渗流模型构建方法,包括:
S1,根据岩心尺度孔隙网络结合非定常流动模型分析流体在孔喉通道中的非定常流动,获得圆形管束中流体轴向速度vr分布,
S2,根据单相液体渗流过程满足的假设条件,采用有限体积法和雷诺运输方程来构建非定常单相液体流固耦合渗流数学模型;
S3,根据单相气体渗流过程中的气体的密度ρg和气体压缩系数Cg,结合雷诺运输方程构建同时满足低压和高压条件的非定常单相气体流固耦合渗流数学模型;
S4,利用非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型,结合孔隙网络模型中混合流体的特性参数构建非混相驱替过程下的流固耦合多相渗流数学模型;
S5,将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合,进行岩心流固耦合渗流模拟。
本实施例中,岩心流固耦合渗流模拟是基于岩心数字化孔隙网络模型进行的,因此本实施例基于核磁共振MRI/T2数据构建岩心数字化孔隙网络模型,该模型具体构建过程如下:
①根据核磁共振MRI/T2数据获取岩心二维切片图像
本发明核磁共振测量采用MesoMR23-60H-I型核磁共振仪,采用倒置脉冲序列和Carr-Purcell-Meiboom-Gill脉冲序列对核磁共振信号进行测量。利用SIRT(同时迭代重建技术)反演算法生成T2分布。采用磁场强度为0.5t的低场核磁共振岩心分析系统(MesoMR-060H-HTHP-I)进行测量,主要测试参数包括主导频率(21.326MHz)、回波间距(TE=0.2ms)、极化时间(TW=3000ms)和回波个数(NECH=8000)。扫描岩心过程中,将长为5cm直径为2.5cm的标准岩心置于恒定的磁场内,并在x,y,z三个方向均施加梯度场,梯度场强=梯度场两端的磁场强度差值/梯度场的长度。采集样本时,起初层面内共振一致,对磁场施加相位编码梯度,撤掉相位编码梯度,然后施加频率编码梯度,对每一个体素标定一个记号,这个过程称为编码或空间定位。对某一层面施加射频脉冲后,接收该层面的MR信号。再进行解码,得到该层面各个体素MR信号的大小,把体素信号的大小显示在荧光屏对应的像素上。具体核磁共振成像扫描流程如图2所示。信号大小用不同的灰度等级表示,信号大,像素亮度大,信号小,像素亮度小。对岩心不同部位进行切片处理,可以得到不同位置的岩心MRI二维图像。根据仪器扫描精度,在仪器可处理范围内选择合适位置和断面数(例如断面数为6,图3),获取如图4所示的岩心不同位置的核磁共振MRI图像。将切片位置坐标与像素体数据保存于TXT文本中。单张二维图片之内的像素点间隔以及两张二维图片之间的间隔取决于扫描设备的分辨率。
②获取岩心的核磁共振MRI/T2三维张量数据体
通常情况下,实际储层岩石孔喉长度l变化范围为50~300微米,平均孔喉长度l约100~150微米(Tang,Y.B.,Li,M.,Bernabé,Y.,&Zhao,J.Z.(2020).Viscous fingeringand preferential flow paths in heterogeneous porous media.Journal ofGeophysical Research:Solid Earth,125(3),e2019JB019306.)。根据岩心孔喉长度的大致范围对核磁共振成像数据体进行插值,使数据体规模满足反应岩石微观孔喉特征的要求。插值算法可采用三线性插值、克里金插值等算法。插值的基本参数根据岩心的实际长度和二维图像切片的空间位置确定,插值后得到的MRI/T2三维张量数据体基本满足建立岩石数字化孔隙网络模型的数据规模,如图4所示。
③空间无序结构孔隙网络模型的构建
根据岩心尺度设置模型大小、配位数、平均孔喉长度等参数。首先采用C++语言和矩阵计算库Eigen构造规则立方体网络结构,生成一个X×Y×Z的三维规则立方体网络,设置网络模型的总节点数为(X-1×(Y-1)×(Z-1),每个节点代表一个孔隙,节点与节点之间由喉道(均匀圆管道)相连(模型二维剖面如图6),其余部分由固体颗粒物质充满。由此建立起来的网络中每个代表孔隙的节点周围都有六个喉道相连,喉道长度设置为常数l=150μm,喉道半径设置为单位1,模型的x,y,z方向边长为分别Lx=(X-1)×l,Ly=(Y-1)×l,Lz=(Z-1)×l。记录网络模型中各节点的坐标。模型中所有网格节点之间均通过圆管实现全连接(此时任意节点配位数z=6),设置孔隙与喉道的半径之比为1。设置Ly和Lz,并将每一层的yoz平面中,距离中心点大于0.5Ly的的点全部移除,从而将模型设置为与真实岩心形状一致的柱塞状模型。将核磁共振MRI/T2三维张量数据体中的数值赋予三维规则立方体网络的节点上,2个节点之间连线的值即为孔喉半径。网络中所有孔喉半径的值可取相邻两节点上数值的平均值。然后,对网络各个节点坐标在球面内空间内进行随机移动,从而生成无序网络空间结构。如图5所示,随机从网络结构中移除部分连接,可得到不同连通性(配位数)特征的孔隙网络模型。
④岩心数字化孔隙网络模型
通过以上方法可得到基于核磁共振MRI/T2的岩心数字化孔隙网络模型(图6)。岩心样本的MRI数据体反应的是岩石内部孔隙空间的相对大小值,而不直接反应岩石的孔喉半径大小。因此,这里可采用试错法得到岩石的实际孔喉半径与MRI数据体之间的转换系数,估算岩石的孔喉半径值并用于构建孔隙网络模型:假定转换系数的初始值,通过转换系数计算孔喉半径Ri得到岩心数字化孔隙网络模型,采用单相稳定渗流孔隙网络模拟算法计算构建的岩心数字化孔隙网络模型的渗透率,检验模拟得到的渗透率与岩心测量的渗透率是否一致;若不一致则调整转换系数,重新计算孔隙网络孔喉半径Ri重新建立岩心数字化孔隙网络模型并计算其渗透率值,直到孔隙网络模型的渗透率与真实岩心的渗透率测量值基本一致,则得到与该实际岩心样本对应的岩心数字化孔隙网络模型。
孔隙网络模型的渗透率计算方法(Bernabé,Y.,Li,M.,Tang,Y.B.,&Evans,B.(2016).Pore space connectivity and the transport properties of rocks.Oil&GasScience and Technology–Revue d’IFP Energies nouvelles,71(4),50.):根据基尔霍夫定律,节点内流入和流出流体的流量之和为零。则流体在孔隙网络中流动满足拉普拉斯方程:g为水力传导率,p为压力。将拉普拉斯方程应用到孔隙网络模型中,即得到孔隙网络模型中流体稳态渗流时满足质量守恒定律:Σjqij=0。根据这一关系,遍历孔隙网络中的所有节点可构建出线性方程组或稀疏矩阵方程[A][X]=[B]。本研究采用共轭梯度求解上述矩阵方程,即可得到流体在网络模型中流动的流体压力场,然后通过入口端和出口端的压差计算模型中流体的流入和流出的流量,以及模型的渗透率。
本实施例中,基于上述岩心数字化孔隙网络模型和种非定常流固耦合多相渗流模型进行流固耦合渗流数值模拟,模拟过程中岩心数字化孔隙网络模型中饱和油,将水以恒定的速率Q从岩心左端面中心入口处注入模型中驱替模型中的油。设置模型右端出口端的压力为大气压0.1MPa,将上述条件带入式27中,采用本发明提出的非定常流固耦合多相渗流数学模型和数值模拟方法,进行岩心的水驱油模拟分析。图7~图11是本发明建立的数字化岩心模型中的两相驱替过程示例,表示的是t1~t5不同时刻注入水侵入占据孔喉空间图像,其中红色是驱替相流体侵入所占据的管道,被驱替相流体未显示。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解,本发明不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明的原理,在不脱离本发明精神和范围的前提下,本发明还会有各种变化和改进,这些变化和改进都落入要求保护的本发明范围内。本发明要求保护的范围由所附的权利要求书及其等效物界定。
Claims (3)
1.一种非定常流固耦合多相渗流模型构建方法,其特征在于,包括:
S1,根据岩心尺度孔隙网络结合非定常流动模型分析流体在孔喉通道中的非定常流动,获得圆形管束中流体轴向速度vr分布,
S2,根据单相液体渗流过程满足的假设条件,采用有限体积法和雷诺运输方程来构建非定常单相液体流固耦合渗流数学模型;
非定常单相液体流固耦合渗流数学模型如下列表达式所示:
式中,Ct=Cρ+Cp,Ct为综合压缩系数,Cp为孔隙压缩系数,Pa-1;Cρ为液体压缩系数,Pa-1;Δt为时间步长;gij为管束内液体的水力传导率;<p>ij为相邻节点i和j之间的平均压力;Δpij为相邻节点i和j之间管束的流体压力差;Vpi0为节点i在初始时刻的孔隙体积,R0ij为初始时刻节点i和j之间孔喉管束的半径;lij为相邻节点i和节点j之间的孔喉长度;n为与控制体中心节点i相连通的节点数;μ为流体粘度;
S3,根据单相气体渗流过程中的气体的密度ρg和气体压缩系数Cg,结合雷诺运输方程构建同时满足低压和高压条件的非定常单相气体流固耦合渗流数学模型;
非定常单相气体流固耦合渗流数学模型如下列表达式所示:
式中,ggij为管束中气体流动的水力传导率;μg为气体黏度;p为控制体中心节点处的孔隙流体压力,Rg为气体常数;
S4,利用非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型,结合孔隙网络模型中混合流体的特性参数构建非混相驱替过程下的流固耦合多相渗流数学模型;
流固耦合多相渗流数学模型如下列表达式所示:
Ct=Cp+CISI+CDSD
式中,gij为管束内流体的水力传导率,pcij为节点i和j之间孔喉通道的毛管力;μeff为节点i和j之间混合流体的有效粘度;CI和CD分别为注入流体和被驱替流体的压缩系数;SI和SD分别为以节点i为中心的控制体内注入流体和被驱替流体的饱和度;
S5,将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合,进行岩心流固耦合渗流模拟。
2.根据权利要求1所述的一种非定常流固耦合多相渗流模型构建方法,其特征在于,所述非定常流固耦合多相渗流数值模拟方法具体包括以下步骤:
S201,首先,对流固耦合多相渗流数学模型中的非线性毛管压力项进行线性化处理:
其中,gij *为线性化处理后的传导率;
从而获得线性化后的流固耦合多相渗流数学模型:
S202,采用隐式数值模拟方法,引入源汇项Qi,对线性化后的流固耦合多相渗流数学模型进行离散处理,得到:
式中,上标t表示当前时刻下的流动状态,上标t+Δt表示下一时刻的流动状态;
S203,对步骤7023公式中不同时刻的流动状态进行分离与合并,得到:
S204,遍历孔隙网络模型中的所有网络节点,孔隙网络模型中所有控制体均可将当前时刻和下一时刻的流动状态带入步骤S703的公式中,整理后形成如下矩阵:
[A]t+Δt[X]t+Δt=[B]t,
式中,[A]t+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵,N为孔隙网络模型的节点数,[X]t+Δt和[B]t是长度为N的两个向量,[X]t+Δt为下一时刻的压力场向量,[B]t为与上一时刻压力场和边界条件有关的向量,利用GPU代数多重网格广义最小残差算法求解以上矩阵获得当前时刻孔隙网络模型中的流体压力场分布;
S205,非定常流固耦合多相渗流数值模拟过程中,采用固定时间步长或变时间步长的方式使该步长内流体界面发生位移,每一时间步下,计算出界面移动后新的流体界面位置,并更新整个孔隙网络模型中所有孔喉通道的水力传导率和模型中各相流体的饱和度,再进行流体压力场分布求解,直至整个孔隙网络空间被侵入流体占满或达到某一饱和度数值时停止。
3.根据权利要求1所述的一种非定常流固耦合多相渗流模型构建方法,其特征在于,所述步骤S5还包括岩心数字化孔隙网络模型的构建过程,包括:
S301,选取岩心样本烘干后测量岩心孔喉长度、孔隙度和渗透率,并对岩心样本进行抽真空饱和度模拟地层水或模拟地层原油;
S302,对岩心样本进行核磁共振MRI/T2扫描测量,获取岩心在不同断面上的二维图像;
S303,对二维图像进行插值,获取岩心MRI/T2三维数据体;
S304,建立三维无序孔隙网络模型,将MRI/T2三维数据体赋值到三维无序孔隙网络模型的节点中;
S305,通过转换系数α计算三维无序孔隙网络模型中各相邻节点的孔喉半径,模拟计算岩心数字化孔隙网络模型的渗透率;调整转换系数α大小,确保孔隙网络渗透率与岩心实测渗透率相近,建立岩心数字化孔隙网络模型。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111662690.3A CN114386302B (zh) | 2021-12-31 | 2021-12-31 | 一种非定常流固耦合多相渗流模型构建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111662690.3A CN114386302B (zh) | 2021-12-31 | 2021-12-31 | 一种非定常流固耦合多相渗流模型构建方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114386302A CN114386302A (zh) | 2022-04-22 |
CN114386302B true CN114386302B (zh) | 2023-02-10 |
Family
ID=81199878
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111662690.3A Active CN114386302B (zh) | 2021-12-31 | 2021-12-31 | 一种非定常流固耦合多相渗流模型构建方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114386302B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117592387B (zh) * | 2023-05-25 | 2024-06-25 | 中国石油大学(北京) | 低渗致密油藏浸润调控渗流规律表征方法、装置及设备 |
CN117113873B (zh) * | 2023-08-15 | 2024-04-09 | 西南石油大学 | 一种多相流体地层渗流的数值模拟方法及应用 |
CN117057271B (zh) * | 2023-08-15 | 2024-03-01 | 西南石油大学 | 一种基于vof的多相流体渗流过程模拟方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109063334A (zh) * | 2018-08-03 | 2018-12-21 | 河海大学 | 一种低渗透多孔介质材料的流固耦合模型构建方法 |
CN113468829A (zh) * | 2021-06-24 | 2021-10-01 | 西南石油大学 | 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法 |
WO2021243765A1 (zh) * | 2020-06-05 | 2021-12-09 | 南京大学 | 基于gpu矩阵的离散元流固耦合数值模拟方法及系统 |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105574320B (zh) * | 2014-11-07 | 2018-02-13 | 中国石油化工股份有限公司 | 低渗砂岩储层有效渗流能力的评价方法 |
US11846183B2 (en) * | 2020-02-11 | 2023-12-19 | Halliburton Energy Services, Inc. | Determination of perforation crush zone permeability |
CN112098293B (zh) * | 2020-08-03 | 2021-06-01 | 西南石油大学 | 基于孔隙裂缝双重介质气藏非稳态气水两相渗流模拟方法 |
CN112082917B (zh) * | 2020-08-03 | 2021-04-30 | 西南石油大学 | 一种基于动态网络模拟的气水非稳态两相渗流模拟方法 |
CN112131803B (zh) * | 2020-08-15 | 2024-10-11 | 唐雁冰 | 基于gpu分布式并行计算的多尺度孔隙网络多相渗流模拟方法 |
CN113392528B (zh) * | 2021-06-17 | 2022-07-29 | 西南石油大学 | 基于数字岩心的气体-液硫两相流动规律模拟方法 |
-
2021
- 2021-12-31 CN CN202111662690.3A patent/CN114386302B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109063334A (zh) * | 2018-08-03 | 2018-12-21 | 河海大学 | 一种低渗透多孔介质材料的流固耦合模型构建方法 |
WO2021243765A1 (zh) * | 2020-06-05 | 2021-12-09 | 南京大学 | 基于gpu矩阵的离散元流固耦合数值模拟方法及系统 |
CN113468829A (zh) * | 2021-06-24 | 2021-10-01 | 西南石油大学 | 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法 |
Non-Patent Citations (1)
Title |
---|
数字岩心技术在致密储层微观渗流特征研究中的应用;盛军等;《西安石油大学学报(自然科学版)》;20180925(第05期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN114386302A (zh) | 2022-04-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114239367B (zh) | 一种室内岩心的数字化多相流固耦合渗流数值模拟方法 | |
CN114386302B (zh) | 一种非定常流固耦合多相渗流模型构建方法 | |
Xu et al. | Lattice B oltzmann simulation of immiscible two‐phase flow with capillary valve effect in porous media | |
Zeng et al. | Numerical simulation of proppant transport in propagating fractures with the multi-phase particle-in-cell method | |
Blunt et al. | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow | |
Bijeljic et al. | Insights into non‐Fickian solute transport in carbonates | |
Blunt | Flow in porous media—pore-network models and multiphase flow | |
Varloteaux et al. | Pore network modelling to determine the transport properties in presence of a reactive fluid: From pore to reservoir scale | |
Pasquier et al. | Modeling two-phase flow of immiscible fluids in porous media: Buckley-Leverett theory with explicit coupling terms | |
CN113468829B (zh) | 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法 | |
CN114283254B (zh) | 基于核磁共振数据的岩心数字化孔隙网络模型构建方法 | |
Fan et al. | Comprehensive study of the interactions between the critical dimensionless numbers associated with multiphase flow in 3D porous media | |
Zhao et al. | Spontaneous imbibition in a square tube with corner films: Theoretical model and numerical simulation | |
Muccino et al. | Toward an improved understanding of multiphase flow in porous media | |
Zhao et al. | Viscous dissipation and apparent permeability of gas flow in nanoporous media | |
Zhang et al. | High-order streamline simulation and macro-scale visualization experimental studies on waterflooding under given pressure boundaries | |
Pei et al. | Simulation of multiphase flow pattern, effective distance and filling ratio in hydraulic fracture | |
Fagbemi et al. | Numerical modeling of strongly coupled microscale multiphase flow and solid deformation | |
Ghasemi et al. | Near miscible relative permeability curves in layered porous media-investigations via diffuse interface Lattice Boltzmann method | |
Chen et al. | Relative permeability of porous media with nonuniform pores | |
Silin et al. | Predicting relative-permeability curves directly from rock images | |
Li et al. | Time-dependent shape factor for scaling of multi-dimensional counter-current imbibition with different boundary conditions | |
Decroux et al. | Computation of effective dynamic properties of naturally fractured reservoirs: Comparison and validation of methods | |
Yang | Multi-scale simulation of multiphase multi-component flow in porous media using the Lattice Boltzmann Method | |
De Souza et al. | Well-reservoir coupling on the numerical simulation of horizontal wells in gas reservoirs |
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 |