本发明涉及一种基于动态网络模拟的气水非稳态两相渗流模拟方法，包括获取通过核磁共振实验得到的岩心T2谱，转换得到孔喉半径分布频率，并根据孔喉半径分布频率拟合得到孔喉分布函数；获取通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数，建立符合储层物性的孔隙网络模型；在建立的孔隙网络模组中引入非稳态渗流理论，根据流体流动、界面移动和孔隙压力扩散过程结合动态网络模拟算法和非稳态渗流理论对气水两相流体非稳态渗流和压力传播过程进行模拟。本发明在传统气水两相渗流过程中考虑了气水两相非稳态渗流这一特征，在模拟过程中考虑了流体的压缩性，且考虑了非稳态过程的动态网络模型能够更精确描述孔隙级气水两相流动过程。
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为气体
根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组并采用梯度下降法对其求解。
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
一种基于动态网络模拟的气水非稳态两相渗流模拟方法
技术领域
[0001] 本发明涉及油气田开发领域，尤其涉及一种基于动态网络模拟的气水非稳态两相渗流模拟方法。
背景技术
[0002] 天然气时世界上清洁、用途广泛且地位十分重要的能源，研究岩石多孔介质孔隙空间内气水两相流动，在石油工程和化学工程等领域都具有重要的意义；因此，通过建立孔
隙网络模型来对多孔介质内部孔隙空间的发育规模、空间分布对流体渗流的影响、流体在
其中的分布规律、相互作用机理等一些决定流体在多孔介质中流动的宏观现象的本质问题
展开研究显得十分重要。
[0003] 目前，国内外许多学者采用连续介质理论对多孔介质多相渗漏进行了研究，但是由于毛管力、粘滞力在孔隙尺度具有非连续性，对孔隙尺度气水两相渗流机理缺乏认识；多
孔介质岩石的气水两相渗流特性在很大程度上取决于岩石的空隙尺度的非均质性，但空隙
尺度非均质性对气水两相渗流的影响尚不明确。传统的空隙网络模拟认为流入和流出节点
的流量之和为零，认为流体能够瞬间从入口传递到出口，即稳态渗流模型，但是实际的流体
和演示通常具有可压缩性，气体可压缩性更强，且流体压力无法瞬间传播到无穷远处，当多
空介质尺度较大时，传统的动态网络模拟技术便无法准确用以描述压力传播的过程，也不
能准确地描述两相流体渗流过程；因此，在孔隙网络模型中引入非稳态两相渗流理论，利用
动态网络深入研究气水两相渗流规律就显得十分必要。
发明内容
[0004] 本发明的目的在于克服现有技术的缺点，提供了一种基于动态网络模拟的气水非稳态两相渗流模拟方法，解决了现有动态网络模拟存在的不足。
[0005] 本发明的目的通过以下技术方案来实现：一种基于动态网络模拟的气水非稳态两相渗流模拟方法，所述模拟方法包括：
[0006] 获取通过核磁共振实验得到的岩心T2谱，根据预设定量关系将T2谱转换得到孔喉半径分布频率，并根据孔喉半径分布频率拟合得到孔喉分布函数；
[0007] 获取通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数，根据所述孔喉分布函数、配位数和无序网络模型建立符合储层物性的孔隙网络模型；
[0008] 在建立的孔隙网络模组中引入非稳态渗流理论，根据流体流动、界面移动和孔隙压力扩散过程结合动态网络模拟算法和非稳态渗流理论对气水两相流体非稳态渗流和压
力传播过程进行模拟。
[0009] 进一步地，所述通过微CT扫描得到用于建立无序网络模型的岩心孔喉长度和配位数包括：
[0010] 利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品、射线源和探测器进行360°的相对旋转，采集样品各个角度的数据；利用计算机断层扫描成像
重构方法进行3D重构，得到样品内外结构的高分辨率3D数据和影像；
[0011] 根据图像不同灰度进行物质区分实现CT数据分析，灰度值低的区域代表物质密度低，通过参考灰度曲线中孔隙的灰度值进行阈值划分，从而在图像中将孔隙进行单独；
[0012] 在样品扫描模型中截取一定像素体积的研究区域，通过二值化分割对孔隙进行提取，计算出在当前分辨率下的孔隙占扫描样品总体的体积百分比，从而通过与物理实验对
比得到建模所需要的孔隙度；通过对大数据量孔隙的连通性进行连通模拟，对连通孔隙进
行识别和提取，剩余的孔隙为孤立的孔隙团，利用等效球直接统计非连通性孔隙；
[0013] 利用最大球算法区分数字岩心三维图像中的孔隙、喉道所占空间及连通性提取相应的孔隙、喉道结构网络模型，同时通过数理统计方法对孔隙结构参数进行定量提取，得到
研究岩石孔喉表征的参数；
[0014] 通过球棒模型建立孔喉网络模型，统计表征参数并从中提取后续建模所需的平均孔喉长度和配位数。
[0015] 进一步地，所述孔隙结构参数包括孔喉尺寸、孔喉体积、孔喉比、配位数和形状因子：所述表征参数包括孔喉半径、孔喉体积、形状因子、连通性或者配位数、以及每个与其连
通的喉道特征。
[0016] 进一步地，所述获取通过核磁共振实验得到的岩心T2谱，根据预设定量关系将T2谱转换得到孔喉半径分布频率包括：
[0017] 对岩心进行洗油洗盐后在一定温度条件下烘干至重量不变化，并利用真空加压饱和仪以KCl2，盐水为介质将岩心饱和一定时间后进行核磁共振测量实验；
[0018] 将准备好的岩心放入到磁体探头中调节设置好参数后，通过T2图像脉冲序列获取不同回波时间的系列T2图像，再将核磁共振T2谱转换为岩石孔喉半径频率分布曲线。
[0019] 进一步地，所述预设定量关系为rm＝cT2m，其中，rm为第m个孔喉半径，T2m为T2谱的第m个幅度值，c为预设的转换系数，m为正整数；所述根据孔喉半径分布频率拟合得到孔喉
分布函数包括：根据所述孔喉半径频率分布通过随机函数建立三维稳定随机场形成初始三
维张量数据体，所述随机函数为对数正态分布随机函数 其
中，通过拟合所述孔喉半径频率分布得到数学期望μ和标准偏差σ，x表示孔喉半径，也就是
将这个分布函数的随机值赋值到网络模型中。
[0020] 进一步地，所述对气水两相流体非稳态渗流和压力传播过程进行模拟包括气驱水非稳态两相渗流模拟和水驱气非稳态两相渗流模拟；所述气驱水非稳态两相渗流模拟包
括：
[0021] 根据基尔霍夫定律引入无流动边界条件得到各节点的流动压力，进而求得各截面的平均流速；在模拟过程中设定整体流动方向为水平方向，根据模型中对于单独的节点i和
j，设两点压力分别为pi和pj，两点间连接孔道的半径和长度分别为rij和lij，气体粘度为μg，
得到两节点间的流体总体积流量
[0022] 在气体沾度为μg非湿相侵入前孔隙网络被粘度为μw的湿相流体占满，模拟驱替过程开始后侵入流体以一定的速率字孔隙网络左端注入，通过杨‑拉普拉斯方程求解毛管压
力pcij＝2γcosθ/rij；
[0023] 选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，选取最小的Δti作为这一步运算总体的时间步长，此时除
到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间
点的水力传导系数gij和两相在孔隙网络中的分布；
[0024] 在实际渗流过程中，考虑流体和岩石的可压缩性，得到非稳态渗流方程
[0025] 通过泰勒展开以及隐式有限差分法将非稳态渗流方程转换为矩阵方程
[0026] 根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组并采用梯度下降法对其
求解。
[0027] 进一步地，所述水驱气非稳态两相渗流模拟包括：
[0028] 在水驱气非稳态两相渗流中此时毛管力pcij成为动力，两节点间的流体总体积流量qij变为
[0029] 在气体粘度为μg非湿相侵入前孔隙网络被粘度为μw的湿相流体占满，模拟驱替过程开始后侵入流体以一定的速率字孔隙网络左端注入，通过杨‑拉普拉斯方程求解毛管动
力pcij＝2γcosθ/rij；
[0030] 选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，选取最小的Δti作为这一步运算总体的时间步长，此时除
到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间
点的水力传导系数gij和两相在孔隙网络中的分布；
[0031] 在实际渗流过程中，考虑流体和岩石的可压缩性，得到非稳态渗流方程
[0032] 通过泰勒展开以及隐式有限差分法将非稳态渗流方程转换为矩阵方程
[0033] 根据两相渗流中每个节点的总体积流量满足守恒规律，构建线性方程组并采用梯度下降法对其
求解。
[0034] 进一步地，所述模拟方法还包括在获取通过核磁共振试验得到的岩心T2谱之前在SC模型的基础上建立无序网络模型的步骤。
[0035] 进一步地，所述在SC模型的基础上建立无序网络模型的步骤包括：
[0036] 指定模型的节点数，构建一个X×Y×Z的三维简单立方体网格，网格中每个节点代表一个孔隙，节点之间通过喉道相连接；
[0037] 确定网络模型的大小和节点数，计算出网格中个节点的坐标；
[0038] 设定概率为为p的概率函数，通过随机数发生器确定x、y和z方向各相邻节点之间是否有管束连通；
[0039] 通过中心点位移修正中心节点坐标生成随机网络
[0040] 进一步地，所述通过随机数发生器建立连通概率函数确定x、y和z方向各相邻节点之间是否有管束连通包括：当渗透概率p为50％时，rand()函数随机生成的整数中有50％
的概率小于50，有另外50％的概率大于50，可实现概率p＝50％的管束连接，即当产生了一
个小于50的数时，表达式if(rand()％100＜p×100)为真，执行管束半径的分配任务；否则
为假，不执行任何操作。
[0041] 本发明的有益效果为：
[0042] 1.本发明中涉及到的无序网络模型建立方法，在传统的SC模型中随机位移中心点坐标，结合了岩心分析技术中的CT扫描实验数据分析与核磁实验分析，考虑了孔喉结构等
物性参数，建立起的模型更符合真是储层特征；
[0043] 2.本发明中涉及到的模拟方法，在传统气水两相渗流过程中考虑了气水两相非稳态渗流这一特征，在模拟过程中考虑了流体的压缩性，且考虑了非稳态过程的动态网络模
型能够更精确描述孔隙级气水两相流动过程；
[0044] 3.本发明涉及到的孔隙尺度非稳态气水两相渗流模拟方法，可分析孔隙流体压力传播影响压力波与区域内流体的黏滞力和毛管力之间的平衡，可以辅助研究影响两相流体
界面移动规律。
附图说明
[0045] 图1为本发明方法的流程示意图；
[0046] 图2为本发明方法的利用核磁T2谱获取孔喉分布示意图；
[0047] 图3为气驱水示意图；
[0048] 图4为水驱气示意图；
[0049] 图5为考虑了岩心物性参数的无序结构网络模型图；
[0050] 图6为气驱水压力波分布图；
[0051] 图7为水驱气压力波分布图。
具体实施方式
[0052] 为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅
是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实
施例的组件可以以各种不同的配置来布置和设计。因此，以下对附图中提供的本申请的实
施例的详细描述并非旨在限制要求保护的本申请的保护范围，而是仅仅表示本申请的选定
实施例。基于本申请的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的
所有其他实施例，都属于本申请保护的范围。下面结合附图对本发明做进一步的描述。
[0053] 如图1所示，本发明涉及一种基于动态网络模型的气水非稳态两相渗流模拟方法，其包括建立无序网络模型步骤、微CT扫描实验统计配位数步骤、核磁共振岩心分析测试统
计孔喉半径分布步骤和气水非稳态两相渗流模拟步骤。
[0054] S1、建立无序网络模型步骤
[0055] (1)指定模型的节点数，构建一个X×Y×Z的三维简单立方体网格。每个节点代表一个孔隙，节点与节点之间由喉道相连。由此建立起来的网络中每个代表孔隙的节点周围
都有六个喉道相连；同理，每个喉道与相邻的两个孔隙相连。其各方向(即x方向、y方向和z
方向)各节点之间的间隔距离设置为l，节点数设置为d，模型的边长为(d‑1)×l。
[0056] (2)计算出网络模型中各节点的坐标。计算式为：(x，y，z)＝[(i‑1)l，(j‑1)l，(k‑1)l]，式中i、j、k分别为x方向、y方向和z方向的节点序号，取值分别为1，2，3…。
[0057] (3)在程序中设定概率为p的概率函数，通过(伪)随机数发生器确定x方向各相邻节点之间是否有管束连通。采用rand()函数产生随机数，从而产生随机概率。具体的C/C++
代码表达式为：if(rand()％100＜p×100)，式中，rand()％100—计算机随机生成0‑99范
围内的任意整数。
[0058] (4)通过(伪)随机数发生器建立连通概率函数，确定y、z方向各相邻节点之间是否有管束连接。方法与x方向的管束分配过程相同。当渗透概率p为50％时，rand()函数随机
生成的整数中有50％的可能性(概率)小于50，有另外50％的可能性(概率)大于50。因此，该
表达式可实现概率p＝50％的管束连接，即当产生了一个小于50的数时，表达式为真，执行
管束半径的分配任务(分配管束半径为r)；否则为假，不执行任何操作。
[0059] (5)通过中心点位移中心节点坐标生成随机网络。如图5所示，为考虑了岩心物性参数的模型模拟图。
[0060] (6)本次建模储层孔喉半径呈对数正态分布(由核磁T2谱得，即步骤3)，通过Matlab拟合频率分布曲线，得到分布的平均值μ和标准偏差σ，得到对数正态分布随机函数：
[0061]
[0062] 其中，x>0。
[0063] S2、微CT扫描实验统计配位数
[0064] 利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品和射线源及探测器进行360°的相对旋转，采集上千帧角度的数据，然后利用计算机断层扫描
成像重构方法进行3D重构，从而得到样品内外结构的高分辨3D数据及影像，CT数据分析是
通过图像不同灰度进行物质区分，灰度值低的区域代表物质密度低，参考灰度曲线中孔隙
的灰度值进行阈值划分，在图像中将孔隙进行单独的提取分离。
[0065] 在样品扫描模型中截取1000×1000×1000像素体积的研究区域，通过二值化分割对孔隙的提取，计算出在当前分辨率下的孔隙占扫描样品总体积的体积百分比，从而通过
与物理实验对比得到建模所需要的孔隙度，通过计算机对大数据量孔隙的连通性进行连通
模拟，对连通孔隙进行识别和提取，剩余的孔隙为孤立的孔隙团，利用等效球直径统计非连
通性孔隙。
[0066] 利用最大球算法，区分数字岩心三维图像中的孔隙、喉道所占空间及连通性，提取相应的孔隙、喉道结构网络模型，同时运用数理统计方法可实现对孔喉尺寸、孔喉体积、孔
喉比、配位数、形状因子等孔隙结构的定量提取，得到研究岩石孔喉表征的参数。
[0067] 再通过球棒模型建立孔喉网络模型，统计半径、体积、形状因子、连通性(配位数)及每个与其连通的喉道特征(喉道长度、形状因子)等表征参数，从中提取后续建模所需的
平均孔喉长度和配位数。
[0068] 具体为，基于计算机高分辨断层扫描成像技术(MicroCT对样品进行了扫描及数字岩心3D重构分别用等效球法和最大球法进行了孔网构建，并对储层进行了结构特征分析。
其中，吼道长度可由下式计算：
[0069] L＝D‑R1‑R2
[0070] 式中，R1，R2分别为该吼道所连接两个孔隙的半径，μm；D为两孔隙中心点的实际坐标距离，μm。配位数由软件自动统计得到。
[0071] S3、核磁共振岩心分布测试统计孔喉半径分布步骤
[0072] 将从地层中采取的岩心进行洗油洗盐后，在温度为80℃的条件下充分烘干至重量不变化，利用真空加压饱和仪，以KCl2盐水为介质，将岩心饱和48小时后进行核磁共振测量
实验，将准备好的岩心放入到磁体探头中，调节共振频率，选择T2 Image脉冲序列，设定系
统参数、采集参数后，用T2图像脉冲序列获取不同回波时间系列T2图像，最后将核磁共振T2
谱转换为岩石孔喉半径频率分布曲线。
[0073] 具体为，如图2所示，核磁共振信号强度与饱和岩样内部的流体氢核数量成正相关。测量方法为：测量一组样品定标，定标样品的孔隙度值分别为2％、4％、8％、15％、30％，
在得到每个定标样品核磁不同回波时间的一系列T2图像，以定标样品的单位体积质子密度
图像信号为横坐标，以定标样的孔隙度为纵坐标，对数据做曲线拟合，获得核磁孔隙度刻度
关系式：
[0074]
[0075] 其中，a，b为拟合参数、 为孔隙度、S为核磁图像信号大小、V为样品体积。然后测量目标岩心，将其测量图像信号值代入刻度关系式中，即可得到岩心孔隙度和整个孔隙度
分布值。
[0076] 同时，将核磁共振弛豫时间分布与常规岩石压汞孔径分布相结合，可以求得一个换算系数c(由压汞实验获取，其值具有地区经验性)值，将T2谱的横坐标乘以c即可获得孔
径分布。
[0077] S4、气水非稳态两相渗流模拟方法
[0078] 在真实岩心中，由于流体的粘性作用，流体质点粘附在物体表面上，形成流体不滑移现象(即相对速度为零)，因而产生摩擦阻力和能量耗散。因此假设孔隙网络中的流体流
动遵循能量耗散最低原则，流动过程中孔隙网络模型遵循的质量守恒定律通过基尔霍夫定
律描述，即流入的流体体积等于流出的流体体积，由此将真实岩心基质流动简化为孔隙网
络模型流动，则可以在无序结构网络模型中进行气水两相非稳态渗流模拟。
[0079] (1)气驱水模拟
[0080] 如图3和图6所示，根据基尔霍夫定律，引入无流动边界条件，解得各节点的流动压力，由此求得各截面的平均流速。模拟过程中整体流动方向设定为水平方向。该模型中对于
单独的节点i和j，设两点压力分别为pi和pj，两点间连接孔道的半径和长度分别为rij和lij，
气体粘度为μg；则两节点间的流体总体积流量qij为：
[0081]
[0082]
[0083]
[0084] 式中，gij为节点i和j之间管束的气体传导率，Bg为气体体积系数，Z和Zsc分别为地下和地面气体偏差因子，T和Tsc分别为地下和地面温度，
地下气体压力，psc为地面大气压(天然气物性参数详细计算过程见第二节)；μeff是单个管道中两相有效粘度(Pa·s)，在
管道中只存在一个两相凹液面的情况下，有效粘度μeff可以用下式计算。
[0085] μeff＝Bgμgxij+μw(1‑xij)
[0086] 这里xij是与凹液面位置有关的无量纲数(0≤xij≤1)，即凹液面所在位置横坐标除以整个孔隙网络的长度，当管道中只有单相流体时，pc＝0。
[0087] 在粘度为μg(Pa·s)非湿相侵入前，孔隙网络被粘度为μw(Pa·s)的湿相流体占满。模拟的驱替过程开始后，侵入流体以一定的速率自孔隙网络左端注入，毛管压力pc(MPa)使
用杨‑拉普拉斯方程求解：
[0088] pcij＝2γcosθ/rij
[0089] 其中，γ是界面张力(N/m)，θ为润湿接触角。
[0090] 在模拟过程中，选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，为此引入最小时间步长和修正时间步长，最小时间步长Δti(i＝1,2,…)是指在所有
管道中的凹液面到达下一节点的时间中，选取最小的Δti作为这一步运算总体的时间步
长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δxij＝vij·Δtmin，由此可以
求得该时间点的水力传导系数gij和两相在孔隙网络中的分布。显然最小时间步长法中时间
步长Δti与压力降Δp和凹液面位置有关，因此每次迭代的步长取决于该次计算的具体情
况，而不是都相等。最小时间步长的引入使得此模型可以使用尽量少的迭代次数得出模拟
结果，保证精度的同时大大提高了计算效率。
[0091] 实际渗流过程中，气体为可压缩流体，岩石骨架具有微可压缩性。若考虑流体和岩石可压缩性，可得到非稳态渗流方程：
[0092]
[0093] 通过泰勒展开以及隐式有限差分法技术，可将上式转换为矩阵方程：
[0094]
[0095]
[0096]
[0097] 式中，Ct为综合压缩系数(以气体压缩系数为主)，Q为气体采出速度，Δt为时间步长。
[0098] 在两相流中，每个节点的总体积流量依然满足守恒规律，即Σjqij＝0，据此可以构建线性方程组：
[0099]
[0100] 上述方程可以简化为矩阵求解式Ap＝B，其中A为主对角占优的对称稀疏矩阵，B为一组向量，p即为需要求解得到的全局压力场向量。上述矩阵可采用梯度下降法进行求解，
梯度下降法如下：
[0101] f(p)＝Ap‑B
[0102] f′(p)＝A
[0103]
[0104] (2)水驱气模拟
[0105] 如图4和图7所示，采用水驱气的方式进行模拟实验时，毛管力pcij成为动力，两节点间的流体总体积流量qij变为：
[0106]
[0107] 其后续处理方式与气驱水模拟一样。
[0108] 在天然气关键物性参数计算中包括气藏温度、压力与相对密度，偏差因子，体积系数，等温压缩系数的计算。
[0109] 1、气藏温度、压力与相对密度
[0110] 地下的天然气是多种气体组分的混合产物，其温度与压力通常采用拟临界参数来处理：
[0111] ppc＝∑yipci，Tpc＝∑yiTci
[0112] 式中，ppc,Tpc为天然气拟临界压力和拟临界温度；pci,Tci为气体组分i的临界压力和临界温度；yi为组分i的摩尔分数。
[0113] 天然气的相对密度γh表示天然气密度ρg；与空气密度ρair之比。
[0114]
[0115] 因此，其拟临界压力与拟临界温度可以通过其相对密度求得：
[0116]
[0117]
[0118] 通过当前天然气的压力p和温度T可以求得其视对比压力ppr与视对比温度Tpr：
[0119]
[0120]
[0121] 2、偏差因子
[0122] 天然气偏差因子z是定量描述真实气体(天然气)与理想气体偏差程度大小的系数,是天然气其他物性计算、天然气藏地质储量计算以及管道天然气产量设计的一项重要
的参数。天然气偏差因子计算的方法很多，本发明采用Dranchuk和Abou－Kassem‑11参数法
计算偏差因子。
[0123] z＝0.27ppr/(ρprTpr)
[0124] 且
[0125]
[0126] 其中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为中间参数，可用牛顿迭代法求取：
[0127] 令原函数为：
[0128]
[0129] 其一阶导数为：
[0130]
[0131] 其K阶导数与K+1阶导数关系为：
[0132]
[0133] 设置迭代精度(本次为误差小于0.05％)满足需求即可求得偏差因子z。
[0134] 3、体积系数
[0135] 天然气的体积系数是天然气地下体积量V与天然气地面标准条件下体积量Vsc的比值，其公式为：
[0136] Bg＝V/Vsc
[0137] 在油气藏条件下，压力为p，温度为T，将天然气状态方程、地面条件带入上式，可得到天然气体积系数计算公式：
[0138] Bg＝3.458×10‑4zT/p
[0139] 式中p、T为地层压力与温度，z为偏差因子。
[0140] 4、等温压缩系数
[0141] 天然气等温压缩系数Cg是指在等温条件下体积随压力变化的变化率，其数学表达式为：
[0142]
[0143] 考虑视对比压力后：
[0144]
[0145] 式中ppc为拟对比压力，ppr为视对比压力，z为偏差因子。
[0146] 以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本
文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进
行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围
内。
