CN112082917B - 一种基于动态网络模拟的气水非稳态两相渗流模拟方法 - Google Patents
一种基于动态网络模拟的气水非稳态两相渗流模拟方法 Download PDFInfo
- Publication number
- CN112082917B CN112082917B CN202010769422.0A CN202010769422A CN112082917B CN 112082917 B CN112082917 B CN 112082917B CN 202010769422 A CN202010769422 A CN 202010769422A CN 112082917 B CN112082917 B CN 112082917B
- Authority
- CN
- China
- Prior art keywords
- pore
- gas
- simulation
- seepage
- phase
- 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
- 238000000034 method Methods 0.000 title claims abstract description 76
- 238000004088 simulation Methods 0.000 title claims abstract description 70
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 49
- 239000011148 porous material Substances 0.000 claims abstract description 139
- 239000012530 fluid Substances 0.000 claims abstract description 49
- 239000011435 rock Substances 0.000 claims abstract description 36
- 230000008569 process Effects 0.000 claims abstract description 30
- 238000001228 spectrum Methods 0.000 claims abstract description 18
- 238000005481 NMR spectroscopy Methods 0.000 claims abstract description 17
- 238000002474 experimental method Methods 0.000 claims abstract description 16
- 238000005315 distribution function Methods 0.000 claims abstract description 9
- 238000010603 microCT Methods 0.000 claims abstract description 6
- 230000000704 physical effect Effects 0.000 claims abstract description 6
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 5
- 230000033001 locomotion Effects 0.000 claims abstract description 4
- 238000009792 diffusion process Methods 0.000 claims abstract description 3
- 239000000523 sample Substances 0.000 claims description 27
- 238000006073 displacement reaction Methods 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 12
- 230000005499 meniscus Effects 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 9
- 238000006243 chemical reaction Methods 0.000 claims description 7
- 230000006835 compression Effects 0.000 claims description 6
- 238000007906 compression Methods 0.000 claims description 6
- 238000011478 gradient descent method Methods 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 5
- 238000001208 nuclear magnetic resonance pulse sequence Methods 0.000 claims description 4
- 229920006395 saturated elastomer Polymers 0.000 claims description 4
- 239000000126 substance Substances 0.000 claims description 4
- 238000012360 testing method Methods 0.000 claims description 4
- FAPWRFPIFSIZLT-UHFFFAOYSA-M Sodium chloride Chemical compound [Na+].[Cl-] FAPWRFPIFSIZLT-UHFFFAOYSA-M 0.000 claims description 3
- 238000004891 communication Methods 0.000 claims description 3
- 238000013170 computed tomography imaging Methods 0.000 claims description 3
- 238000007405 data analysis Methods 0.000 claims description 3
- 238000001035 drying Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 230000000149 penetrating effect Effects 0.000 claims description 3
- 150000003839 salts Chemical class 0.000 claims description 3
- 230000011218 segmentation Effects 0.000 claims description 3
- 239000011780 sodium chloride Substances 0.000 claims description 3
- 238000005406 washing Methods 0.000 claims description 3
- 238000012512 characterization method Methods 0.000 claims description 2
- 230000004069 differentiation Effects 0.000 claims description 2
- 238000000605 extraction Methods 0.000 claims description 2
- 238000009736 wetting Methods 0.000 claims description 2
- 238000005325 percolation Methods 0.000 claims 1
- 230000005514 two-phase flow Effects 0.000 abstract description 3
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 36
- 239000007789 gas Substances 0.000 description 32
- 239000003345 natural gas Substances 0.000 description 20
- 238000004458 analytical method Methods 0.000 description 4
- 238000002591 computed tomography Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- QSHDDOUJBYECFT-UHFFFAOYSA-N mercury Chemical compound [Hg] QSHDDOUJBYECFT-UHFFFAOYSA-N 0.000 description 2
- 229910052753 mercury Inorganic materials 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000035515 penetration Effects 0.000 description 2
- 230000021715 photosynthesis, light harvesting Effects 0.000 description 2
- 239000013598 vector Substances 0.000 description 2
- 239000011800 void material Substances 0.000 description 2
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000003889 chemical engineering Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000000875 corresponding effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 239000001257 hydrogen Substances 0.000 description 1
- 229910052739 hydrogen Inorganic materials 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000009545 invasion Effects 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000005311 nuclear magnetism Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
- G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
- G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
- G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N24/00—Investigating or analyzing materials by the use of nuclear magnetic resonance, electron paramagnetic resonance or other spin effects
- G01N24/08—Investigating or analyzing materials by the use of nuclear magnetic resonance, electron paramagnetic resonance or other spin effects by using nuclear magnetic resonance
- G01N24/081—Making measurements of geologic samples, e.g. measurements of moisture, pH, porosity, permeability, tortuosity or viscosity
Landscapes
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- General Physics & Mathematics (AREA)
- Pathology (AREA)
- Immunology (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- High Energy & Nuclear Physics (AREA)
- Theoretical Computer Science (AREA)
- Radiology & Medical Imaging (AREA)
- Pulmonology (AREA)
- Environmental & Geological Engineering (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geochemistry & Mineralogy (AREA)
- Geology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Dispersion Chemistry (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种基于动态网络模拟的气水非稳态两相渗流模拟方法，包括获取通过核磁共振实验得到的岩心T2谱，转换得到孔喉半径分布频率，并根据孔喉半径分布频率拟合得到孔喉分布函数；获取通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数，建立符合储层物性的孔隙网络模型；在建立的孔隙网络模组中引入非稳态渗流理论，根据流体流动、界面移动和孔隙压力扩散过程结合动态网络模拟算法和非稳态渗流理论对气水两相流体非稳态渗流和压力传播过程进行模拟。本发明在传统气水两相渗流过程中考虑了气水两相非稳态渗流这一特征，在模拟过程中考虑了流体的压缩性，且考虑了非稳态过程的动态网络模型能够更精确描述孔隙级气水两相流动过程。
Description
技术领域
本发明涉及油气田开发领域，尤其涉及一种基于动态网络模拟的气水非稳态两相渗流模拟方法。
背景技术
天然气时世界上清洁、用途广泛且地位十分重要的能源，研究岩石多孔介质孔隙空间内气水两相流动，在石油工程和化学工程等领域都具有重要的意义；因此，通过建立孔隙网络模型来对多孔介质内部孔隙空间的发育规模、空间分布对流体渗流的影响、流体在其中的分布规律、相互作用机理等一些决定流体在多孔介质中流动的宏观现象的本质问题展开研究显得十分重要。
目前，国内外许多学者采用连续介质理论对多孔介质多相渗漏进行了研究，但是由于毛管力、粘滞力在孔隙尺度具有非连续性，对孔隙尺度气水两相渗流机理缺乏认识；多孔介质岩石的气水两相渗流特性在很大程度上取决于岩石的空隙尺度的非均质性，但空隙尺度非均质性对气水两相渗流的影响尚不明确。传统的空隙网络模拟认为流入和流出节点的流量之和为零，认为流体能够瞬间从入口传递到出口，即稳态渗流模型，但是实际的流体和演示通常具有可压缩性，气体可压缩性更强，且流体压力无法瞬间传播到无穷远处，当多空介质尺度较大时，传统的动态网络模拟技术便无法准确用以描述压力传播的过程，也不能准确地描述两相流体渗流过程；因此，在孔隙网络模型中引入非稳态两相渗流理论，利用动态网络深入研究气水两相渗流规律就显得十分必要。
发明内容
本发明的目的在于克服现有技术的缺点，提供了一种基于动态网络模拟的气水非稳态两相渗流模拟方法，解决了现有动态网络模拟存在的不足。
本发明的目的通过以下技术方案来实现：一种基于动态网络模拟的气水非稳态两相渗流模拟方法，所述模拟方法包括：
获取通过核磁共振实验得到的岩心T2谱，根据预设定量关系将T2谱转换得到孔喉半径分布频率，并根据孔喉半径分布频率拟合得到孔喉分布函数；
获取通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数，根据所述孔喉分布函数、配位数和无序网络模型建立符合储层物性的孔隙网络模型；
在建立的孔隙网络模组中引入非稳态渗流理论，根据流体流动、界面移动和孔隙压力扩散过程结合动态网络模拟算法和非稳态渗流理论对气水两相流体非稳态渗流和压力传播过程进行模拟。
进一步地，所述通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数包括：
利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品、射线源和探测器进行360°的相对旋转，采集样品各个角度的数据；利用计算机断层扫描成像重构方法进行3D重构，得到样品内外结构的高分辨率3D数据和影像；
根据图像不同灰度进行物质区分实现CT数据分析，灰度值低的区域代表物质密度低，通过参考灰度曲线中孔隙的灰度值进行阈值划分，从而在图像中将孔隙进行单独；
在样品扫描模型中截取一定像素体积的研究区域，通过二值化分割对孔隙进行提取，计算出在当前分辨率下的孔隙占扫描样品总体的体积百分比，从而通过与物理实验对比得到建模所需要的孔隙度；通过对大数据量孔隙的连通性进行连通模拟，对连通孔隙进行识别和提取，剩余的孔隙为孤立的孔隙团，利用等效球直接统计非连通性孔隙；
利用最大球算法区分数字岩心三维图像中的孔隙、喉道所占空间及连通性提取相应的孔隙、喉道结构网络模型，同时通过数理统计方法对孔隙结构参数进行定量提取，得到研究岩石孔喉表征的参数；
通过球棒模型建立孔喉网络模型，统计表征参数并从中提取后续建模所需的平均孔喉长度和配位数。
进一步地，所述孔隙结构参数包括孔喉尺寸、孔喉体积、孔喉比、配位数和形状因子：所述表征参数包括孔喉半径、孔喉体积、形状因子、连通性或者配位数、以及每个与其连通的喉道特征。
进一步地，所述获取通过核磁共振实验得到的岩心T2谱，根据预设定量关系将T2谱转换得到孔喉半径分布频率包括：
对岩心进行洗油洗盐后在一定温度条件下烘干至重量不变化，并利用真空加压饱和仪以KCl2，盐水为介质将岩心饱和一定时间后进行核磁共振测量实验；
将准备好的岩心放入到磁体探头中调节设置好参数后，通过T2图像脉冲序列获取不同回波时间的系列T2图像，再将核磁共振T2谱转换为岩石孔喉半径频率分布曲线。
进一步地，所述预设定量关系为rm＝cT2m，其中，rm为第m个孔喉半径，T2m为T2谱的第m个幅度值，c为预设的转换系数，m为正整数；所述根据孔喉半径分布频率拟合得到孔喉分布函数包括：根据所述孔喉半径频率分布通过随机函数建立三维稳定随机场形成初始三维张量数据体，所述随机函数为对数正态分布随机函数
其中，通过拟合所述孔喉半径频率分布得到数学期望μ和标准偏差σ，x表示孔喉半径，也就是将这个分布函数的随机值赋值到网络模型中。
进一步地，所述对气水两相流体非稳态渗流和压力传播过程进行模拟包括气驱水非稳态两相渗流模拟和水驱气非稳态两相渗流模拟；所述气驱水非稳态两相渗流模拟包括：
根据基尔霍夫定律引入无流动边界条件得到各节点的流动压力，进而求得各截面的平均流速；在模拟过程中设定整体流动方向为水平方向，根据模型中对于单独的节点i和j，设两点压力分别为pi和pj，两点间连接孔道的半径和长度分别为rij和lij，气体粘度为μg，得到两节点间的流体总体积流量
在气体沾度为μg非湿相侵入前孔隙网络被粘度为μw的湿相流体占满，模拟驱替过程开始后侵入流体以一定的速率字孔隙网络左端注入，通过杨-拉普拉斯方程求解毛管压力pcij＝2γcosθ/rij；
选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，选取最小的Δti作为这一步运算总体的时间步长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间点的水力传导系数gij和两相在孔隙网络中的分布；
在实际渗流过程中，考虑流体和岩石的可压缩性，得到非稳态渗流方程
通过泰勒展开以及隐式有限差分法将非稳态渗流方程转换为矩阵方程
根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组
并采用梯度下降法对其求解。
进一步地，所述水驱气非稳态两相渗流模拟包括：
在水驱气非稳态两相渗流中此时毛管力pcij成为动力，两节点间的流体总体积流量qij变为
在气体粘度为μg非湿相侵入前孔隙网络被粘度为μw的湿相流体占满，模拟驱替过程开始后侵入流体以一定的速率字孔隙网络左端注入，通过杨-拉普拉斯方程求解毛管动力pcij＝2γcosθ/rij；
选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，选取最小的Δti作为这一步运算总体的时间步长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间点的水力传导系数gij和两相在孔隙网络中的分布；
在实际渗流过程中，考虑流体和岩石的可压缩性，得到非稳态渗流方程
通过泰勒展开以及隐式有限差分法将非稳态渗流方程转换为矩阵方程
根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组
并采用梯度下降法对其求解。
进一步地，所述模拟方法还包括在获取通过核磁共振试验得到的岩心T2谱之前在SC模型的基础上建立无序网络模型的步骤。
进一步地，所述在SC模型的基础上建立无序网络模型的步骤包括：
指定模型的节点数，构建一个X×Y×Z的三维简单立方体网格，网格中每个节点代表一个孔隙，节点之间通过喉道相连接；
确定网络模型的大小和节点数，计算出网格中个节点的坐标；
设定概率为为p的概率函数，通过随机数发生器确定x、y和z方向各相邻节点之间是否有管束连通；
通过中心点位移修正中心节点坐标生成随机网络
进一步地，所述通过随机数发生器建立连通概率函数确定x、y和z方向各相邻节点之间是否有管束连通包括：当渗透概率p为50％时，rand()函数随机生成的整数中有50％的概率小于50，有另外50％的概率大于50，可实现概率p＝50％的管束连接，即当产生了一个小于50的数时，表达式if(rand()％100＜p×100)为真，执行管束半径的分配任务；否则为假，不执行任何操作。
本发明的有益效果为：
1.本发明中涉及到的无序网络模型建立方法，在传统的SC模型中随机位移中心点坐标，结合了岩心分析技术中的CT扫描实验数据分析与核磁实验分析，考虑了孔喉结构等物性参数，建立起的模型更符合真是储层特征；
2.本发明中涉及到的模拟方法，在传统气水两相渗流过程中考虑了气水两相非稳态渗流这一特征，在模拟过程中考虑了流体的压缩性，且考虑了非稳态过程的动态网络模型能够更精确描述孔隙级气水两相流动过程；
3.本发明涉及到的孔隙尺度非稳态气水两相渗流模拟方法，可分析孔隙流体压力传播影响压力波与区域内流体的黏滞力和毛管力之间的平衡，可以辅助研究影响两相流体界面移动规律。
附图说明
图1为本发明方法的流程示意图；
图2为本发明方法的利用核磁T2谱获取孔喉分布示意图；
图3为气驱水示意图；
图4为水驱气示意图；
图5为考虑了岩心物性参数的无序结构网络模型图；
图6为气驱水压力波分布图；
图7为水驱气压力波分布图。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。因此，以下对附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的保护范围，而是仅仅表示本申请的选定实施例。基于本申请的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。下面结合附图对本发明做进一步的描述。
如图1所示，本发明涉及一种基于动态网络模型的气水非稳态两相渗流模拟方法，其包括建立无序网络模型步骤、微CT扫描实验统计配位数步骤、核磁共振岩心分析测试统计孔喉半径分布步骤和气水非稳态两相渗流模拟步骤。
S1、建立无序网络模型步骤
(1)指定模型的节点数，构建一个X×Y×Z的三维简单立方体网格。每个节点代表一个孔隙，节点与节点之间由喉道相连。由此建立起来的网络中每个代表孔隙的节点周围都有六个喉道相连；同理，每个喉道与相邻的两个孔隙相连。其各方向(即x方向、y方向和z方向)各节点之间的间隔距离设置为l，节点数设置为d，模型的边长为(d-1)×l。
(2)计算出网络模型中各节点的坐标。计算式为：(x，y，z)＝[(i-1)l，(j-1)l，(k-1)l]，式中i、j、k分别为x方向、y方向和z方向的节点序号，取值分别为1，2，3…。
(3)在程序中设定概率为p的概率函数，通过(伪)随机数发生器确定x方向各相邻节点之间是否有管束连通。采用rand()函数产生随机数，从而产生随机概率。具体的C/C++代码表达式为：if(rand()％100＜p×100)，式中，rand()％100—计算机随机生成0-99范围内的任意整数。
(4)通过(伪)随机数发生器建立连通概率函数，确定y、z方向各相邻节点之间是否有管束连接。方法与x方向的管束分配过程相同。当渗透概率p为50％时，rand()函数随机生成的整数中有50％的可能性(概率)小于50，有另外50％的可能性(概率)大于50。因此，该表达式可实现概率p＝50％的管束连接，即当产生了一个小于50的数时，表达式为真，执行管束半径的分配任务(分配管束半径为r)；否则为假，不执行任何操作。
(5)通过中心点位移中心节点坐标生成随机网络。如图5所示，为考虑了岩心物性参数的模型模拟图。
(6)本次建模储层孔喉半径呈对数正态分布(由核磁T2谱得，即步骤3)，通过Matlab拟合频率分布曲线，得到分布的平均值μ和标准偏差σ，得到对数正态分布随机函数：
其中，x>0。
S2、微CT扫描实验统计配位数
利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品和射线源及探测器进行360°的相对旋转，采集上千帧角度的数据，然后利用计算机断层扫描成像重构方法进行3D重构，从而得到样品内外结构的高分辨3D数据及影像，CT数据分析是通过图像不同灰度进行物质区分，灰度值低的区域代表物质密度低，参考灰度曲线中孔隙的灰度值进行阈值划分，在图像中将孔隙进行单独的提取分离。
在样品扫描模型中截取1000×1000×1000像素体积的研究区域，通过二值化分割对孔隙的提取，计算出在当前分辨率下的孔隙占扫描样品总体积的体积百分比，从而通过与物理实验对比得到建模所需要的孔隙度，通过计算机对大数据量孔隙的连通性进行连通模拟，对连通孔隙进行识别和提取，剩余的孔隙为孤立的孔隙团，利用等效球直径统计非连通性孔隙。
利用最大球算法，区分数字岩心三维图像中的孔隙、喉道所占空间及连通性，提取相应的孔隙、喉道结构网络模型，同时运用数理统计方法可实现对孔喉尺寸、孔喉体积、孔喉比、配位数、形状因子等孔隙结构的定量提取，得到研究岩石孔喉表征的参数。
再通过球棒模型建立孔喉网络模型，统计半径、体积、形状因子、连通性(配位数)及每个与其连通的喉道特征(喉道长度、形状因子)等表征参数，从中提取后续建模所需的平均孔喉长度和配位数。
具体为，基于计算机高分辨断层扫描成像技术(MicroCT对样品进行了扫描及数字岩心3D重构分别用等效球法和最大球法进行了孔网构建，并对储层进行了结构特征分析。其中，吼道长度可由下式计算：
L＝D-R1-R2
式中，R1，R2分别为该吼道所连接两个孔隙的半径，μm；D为两孔隙中心点的实际坐标距离，μm。配位数由软件自动统计得到。
S3、核磁共振岩心分布测试统计孔喉半径分布步骤
将从地层中采取的岩心进行洗油洗盐后，在温度为80℃的条件下充分烘干至重量不变化，利用真空加压饱和仪，以KCl2盐水为介质，将岩心饱和48小时后进行核磁共振测量实验，将准备好的岩心放入到磁体探头中，调节共振频率，选择T2 Image脉冲序列，设定系统参数、采集参数后，用T2图像脉冲序列获取不同回波时间系列T2图像，最后将核磁共振T2谱转换为岩石孔喉半径频率分布曲线。
具体为，如图2所示，核磁共振信号强度与饱和岩样内部的流体氢核数量成正相关。测量方法为：测量一组样品定标，定标样品的孔隙度值分别为2％、4％、8％、15％、30％，在得到每个定标样品核磁不同回波时间的一系列T2图像，以定标样品的单位体积质子密度图像信号为横坐标，以定标样的孔隙度为纵坐标，对数据做曲线拟合，获得核磁孔隙度刻度关系式：
其中，a，b为拟合参数、
为孔隙度、S为核磁图像信号大小、V为样品体积。然后测量目标岩心，将其测量图像信号值代入刻度关系式中，即可得到岩心孔隙度和整个孔隙度分布值。
同时，将核磁共振弛豫时间分布与常规岩石压汞孔径分布相结合，可以求得一个换算系数c(由压汞实验获取，其值具有地区经验性)值，将T2谱的横坐标乘以c即可获得孔径分布。
S4、气水非稳态两相渗流模拟方法
在真实岩心中，由于流体的粘性作用，流体质点粘附在物体表面上，形成流体不滑移现象(即相对速度为零)，因而产生摩擦阻力和能量耗散。因此假设孔隙网络中的流体流动遵循能量耗散最低原则，流动过程中孔隙网络模型遵循的质量守恒定律通过基尔霍夫定律描述，即流入的流体体积等于流出的流体体积，由此将真实岩心基质流动简化为孔隙网络模型流动，则可以在无序结构网络模型中进行气水两相非稳态渗流模拟。
(1)气驱水模拟
如图3和图6所示，根据基尔霍夫定律，引入无流动边界条件，解得各节点的流动压力，由此求得各截面的平均流速。模拟过程中整体流动方向设定为水平方向。该模型中对于单独的节点i和j，设两点压力分别为pi和pj，两点间连接孔道的半径和长度分别为rij和lij，气体粘度为μg；则两节点间的流体总体积流量qij为：
式中，gij为节点i和j之间管束的气体传导率，Bg为气体体积系数，Z和Zsc分别为地下和地面气体偏差因子，T和Tsc分别为地下和地面温度，<p>地下气体压力，psc为地面大气压(天然气物性参数详细计算过程见第二节)；μeff是单个管道中两相有效粘度(Pa·s)，在管道中只存在一个两相凹液面的情况下，有效粘度μeff可以用下式计算。
μeff＝Bgμgxij+μw(1-xij)
这里xij是与凹液面位置有关的无量纲数(0≤xij≤1)，即凹液面所在位置横坐标除以整个孔隙网络的长度，当管道中只有单相流体时，pc＝0。
在粘度为μg(Pa·s)非湿相侵入前，孔隙网络被粘度为μw(Pa·s)的湿相流体占满。模拟的驱替过程开始后，侵入流体以一定的速率自孔隙网络左端注入，毛管压力pc(MPa)使用杨-拉普拉斯方程求解：
pcij＝2γcosθ/rij
其中，γ是界面张力(N/m)，θ为润湿接触角。
在模拟过程中，选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，最小时间步长Δti(i＝1,2,…)是指在所有管道中的凹液面到达下一节点的时间中，选取最小的Δti作为这一步运算总体的时间步长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间点的水力传导系数gij和两相在孔隙网络中的分布。显然最小时间步长法中时间步长Δti与压力降Δp和凹液面位置有关，因此每次迭代的步长取决于该次计算的具体情况，而不是都相等。最小时间步长的引入使得此模型可以使用尽量少的迭代次数得出模拟结果，保证精度的同时大大提高了计算效率。
实际渗流过程中，气体为可压缩流体，岩石骨架具有微可压缩性。若考虑流体和岩石可压缩性，可得到非稳态渗流方程：
通过泰勒展开以及隐式有限差分法技术，可将上式转换为矩阵方程：
式中，Ct为综合压缩系数(以气体压缩系数为主)，Q为气体采出速度，Δt为时间步长。
在两相流中，每个节点的总体积流量依然满足守恒规律，即Σjqij＝0，据此可以构建线性方程组：
上述方程可以简化为矩阵求解式Ap＝B，其中A为主对角占优的对称稀疏矩阵，B为一组向量，p即为需要求解得到的全局压力场向量。上述矩阵可采用梯度下降法进行求解，梯度下降法如下：
f(p)＝Ap-B
f′(p)＝A
(2)水驱气模拟
如图4和图7所示，采用水驱气的方式进行模拟实验时，毛管力pcij成为动力，两节点间的流体总体积流量qij变为：
其后续处理方式与气驱水模拟一样。
在天然气关键物性参数计算中包括气藏温度、压力与相对密度，偏差因子，体积系数，等温压缩系数的计算。
1、气藏温度、压力与相对密度
地下的天然气是多种气体组分的混合产物，其温度与压力通常采用拟临界参数来处理：
ppc＝∑yipci，Tpc＝∑yiTci
式中，ppc,Tpc为天然气拟临界压力和拟临界温度；pci,Tci为气体组分i的临界压力和临界温度；yi为组分i的摩尔分数。
天然气的相对密度γh表示天然气密度ρg；与空气密度ρair之比。
因此，其拟临界压力与拟临界温度可以通过其相对密度求得：
通过当前天然气的压力p和温度T可以求得其视对比压力ppr与视对比温度Tpr：
2、偏差因子
天然气偏差因子z是定量描述真实气体(天然气)与理想气体偏差程度大小的系数,是天然气其他物性计算、天然气藏地质储量计算以及管道天然气产量设计的一项重要的参数。天然气偏差因子计算的方法很多，本发明采用Dranchuk和Abou－Kassem-11参数法计算偏差因子。
z＝0.27ppr/(ρprTpr)
且
其中A1＝0.3265；A2＝-1.07；A3＝-0.5339；A4＝0.01569；A5＝-0.05165；A6＝0.5475；A7＝-0.7361；A8＝0.1844；A9＝0.1056；A10＝0.6134；A11＝0.721，Tpr为给定条件下的视对比温度；Ppr为给定条件下的视对比压力；ρpr为中间参数，可用牛顿迭代法求取：
令原函数为：
其一阶导数为：
其K阶导数与K+1阶导数关系为：
设置迭代精度(本次为误差小于0.05％)满足需求即可求得偏差因子z。
3、体积系数
天然气的体积系数是天然气地下体积量V与天然气地面标准条件下体积量Vsc的比值，其公式为：
Bg＝V/Vsc
在油气藏条件下，压力为p，温度为T，将天然气状态方程、地面条件带入上式，可得到天然气体积系数计算公式：
Bg＝3.458×10-4zT/p
式中p、T为地层压力与温度，z为偏差因子。
4、等温压缩系数
天然气等温压缩系数Cg是指在等温条件下体积随压力变化的变化率，其数学表达式为：
考虑视对比压力后：
式中ppc为拟对比压力，ppr为视对比压力，z为偏差因子。
以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。
Claims (9)
1.一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述模拟方法包括：
获取通过核磁共振实验得到的岩心T2谱，根据预设定量关系将T2谱转换得到孔喉半径分布频率，并根据孔喉半径分布频率拟合得到孔喉分布函数；
获取通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数，根据所述孔喉分布函数、配位数和无序网络模型建立符合储层物性的孔隙网络模型；
在建立的孔隙网络模组中引入非稳态渗流理论，根据流体流动、界面移动和孔隙压力扩散过程结合动态网络模拟算法和非稳态渗流理论对气水两相流体非稳态渗流和压力传播过程进行模拟；
其中，所述对气水两相流体非稳态渗流和压力传播过程进行模拟包括气驱水非稳态两相渗流模拟和水驱气非稳态两相渗流模拟；所述气驱水非稳态两相渗流模拟包括：
根据基尔霍夫定律引入无流动边界条件得到各节点的流动压力，进而求得各截面的平均流速；在模拟过程中设定整体流动方向为水平方向，根据模型中对于单独的节点i和j，设两点压力分别为pi和pj，两点间连接孔道的半径和长度分别为rij和lij，气体粘度为μg，得到两节点间的流体总体积流量
其中，Δpij表示pi和pj的差值，μeff是单个管道中两相有效粘度；
在气体沾度为μg非湿相侵入前孔隙网络被粘度为μw的湿相流体占满，模拟驱替过程开始后侵入流体以一定的速率在孔隙网络左端注入，通过杨-拉普拉斯方程求解毛管压力pcij＝2γcosθ/rij，γ是界面张力，θ为润湿接触角；
选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，选取最小的Δti作为这一步运算总体的时间步长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，Δtmin表示最小的时间步长Δti，vij表示节点i到j的速度，由此可以求得该时间点的水力传导系数gij和两相在孔隙网络中的分布；
在实际渗流过程中，考虑流体和岩石的可压缩性，得到非稳态渗流方程
Ct为综合压缩系数；
通过泰勒展开以及隐式有限差分法将非稳态渗流方程转换为矩阵方程
Q为气体采出速度，Δt为时间步长；
根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组
并采用梯度下降法对其求解。
2.根据权利要求1所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数包括：
利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品、射线源和探测器进行360°的相对旋转，采集样品各个角度的数据；利用计算机断层扫描成像重构方法进行3D重构，得到样品内外结构的高分辨率3D数据和影像；
根据图像不同灰度进行物质区分实现CT数据分析，灰度值低的区域代表物质密度低，通过参考灰度曲线中孔隙的灰度值进行阈值划分，从而在图像中将孔隙进行单独的提取分离；
在样品扫描模型中截取一定像素体积的研究区域，通过二值化分割对孔隙进行提取，计算出在当前分辨率下的孔隙占扫描样品总体的体积百分比，从而通过与物理实验对比得到建模所需要的孔隙度；通过对大数据量孔隙的连通性进行连通模拟，对连通孔隙进行识别和提取，剩余的孔隙为孤立的孔隙团，利用等效球直接统计非连通性孔隙；
利用最大球算法区分数字岩心三维图像中的孔隙、喉道所占空间及连通性提取相应的孔隙、喉道结构网络模型，同时通过数理统计方法对孔隙结构参数进行定量提取，得到研究岩石孔喉表征的参数；
通过球棒模型建立孔喉网络模型，统计表征参数并从中提取后续建模所需的平均孔喉长度和配位数。
3.根据权利要求2所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述孔隙结构参数包括孔喉尺寸、孔喉体积、孔喉比、配位数和形状因子：所述表征参数包括孔喉半径、孔喉体积、形状因子、连通性或者配位数、每个与其连通的喉道特征。
4.根据权利要求1所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述获取通过核磁共振实验得到的岩心T2谱，根据预设定量关系将T2谱转换得到孔喉半径分布频率包括：
对岩心进行洗油洗盐后在一定温度条件下烘干至重量不变化，并利用真空加压饱和仪以KCl2，盐水为介质将岩心饱和一定时间后进行核磁共振测量实验；
将准备好的岩心放入到磁体探头中调节设置好参数后，通过T2图像脉冲序列获取不同回波时间的系列T2图像，再将核磁共振T2谱转换为岩石孔喉半径频率分布曲线。
5.根据权利要求4所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述预设定量关系为rm＝cT2m，其中，rm为第m个孔喉半径，T2m为T2谱的第m个幅度值，c为预设的转换系数，m为正整数；所述根据孔喉半径分布频率拟合得到孔喉分布函数包括：根据所述孔喉半径频率分布通过随机函数建立三维稳定随机场形成初始三维张量数据体，所述随机函数为对数正态分布随机函数
其中，通过拟合所述孔喉半径频率分布得到数学期望μ和标准偏差σ，x表示孔喉半径，也就是将这个分布函数的随机值赋值到网络模型中。
6.根据权利要求1所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述水驱气非稳态两相渗流模拟包括：
在水驱气非稳态两相渗流中此时毛管力pcij成为动力，两节点间的流体总体积流量qij变为
在气体粘度为μg非湿相侵入前，孔隙网络被粘度为μw的湿相流体占满，模拟驱替过程开始后侵入流体以一定的速率从孔隙网络左端注入，通过杨-拉普拉斯方程求解毛管动力pcij＝2γcosθ/rij；
选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，选取最小的Δti作为这一步运算总体的时间步长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间点的水力传导系数gij和两相在孔隙网络中的分布；
在实际渗流过程中，考虑流体和岩石的可压缩性，得到非稳态渗流方程
通过泰勒展开以及隐式有限差分法将非稳态渗流方程转换为矩阵方程
根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组
并采用梯度下降法对其求解。
7.根据权利要求1-6中任意一项所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述模拟方法还包括在获取通过核磁共振试验得到的岩心T2谱之前，在SC模型的基础上建立无序网络模型的步骤。
8.根据权利要求7所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述在SC模型的基础上建立无序网络模型的步骤包括：
指定模型的节点数，构建一个X×Y×Z的三维简单立方体网格，网格中每个节点代表一个孔隙，节点之间通过喉道相连接；
计算出网格中个节点的坐标；
设定概率为为p的概率函数，通过随机数发生器确定x、y和z方向各相邻节点之间是否有管束连通；
通过中心点位移修正中心节点坐标生成随机网络。
9.根据权利要求8所述的一种基于动态网络模拟的气水非稳态两相渗流模拟方法，其特征在于：所述通过随机数发生器确定x、y和z方向各相邻节点之间是否有管束连通包括：当逾渗概率p为50％时，rand()函数随机生成的整数中有50％的概率小于50，有另外50％的概率大于50，可实现概率p＝50％的管束连接，即当产生了一个小于50的数时，表达式if(rand()％100<p×100)为真，执行管束半径的分配任务；否则为假，不执行任何操作。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010769422.0A CN112082917B (zh) | 2020-08-03 | 2020-08-03 | 一种基于动态网络模拟的气水非稳态两相渗流模拟方法 |
| AU2021104836A AU2021104836A4 (en) | 2020-08-03 | 2021-08-02 | Method for simulation of gas-water unsteady two-phase seepage flow based on dynamic network simulation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010769422.0A CN112082917B (zh) | 2020-08-03 | 2020-08-03 | 一种基于动态网络模拟的气水非稳态两相渗流模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112082917A CN112082917A (zh) | 2020-12-15 |
| CN112082917B true CN112082917B (zh) | 2021-04-30 |
Family
ID=73736051
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010769422.0A Active CN112082917B (zh) | 2020-08-03 | 2020-08-03 | 一种基于动态网络模拟的气水非稳态两相渗流模拟方法 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN112082917B (zh) |
| AU (1) | AU2021104836A4 (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112881259A (zh) * | 2021-01-18 | 2021-06-01 | 山东科技大学 | 一种基于稳态法测节理网络气-水相对渗透率的可视化装置及方法 |
| CN113433157B (zh) * | 2021-06-24 | 2023-12-15 | 西南石油大学 | 基于核磁共振t2谱建立随机单元等效岩心模型的方法 |
| CN114386302B (zh) * | 2021-12-31 | 2023-02-10 | 西南石油大学 | 一种非定常流固耦合多相渗流模型构建方法 |
| CN114283254B (zh) * | 2021-12-31 | 2022-09-16 | 西南石油大学 | 基于核磁共振数据的岩心数字化孔隙网络模型构建方法 |
Family Cites Families (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8311788B2 (en) * | 2009-07-01 | 2012-11-13 | Schlumberger Technology Corporation | Method to quantify discrete pore shapes, volumes, and surface areas using confocal profilometry |
| US10551520B2 (en) * | 2014-11-13 | 2020-02-04 | Colorado School Of Mines | Surface relaxivity calculation using nuclear magnetic resonance (NMR) measurement, three dimensional (3D) rock model and NMR response simulation |
| CN106574981B (zh) * | 2015-08-17 | 2019-03-08 | 数岩科技（厦门）股份有限公司 | 针对多孔介质的核磁共振分析系统和方法 |
| CN106547938B (zh) * | 2015-11-09 | 2019-10-01 | 中国地质大学（北京） | 裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法 |
| CN106202695B (zh) * | 2016-07-07 | 2018-04-20 | 清能艾科（深圳）能源技术有限公司 | 一种采用数字岩心模拟计算岩心渗透率的方法 |
| CN107449707B (zh) * | 2017-07-03 | 2020-01-07 | 中国石油天然气股份有限公司 | 页岩储层中不同尺度孔隙定量的三维表征确定方法和装置 |
| CN110362842A (zh) * | 2018-04-09 | 2019-10-22 | 长江大学 | 基于多种形状孔喉的随机孔隙网络模型建模方法 |
| CN108876923A (zh) * | 2018-06-17 | 2018-11-23 | 西南石油大学 | 一种基于岩石微ct图像的三维孔隙尺度模型重建方法 |
| CN108918829B (zh) * | 2018-07-11 | 2021-11-02 | 中国石油天然气股份有限公司 | 一种基于形态学的模拟数字岩心微观变形方法及装置 |
| CN109242985B (zh) * | 2018-10-29 | 2020-06-05 | 中国科学院力学研究所 | 一种从三维图像确定孔隙结构关键参数的方法 |
| CN110853138B (zh) * | 2019-11-21 | 2023-08-18 | 科吉思石油技术咨询(北京)有限公司 | 双重介质碳酸盐岩孔隙-裂缝双重网络模型构建方法 |
-
2020
- 2020-08-03 CN CN202010769422.0A patent/CN112082917B/zh active Active
-
2021
- 2021-08-02 AU AU2021104836A patent/AU2021104836A4/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| AU2021104836A4 (en) | 2021-09-30 |
| CN112082917A (zh) | 2020-12-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112098293B (zh) | 基于孔隙裂缝双重介质气藏非稳态气水两相渗流模拟方法 | |
| CN112082917B (zh) | 一种基于动态网络模拟的气水非稳态两相渗流模拟方法 | |
| Prasser et al. | Bubble size measurement using wire-mesh sensors | |
| Prasser et al. | Evolution of the two-phase flow in a vertical tube—decomposition of gas fraction profiles according to bubble size classes using wire-mesh sensors | |
| Vogel et al. | A new approach for determining effective soil hydraulic functions | |
| CN113468829B (zh) | 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法 | |
| Zhang et al. | Pore scale simulation of liquid and gas two-phase flow based on digital core technology | |
| Ganapathisubramani et al. | Investigation of three-dimensional structure of fine scales in a turbulent jet by using cinematographic stereoscopic particle image velocimetry | |
| Suh et al. | Capillary pressure at irregularly shaped pore throats: Implications for water retention characteristics | |
| CN111624147A (zh) | 岩心的相对渗透率测定方法及装置 | |
| JP2023099406A (ja) | 室内岩石コア実験における、デジタル化された多相流体－構造連成（ｆｓｉ）浸透数値シミュレーション手法 | |
| CN114283254B (zh) | 基于核磁共振数据的岩心数字化孔隙网络模型构建方法 | |
| Bordoloi et al. | Rotational kinematics of large cylindrical particles in turbulence | |
| CN112577979B (zh) | 一种岩石内部流体饱和度空间分布的定量分析装置及方法 | |
| CN114386302B (zh) | 一种非定常流固耦合多相渗流模型构建方法 | |
| Dou et al. | Lattice Boltzmann simulation of solute transport in a single rough fracture | |
| CN110441209A (zh) | 一种基于致密储层数字岩心计算岩石渗透率的方法 | |
| Taylor et al. | Sub-particle-scale investigation of seepage in sands | |
| CN112179815B (zh) | 一种基于孔隙网络模型的单相非稳态渗流模型建立方法 | |
| Ju et al. | Prediction of preferential fluid flow in porous structures based on topological network models: Algorithm and experimental validation | |
| Chen et al. | Numerical study of a round buoyant jet under the effect of JONSWAP random waves | |
| Hu et al. | Analysing flow in rocks by combined positron emission tomography and computed tomography imaging | |
| Dakhelpour-Ghoveifel et al. | Prediction of gas-oil capillary pressure of carbonate rock using pore network modeling | |
| Cui et al. | Adaptive error modelling in MCMC sampling for large scale inverse problems | |
| Santos et al. | Bubble shape estimation in gas-liquid slug flow using wire-mesh sensor and advanced data processing |
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 | ||
| OL01 | Intention to license declared | ||
| OL01 | Intention to license declared |