基于CT扫描的预测应力加载条件下煤体渗透率的方法与流程

本发明涉及煤层注水渗流技术领域，尤其是一种基于ct扫描的预测并计算应力加载条件下煤体渗透率的方法。

背景技术：

煤矿生产过程中，需要通过煤层注水实现对工作面粉尘污染、煤与瓦斯突出、冲击地压、自然发火、煤体软化等矿井实际生产问题进行综合治理，而注水后的渗流是自然界最常见的物理现象。在一定压差下，岩石允许流体通过的能力叫渗透率，渗透率是表征土或岩石本身传导液体能力的参数，其大小与孔隙度、液体渗透方向上空隙的几何形状、颗粒大小以及排列方向等因素有关，渗透率用来表示渗透性的大小。

传统的渗透率研究方法包括：(1)通过岩石试验系统和真三轴试验系统等对稳定渗流状态下的煤体应力-渗流耦合作用进行研究，在轴压和侧压稳定的情况下，通过测取流量、流动长度、压力梯度等参数，并结合达西渗流公式也可以实现对煤体渗透率的计算。但是由于煤体结构的非均质不透明等特性，所以该方法不能对煤体内部微细观孔裂隙结构的渗流过程进行解释和描述，而复杂的外界环境会影响流量等参数的确定，从而降低渗透率计算结果的准确性。(2)对于非稳态渗流过程中的渗透率变化规律，常采用脉冲衰减法进行研究，其基本工作原理为通过测试岩样一维非稳态渗流过程中孔隙压力随时间的衰减数据，在此基础上，结合数学模型对煤样脉冲渗透率进行测取。但是该方法在测试高渗透煤样时误差很大，这限制了此方法的测试范围。

随着ct扫描技术在岩石领域的大量应用，其具有的高精度、无损伤等特点，可以帮助人们捕获各种储层岩石的孔裂隙空间几何形状，实现对于孔隙率的测定以及孔裂隙结构特征、矿物形态的分析；同时，借助载荷加载装置还可以实现对于孔裂隙结构演化规律的研究。在气体渗流方面，ct扫描技术主要侧重于气体吸附和解吸作用研究；而在注水渗流方面，利用该技术则可以实现渗流通道的构建以及揭示注水渗流过程中微细观结构的损伤破坏机理。但是受到监测手段的制约，目前的ct扫描实验仍然无法实现对渗流通道中流速、压力梯度变化情况进行详细的测量统计，同时复杂的实验测试也增大了研究成本。因此，为进一步的提高煤体渗透率计算的准确性，实现对煤体内部渗流分布情况的描述，有效预测不同应力条件下渗透率的变化情况，需要对现有的技术有待于更进一步的改善和发展。

技术实现要素：

为解决实现直观描述变形煤样内部渗流分布情况，并有效预测不同应力条件下煤体渗透率变化的技术问题，本发明提供了一种基于ct扫描的预测应力加载条件下煤体渗透率的方法，具体技术方案如下：

基于ct扫描的预测应力加载条件下煤体渗透率的方法，其特征在于，包括以下步骤：

a.制作煤样试件，进行单轴压缩实验，同时对煤样试件进行ct扫描；

b.采用数字地形模型法进行阀值分割，并得到煤样试件的阀值，将ct扫描图像导入simpleware软件中建立三维数值模型；

c.将三维数值模型导入hypermesh软件中并设置渗流条件；

d.对三维数值模型施加多个压力梯度值▽p，设置初始流动速度v0，调整▽p和v0参数分别模拟，得到渗流模型；

e.将渗流模型导入ls-dyna软件中模拟计算，得到计算结果；

f.将计算结果导入hyperview软件中，沿渗流方向等距提取渗流速度，得到渗流速度和压力梯度的关系，计算渗透率。

优选的是，煤样试件为圆柱体，所述单轴压缩实验使用原位拉伸、压缩及温度控制实验装置；所述单轴压缩试验中分别进行加载速度不同的多组煤样试件的单轴压缩，同一煤样试件的加载速度恒定。

进一步优选的是，步骤a中，在ct扫描前首先根据x射线稳定性、煤样试件尺寸、煤样试件x射线衰减分数和曝光时间来确定扫描电压、扫描功率和视场大小；所述ct扫描是以恒定的速度旋转扫描，检测器捕捉由x射线源发出穿过煤样试件的x射线，以电信号的形式储存ct扫描图像。

还优选的是，步骤b中确定煤样试件的阀值包括：将ct扫描图像转化为数字地形模型，表示孔裂隙体积和煤体总体积，孔裂隙体积与煤体总体积的比值作为孔隙率；建立孔隙率和图像灰度值的函数关系，计算孔隙率函数所有极小值中的最大值，作为ct扫描图像的阀值。

还进一步优选的是，步骤b中在simpleware软件建立三维数值模型后进行网格划分。

还优选的是，步骤c中的渗流条件包括水场模拟、空气场模拟和煤体模拟；所述水场模拟选用mat-null本构模型，设置流体沿一个方向流动，在流动方向上的两个侧面施加垂直于流动方向的两个方向的约束，在另外四个侧面施加固定约束；通过布置线性多项式状态方程为水头处流体的运动提供压力。

还进一步优选的是，空气场模拟中选用mat-vacuum空隙材料模型，并选择ale算法作为运算方法。

还进一步优选的是，煤体模拟选用拉格朗日算法作为运算方法。

还优选的是，压力梯度值▽p的选择范围为0～300pa/mm；所述初始流动速度v0的选择范围为0.01～0.09mm/s。

还进一步优选的是，骤d中渗流模型施加的压力梯度值▽p分别为0pa/mm、10pa/mm、24pa/mm和213pa/mm；初始流动速度v0设置为0.02mm/s。

本发明的有益效果包括：

(1)本方法通过室内试验和ct扫描获得煤体的结构特征，并采用数字地形模型法进行阀值分割，得到煤样阀值并建立三维数值模型，从而降低了试验和测量成本，提高了计算效率，保证了数值模拟的准确性。

(2)通过设置渗流条件得到渗流模型，在煤体模拟中选用了拉格朗日算法能够精确地描述结构边界运动；空气场模拟中选用mat-vacuum空隙材料模型，并选择ale算法作为运算方法，能够克服因流体穿过煤体骨架结构而引起的自身网格的严重畸变。

(3)利用数值模拟的方法对渗流进行模拟，根据模拟结果能够实现对煤体内部渗流分布情况的直观描述，并根据变化过对对渗流变化进行有效的预测。

(4)本方法中利用达西定律拟合压力梯度和渗流速度的关系，得到形变煤体结构的渗透率，预测结果准确性高，测试方法的适用范围广。

附图说明

图1是煤样试件形状结构示意图；

图2是ct扫描结构示意图；

图3是原位拉伸、压缩及温度控制实验装置结构示意图；

图4是煤样试件微观结构示意图；

图5是渗流模型和渗流条件示意图；

图6是裂隙渗流速度曲线图；

图7是孔隙渗流速度曲线图；

图中：1-煤样试件；2-x射线源；3-载物台；4-ccd探测器；5-x射线；6-刚性压头；7-上夹持装置；8-上固定螺栓；9-下加持装置；10-下固定螺栓；11-煤样试件微观结构；12-水场；13-空气场。

具体实施方式

结合图1至图7所示，本发明提供的基于ct扫描的预测应力加载条件下煤体渗透率的方法具体实施方式如下：

实施例1

基于ct扫描的预测应力加载条件下煤体渗透率的方法，包括以下步骤：

a.制作煤样试件，进行单轴压缩实验，同时对煤样试件进行ct扫描；

其中，煤样试件为圆柱体，所述单轴压缩实验使用原位拉伸、压缩及温度控制实验装置；所述单轴压缩试验中分别进行加载速度不同的多组煤样试件的单轴压缩，同一煤样试件的加载速度恒定。

在本步骤中，ct扫描前首先根据x射线稳定性、煤样试件尺寸、煤样试件x射线衰减分数和曝光时间来确定扫描电压、扫描功率和视场大小；所述ct扫描是以恒定的速度旋转扫描，检测器捕捉由x射线源发出穿过煤样试件的x射线，以电信号的形式储存ct扫描图像。

b.采用dtm数字地形模型法进行阀值分割，并得到煤样试件的阀值，将ct扫描图像导入simpleware软件建立三维数值模型；在simpleware软件建立三维数值模型后进行网格划分。

确定煤样试件的阀值包括：将ct扫描图像转化为数字地形模型，表示孔裂隙体积和煤体总体积，孔裂隙体积与煤体总体积的比值作为孔隙率；建立孔隙率和图像灰度值的函数关系，计算孔隙率函数所有极小值中的最大值，作为ct扫描图像的阀值。

其中simpleware软件是基于核心图形处理平台scanip的图像处理模块，可选择的网格生成femodule有限元模块，以及cad集成+cad模块和nurbsmodule曲面建模模块，在处理与整合图像、cad与仿真技术领域有显著效果。

c.将三维数值模型导入hypermesh软件中并设置渗流条件；其中hypermesh软件具有强大的有限元网格划分前处理功能。

渗流条件包括水场模拟、空气场模拟和煤体模拟；所述水场模拟选用mat-null本构模型，设置流体沿一个方向流动，在流动方向上的两个侧面施加垂直于流动方向的两个方向的约束，在另外四个侧面施加固定约束；通过布置线性多项式状态方程为水头处流体的运动提供压力。空气场模拟中选用mat-vacuum空隙材料模型，并选择ale算法作为运算方法。煤体模拟选用拉格朗日算法作为运算方法。

d.对三维数值模型施加多个压力梯度值▽p，设置初始流动速度v0，调整▽p和v0参数分别模拟，得到渗流模型；

压力梯度值▽p的选择范围为0～300pa/mm；所述初始流动速度v0的选择范围为0.01～0.09mm/s。在试验中，优选的渗流模型施加的压力梯度值▽p分别为0pa/mm、10pa/mm、24pa/mm和213pa/mm；初始流动速度v0设置为0.02mm/s。

e.将渗流模型导入ls-dyna软件中模拟计算，得到计算结果；其中渗流模型为“k”文件格式，计算后得到“d3plot”文件。ls-dyna软件是功能齐全的几何非线性、材料非线性和接触非线性软件，它以lagrange算法为主，兼有ale和euler算法；以结构分析为主，兼有热分析、流体-结构耦合功能；以非线性动力分析为主，兼有静力分析功能的非线性有限元程序。

f.将计算结果导入hyperview软件中，沿渗流方向等距提取渗流速度，得到渗流速度和压力梯度的关系，计算渗透率。将“d3plot”文件导入hyperview软件中，进行后处理，沿模型等距提取渗流模型内部的多个稳定渗流速度。

实施例2

本发明提供的基于ct扫描的预测应力加载条件下煤体渗透率的方法实施步骤包括：使用nanovoxel-3502e型扫描仪和原位拉伸、压缩控制实验装置，对多组煤样试件进行不同应力加载条件下的ct扫描，获取煤样试件压缩过程中的结构ct扫描图像。将获取的ct扫描图像通过simpleware软件建立三维数值模型，从而将形变煤体结构的ct扫描图像转化为具有可视化特征的数值模型，以方便对煤体内部渗流分布情况的直观描述。通过hypermesh软件进行煤体结构的渗流条件设置，主要包括空气场模拟和水场模拟。通过ls-dyna软件进行模拟计算，将计算结果导入到hyperview软件中提取渗流速度，通过拟合压力梯度和渗流速度记录数据，得到形变煤体结构的渗透率。

具体实施步骤包括：

步骤a.首先将煤样试件制作成截面直径9mm的煤柱，通过原位拉伸、压缩及温度控制实验装置和ct扫描试验装置，对煤样试件在进行单轴压缩过程进行ct扫描。其中ct扫描试验装置包括x射线源2、载物台3、ccd探测器4和x射线5；原位拉伸、压缩及温度控制实验装置包括刚性压头6、上夹持装置7、上固定螺栓8、下加持装置9和下固定螺栓10。

具体操作过程为：将煤样试件安放在载物台3上，通过下固定螺栓10调节下加持装置9，完成对煤样试件的调平；通过上加持装置7和上固定螺栓8对煤样试件的顶端进行固定。根据x射线稳定性、样品尺寸、样品的x射线衰减分数和曝光时间等因素的影响，对x射线源2和ccd探测器4进行调节，选择适当的扫描电压、扫描功率和视场大小等实验条件。利用刚性压头6对煤样试件分别进行单轴压缩加载，加载速度分别选取0mm/s、0.001mm/s、0.002mm/s和0.003mm/s；同时以恒定的扫描速度对载物台3进行旋转，x射线5穿过煤柱1后由ccd探测器4进行捕捉，以电信号的形式进行储存并得到ct扫描图像。

步骤b.使用dtm法进行阈值分割，确定各煤样的阈值情况。其原理是通过将ct的灰度图像转化成数字地形模型，基于此对孔裂隙的体积ve、煤体的总体积vt进行表征，进而构建孔隙率φ(x)与灰度值x之间的变化函数，最终通过求取该函数的极小值中的最大值确定作为ct扫描图像的阈值。其中孔隙率φ(x)与灰度值x的函数关系式为：

式中：ri为各像素点的灰度值，ri∈[rmin，rmax]；h(ri)为[rmin，rmax]范围内的灰度直方图，其数值大小等于图像中灰度值为ri的像素的个数与图像像素个数总和的比值。

对不同加载速度下获得的4组ct扫描图像进行适当裁剪，并将其分别导入到simpleware软件建立三维数值模型，选取适当网格类型对该模型进行网格划分，构建出4组煤体微细观结构，如图4所示，煤样试件微观结构示意图，是原始煤体的微细观结构，此时单轴加载速度为0mm/s。

步骤c.将构建的三维数值模型煤体微观结构导入hypermesh软件中并设置渗流条件。模拟非稳态渗流过程，其中渗流条件模拟包括水场模拟、空气场模拟和煤体模拟。其中水场模拟设置流体并提供动能，空气场模拟和煤体结构耦合接触，空气场是流固耦合的主要场所。

在水场模拟中，规定流体只沿y方向进行流动，而其余侧面不允许流体流出。在设置时，对煤体流动方向上的两个侧面施加x、z方向上的约束，其它4个侧面施加x、y、z方向的约束。煤体模拟选用拉格朗日算法作为运算方法，在渗流模拟设置中lagrange算法能够精确地描述结构边界运动。由于流体的流动会引起水头和空气场网格的严重变形，所以选择能够克服因流体穿过煤体骨架结构而引起的自身网格严重畸变的ale算法作为水头和空气场的运算方法，其控制流体问题的方程是navier-stokes方程的ale描述：

式中：μ是速度矢量；σ＝(-pi+τ(μ))应力张量，p是压力，i是二阶恒等张量；对于牛顿流体，粘性应力张量由τ(u)＝μ(▽u+(▽u)^t)给出，其中μ是动态粘度；f是质量体力；ρ是密度。

在材料属性设置方面，选择能够使自身结构不发生变形的刚体属性对煤骨架结构进行定义。水场模拟采用适合于流体属性的“mat-null”本构模型；空气场模拟作为流固耦合发生的主要场所，在属性定义时需要同时兼顾空气物理性质以及煤的多孔结构特点，而“mat-vacuum”是一种结合ale算法的孔隙材料模型，它能够满足模拟的设置要求，因此选其作为空气场的本构模型；通过布置线性多项式状态方程为水头处流体的运动提供压力，其中线性多项式状态方程为：

p＝(c0+c1μ+c2μ²+c3μ³)+(c4+c5μ+c6μ²)ev

式中：p为流体压力；ev为内能与初始体积之比；μ为比体积(μ＝ρ/ρ0-1，其中ρ0为原始流体密度，ρ为当前流体密度)；c0～c6为系数；方程参数按照一般压力方程进行取值：c0＝p0＝1.01×10⁵pa，c2＝c3＝c4＝c5＝c6＝0，c1＝2.25×10⁹pa为水的体积模量。

步骤d.对4组煤体微细观结构分别进行渗流模拟设置，并施加0.02mm/s的初始流动速度v01，在此基础上设置▽p1、▽p2、▽p3的压力梯度，其中▽p1为10pa/mm，▽p2为24pa/mm，▽p3为213pa/mm。

步骤e.将建立好的渗流模拟以“k”文件的格式导入到ls-dyna软件中进行模拟计算，得到计算结果。

步骤f.将计算结果导入hyperview软件中，具体是将生成的“d3plot”文件导入到hyperview软件中进行后处理。沿模型等距提取内部流体的稳定渗流速度，通过拟合压力梯度与渗流速度的关系，获取形变煤体结构的渗透率。

实施例3

以新疆金塔大黄山煤矿2zw11工作面气煤煤样为例，对本发明的方法应用及原理做进一步的说明。

步骤a.首先将煤样试件制作成截面直径9mm的煤柱，通过原位拉伸、压缩及温度控制实验装置和ct扫描试验装置，对煤样试件在进行单轴压缩过程进行ct扫描。其中ct扫描试验装置包括x射线源2、载物台3、ccd探测器4和x射线5；原位拉伸、压缩及温度控制实验装置包括刚性压头6、上夹持装置7、上固定螺栓8、下加持装置9和下固定螺栓10。

具体操作过程为：将煤样试件安放在载物台3上，通过下固定螺栓10调节下加持装置9，完成对煤样试件的调平；通过上加持装置7和上固定螺栓8对煤样试件的顶端进行固定。根据x射线稳定性、样品尺寸、样品的x射线衰减分数和曝光时间等因素的影响，对x射线源2和ccd探测器4进行调节，选择扫描电压为60kv，功率为5w，视场大小为9.5×9.5mm²的实验条件。利用刚性压头6对煤样试件分别进行单轴压缩加载，加载速度分别选取0mm/s、0.001mm/s、0.002mm/s和0.003mm/s；同时以恒定的扫描速度对载物台3进行旋转，x射线5穿过煤柱1后由ccd探测器4进行捕捉，以电信号的形式进行储存并得到ct扫描图像，经过4.3h完成扫描。

步骤b.使用dtm法进行阈值分割，确定各煤样的阈值为162pixel。将获得的4组ct图片进行裁剪，得到直径为0.72mm的ct扫描图像，将裁剪后的ct扫描图像分别导入到simpleware软件建立三维数值模型，选取适当网格类型对该模型进行网格划分，构建出4组煤体微细观结构，如图4所示，煤样试件微观结构示意图，是原始煤体的微细观结构，此时单轴加载速度为0mm/s。

步骤c.将构建的三维数值模型煤体微观结构导入hypermesh软件中并设置渗流条件。模拟非稳态渗流过程，其中渗流条件模拟包括水场模拟、空气场模拟和煤体模拟。其中水场模拟设置流体并提供动能，空气场模拟和煤体结构耦合接触，空气场是流固耦合的主要场所。

在水场模拟中，规定流体只沿y方向进行流动，而其余侧面不允许流体流出。在设置时，对煤体流动方向上的两个侧面施加x、z方向上的约束，其它4个侧面施加x、y、z方向的约束。煤体模拟选用拉格朗日算法作为运算方法，在渗流模拟设置中lagrange算法能够精确地描述结构边界运动。由于流体的流动会引起水头和空气场网格的严重变形，所以选择能够克服因流体穿过煤体骨架结构而引起的自身网格严重畸变的ale算法作为水头和空气场的运算方法，其控制流体问题的方程是navier-stokes方程的ale描述。

煤体模拟选择能够使自身结构不发生变形的刚体属性对煤骨架结构进行定义。水场模拟采用适合于流体属性的“mat-null”本构模型。空气场模拟选择“mat-vacuum”模型。通过布置线性多项式状态方程为水头处流体的运动提供压力。

步骤d.对4组煤体微细观结构进行渗流模拟设置，并施加0.02mm/s的初始流动速度v01，在此基础上设置▽p1、▽p2、▽p3的压力梯度，其中▽p1为10pa/mm，▽p2为24pa/mm，▽p3为213pa/mm。

步骤e.将建立好的渗流模拟以“k”文件的格式导入到ls-dyna软件中进行模拟计算，得到计算结果。

步骤f.通过ls-dyna软件进行模拟得到计算结果，将计算结果导入hyperview软件中，也就是将生成的“d3plot”文件导入到hyperview软件中进行后处理。利用hyperview软件沿模型等距提取3个截面的内部孔隙、裂隙结构中流体速度随时间的变化情况，得到裂隙渗流速度曲线图和孔隙渗流速度曲线图，如图6和图7所示。并对4种渗流模型在压力梯度作用下的稳定渗流速度进行统计，结果如表1所示。

表1压力梯度作用下4组渗流模型的稳定渗流速度

将上述压力梯度、渗流速度数据导入到origin软件中，分别采用内置的“线性函数”、“非线性函数”进行拟合求解，只有“线性函数”的相关系数达到0.99。由此可以看出，在低压力梯度10～213pa/mm的条件下，随着压力梯度的增加，渗流速度呈线性增大，其关系符合达西定律的描述。通过达西定律对“稳定渗流值-压力梯度”关系进行拟合：

式中：v为模型的稳定速度值；▽p为压力梯度；k为渗透率；μ为动力粘滞系数。

通过计算，分别获得原始模型、0.001mm/s、0.002mm/s、0.003mm/s形变模型结构的渗透率为37.4md、36.8md、37.1md和37.2md。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。