CN114239367B - 一种室内岩心的数字化多相流固耦合渗流数值模拟方法 - Google Patents

一种室内岩心的数字化多相流固耦合渗流数值模拟方法 Download PDF

Info

Publication number
CN114239367B
CN114239367B CN202111664901.7A CN202111664901A CN114239367B CN 114239367 B CN114239367 B CN 114239367B CN 202111664901 A CN202111664901 A CN 202111664901A CN 114239367 B CN114239367 B CN 114239367B
Authority
CN
China
Prior art keywords
fluid
pore
core
seepage
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
CN202111664901.7A
Other languages
English (en)
Other versions
CN114239367A (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 CN202111664901.7A priority Critical patent/CN114239367B/zh
Priority to JP2022022016A priority patent/JP2023099406A/ja
Publication of CN114239367A publication Critical patent/CN114239367A/zh
Application granted granted Critical
Publication of CN114239367B publication Critical patent/CN114239367B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/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
    • 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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/12—Symbolic schematics
    • 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
    • Y—GENERAL 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
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

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地质封存技术。
背景技术
石油和天然气是保持国民经济高速发展的重要能源之一，如何合理开采石油和天然气并提高其采收率一直是油气田发发过程中的重要问题。实际地下储层岩石内部孔喉结构复杂，很难通过实验手段厘清流体在其中的渗流规律。许多研究人员利用多孔介质模型来模拟不同类型流体在岩石内部流动，从而找到有利于提高石油采收率的方法。数字岩心技术作为多孔介质模型的一个分支，可以用于石油与天然气行业的地质、地震、测井和开发以及提高采收率等各个领域。数字岩心能有效地保留岩心微观物理特征，能确保岩心可以无限次地被使用，是岩石物理实验数值模拟的重要平台，能定量研究岩心内部各种微观因素(如孔隙连通性、润湿性等)对储层渗流过程的影响，能够计算传统物理实验无法直接测量的物理性质，如油气水三相相对渗透率。鉴于其广泛的应用性，开展数字岩心的研究对提高石油和天然气采收率具有重要意义。
应用数字岩心进行各种岩石物理实验研究时最基础的工作就是建立一个准确的、与实际岩石相符的三维数字岩心模型。在过去的15年里，随着实验仪器的创新和新理论的突破，国内外研究团队不断提出新的构建数字岩心模型的方法。经过数年的研究，把构建数字岩心的方法分为三大类，分别是数值重建法、物理实验法和混合法。数值重建法是以少量的二维薄片图像为基础，利用二维图片中包含的信息，通过随机模拟法或沉积岩的过程模拟法来重建三维数字岩心的方法。该方法所建模型的准确度和建模效率较低，并且建模方法中约束条件的选择会使得模拟结果存在偶然性，难以还原储层真实岩心特征。物理实验法是指利用实验仪器对岩心样品拍摄或扫描以获取大量的岩心二维图片，然后通过建模程序或软件把二维图片叠加重构成三维数字岩心。但该方法受限于实验仪器(如CT扫描仪)的分辨率与精度，建立的模型尺度较小(一般为毫米尺度)，模型的代表性和工程应用受到较大限制，对于提取和分析具有溶洞裂缝特征岩心的微观参数较为困难，且物理实验成本高、周期长。混合法通过将多种建模方法结合，各取所长，能够建立较准确的三维模型，但所建立的模型尺度仍然与真实岩心存在差异。同时，受限于孔隙尺度渗流理论体系和数学模型的不完善，通过上述方法建立的模型中模拟的结果与实际实验结果存在不同程度的差异。近年来，随核磁共振成像技术(MRI)的快速进步和计算机GPU芯片算力的显著提升，为构建与室内实际岩石样品相对应的三维岩心数字化孔隙网络模型，完善孔隙网络模型多相渗流理论、数学模型和数值模拟方法，以及结合岩心数字化孔隙网络模型和孔隙网络模型多相渗流数值模拟方法开展岩心渗流模拟研究奠定了重要的物质基础。同时，也催生了岩心数字化孔隙网络模型和孔隙网络模型多相渗流数值模拟方法这一技术方法，对石油和天然气相关行业的发展具有重要意义。
发明内容
本发明的目的在于克服现有技术的不足，提供一种室内岩心的数字化多相流固耦合渗流数值模拟方法，将岩心数字化孔隙网络模型与流固耦合渗流数值模拟方法相结合，即可模拟再现室内岩心流动实验过程，实现室内岩心的数字化单相和多相(油-水，气-水，油-气)流固耦合渗流模拟分析与测试。
本发明的目的是通过以下技术方案来实现的：
一种室内岩心的数字化多相流固耦合渗流数值模拟方法，包括：
步骤一：选取岩心样本烘干后测量岩心孔喉长度、孔隙度和渗透率，并对岩心样本进行抽真空饱和度模拟地层水或模拟地层原油；
步骤二：对岩心样本进行核磁共振MRI/T2扫描测量，获取岩心在不同断面上的二维图像；
步骤三：对二维图像进行插值，获取岩心MRI/T2三维数据体；
步骤四：建立三维无序孔隙网络模型，将MRI/T2三维数据体赋值到三维无序孔隙网络模型的节点中；
步骤五：通过转换系数α计算三维无序孔隙网络模型中各相邻节点的孔喉半径，模拟计算岩心数字化孔隙网络模型的渗透率；调整转换系数α大小，确保孔隙网络渗透率与岩心实测渗透率相近，建立岩心数字化孔隙网络模型；
步骤六：将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合，进行岩心流固耦合渗流模拟。
具体的，步骤二具体包括：利用核磁共振仪器对抽真空饱和度模拟后的岩心样本进行核磁共振扫描，利用SIRT反演算法生成T2分布；对岩心不同部位进行切片处理，根据核磁共振仪器的扫描精度，在核磁共振仪器可处理范围内选择合适切片位置和断面数对岩心进行核磁共振扫描，获取岩心不同切片位置的岩心核磁共振MRI二维图像，将切片位置坐标与二维图像的像素数据保存于二维图像TXT文本中。
具体的，步骤三具体包括：利用插值算法结合步骤二中获取的二维图像TXT文本和步骤一中的岩心孔喉长度，对核磁共振成像数据体，即岩心核磁共振MRI二维图像进行插值处理，使数据体规模满足反应岩石微观孔喉特征的要求，得到关于岩心的MRI/T2三维张量数据体。
具体的，步骤四具体包括以下子步骤：
S401，根据岩心尺度设置模型大小、配位数和平均孔喉长度；
S402，采用计算机编程语言和矩阵计算库Eigen构造三维规则立方体网络结构，生成一个X×Y×Z的三维规则立方体网络；
S403，设置立方体网络的总节点数为(X-1×(Y-1)×(Z-1)，每个节点代表一个孔隙，节点与节点之间由喉道相连，其余部分为岩石骨架；
S404，建立的立方体网络模型中每个代表孔隙的节点周围都有六个喉道相连，喉道长度为岩石平均孔喉长度<l>；立方体网络模型的x，y，z方向边长为分别Lx＝(X-1)<l>，Ly＝(Y-1)<l>，Lz＝(Z-1)<l>；立方体网络模型中所有网格节点之间均通过圆管实现全连接，设置孔隙与喉道的半径之比为1；设置Ly和Lz为实际岩心的直径，并将每一层的yoz平面中，距离中心点大于Ly的点全部移除，使得立方体网络模型成为与真实岩心形状一致的柱塞状模型；
S405，将核磁共振MRI/T2三维张量数据体中的数值赋予三维规则立方体网络模型的各个节点，2个节点之间连线的值即为孔喉半径R；网络中所有孔喉半径R均取相邻2个节点上MRI/T2值的平均值；对模型中各个节点坐标在球面内空间内进行随机移动，生成无序网络空间结构并产生孔喉长度随机变化；随机从网络结构中移除部分连接，得到不同连通性特征的孔隙网络模型。
具体的，步骤五具体包括假定转换系数的初始值为α，根据步骤四的方法并通过转换系数α计算孔喉半径Ri得到岩心数字化孔隙网络模型，采用单相稳定渗流孔隙网络模拟算法计算构建的岩心数字化孔隙网络模型的渗透率，检验模拟得到的渗透率与岩心测量的渗透率是否一致；若不一致则调整转换系数α，重新计算孔隙网络孔喉半径Ri，重新建立岩心数字化孔隙网络模型并计算其渗透率值，直到岩心数字化孔隙网络模型的渗透率与真实岩心的渗透率测量值基本一致，则得到与该实际岩心样本对应的岩心数字化孔隙网络模型。
具体的，步骤六具体包括以下子步骤：
S601，根据岩心尺度孔隙网络结合非定常流动模型分析流体在孔喉通道中的非定常流动，获得圆形管束中流体轴向速度vr分布，
S602，根据单相液体渗流过程满足的假设条件，采用有限体积法和雷诺运输方程来构建非定常单相液体流固耦合渗流数学模型；
S603，根据单相气体渗流过程中的气体的密度ρg和气体压缩系数Cg，结合雷诺运输方程构建同时满足低压和高压条件的非定常单相气体流固耦合渗流数学模型；
S604，利用非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型，结合孔隙网络模型中混合流体的特性参数构建非混相驱替过程下的流固耦合多相渗流数学模型：
S605，将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合，进行岩心流固耦合渗流模拟。
具体的，非定常单相液体流固耦合渗流数学模型如下列表达式所示：
Figure BDA0003450821430000041
Figure BDA0003450821430000042
式中，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相连通的节点数；μ为流体粘度。
具体的，非定常单相气体流固耦合渗流数学模型如下列表达式所示：
Figure BDA0003450821430000043
Figure BDA0003450821430000044
式中，ggij为管束中气体流动的水力传导率；μg为气体黏度；p/Zμg为非线性项；p为控制体中心节点处的孔隙流体压力，Rg为气体常数。
具体的，流固耦合多相渗流数学模型如下列表达式所示：
Figure BDA0003450821430000045
Figure BDA0003450821430000046
Ct＝Cp+CISI+CDSD
式中，gij为管束内流体的水力传导率，pcij为节点i和j之间孔喉通道的毛管力；μeff为节点i和j之间混合流体的有效粘度；CI和CD分别为注入流体和被驱替流体的压缩系数；SI和SD分别为以节点i为中心的控制体内注入流体和被驱替流体的饱和度。
具体的，非定常流固耦合多相渗流数值模拟方法具体包括以下步骤：
S701，首先，对流固耦合多相渗流数学模型中的非线性毛管压力项进行线性化处理：
Figure BDA0003450821430000051
Figure BDA0003450821430000052
从而获得线性化后的流固耦合多相渗流数学模型：
Figure BDA0003450821430000053
S702，采用隐式数值模拟方法，引入源汇项Qi，对线性化后的流固耦合多相渗流数学模型进行离散处理，得到：
Figure BDA0003450821430000054
式中，上标t表示当前时刻下的流动状态，上标t+Δt表示下一时刻的流动状态；
S703，对步骤7023公式中不同时刻的流动状态进行分离与合并，得到：
Figure BDA0003450821430000055
S704，遍历孔隙网络模型中的所有网络节点，孔隙网络模型中所有控制体均可将当前时刻和下一时刻的流动状态带入步骤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代数多重网格广义最小残差算法求解以上矩阵获得当前时刻孔隙网络模型中的流体压力场分布；
S705，非定常流固耦合多相渗流数值模拟过程中，采用固定时间步长或变时间步长的方式使该步长内流体界面发生位移，每一时间步下，计算出界面移动后新的流体界面位置，并更新整个孔隙网络模型中所有孔喉通道的水力传导率和模型中各相流体的饱和度，再进行流体压力场分布求解，直至整个孔隙网络空间被侵入流体占满或达到某一饱和度数值时停止。
本发明涉及基于核磁共振MRI/T2的岩心数字化孔隙网络模型，非定常流固耦合渗流数学模型和模拟方法，因此本发明所提供的方法详细实现过程可具体分为以下两大部分流程：
1.基于核磁共振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算法进行计算和求解。
上述方法建立的岩心数字化孔隙网络模型能够将微观渗流机理应用于岩心尺度孔隙网络进行数值模拟研究，并能直接与室内岩心实验分析得到的宏观实验结果进行对比验证，还能够在室内岩心分析的基础上进行更大规模的升尺度分析，是室内岩心分析和宏观尺度油气藏数值模拟的重要补充，在室内岩心分析和宏观油气藏数值模拟之间架起一座桥梁。
2.适用于岩心数字化孔隙网络模型的非定常流固耦合渗流数学模型与数值模拟方法
在构建的岩心数字化孔隙网络模型基础上，为了准确描述多相流体流动过程以及流体-岩石骨架相互作用关系，通过流体力学雷诺输运方程、有限体积法和流固耦合渗流机理，推导出适用于岩心数字化孔隙网络模型的非定常流固耦合多相渗流数学模型，并提出相应的隐式数值模拟方法，形成一套基于核磁共振MRI/T2的岩心数字化建模和非定常流固耦合渗流模拟方法。
①流体在孔喉通道中的非定常流动
孔隙尺度的非定常流动模型是基于假设压力下降是由惯性和流体粘度引起的。假定流动为层流状态，由Navier-Stokes方程简化得到圆形管束内流体的轴向速度vr为:
Figure BDA0003450821430000071
根据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 ofMechanical Engineering Science,2010,224(6),1231-1243.)的研究，圆形管束中流体轴向速度vr分布为：
Figure BDA0003450821430000072
式中，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)液体和岩石骨架具有一定的可压缩性。实际渗流过程中，由于岩石和流体具有微可压缩性，使得流动过程压力的传播过程不稳定，产生压力波传播。以下是非定常单相液体流固耦合渗流数学模型的推导过程。
流体密度和岩石孔隙空间体积随压力变化满足(忽略孔喉长度随压力变化而变化)：
Figure BDA0003450821430000081
Figure BDA0003450821430000082
式中，ρ为流体密度，kg/m3；p为孔隙流体压力，Pa；R为管道半径；Vp为管道体积，m3；Cp为孔隙压缩系数，Pa-1；Cρ为液体压缩系数，Pa-1。上述式子分离变量后积分得到：
Figure BDA0003450821430000083
Figure BDA0003450821430000084
同理，孔喉半径R满足：
Figure BDA0003450821430000085
式中：Vp0，R0，ρ0分别为初始时刻(p0压力下)的孔隙体积、孔喉半径和液体密度。将式4按麦克劳林级数展开并取前两项可得：
V≈V0[1+Cp(p-p0)] (5c)
Figure BDA0003450821430000086
ρ≈ρ0[1+Cρ(p-p0)] (5c)
孔隙网络中，任意控制体内流体流动的质量守恒定律采用雷诺输运方程进行表示：
Figure BDA0003450821430000091
式中，V表示与任意节点为中心的控制体的体积，ρvdA表示任意时刻流体流入或流出控制体的质量流量。雷诺输运方程描述了在没有汇或源的情况下，通过控制面上净输出的该物理量的通量与控制体内流体质量的当地时间变化率。任意控制体内，
Figure BDA0003450821430000092
式中，根据式4a和4b，
ρVp≈ρ0Vp0+ρ0Vp0(Cρ+Cp)(p-p0)+ρ0Vp0CρCp(p-p0)2 (7b)由于Cp和Cρ均为很小的数，忽略式7b中的高阶项可得，
Figure BDA0003450821430000093
Figure BDA0003450821430000094
式中，lij为相邻节点i和节点j之间的孔喉长度。
对于控制体中任意方向孔喉通道内流体的流入和流出的净流量表示为：
Figure BDA0003450821430000095
式中，Aij为相邻节点i和节点j之间孔喉的横截面积。管道变形和压力波动会破坏稳态层流速度分布，但根据式2的分析结果可以发现，微观孔喉通道中流体流动的速度剖面会快速恢复至稳态抛物线速度分布。因此，任意孔喉通道的流量满足：
Figure BDA0003450821430000096
由式9a和9b可得，
Figure BDA0003450821430000101
式中，<p>ij为相邻节点i和j之间的平均压力,Δpij为相邻节点i和j之间管束的流体压力差(压差)。根据式6、8和10，可得到孔隙网络模型的非定常单相液体流动数学模型：
Figure BDA0003450821430000102
Figure BDA0003450821430000103
式中，Ct＝Cρ+Cp(Ct为综合压缩系数)，方程两端的流体密度被约去，Δt为时间步长。式11a中有单管流动方程qij＝gijΔpij(式中qij为管束内流体的体积流量，gij为管束内液体的水力传导率，)。式11是扩散方程，该方程描述了孔隙网络模型中流体流动过程中的压力扩散(传播)。当Ct＝Cρ+Cp＝0时，式11退化为基尔霍夫定律和拉普拉斯方程(稳态单相渗流数学模型)。
③非定常单相气体流固耦合渗流数学模型
对于真实气体，气体的密度ρg和气体压缩系数Cg计算式如下：
Figure BDA0003450821430000104
Figure BDA0003450821430000105
式中，M为气体分子摩尔质量，Z为气体偏差因子，p为控制体中心节点的压力，Rg为气体常数，T为开氏温度。气体在孔隙网络模型中的流动过程仍然满足雷诺输运方程(式6)。此时，式6的左边第一项可写成如下关系：
Figure BDA0003450821430000106
由式4、5和12，
Figure BDA0003450821430000107
其中，
Figure BDA0003450821430000108
同样地，相较于Cp，CpCρ是很小的数，忽略该项后，根据式14可得，
Figure BDA0003450821430000111
对于气体流动，雷诺输运方程(式6)左边第二项满足如下关系：
Figure BDA0003450821430000112
此时，气体流动穿过控制面的净流量方程如下：
Figure BDA0003450821430000113
联立式6、式15和式17，得到适用于孔隙网络的气体单相非定常流固耦合渗流数学模型通式：
Figure BDA0003450821430000114
对式18中非线性项p/Zμg的处理如下：低压条件(气体压力小于10MPa)下，假定气体为理想气体，此时气体偏差因子Z为1；气体黏度μg受压力影响较小，可近似认为μg为常数；采用数值解法时，式18中的压力p可将上一时间步的压力值直接带入计算下一时间步的压力场，此时式18两边的压力可约去。高压条件(气体压力大于10MPa)下，p/Zμg近似为常数。当满足以上两种条件时，时18可简化为：
Figure BDA0003450821430000115
Figure BDA0003450821430000116
式中，ggij为管束中气体流动的水力传导率。此时，气体压缩系数Cg远大于岩石的孔隙压缩系数Cp，因此上式右端项忽略了孔隙压缩系数。对比式11和式19可以发现，非定常单相液体流固耦合渗流(式11)和非定常单相气体流固耦合渗流(式18、19)在不同的约束条件下，具有相同的数学表达式。
④非定常流固耦合多相渗流数学模型
这里的多相是指，油-水两相，气-水两相和油-气两相。以下以油水两相驱替过程为例进行介绍。水驱油过程满足如下假设：1)任意孔喉中不同流体之间最多只有一个流体界面，2)流动过程中不发生混相，3)孔喉中发生活塞式驱替，4)忽略重力的影响。
驱替过程中，任意两相邻节点i和j之间的孔喉通道内，混合流体的密度计算式为：
Figure BDA0003450821430000121
任意以节点i为中心的控制体内，混合流体的密度计算式为：
Figure BDA0003450821430000122
式中，ρ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。孔隙网络模型中，任意控制体内混合流体的压缩性满足如下关系：
Figure BDA0003450821430000123
Figure BDA0003450821430000124
Ceffρ＝CISI+CDSD (22c)
式中，Ceffρ为混合流体的有效压缩系数，ρeff0为初始条件(压力p0)下混合流体的密度，CI和CD分别为注入流体和被驱替流体的压缩系数。非混相驱替过程中，任意时刻多相流体在孔隙网络模型中的流动过程均服从雷诺输运方程。与非定常单相流体(气体和液体)流固耦合渗流数学模型推导过程中采用的方法类似，可以得到非混相驱替过程下的流固耦合多相渗流数学模型：
Figure BDA0003450821430000131
Figure BDA0003450821430000132
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中非线性毛管压力项进行线性化处理：
Figure BDA0003450821430000133
Figure BDA0003450821430000134
从而将式23a改写为：
Figure BDA0003450821430000135
式25与式23完全等价，但数学模型的形式与非定常单相渗流数学模型(式11和19)更为接近。采用隐式数值模拟方法，引入源汇项Qi(表示某一节点内可能有外部流体注入，或有流体从孔隙网络模型中流出；否则Qi＝0)，对式25进行离散处理：
Figure BDA0003450821430000136
式中，上标t表示当前时刻下的流动状态(流体水力传导率和压力)，上标t+Δt表示下一时刻的流动状态。对式26不同时刻的流动状态进行分离与合并(将当前时刻的流动状态量放到方程右边，下一时刻的流动状态量放到左边)，得到：
Figure BDA0003450821430000141
遍历孔隙网络模型中的所有网络节点，孔隙网络模型中所有控制体均可将当前时刻和下一时刻的流动状态带入式26b中，从而可得到由N个与式26b相似的方程组成的方程组。根据方程组的下标，可将方程组整理为矩阵形式[A]t+Δt[X]t+Δt＝[B]t：
Figure BDA0003450821430000142
Figure BDA0003450821430000143
Figure BDA0003450821430000144
其中，
Figure BDA0003450821430000145
式中，[A]t+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵(N为孔隙网络模型的节点数)，[X]t+Δt和[B]t是长度为N的两个向量，[X]t+Δt为下一时刻的压力场向量，[B]t为与上一时刻压力场和边界条件有关的向量。求解以上矩阵方程即得到当前时刻孔隙网络模型中的流体压力场。
由于流动过程中两相界面随时间推进，必须选定一个时间步长使该步长内流体界面均发生了适量的位移Δx＝Δt(qij/Aij)(Aij为孔喉通道的横截面积)，这里的“适量”是指必须在保证精度的前提下尽量减少运算次数，这里可采用固定时间步长。固定时间步长是指将渗流瞬态划分为长度相等的时间间隔Δtk(k＝1,2,…)，当Δtk足够小时，相邻时间步长的压力场变化可以认为是稳定、线性的，此时可以使用与单相流部分相同的方法对压力场进行求解。
模拟过程也可以采用变时间步长的方法以加快数值模拟计算速度。驱替过程一直进行直到。每一时间步下，计算出界面移动后新的流体界面位置，并更新整个网络模整个孔隙网络空间被侵入流体占满或达到某一饱和度数值时停止型中所有孔喉通道的水力传导率和模型中各相流体的饱和度，再进行压力场求解。模拟过程中，流体界面发生移动后需要更新水力传导率，甚至发生流体界面“锁死”的情况，因此求解过程中会产生一个严重病态的矩阵(式27)，大大增加计算难度。模拟过程中使用AMGCL算法库中的GPU代数多重网格广义最小残差算法(AMG-GMRES)求解该矩阵。整个模拟过程保持注入流量Q不变。
本发明将岩心数字化孔隙网络模型与流固耦合渗流数值模拟方法相结合，即可模拟再现室内岩心流动实验过程，实现室内岩心的数字化单相和多相(油-水，气-水，油-气)流固耦合渗流模拟分析与测试。本发明提出的流固耦合渗流数值模拟方法也可以与基于CT扫描的数字岩心/孔隙网络模型等相关方法相结合，进行多相渗流模拟和分析。
本发明提出的模拟技术通过隐式方法计算孔隙网络模型当前时刻的压力场分布，然后显式计算两相流体界面移动。本发明提出的岩心数字化孔隙网络模型的构建和多相渗流模拟方法可通过GPU加速计算技术提升模型尺度、精度和计算效率，实现基于同一渗流数学模型，从孔隙网络、室内岩心和物理模型实验、到井筒-单井油藏模型尺度的多尺度跨越。
本发明的有益效果：
1.采用高精度核磁共振成像MRI和T2谱数据，通过插值算法与随机无序孔隙网络模型相结合，建立了室内岩石样品的数字化孔隙网络模型。该模型克服了早期数字岩心尺度较小、不利于开展多相渗流分析的问题，具有尺度大、精度高、且能与实际岩石样品大小和特征直接对应的优点，能更全面地分析岩石的微观孔喉非均质性和宏观非均质性对渗流过程的影响，辅助提高石油与天然气采收率的方案研究；
2.分析了单个微观圆形通道中流体的非定常流动特征，证明泊肃叶公式和抛物线形态的速度剖面仍然适用于孔隙尺度非定常渗流过程；
3.单相液体渗流模型中，同时考虑了液体和岩石骨架的微可压缩性特征及其耦合作用关系，得到单相液体流固耦合渗流数学模型；单相气体渗流模型中，同时考虑岩石骨架的微可压缩性和气体较强的可压缩性，根据气体压缩和固体岩石骨架形变的耦合作用关系，得到了单相气体流固耦合渗流数学模型。进一步完善了孔隙尺度流固耦合渗流理论和方法；
4.孔隙网络模型中，单相气体非定常渗流与单相液体非定常渗流具有相似的数学表达形式，因此可进一步推导多相(油-水，气-水，油-气)渗流数学模型；在孔隙网络多相渗流数学模型中引入了流固耦合(考虑了液体和岩石骨架的微可压缩性，以及气体的可压缩性)，得到非定常流固耦合多相渗流数学模型，并导出与数学模型对应的隐式数值模拟方法，完善了孔隙尺度流固耦合多相渗流理论和数值模拟方法；
5.将岩心的数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法相结合，可模拟再现室内岩心流动(驱替、自吸等)真实实验过程和结果，且数值模拟能跟岩心流动实验的磁磁共振MRI图像分析结果进行直接对比和分析，确保模拟结果的真实可靠；
6.通过GPU算法可大幅度提高模型的构建和多相渗流模拟速度。因此，可在同一岩石样品上开展不同驱替速度(或压力)下的渗流模拟研究，相比室内岩心实验具有周期短、速度快、操作方便、结果准确等优势，是室内岩心流动实验的重要补充；
7.部分室内岩心流动实验分析需要在高温高压下进行，操作较为复杂；采用本发明提出的建模和模拟方法只需根据温度和压力调整岩石骨架和流体相关参数，即可实现高温高压下的多相渗流模拟研究和分析，能够方便快速地获得岩心流动分析结果；
8.本发明提出的非定常流固耦合渗流数值模拟方法同样能与基于CT扫描的数字岩心/孔隙网络模型等相关方法相结合，进行多相渗流模拟和分析。
附图说明
图1是本发明的技术流程图；
图2是核磁共振MRI扫描流程图；
图3是岩心不同位置的MRI二维图像示例图；
图4是基于核磁共振MRI/T2的三维张量数据可视化图像示意图；
图5是空间无序网络结构的构建方法示意图；
图6是岩心数字化孔隙网络模型示意图；
图7是数字化岩心模型中的两相驱替过程示例图一；
图8是数字化岩心模型中的两相驱替过程示例图二；
图9是数字化岩心模型中的两相驱替过程示例图三；
图10是数字化岩心模型中的两相驱替过程示例图四；
图11是数字化岩心模型中的两相驱替过程示例图五。
具体实施方式
为了对本发明的技术特征、目的和有益效果有更加清楚的理解，现对本发明的技术方案精选以下详细说明。显然，所描述的实施案例是本发明一部分实施例，而不是全部实施例，不能理解为对本发明可实施范围的限定。基于本发明的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的其他所有实施例，都属于本发明的保护范围。
实施例一：
本实施例中，如图1所示，一种室内岩心的数字化多相流固耦合渗流数值模拟方法，技术流程包括：选取岩心样本烘干后测量岩心孔喉长度、孔隙度和渗透率，并对岩心样本进行抽真空饱和度模拟地层水或模拟地层原油；对岩心样本进行核磁共振MRI/T2扫描测量，获取岩心在不同断面上的二维图像；对二维图像进行插值，获取岩心MRI/T2三维数据体；建立三维无序孔隙网络模型，将MRI/T2三维数据体赋值到三维无序孔隙网络模型的节点中；通过转换系数α计算三维无序孔隙网络模型中各相邻节点的孔喉半径，模拟计算岩心数字化孔隙网络模型的渗透率；调整转换系数α大小，确保孔隙网络渗透率与岩心实测渗透率相近，建立岩心数字化孔隙网络模型；将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合，进行岩心流固耦合渗流模拟。
本实施例中，方法的具体实现过程如下：
1.基于核磁共振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微米(Bernabé,Y.,Li,M.,Tang,Y.B.,&Evans,B.,Pore space connectivity andthe transport properties ofrocks.Oil&Gas Science and Technology,2016,71(4),50)。根据岩心孔喉长度的大致范围对核磁共振成像数据体进行插值，使数据体规模满足反应岩石微观孔喉特征的要求。插值算法可采用三线性插值、克里金插值等算法。插值的基本参数根据岩心的实际长度和二维图像切片的空间位置确定，插值后得到的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重新建立岩心数字化孔隙网络模型并计算其渗透率值，直到孔隙网络模型的渗透率与真实岩心的渗透率测量值基本一致，则得到与该实际岩心样本对应的岩心数字化孔隙网络模型。
孔隙网络模型的渗透率计算方法(Tang,Y.B.,Li,M.,Bernabé,Y.,&Zhao,J.Z.,Viscous fingering and preferential flow paths in heterogeneous porousmedia.Journal of Geophysical Research:SolidEarth,2020,125(3),e2019JB019306)：根据基尔霍夫定律，节点内流入和流出流体的流量之和为零。则流体在孔隙网络中流动满足拉普拉斯方程：
Figure BDA0003450821430000181
g为水力传导率，p为压力。将拉普拉斯方程应用到孔隙网络模型中，即得到孔隙网络模型中流体稳态渗流时满足质量守恒定律：Σjqij＝0。根据这一关系，遍历孔隙网络中的所有节点可构建出线性方程组或稀疏矩阵方程[A][X]＝[B]。本研究采用共轭梯度求解上述矩阵方程，即可得到流体在网络模型中流动的流体压力场，然后通过入口端和出口端的压差计算模型中流体的流入和流出的流量，以及模型的渗透率。
2.基于岩心数字化孔隙网络模型的流固耦合渗流数值模拟
模拟过程中岩心数字化孔隙网络模型中饱和油，将水以恒定的速率Q从岩心左端面中心入口处注入模型中驱替模型中的油。设置模型右端出口端的压力为大气压0.1MPa，将上述条件带入式27中，采用本发明提出的非定常流固耦合多相渗流数学模型和数值模拟方法，进行岩心的水驱油模拟分析。图7～图11是本发明建立的数字化岩心模型中的两相驱替过程示例，表示的是t1～t5不同时刻注入水侵入占据孔喉空间图像，其中红色是驱替相流体侵入所占据的管道，被驱替相流体未显示。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护的范围由所附的权利要求书及其等效物界定。

Claims (9)

1.一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，包括：
步骤一：选取岩心样本烘干后测量岩心孔喉长度、孔隙度和渗透率，并对岩心样本进行抽真空饱和度模拟地层水或模拟地层原油；
步骤二：对岩心样本进行核磁共振
Figure 271465DEST_PATH_IMAGE002
扫描测量，获取岩心在不同断面上的二维图像；
步骤三：对二维图像进行插值，获取岩心
Figure 641136DEST_PATH_IMAGE002
三维数据体；
步骤四：建立三维无序孔隙网络模型，将
Figure 92977DEST_PATH_IMAGE002
三维数据体赋值到三维无序孔隙网络模型的节点中；
步骤五：通过转换系数
Figure 471874DEST_PATH_IMAGE004
计算三维无序孔隙网络模型中各相邻节点的孔喉半径，模拟计算岩心数字化孔隙网络模型的渗透率；调整转换系数
Figure 958351DEST_PATH_IMAGE004
大小，确保孔隙网络渗透率与岩心实测渗透率一致，建立岩心数字化孔隙网络模型；
步骤六：将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合，进行岩心流固耦合渗流模拟；
步骤六具体包括以下子步骤：
S601，根据岩心尺度孔隙网络结合非定常流动模型分析流体在孔喉通道中的非定常流动，获得圆形管束中流体轴向速度
Figure 954382DEST_PATH_IMAGE006
分布，
S602，根据单相液体渗流过程满足的假设条件，采用有限体积法和雷诺运输方程来构建非定常单相液体流固耦合渗流数学模型；
S603，根据单相气体渗流过程中的气体的密度
Figure 893519DEST_PATH_IMAGE008
和气体压缩系数
Figure 482632DEST_PATH_IMAGE010
，结合雷诺运输方程构建同时满足低压和高压条件的非定常单相气体流固耦合渗流数学模型；低压条件为气体压力小于 10MPa，所述高压条件为气体压力大于 10MPa；
S604，利用非定常单相液体流固耦合渗流数学模型和非定常单相气体流固耦合渗流数学模型，结合孔隙网络模型中混合流体的特性参数构建非混相驱替过程下的流固耦合多相渗流数学模型：
S605，将岩心数字化孔隙网络模型与非定常流固耦合多相渗流数值模拟方法结合，进行岩心流固耦合渗流模拟。
2.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述步骤二具体包括：利用核磁共振仪器对抽真空饱和度模拟后的岩心样本进行核磁共振扫描，利用SIRT反演算法生成
Figure 72882DEST_PATH_IMAGE012
分布；对岩心不同部位进行切片处理，根据核磁共振仪器的扫描精度，在核磁共振仪器可处理范围内选择合适切片位置和断面数对岩心进行核磁共振扫描，获取岩心不同切片位置的岩心核磁共振MRI二维图像，将切片位置坐标与二维图像的像素数据保存于二维图像TXT文本中。
3.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述步骤三具体包括：利用插值算法结合步骤二中获取的二维图像TXT文本和步骤一中的岩心孔喉长度，对核磁共振成像数据体，即岩心核磁共振MRI二维图像进行插值处理，使数据体规模满足反应岩石微观孔喉特征的要求，得到关于岩心的
Figure 719109DEST_PATH_IMAGE002
三维张量数据体。
4.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述步骤四具体包括以下子步骤：
S401，根据岩心尺度设置模型大小、配位数和平均孔喉长度；
S402，采用计算机编程语言和矩阵计算库Eigen构造三维规则立方体网络结构，生成一个X×Y×Z的三维规则立方体网络；
S403，设置立方体网络的总节点数为(X-1)×(Y-1)×(Z-1)，每个节点代表一个孔隙，节点与节点之间由喉道相连，其余部分为岩石骨架；
S404，建立的立方体网络模型中每个代表孔隙的节点周围都有六个喉道相连，喉道长度为岩石平均孔喉长度<l>；立方体网络模型的x，y，z方向边长为分别
Figure 207859DEST_PATH_IMAGE014
，
Figure 413712DEST_PATH_IMAGE016
，
Figure 405939DEST_PATH_IMAGE018
；立方体网络模型中所有网格节点之间均通过圆管实现全连接，设置孔隙与喉道的半径之比为1；设置
Figure 22734DEST_PATH_IMAGE020
和
Figure 733201DEST_PATH_IMAGE022
为实际岩心的直径，并将每一层的yoz平面中，距离中心点大于
Figure 477166DEST_PATH_IMAGE024
的点全部移除，使得立方体网络模型成为与真实岩心形状一致的柱塞状模型；
S405，将核磁共振
Figure 42009DEST_PATH_IMAGE002
三维张量数据体中的数值赋予三维规则立方体网络模型的各个节点，2个节点之间连线的值即为孔喉半径R；网络中所有孔喉半径R均取相邻2个节点上
Figure 314858DEST_PATH_IMAGE002
值的平均值；对模型中各个节点坐标在球面内空间内进行随机移动，生成无序网络空间结构并产生孔喉长度随机变化；随机从网络结构中移除部分连接，得到不同连通性特征的孔隙网络模型。
5.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述步骤五具体包括假定转换系数的初始值为
Figure 761889DEST_PATH_IMAGE004
，根据步骤四的方法并通过转换系数
Figure 309545DEST_PATH_IMAGE004
计算孔喉半径得到岩心数字化孔隙网络模型，采用单相稳定渗流孔隙网络模拟算法计算构建的岩心数字化孔隙网络模型的渗透率，检验模拟得到的渗透率与岩心测量的渗透率是否一致；若不一致则调整转换系数
Figure 541943DEST_PATH_IMAGE004
，重新计算孔隙网络孔喉半径
Figure 175574DEST_PATH_IMAGE026
，重新建立岩心数字化孔隙网络模型并计算其渗透率值，直到岩心数字化孔隙网络模型的渗透率与真实岩心的渗透率测量值一致，则得到与该实际岩心样本对应的岩心数字化孔隙网络模型。
6.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述非定常单相液体流固耦合渗流数学模型如下列表达式所示：
Figure 860633DEST_PATH_IMAGE028
Figure 461248DEST_PATH_IMAGE030
式中，
Figure 282573DEST_PATH_IMAGE032
，
Figure 897225DEST_PATH_IMAGE034
为综合压缩系数；
Figure 318848DEST_PATH_IMAGE036
为孔隙压缩系数，
Figure 411569DEST_PATH_IMAGE038
；
Figure 618560DEST_PATH_IMAGE040
为液体压缩系数，
Figure 450118DEST_PATH_IMAGE038
；Δt为时间步长；
Figure 844191DEST_PATH_IMAGE042
为管束内流体的体积流量，
Figure 537340DEST_PATH_IMAGE044
为管束内液体的水力传导率；
Figure DEST_PATH_IMAGE046
为相邻节点i和j之间的平均压力；
Figure DEST_PATH_IMAGE048
为相邻节点i和j之间管束的流体压力差；
Figure DEST_PATH_IMAGE050
为节点i在初始时刻的孔隙体积，
Figure DEST_PATH_IMAGE052
为初始时刻节点i和j之间孔喉管束的半径；
Figure DEST_PATH_IMAGE054
为相邻节点i和节点j之间的孔喉长度；n为与控制体中心节点i相连通的节点数；
Figure DEST_PATH_IMAGE056
为流体粘度。
7.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述非定常单相气体流固耦合渗流数学模型如下列表达式所示：
Figure DEST_PATH_IMAGE058
Figure DEST_PATH_IMAGE060
式中，
Figure DEST_PATH_IMAGE062
为管束中气体流动的水力传导率；
Figure DEST_PATH_IMAGE064
为气体黏度；p为控制体中心节点处的孔隙流体压力，
Figure DEST_PATH_IMAGE066
为气体常数，
Figure DEST_PATH_IMAGE068
为相邻节点i和j之间管束的流体压力差；
Figure DEST_PATH_IMAGE070
为节点i在初始时刻的孔隙体积，
Figure DEST_PATH_IMAGE072
为相邻节点i和j之间的平均压力，
Figure DEST_PATH_IMAGE074
为初始时刻节点i和j之间孔喉管束的半径；
Figure DEST_PATH_IMAGE076
为孔隙压缩系数，
Figure 956427DEST_PATH_IMAGE054
为相邻节点i和节点j之间的孔喉长度。
8.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述流固耦合多相渗流数学模型如下列表达式所示：
Figure DEST_PATH_IMAGE078
Figure DEST_PATH_IMAGE080
Figure DEST_PATH_IMAGE082
式中，
Figure DEST_PATH_IMAGE084
为管束内流体的水力传导率，
Figure DEST_PATH_IMAGE086
为节点i和j之间孔喉通道的毛管力；
Figure DEST_PATH_IMAGE088
为节点i和j之间混合流体的有效粘度；
Figure DEST_PATH_IMAGE090
和
Figure DEST_PATH_IMAGE092
分别为注入流体和被驱替流体的压缩系数；
Figure DEST_PATH_IMAGE094
和
Figure DEST_PATH_IMAGE096
分别为以节点i为中心的控制体内注入流体和被驱替流体的饱和度；
Figure 83521DEST_PATH_IMAGE072
为相邻节点i和j之间的平均压力；
Figure 500243DEST_PATH_IMAGE068
为相邻节点i和j之间管束的流体压力差；
Figure 997084DEST_PATH_IMAGE070
为节点i在初始时刻的孔隙体积，Δt为时间步长；
Figure 647508DEST_PATH_IMAGE074
为初始时刻节点i和j之间孔喉管束的半径；
Figure 758552DEST_PATH_IMAGE076
为孔隙压缩系数，
Figure 127217DEST_PATH_IMAGE054
为相邻节点i和节点j之间的孔喉长度。
9.根据权利要求1所述的一种室内岩心的数字化多相流固耦合渗流数值模拟方法，其特征在于，所述非定常流固耦合多相渗流数值模拟方法具体包括以下步骤：
S701，首先，对流固耦合多相渗流数学模型中的非线性毛管压力项进行线性化处理：
Figure DEST_PATH_IMAGE098
Figure DEST_PATH_IMAGE100
从而获得线性化后的流固耦合多相渗流数学模型：
Figure DEST_PATH_IMAGE102
式中，
Figure 83540DEST_PATH_IMAGE084
为管束内液体的水力传导率；
Figure 57313DEST_PATH_IMAGE068
为相邻节点 i 和 j 之间管束的流体压力差；
Figure 135996DEST_PATH_IMAGE086
为节点 i 和 j 之间孔喉通道的毛管力；
Figure 991957DEST_PATH_IMAGE034
为综合压缩系数；
Figure 299441DEST_PATH_IMAGE070
为节点 i 在初始时刻的孔隙体积；Δt为时间步长；
S702，采用隐式数值模拟方法，引入源汇项
Figure DEST_PATH_IMAGE104
，对线性化后的流固耦合多相渗流数学模型进行离散处理，得到：
Figure DEST_PATH_IMAGE106
式中，上标t表示当前时刻下的流动状态，上标t+Δt表示下一时刻的流动状态；
S703，对步骤7023公式中不同时刻的流动状态进行分离与合并，得到：
Figure DEST_PATH_IMAGE108
S704，遍历孔隙网络模型中的所有网络节点，孔隙网络模型中所有控制体均可将当前时刻和下一时刻的流动状态带入步骤S703的公式中，整理后形成如下矩阵：
Figure DEST_PATH_IMAGE110
，
式中，
Figure DEST_PATH_IMAGE112
为N×N大小的与流体水力传导率有关的稀疏矩阵，N为孔隙网络模型的节点数，
Figure DEST_PATH_IMAGE114
和
Figure DEST_PATH_IMAGE116
是长度为N的两个向量，
Figure 425922DEST_PATH_IMAGE114
为下一时刻的压力场向量，
Figure 691819DEST_PATH_IMAGE116
为与上一时刻压力场和边界条件有关的向量，利用GPU代数多重网格广义最小残差算法求解以上矩阵获得当前时刻孔隙网络模型中的流体压力场分布；
S705，非定常流固耦合多相渗流数值模拟过程中，采用固定时间步长或变时间步长的方式使该步长内流体界面发生位移，每一时间步下，计算出界面移动后新的流体界面位置，并更新整个孔隙网络模型中所有孔喉通道的水力传导率和模型中各相流体的饱和度，再进行流体压力场分布求解，直至整个孔隙网络空间被侵入流体占满或达到预定饱和度数值时停止。
CN202111664901.7A 2021-12-31 2021-12-31 一种室内岩心的数字化多相流固耦合渗流数值模拟方法 Active CN114239367B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202111664901.7A CN114239367B (zh) 2021-12-31 2021-12-31 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
JP2022022016A JP2023099406A (ja) 2021-12-31 2022-02-16 室内岩石コア実験における、デジタル化された多相流体－構造連成（ｆｓｉ）浸透数値シミュレーション手法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111664901.7A CN114239367B (zh) 2021-12-31 2021-12-31 一种室内岩心的数字化多相流固耦合渗流数值模拟方法

Publications (2)

Publication Number Publication Date
CN114239367A CN114239367A (zh) 2022-03-25
CN114239367B true CN114239367B (zh) 2022-12-30

Family

ID=80745068

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111664901.7A Active CN114239367B (zh) 2021-12-31 2021-12-31 一种室内岩心的数字化多相流固耦合渗流数值模拟方法

Country Status (2)

Country Link
JP (1) JP2023099406A (zh)
CN (1) CN114239367B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115221735B (zh) * 2022-09-20 2022-12-13 中国石油大学(华东) 一种电加热油页岩原位转化数值模拟的升尺度方法
CN116413187B (zh) * 2023-04-14 2024-01-05 西南石油大学 一种基于毛管压力资料的储层渗透率预测方法及系统
CN117592387B (zh) * 2023-05-25 2024-06-25 中国石油大学(北京) 低渗致密油藏浸润调控渗流规律表征方法、装置及设备
CN117115370B (zh) * 2023-08-15 2024-03-19 西南石油大学 一种高精度数字岩心模型构建方法
CN117057271B (zh) * 2023-08-15 2024-03-01 西南石油大学 一种基于vof的多相流体渗流过程模拟方法
CN117113873B (zh) * 2023-08-15 2024-04-09 西南石油大学 一种多相流体地层渗流的数值模拟方法及应用
CN118010774B (zh) * 2024-04-09 2024-05-31 北京科技大学 一种基于ct原位实验的页岩油加热改质流固界面作用表征的装置和方法
CN118150441B (zh) * 2024-05-11 2024-07-30 中国地质大学（北京） 一种基于流体流动特征评价岩石微观孔隙结构的方法
CN118410684B (zh) * 2024-06-28 2024-09-24 长安大学 用于岩土类材料渗透性评估模型的构建方法及应用

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103926267A (zh) * 2014-04-28 2014-07-16 西安石油大学 一种定量评价应力敏感过程中孔喉变化程度的方法
CN112326524A (zh) * 2020-10-22 2021-02-05 中国石油大学(华东) 一种基于ct扫描图像的岩石孔渗测量方法
WO2021243765A1 (zh) * 2020-06-05 2021-12-09 南京大学 基于gpu矩阵的离散元流固耦合数值模拟方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103926267A (zh) * 2014-04-28 2014-07-16 西安石油大学 一种定量评价应力敏感过程中孔喉变化程度的方法
WO2021243765A1 (zh) * 2020-06-05 2021-12-09 南京大学 基于gpu矩阵的离散元流固耦合数值模拟方法及系统
CN112326524A (zh) * 2020-10-22 2021-02-05 中国石油大学(华东) 一种基于ct扫描图像的岩石孔渗测量方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Modelling Fluid-Structure Interaction Problems with Coupled DEM-LBM;Rodrigo Guadarrama Lara;《https://etheses.whiterose.ac.uk/17444/》;20180725;全文 *
数字岩心建模方法研究现状及展望;林承焰 等;《地球物理学进》;20180430;第33卷(第2期);全文 *
流固化学耦合作用下裂隙岩体渗透特性研究进展;盛金昌 等;《岩土工程学报》;20110731;第33卷(第7期);全文 *

Also Published As

Publication number Publication date
JP2023099406A (ja) 2023-07-13
CN114239367A (zh) 2022-03-25

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
Abd et al. A review of the phenomenon of counter-current spontaneous imbibition: Analysis and data interpretation
Ferrari et al. Challenges in modeling unstable two‐phase flow experiments in porous micromodels
Blunt Flow in porous media—pore-network models and multiphase flow
Blunt et al. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow
Varloteaux et al. Pore network modelling to determine the transport properties in presence of a reactive fluid: From pore to reservoir scale
Bhattad et al. Effect of network structure on characterization and flow modeling using X-ray micro-tomography images of granular and fibrous porous media
CN114283254B (zh) 基于核磁共振数据的岩心数字化孔隙网络模型构建方法
Fan et al. Comprehensive study of the interactions between the critical dimensionless numbers associated with multiphase flow in 3D porous media
CN113468829A (zh) 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法
Liu et al. Improvement of predictions of petrophysical transport behavior using three-dimensional finite volume element model with micro-CT images
Zhao et al. Viscous dissipation and apparent permeability of gas flow in nanoporous media
Zhang et al. Evolution of the effective permeability for transient and pore-scale two-phase flow in real porous media
Mei et al. Fractal analysis of shape factor for matrix-fracture transfer function in fractured reservoirs
Liu et al. Investigation on Nonlinear Flow Behavior through Rock Rough Fractures Based on Experiments and Proposed 3‐Dimensional Numerical Simulation
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
Yang Multi-scale simulation of multiphase multi-component flow in porous media using the Lattice Boltzmann Method
Liang et al. Prediction of permeability from the skeleton of three-dimensional pore structure
Evseev et al. Coupling multiphase hydrodynamic and nmr pore-scale modeling for advanced characterization of saturated rocks
Lampe Modelling fluid flow and heat transport in fractured porous media
Sinha et al. A dynamic network simulator for immiscible two-phase flow in porous media
Pei et al. Anisotropic Characteristics of Relative Permeability in Sedimentary 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