基于空间相关性不同组构干热岩裂隙发育

2021-12-02杜广盛陈世江马万里

科学技术与工程 2021年33期

杜广盛，陈世江*，马万里

(1.内蒙古科技大学矿业与煤炭学院，包头 014010；2.内蒙古自治区矿业工程重点实验室，包头 014010)

随着能源需求的清洁性和可再生性不断攀升，地热能、天然气等能源被进一步开发出来。早在1974年，美、英、法、德、日等国家就联合建立了研究中心。干热岩(hot dry rock)作为地热能开采最具价值的部分，它的开采方式是通过地下注入高压水使其在温度和地应力作用下致裂其上部低应力岩石，并通过该裂隙通道从而与周围岩石交换能量最终回收[1-3]；以往的研究表明干热岩普遍存在于地下较深部的以花岗岩为主的高温岩石，因此对花岗岩的裂隙发育规律的研究有着十分重要的意义。

在裂隙研究的过程中，通常认为岩石材料的破坏是由岩石内部微裂纹逐渐展开贯通，最终连接形成一条或多条较粗的宏观裂纹[4-5]。众多研究表明岩石的细观结构特征决定其宏观力学响应同时也影响其裂隙发育过程，岩石不同的力学特性就是其组构非均质性的宏观体现[6]。部分学者对于矿物组分、颗粒尺度分布特征与岩石力学性能的关系进行了研究[7-9]，结果表明矿物组分、颗粒尺度分布对岩石强度有明显影响。对于矿物质组分含量及分布的研究，大量学者分别从均质、随机非均质、原岩非均质等不同角度进行分析研究[10-14]，其中采用原岩非均质性研究能够更好地揭示岩石破坏机理和裂隙发展规律。但该方法采用的图像处理过程中存在一定的误差，同时对图像中岩石组构分布状态无法进行准确的定量。唐欣薇等[15]、Tang等[16]、黄文敏等[17]引入了空间相关性函数用来对岩石组分分布和各向异性进行描述，实现了通过改变参数表示不同矿物聚合分布的岩石结构状态。

随着岩石力学研究的不断深入，人们在岩石力学理论分析、试验、数值模拟等方面有了很大的进展，针对岩石结构的研究，其中基于离散元理论的数值模型能够较好地表征岩石间结构破坏时的力学关系[18-21]。现通过引入空间相关性函数，生成可控制空间相关长度进而控制组构分布状态的花岗岩数字图像；然后应用数字图像技术对生成的图像进行区分，一方面消除岩石图片采集的误差，另一方面使分析过程中使用的图像组分分布可以进行有明确的参数标定；最后采用细观颗粒流数值模拟方法对花岗岩进行抗拉强度试验模拟，观察其破坏过程中裂隙发展、结构状态对岩石强度的影响，进而从破坏过程中的细观角度分析不同组分状态与岩石试件强度的关系。为研究不同组分花岗岩破裂行为提供一种新的分析思路。

1 含组构信息花岗岩数值模型的构建

基于岩石的细观层次对岩石内部结构破坏规律与强度的关系进行分析，需要获得数值分析中最重要的岩石力学细观参数，首先采用相应室内试验的方式获得部分岩石物理参数，然后对原岩结构获取并进行模型构建，进而进行力学试验与数值试验匹配细观参数。最终获得该岩石较为可靠的细观参数。

1.1 试件细观物理参数获取

试验采用苏州纽迈电子科技有限公司生产的miniMR-60核磁共振分析测试仪(图1)对试件进行孔隙度检测[22]。5个试样测得孔隙度分别为1.01%、0.974%、0.921%、0.982%、0.881%，平均孔隙度为0.953 6%。

图1 miniMR-60核磁共振分析测试仪Fig.1 Nuclear magnetic resonance analysis tester miniMR-60

为获取花岗岩颗粒尺寸及分布情况，采用电子显微镜对花岗岩表面进行观测，在每个试件上选取八九个观测点进行图像获取，对其中较为明显的颗粒结构进行测量(图2)，观测结果表明石英、长石、云母颗粒直径分别为0.14～0.25、0.20～0.32、0.07～0.13 mm。

图2 岩石颗粒显微成像Fig.2 Microscopic imaging of rock particles

1.2 不同组构分布矢量化模型构建

采用数码相机设置相同参数获取花岗岩表面图像并进行调整。

将拍摄的图像在MATLAB软件中进行数值图像处理，首先进行灰度处理，使得图像变成以0～255表示不同灰度的图像；然后进行以灰度值0～96为云母、96～195为石英、196～255为长石不同阈值组分划分；最后采用中值去噪获得含三种组构的图像。

划分后采用canny边缘检测算法对轮廓进行确定，将获取到的边缘图像导入CAD进行矢量化形成可被识别的矢量边界模型。

数值图像处理过程及组分区分如图3所示。

图3 试件组分位置结构获取过程Fig.3 Acquisition process of component position structure of test piece

1.3 花岗岩抗拉强度的测定

在花岗岩薄板上钻取Φ50 mm×25 mm的圆形试样。选取有明显组构特征的试件，对其表面进行打磨，使平整度控制在±0.05 mm范围内，试验共加工了5个试件。在钻取过程中因钻进速度及薄板平整度等原因，使得加工试件尺寸略有不同。

岩石巴西劈裂试验采用Wance106微机控制电液伺服万能测试机如图4所示，其最大垂直加载力为10 T，试验采用压力控制方式，加载速度设置为200 N/s，在达到极限强度时自动卸载5 s，自动记录应力及加载时间。将待测件圆盘直接放置在巴西劈裂载具之间，通过两端弹簧升起对试件抬升防止试件偏移，保证加载方向。

图4 花岗岩巴西劈裂试验Fig.4 Brazilian splitting test of granite

按照上述试验方案，进行了5个试件的巴西劈裂试验，试验结果见表1，花岗岩试件的抗拉强度采用国际岩石力学协会推荐计算公式[式(1)]进行计算。

表1 花岗岩试件尺寸及实验结果Table 1 Size and experimental results of granite specimen

(1)

式(1)中：σt为岩石的抗拉强度；P为岩石的峰值载荷；D为试样的直径；t为试样的厚度。

1.4 数值模型的建立及参数匹配

颗粒流分析程序(particle flow code，PFC)采用键粒子模型(bond particle model，BPM)平行黏结模型。该模型可以在颗粒与颗粒之间设置黏结组件，当黏结组件生效时可以传递力的作用关系，当到达设置黏结组件的极限强度时该组件失效，模型转变为linear Model线性接触模型。将上述步骤处理后形成花岗岩试件图像矢量边界导入PFC模拟软件中，生成含有岩石特征的空间相似模型。参照以往花岗岩离散元数值模拟的细观参数标定[23-26]，结合试验测得部分物理力学参数，对模型进行参数匹配，经过反复校准，获得与室内试验应力曲线相吻合的数值试验如图5所示。细观参数标定结果如表2、表3所示。

图5 数值模型与室内试验参数匹配曲线Fig.5 Matching curve of numerical model and indoor test parameters

表2 颗粒单元的细观参数Table 2 Mesoscopic parameters of the particle unit

表3 接触黏结模型的细观参数Table 3 Mesoscopic parameters of the contact bond model

2 不同组构花岗岩劈裂数值试验

2.1 不同组构花岗岩图像生成过程

干热岩开采过程中需要对地下深部的花岗岩进行一定程度的致裂，而在成岩过程中，各个矿物组分不是随机的具有一定的相关性特征，即一种矿物材料周围一定区域内的矿物分布状态有一定的聚集性，而不同矿物则保持一定的区分度。由于不同组分的力学性质有所区别，这对岩石破裂时的裂隙发育有着重要影响。

引入空间相关性函数，通过改变相关参数，生成不同组构特征的岩石图像，具体生成步骤如下。

步骤1在MATLAB软件中建立一个二维随机数组A，该数组中各元素ai,j的值均匀分布在(0,1)内。若岩石仅由两种含量不同的物质组成，则可将数组A中元素进行二值化；当0

步骤2为了生成矿物组构聚集性特征明显的岩石图像，须对图6(a)岩石图像进行处理，处理方法就是引入空间相关函数对数组A进行变换，此空间相关函数式见式(2)，叠加方法见式(3)。变换后得数组B，具体计算方法参见文献[16]。由数组B可生成矿物组构聚集特征明显的岩石图像；式(2)中参数L控制岩石矿物组构聚集程度。图6(b)是L=6 时的岩石图像。

(2)

式(2)中:L为空间相关长度；d为两个元素点之间的有效距离。

(3)

式(3)中：bi,j为数组B在(i,j)坐标下的元素值；ai+p,j+q为与元素bi,j点距离为(p,q)时数组A中的元素值；n为计算区域取值尺寸，选取为45。

步骤3将求解出的数组B乘对应系数转换为数值范围在0～255的矩阵，此时该矩阵所含数值范围与图像所呈现的灰度值保持一致，借助灰度值组分划分的方法采用相同的参数进行划分，既0～96为云母；96～195为石英；196～255为长石。图6(c)是经灰度划分后的3种不同组分的分布图。

图6 岩石组构图像数值实现过程Fig.6 Numerical realization of rock fabric image

2.2 不同组构花岗岩试件模型表征

当只考虑岩石结晶形成时细观组分的空间相关性时，在计算过程中选取空间相关长度分别为L=2、4、6、8，介质含量x=0.45，此时不同的空间相关长度计算结果如图7(a)所示。针对两相介质的矿物成分含量比例的情况，当成分含量不同时，考虑其含空间相关性的细观结构表征。以空间相关长度L=6为例，取其中某一介质使其为低强度组分，分别赋予含量x=0.25、0.35、0.45，得到不同矿物成分含量计算结果如图7(b)所示，很明显在选取不同矿物成分含量时，计算结果不同介质比例也有明显的不同。

2.3 试验结果及分析

将表中细观参数赋予模型试件，对生成的含不同空间参数的结构进行巴西劈裂试验，每种空间参数生成了3组随机图像，对同一图像结构在保证颗粒尺寸在一定范围内取值随机变量，得到不同的5个试验结果，如图8所示。

从图8数据和图7(a)中可以看出，空间相关性系数L能够影响矿物组分聚合的排布，L值升高矿物组分更趋于大块聚集。通过对不同L值的巴西劈裂模拟试验，L=4、6、8试件极限载荷平均值分别为13.577、11.482 77、11.200 55 kN。抗拉强度随L的增大而减小。表4中数据表明当L=6，x=0.25、0.35、0.45时矿物组分含量有明显的变化同时其所对应的试件极限载荷平均值为15.208、13.881、11.482 kN。当x所表示的低强度组分含量增大，试件强度明显降低。

图7 不同组构花岗岩试件图像Fig.7 Image of different fabric granite specimens

3 不同组构试件裂隙发育的分形分析

大量试验表明流体通过岩石间裂隙时，裂隙的复杂程度和裂隙总长度分布情况，会影响流体通过裂隙的速度和热交换的效率。因此对受载状态下花岗岩裂隙的发育状态研究有着很大意义。

事实上花岗岩内部不同组构的物理力学性质有所差异，其微裂隙发育状态也有着明显的不同，因此微裂隙发育状态的分析是连接花岗岩组构与强度之间关系的桥梁。

3.1 分形维数计算过程及原理

分形维数是一种能够对微裂隙分布状态描述的手段[27-28]。选取其中组构关系较为明显的图片，从加载过程中破裂点产生的附近，每间隔500步获取试件数值模拟加载过程中内部破裂图像并记录其相应的破裂点数目，利用MATLAB软件对产生的破裂形貌图像去噪和灰度处理并对其分形维数计算，图9为20 000步L=4的计算结果图。在计算分形维数时，采用计盒维数算法，即对图像进行2n(n=1，2，…)不同尺度的划分，并对其中所含裂隙的进行计数，计算公式见式(4)。本次计算图像大小设置为512×512像素，最大网格划分为256×256。各阶段破裂点数据及相应分维数据发展趋势如图10、图11所示。

图9 20 000步L=4试件分形维数计算结果Fig.9 Calculation results of fractal dimension of 20 000 step L =4 specimen

(4)

式(4)中：D为计盒维数；N(C)为在网格尺度为C时所覆盖分形体的格子数；C为网格尺度。

3.2 不同空间相关长度加载阶段破裂特征分形分析

由图10可知，破裂点以及相应分形维数随加载步长增加总体呈现增长趋势。加载破裂初始阶段，当L较小时，组分分布离散程度较高。在加载过程中空间相关长度L高的结构产生的破裂点较多，与其相对应的分维值也相对较高；相同加载步长下L=4的结构比L=8的结构强度高。

图10 不同空间相关系数破裂点及相应分形维数发展趋势Fig.10 Fracture points with different spatial correlation coefficients and development trend of corresponding fractal dimensions

从裂隙的发育状态可以看出，当L=4时，花岗岩中不同的成分混合在一起，不同成分间产生的微小结构面也较为均匀的散布在整个试件中，该状态下加载时产生的裂隙数也相对较小，与其相应的分形维数也较小，不利于应力集中现象的产生；当L=8时，花岗岩中不同成分区分较为明显，不同组分之间有较大、较明显的软弱结构，此时试件加载过程中产生的裂隙数相对较大，其分维数也较高，更加有利于应力集中现象的出现。当加载至17 500步时，L为4、8时，破裂点数为1 587、3 009，此时分维值为1.240 8、1.232 1，说明L=8时破坏点虽多，但其规律性较强，都是在加载轴线附近的较大软弱结构处产生的破坏，应力相对集中。而L=4破坏点虽少，但其分散性较高，在加载到试件即将破坏的后期，有部分裂隙趋向于分散发展，使其应力分散，从而破坏强度上升。

3.3 不同组分含量加载阶段破裂特征分形分析

如图11所示，不同组分含量试件，在加载过程中的破坏点数与破坏图像分维数随加载步长整体呈现增长趋势，当x=0.25时，云母含量较少，试件中大量成分被设置为长石和石英，因此形成的低强度结构面较少，破坏过程中裂隙产生数量降低，初始裂隙分维值也因此较低，强度较高；当x=0.45时，云母含量较多，试件中不同组分间的结构面形成较多，破坏过程中裂隙产生较多，裂隙初始分维值较高，强度较低。当加载至18 000步时，不同组分含量x为0.25、0.45时，破裂点数为1 193、2 222，此时分维值为1.224 9、1.200 5，这说明x=0.25试件中的破裂点更加均匀地分布在整个空间内，裂隙趋向于分散发育，而x=0.45的试件中破裂点由于集中在较弱组构附近使得裂隙发育趋于集中，应力集中现象明显，因此其强度更低。从上述试验数据可以看出花岗岩抗拉强度受到破坏点数和破坏点分维数的影响。