CN107680131B - 一种快速确定多孔介质表征单元体积尺寸的方法 - Google Patents

Abstract

本发明涉及一种快速确定多孔介质三维重建模型的表征单元体积的尺寸方法，属于多孔介质参数研究领域；根据这种多孔介质的生成过程，利用离散单元法对符合特定粒径分布的球形颗粒进行随机堆积，得到多孔介质重建模型。选取足够大的体积单元，选取不同的分辨率对重建模型进行采样并计算渗透率，最终确定合适的分辨率，获得基于重建模型的三维数字图像；计算渗透率、比表面积、渗透率以及各参数变异系数随体积单元尺寸的变化，获取变化关系曲线，确定重建模型的表征单元体积，并利用最小二乘法确定表征单元体积尺寸与平均粒径的关系；能够被用于快速确定这种类型多孔介质重建模型的表征体积单元的大小。

Description

一种快速确定多孔介质表征单元体积尺寸的方法

技术领域

本发明涉及一种快速确定多孔介质三维重建模型的表征单元体积的尺寸方法，属于多孔介质参数研究领域。

背景技术

基于界面双电层效应的液环式角加速度计是一种新型角加速度计，在精度、稳定性、带宽等方面都有良好的性能。固相转换器是角加速度计的核心部件，完成压力信号向电信号的转换。固相转换器是由符合特定粒径分布的玻璃微珠堆积、烧结而成，利用离散单元法可以完成对固相转换器的三维重建，得到其重建模型。对重建模型进行分析，可以指导固相转换器的生产过程。

固相转换器是由固体相与孔隙相组成的多孔介质。它属于“不连续介质”，即其孔隙率、比表面积、渗透率等特性参数随其研究尺寸的不同而变化。因此在对固相转换器重建模型进行研究时，需要研究其特性参数随尺寸的变化规律，找到能够表征重建模型特性的最小尺寸。表征单元体积(Representative Element Volume,REV)是人们在研究岩层等多孔介质时，当其特性参数趋于基本稳定时选取的最小体积。当选取的体积小于REV时，多孔介质的特性参数会随着体积的变化而变化；当选取的体积大于REV时，多孔介质的特性参数基本稳定，可以将多孔介质视为等效连续介质。

离散单元法最初由Cundall提出，它是将系统分成单独的块体单元，实时分析每个块体单元的受力情况并更新其加速度、速度与位移等状态，直到系统达到稳定。被广泛应用于颗粒堆积系统的仿真。

目前确定多孔介质的表征单元体积的常用方法有理论与实验方法。采用理论方法需要的假设较多，同时需要针对不同的多孔介质做出不同的假设并建模，较为复杂；很多学者均采用实验的方法确定多孔介质的表征单元体积，但是实验次数较少，随机性较大。

很多学者对于特定岩体的表征单元体积进行了计算与分析。Clausnitzer.V等人通过实验，得出表征单元体积的边长为球粒径的5.15倍。但是仅考虑了球形颗粒粒径相同的情况。Razavi M R对实际多孔介质进行分析，设置REV的体积为球区域，对于由细长颗粒组成的硅砂，REV的半径为粒径中位数的5～11倍；对于由此球形颗粒组成的渥太华砂，REV的半径为粒径中位数的9～16倍。但是中位数难以完全表征粒径的分布特征，得到的无量纲化表征单元体积的尺寸差异较大。

基于在对符合特定粒径分布的球堆积模型表征单元体积及其无量纲化方面研究较少的现状，有必要研究这种多孔介质重建模型的REV及其无量纲化，便于后续固相转换器这类多孔介质以及数字岩心的研究。

发明内容

针对现有技术的不足，本发明提供了一种多孔介质重建模型的表征单元体积无量纲化的计算方法。重建模型是由符合特定粒径分布的球形颗粒随机堆积而成的系统。利用该方法可以快速确定这种多孔介质重建模型的表征体积单元的尺寸。

一种确定多孔介质表征单元体积尺寸的方法，包括如下步骤：

第一步，选取至少三组对数正态分布的多孔介质粒径分布参数，所述粒径分布参数为对数均值μ与对数标准差σ；

第二步，针对每一组粒径分布参数，建立多孔介质的三维模型；

第三步，确定三维数字图像分辨率；

第四步，基于第三步获得的分辨率，得到每一组粒径分布参数下三维模型的三维数字图像；利用三维数字图像计算得到对应的表征单元体积REV尺寸R_L；则每一组粒径分布参数对应一个表征单元体积REV尺寸R_L；计算每一组粒径分布参数下的表面积平均径D_p：

其中，x表示对数正态分布的随机变量，f(x；μ,σ)表示对数正态分布的随机变量x的概率密度函数；

第五步，针对得到的所有组粒径分布参数下的R_L与D_p，以D_p为自变量，以R_L为因变量，利用最小二乘法做线性拟合，得到R_L与D_p的线性关系式R_L＝kD_p+b；其中，k和b表示系数，通过线性拟合获得；

第六步，针对待确定表征单元体积尺寸的多孔介质，获得该多孔介质的表面积平均径D_p，再根据第五步确定的所述线性关系式得到该多孔介质的表征单元体积尺寸R_L。

较佳的，第三步中确定三维数字图像分辨率的方法为：

在第二步获得的三维模型中选取一个区域，在不同的分辨率下对该区域进行采样，制作不同尺寸大小的三维数字图像并计算沿流动方向的渗透率，得到该区域的“分辨率-渗透率”曲线，选取渗透率值稳定的部分对应的分辨率值作为分辨率。

较佳的，所述第四步中，计算表征单元体积REV尺寸R_L的方法为：

利用第三步确定的分辨率，在第二步建立的三维模型中随机选取至少10组三维图像；计算每一组三维图像的三个特性参数：孔隙率、比表面积以及渗透率；将所有三维图像对应的三个特性参数组成参数序列seq，并计算参数序列的变异系数Cov：

其中，Var(·)表示求方差；

当Cov＜0.05时，此时对应的尺寸即为表征单元体积尺寸。

较佳的，所述第二步中，采用离散单元法建立所述三维模型。

本发明具有如下有益效果：

本发明将基于离散单元法构建的多孔介质重建模型作为研究对象。通过实验选取合适的分辨率对其进行采样，构建三维数字模型；通过计算孔隙率、比表面积与渗透率随三维数字图像尺寸的变化，确定不同粒径分布参数的球形颗粒堆积模型表征单元体积的尺寸；最终获得了重建模型表征单元体积尺寸的取值方法；可以快速确定符合这种粒径分布的多孔介质重建模型表征单元体积的尺寸，而无需重新通过实验获得。

附图说明

图1为不同球形颗粒的粒径分布；

图2为基于离散单元法的球形颗粒堆积模型的堆积过程；其中，左图的iter＝0，中图的iter＝100000，右图的iter＝1000000；

图3为在一定分辨率下采样得到的多孔介质数字图像；

图4为type1堆积模型的渗透率随分辨率的变化关系；

图5为type1堆积模型孔隙率随选取的图像尺寸的变化关系；

图6为type1堆积模型孔隙率变异系数随选取的图像尺寸的变化关系；

图7为type1堆积模型比表面积随选取的图像尺寸的变化关系；

图8为type1堆积模型比表面积变异系数随选取的图像尺寸的变化关系；

图9为type1堆积模型渗透率随选取的图像尺寸的变化关系；

图10为type1堆积模型渗透率变异系数随选取的图像尺寸的变化关系；

图11为本发明的方法流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述，如图11所示：

第一步，基于离散单元法对多孔介质的三维模型进行重建。

1.1确定球形颗粒的粒径分布特征与堆积容器

多孔介质由玻璃微珠堆积、烧结而成,对实际玻璃微珠的粒径分布特征进行拟合,其符合对数正态分布的特征,即确定粒径分布的对数均值μ与对数标准差σ。

对于符合对数正态分布的随机变量x，其概率密度函数为

偏度为期望值和标准差分别为

堆积容器的形状包括但不限于圆柱形、长方体等。

1.2在容器中随机生成球形颗粒并进行受力分析

利用离散单元法对小球进行随机堆积，每个小球被视为单独的块体单元。对于系统中的任意小球i，其受力分析如下

其中m_i是小球i的质量，I_i是其惯量，

分别是球与球之间以及球与容器边界之间的接触力，l_i、l_j分别是

与小球的接触点到球心的距离。设玻璃微珠之间的弹性模量为E，两个小球之间没有发生碰撞时，它们之间没有相互作用力，即

当小球之间发生碰撞时，计算方法为

r_i、r_j分别是小球i、j的半径，d_ij为两个球心的距离，即球与球或者球与容器之间的作用力与它们重叠的距离大小成正比，比例系数为玻璃微珠的弹性模量。

同理。

小球的状态包含其位置、速度与加速度信息，在每一个迭代步长中采用中心差分法更新小球的状态，首先计算小球的受力，再更新其加速度、速度与位移信息。迭代时间步长设置为Δt。

1.3停止迭代

当系统迭代次数达到设置的最大次数Iter_max次之后，停止迭代，记录所有小球的半径及其位置、速度、加速度等状态信息。

第二步，确定重建模型的三维数字图像分辨率：

选取特定的分辨率对重建模型进行采样，生成三维数字图像。需要选取合适的分辨率，避免分辨率太低导致细节信息丢失或者分辨率太高导致计算特性参数耗时的问题。渗透率是表征多孔介质渗透特性最重要的参数。计算不同分辨率下生成的三维图像的渗透率，确定合适的分辨率为后续生成基于重建模型的三维数字图像做准备。

这种多孔介质中流体流速很小，利用达西定律可以计算其渗透率值。设通过多孔介质某一横截面的流量为Q，渗透率为

其中A为横截面的面积，ΔP为流体沿流动方向的压降，μ为流体粘度，L为流体流动方向的距离。对于数字图像，采用格子玻尔兹曼方法计算多孔介质的渗透率。

根据离散单元法得到的球形颗粒堆积结果，选取相同大小的区域，在不同的分辨率下对区域进行采样，制作不同尺寸大小的三维数字图像并计算沿流动方向的渗透率。对选定的区域，选取不同的分辨率进行采样，制作三维数字图像，并计算在不同分辨率下的渗透率，得到该区域的“分辨率-渗透率”曲线，选取渗透率值稳定的部分对应的分辨率值，设为[R_min,R_max]，作为第三步对该类型多孔介质的采样分辨率。

分辨率太大会丢失重建模型的细节信息，可以通过直径最小的球形颗粒来表征。上述实验中重建模型的最小粒径为d_min。对于任意球形颗粒堆积得到的多孔介质重建模型，最小粒径为d′_min，建议分辨率值的范围是

可以快速确定这种多孔介质重建模型的分辨率。

第三步，不同粒径分布参数下的特性参数计算以及表征单元体积尺寸选取。

孔隙率是孔隙部分与多孔介质总体积的比值；比表面积是指单位质量或体积所具有的表面积；渗透率是指在一定压力差下，介质允许流体流过的能力。上述三个指标能够反映多孔介质的流体传输特性。选取这三个参数作为研究指标，当三者均收敛于稳定值时，记录对应的体积单元尺寸作为REV的尺寸。

对于三维数字图像img，固体相为1，孔隙相为0。设图像中任意一点r，img(r)为图像在点r处的值，P{A}为事件A的概率。则孔隙率

p＝P{img(r)＝0}

孔隙相两点相关函数是基于图像中两点距离的函数，其定义为

S(d)＝P{img(r)＝0and img(r+d)＝0}

其中d是图像中两点之间距离。设图像维度为β(2或3)，图像比表面积为

s＝2β[S(0)-S(1)]

对于每个特定的尺寸，随机选取N＝10组三维图像，计算特性参数得到序列seq及其变异系数Cov，计算方法为：

Cov＜0.05时，认为参数逐渐稳定，此时对应的尺寸即为该特性参数的表征单元体积尺寸。

第四步，表征单元体积的计算与拟合。

平均粒径有多种表示方式，包括算术平均径、体积平均径、表面积平均径等。表面积平均径指与给定颗粒群总体积相同、总表面积相同、且粒度均匀的一个假想颗粒群的粒径，计算公式为：

表面积平均径考虑到多孔介质的表面积特征，能够更好地反映多孔介质的渗流特性，选取表面积平均径作为球形颗粒的平均粒径。选取多组对数正态分布的粒径分布参数，并利用这些参数进行三维重建并计算得到对应的REV尺寸R_L，对R_L与D_p利用最小二乘法做线性拟合，得到其线性关系式R_L＝kD_p+b。

实施例：

第一步：基于离散单元法对多孔介质进行三维重建。

角加速度计中固相转换器是由玻璃微珠堆积、烧结而成，对玻璃微珠的粒径拟合得到其粒径分布。它符合对数正态分布，不同类型固相转换器的概率密度分布如图1，具体的分布参数(变量对数的均值与标准差、表面积平均径)如表1。

表1不同粒径类型玻璃微珠的参数

类型	对数均值μ	对数标准差σ	表面积平均径D<sub>p</sub>(um)
				type1	3.2941	0.2768	32.3929
type2	3.7275	0.2732	49.7318
				type3	4.1702	0.2568	75.8470

对于任意类型的固相转换器，设其平均粒径，即表面积平均径为D_p。为避免壁面效应对三维数字图像的影响，设置圆柱体的容器的底面直径为30D_p，高度为60D_p，保证堆积完成后的球形颗粒的最高高度H满足H＞20D_p。根据玻璃微珠的实际参数值，设置仿真系统中球形颗粒的弹性模量E＝46GPa，迭代时间步长设置为Δt＝0.00001s，最大迭代次数为Iter_max＝1000000。在容器内随机生成小球并开始迭代过程。迭代初始状态、迭代过程中的状态以及迭代完成后的状态如图2。

迭代完成后，记录所有球形颗粒的半径与位置信息，方便后续生成基于多孔介质重建模型的三维数字图像。

第二步：重建模型生成三维数字图像分辨率的确定。

选取type1进行分析，选取研究的立方体实际尺寸为500um×500um×500um。在1～3μm/pixel的范围内间隔0.1um/pixel选取分辨率，生成重建模型基于三维多孔介质的三维数字图像如图3。利用格子玻尔兹曼方法计算重建模型沿堆积方向的渗透率，得到渗透率值随分辨率的变化曲线如图4。对于type1，分辨率在1～2μm/pixel时，渗透率计算的结果较为稳定。而分辨率太大会丢失重建模型的细节信息，可以通过直径最小的球形颗粒来表征。type1中最小的粒径为

d_min＝e^μ-3σ＝11.747μm

d_min在分辨率值为1～2μm/pixel时能够较好地保留细节特征。对于任意球形颗粒堆积得到的多孔介质重建模型，建议分辨率值的范围是

制作type1在分辨率1μm/pixel的三维切片图像。type2～type3同理。

第三步：表征单元体积的尺寸确定。

选取三维图像的孔隙率、比表面积以及渗透率作为确定表征单元体积尺寸的特性参数。以type1的重建模型为例。从小到大依次选取不同尺寸的体积单元，每种尺寸随机选取10个体积单元，计算其孔隙率、比表面积、渗透率以及各参数的变异系数。得到图像如图5～10。

设变异系数的阈值为Cov_th＝0.05，当参数的Cov≤Cov_th时，认为参数收敛，此时的体积单元尺寸即为REV的尺寸。计算所有特性参数的REV尺寸。选取最大的REV尺寸作为重建模型的尺寸，得到最终的REV尺寸为466um。type1～type3的REV尺寸见表2。

表2不同类型重建模型特性参数的表征单元体积尺寸

通过最小二乘法拟合得到REV尺寸与平均粒径的关系为R_L＝11.9128D_p+9.7044。

综上所述，以上仅为本发明的几个实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种确定多孔介质表征单元体积尺寸的方法，其特征在于，包括如下步骤：

第一步，选取至少三组对数正态分布的多孔介质粒径分布参数，所述粒径分布参数为对数均值μ与对数标准差σ；

第二步，针对每一组粒径分布参数，建立多孔介质的三维模型；

第三步，确定三维数字图像分辨率；

第四步，基于第三步获得的分辨率，得到每一组粒径分布参数下三维模型的三维数字图像；利用三维数字图像计算得到对应的表征单元体积REV尺寸R_L；则每一组粒径分布参数对应一个表征单元体积REV尺寸R_L，具体方法为：

利用第三步确定的分辨率，在第二步建立的三维模型中随机选取至少10组三维图像；计算每一组三维图像的三个特性参数：孔隙率、比表面积以及渗透率；将所有三维图像对应的三个特性参数组成参数序列seq，并计算参数序列的变异系数Cov：

其中，Var(·)表示求方差；

当Cov<0.05时，此时对应的尺寸即为表征单元体积尺寸；

计算每一组粒径分布参数下的表面积平均径D_p：

其中，x表示对数正态分布的随机变量，f(x；μ,σ)表示对数正态分布的随机变量x的概率密度函数；

第五步，针对得到的所有组粒径分布参数下的R_L与D_p，以D_p为自变量，以R_L为因变量，利用最小二乘法做线性拟合，得到R_L与D_p的线性关系式R_L＝kD_p+b；其中，k和b表示系数，通过线性拟合获得；

第六步，针对待确定表征单元体积尺寸的多孔介质，获得其表面积平均径D_p，再根据第五步确定的所述线性关系式得到表征单元体积尺寸R_L。

2.如权利要求1所述的一种确定多孔介质表征单元体积尺寸的方法，其特征在于，第三步中确定三维数字图像分辨率的方法为：

在第二步获得的三维模型中选取一个区域，在不同的分辨率下对该区域进行采样，制作不同尺寸大小的三维数字图像并计算沿流动方向的渗透率，得到该区域的“分辨率-渗透率”曲线，选取渗透率值稳定的部分对应的分辨率值作为分辨率。

3.如权利要求1所述的一种确定多孔介质表征单元体积尺寸的方法，其特征在于，所述第二步中，采用离散单元法建立所述三维模型。

Images