3D Voronoi等效晶质模型在岩石破坏细观研究中的应用*

韩 森， 张钦礼

(1. 中南大学 资源与安全工程学院，湖南 长沙 410083；2. 湖南有色金属研究院，湖南 长沙 410100)

0 引言

采矿、隧道、水利水电和边坡等工程与岩石力学特性密切相关，岩石是由不同性质矿物组成的集合体，受到岩石内部矿物晶粒形状和组分、晶粒间的结构面及微裂隙等微观结构特性影响[1]导致岩石宏观力学性质千差万别，即使同样的岩石其力学特性也往往表现为各向异性[2]。因此，从岩石微观结构入手，建立精细化模型，开展岩石破坏特征的细观研究，对于岩体工程优化设计，具有重要的工程应用价值。

离散元法是一种非连续介质方法，不需要建立复杂的宏观本构方程，且其将岩石视为由颗粒和不连续面组成的离散体，这与岩石的微观结构较相似[3]。因此，离散元法在岩石细观力学特性等方面的研究取得了较快发展，Itasca公司推出的颗粒流程序(Particle Flow Code,PFC)商业软件，是由Cundall等[4]在离散元理论基础上发展而来的，能够实现具有不同力学参数多相介质的建模，并可以方便的展现细观力学机理。PFC通过建立颗粒之间接触的本构关系对岩石的力学特性进行数值模拟，首先Cundall和Potyondy提出了黏结颗粒接触模型和平行黏结接触模型,并对Lac du Bonnet花岗岩的破坏特征开展了细观研究[5-7]，但应用于单轴试验数值模拟时往往会高估岩石抗拉强度，偏离了岩石的实际情况[8]，在此基础上，簇黏结模型[9]、等效晶质模型[10-11]、平直节理模型[12]以及蒋明镜微观胶结模型[13]等被相继提出，有效解决了岩石单轴抗拉强度过高问题。本文重点建立反映岩石微观结构的精细化模型，开展岩石破坏特征的细观研究，根据郑达等[14]采用电镜对岩石断裂后的断口扫描研究发现微观断裂方式主要为穿晶断裂和沿晶断裂，等效晶质模型采用平行黏结和光滑节理接触模型表征岩石的晶粒和晶粒间结构面，能够展现穿(沿)晶断裂行为。周喻等[15]构建了二维颗粒流的等效晶质模型体，通过计算与试验结果对比验证了等效晶质模型的适宜性。然而，自然界中的岩石为三维形态，岩石中颗粒的形态、组分及颗粒间的结构面等均与二维形态存在较大差异，本文基于3D Voronoi图，采用PFC3D构建了能够反映岩石内部晶粒结构的等效晶质模型，分析了晶粒排列均匀化程度及大小等微观结构的变化对模型宏观力学性质的影响程度，模拟了板岩单轴压缩和巴西试验，并从细观角度揭示了岩石的破坏机理。

1 3D Voronoi晶质模型构建

1.1 生成3D Voronoi晶粒网

Voronoi晶胞的几何形状与多晶粒的微观结构相似，因此在计算机仿真中常采用Voronoi图构建多晶粒，本文采用Rhino6.0生成3D Voronoi晶粒网模拟岩石的多晶粒结构，具体步骤为：在PFC3D中随机生成“粗”颗粒球体，其球体半径与岩石晶粒尺寸相当，导出球体中心坐标点，形成三维空间点集合，在Rhino6.0中构建由点集产生的3D Voronoi晶粒网。

1.2 构建岩石晶粒结构

构建不同矿物组成成分的岩石晶粒结构，需要以下3个步骤：1)根据岩石的晶粒大小生成尺寸相当的“粗”颗粒体模型，依照岩石的矿物组分情况分组，并记录分组信息；2)按照上文表述的方法构建晶粒网，并将“粗”颗粒的分组信息逐一对映至每个晶粒结构单元中；3)建立“细”颗粒体模型，将晶粒网覆盖至该颗粒体模型中，利用向量数量积方法判断每个“细”颗粒所在的晶粒结构单元，并将分组信息映射至对应的“细”颗粒体中，最后“细”颗粒模型被赋予平行接触模型，晶粒网被赋予光滑节理接触模型，平行黏结模型表征岩石晶粒体、光滑节理接触模型表征岩石晶粒间的结构面，形成了与岩石内部结构相似的结构体。图1显示了含晶粒结构的圆柱体岩石模型的构建过程(不同的颜色表示不同的颗粒组分)。

图1 岩石等效晶质模型构建过程Fig.1 Construction process of rock equivalent crystal model

2 3D Voronoi晶粒结构的优化研究

本节分析了3D Voronoi晶粒排列均匀化程度及大小等微观结构的变化对模型宏观力学性质的影响程度，并对模型中的晶粒结构进行优化。

选用粒径比为1∶1，1∶1.5，1∶2，1∶2.5及1∶3(半径均值4.5 mm)呈均匀分布“粗”颗粒集构建长×宽×高为150 mm×150 mm×150 mm的晶粒网，其晶粒个数分别为7 073，7 885，6 322，5 959和5 458个。

按照上述方法，沿晶粒网xy平面不同位置随机构建3个圆柱体(φ50 mm×100 mm)模型，并对平行黏结接触和光滑节理接触赋予表1中的细观参数。

表1 颗粒及接触模型细观参数Table 1 Microscopic parameters of particle and contact model

对圆柱体模型进行单轴压缩(拉伸)试验，得到的宏观力学参数与颗粒粒径比的关系如图2所示。在单轴压缩试验中，模型的抗压强度受粒径比的影响较小，但随着粒径比的增大其离散性随之增大；模型的弹性模量随粒径比增大逐步增长，但增长趋势逐渐减小；泊松比波动性较大，但除了粒径比1∶1时岩样的泊松比显著低于其他水平外，泊松比整体受粒径比的影响并不明显。在单轴拉伸试验中，粒径比1∶1时模型抗拉强度明显低于其他水平，除此在外，模型抗拉强度受粒径比的影响较小，但随着粒径比的增大其离散性随之增大。当1∶2.5时弹性模量趋于稳定，且抗压(拉)强度的离散程度在控制范围之内，因此，本文选取“粗”颗粒粒径比为1∶2.5。

为了研究晶粒体大小对模型宏观力学性质的影响，在上述所构建的晶质等效模型体中(粒径比1∶2.5)，沿模型体中心位置分别取φ30 mm×60 mm,φ40 mm×80 mm,φ50 mm×100 mm和φ60 mm×120 mm这4种尺寸的圆柱体模型进行单轴压缩试验。不同尺寸模型轴向应力-应变曲线如图3所示。随着模型尺寸的增大，应力-应变曲线趋于平稳，晶粒尺寸的影响逐渐减小，结合计算效率，本文选取φ50 mm×100 mm的岩样模型。

图2 宏观力学参数与颗粒粒径比的关系Fig.2 Relationship between macro-mechanical parameters and particle size ratio

图3 不同尺寸模型轴向应力—应变曲线Fig.3 Axial stress-strain curves of models with different sizes

3 岩石破坏特征细观研究

为了验证3D Voronoi晶粒模型适用性和观察其破坏时的细观特征，本文选取某黑钨矿围岩——粉砂质板岩制成尺寸为φ50 mm×100 mm和φ50 mm×40 mm的岩样，通过室内单轴压缩和巴西试验得出宏观物理力学参数值。利用扫描电镜(SEM)对岩样受压破坏后的断口观察，如图4所示,岩石断裂口呈现出穿(沿)晶断裂的组合形态，图4中圈定区域内以穿晶断裂为主，断裂面不平整，多为撕裂痕迹，其他区域以沿晶断裂为主，断裂面相对平整光滑，穿晶断裂与沿晶断裂相互穿插呈无规则分布。

图4 岩样断裂断口断裂特征Fig.4 Characteristics of rock sample’s rupture surface

将岩样切割打磨形成5组岩样，利用x-射线仪对岩样矿物组成成分进行分析。

该板岩主要由石英、黑云母、长石、绿泥岩及少量矿物成分组成，岩样矿物组成成分半定量分析结果，如表2所示。

表2 岩样主要矿物成分分析结果Table 2 Analysis results of main mineral components in rock samples %

根据表2岩样矿物成分分析结果，选择表征石英、黑云母、长石、绿泥石占比为28%，45%，22%和5%，粒径为2.57～6.43 mm(粒径比1∶2.5)均匀分布的“粗”颗粒集构建的晶粒网，并覆盖至尺寸为φ50 mm×100 mm 和φ50 mm×40 mm的细颗粒(最小颗粒半径1.11 mm,粒径比1∶1.66)圆柱体中，组成晶质等效模型。

利用上述构建的粉质板岩晶质等效模型进行单轴压缩和巴西试验模拟。设置伺服加压(拉)速率0.3，采用试错方法对PFC3D细观参数进行标定，标定后的矿物颗粒及平行黏结模型细观参数见表3，光滑节理接触模型细观参数沿用表1中的参数。最终得到了与室内试验结果相近的宏观力学参数，如表4所示。图5～6为应力-轴应变曲线图。从图中可以看出，数值模拟的曲线斜率和峰值强度均与室内试验的曲线接近。峰后均表现出宏观脆性破坏特征,如图7所示。从而验证了3D Voronoi等效晶质模型的适用性。

表3 矿物颗粒及平行黏结接触模型细观参数Table 3 Microscopic parameters of mineral grain and parallel bond contact model

在单轴压缩试验数值模拟过程中，裂纹沿岩样晶粒结构面随机产生，并缓慢增长，该阶段应力-应变呈现线性关系；当应力接近峰值的50%时，沿晶断裂出现快速生长和扩展，同时穿晶断裂点在加载板接触面附近出现和缓慢增长；当应力接近峰值的90%时，沿晶断裂贯通形成主裂纹，穿晶断裂生长速度加快，主要分布在周边结构面发生断裂和错动的晶粒内部，特别是力学性质较弱的黑云母和绿泥石晶粒内部穿晶断裂较明显，这是由局部应力不平衡而引起的，穿晶断裂扩展和贯通导致了局部失稳，该阶段应力与应变的非线性关系明显；峰值后，应力陡然下降，该阶段沿晶和穿晶断裂增长速率最快，局部失稳演化为宏观断裂，断裂面呈现沿晶和穿晶断裂的组合形态。岩样在单轴压缩破坏过程中以拉伸断裂为主,在宏观上表现为脆性断裂，图8为单轴压缩试验不同阶段的断裂数统计图。

表4 粉砂质板岩宏观参数试验值和模拟值Table 4 Tested and simulated macroscopic parameters of silty slate

图5 单轴试验轴向应力(裂纹增长)—应变曲线Fig.5 Axial stress (crack growth)- strain curves of uniaxial tests

图6 巴西试验径向应力(裂纹增长)—应变曲线Fig.6 Radial stress (crack growth)-strain curves of Brazilian tests

图7 岩样宏观破坏后的对照Fig.7 Contrast diagram of of rock samples after macro-failure

图8 单轴试验不同阶段断裂数统计图Fig.8 Statistical chart of fracture number at different stages of uniaxial tests

巴西试验数值模拟过程中，其裂纹的生成同样经历了稳定发展、峰值前的快速增长及峰值后增长速率最快3个阶段，但沿晶断裂首先从加载板接触面附近生成，然后沿加载板垂直方向扩展和贯通形成主裂纹，随着晶粒周边结构面发生断裂和错动，诱发晶粒内部的穿晶断裂，穿晶断裂扩展和贯通导致了局部失稳，应力达到峰值，最终形成与加载板大致呈垂直状的宏观断裂，岩样在巴西劈裂破坏过程中同样以拉伸断裂为主,在宏观上表现为脆性断裂，图9为巴西试验不同阶段的断裂数统计图。

图9 巴西试验不同阶段断裂数统计图Fig.9 Statistical chart of fracture number at different stages of Brazilian tests

4 结论

1)基于3D Voronoi图，结合PFC3D中颗粒平行黏结及光滑节理模型，构建了一种能够反映岩石矿物晶粒结构三维精细化建模方法，为岩石破坏特征的细观研究提供了一种数值分析方法。

2)分析了晶粒排列均匀化程度及大小等微观结构的变化对模型宏观力学性质的影响程度，随晶粒不规则程度的提高，弹性模量逐步增大，单轴抗压(拉)强度离散型越来越明显。晶粒大小不变的情况下，随着模型尺寸的增大，晶粒结构对模型的宏观力学性质影响逐渐减小。

3)通过裂纹生长经历的稳定发展、快速增长及峰值后增长速率最快3个阶段可以大致表征室内试验弹性变形阶段、非弹性变形阶段和峰值后宏观破坏阶段；岩样在单轴压缩和巴西劈裂破坏过程中以拉伸断裂为主,在宏观上表现为脆性断裂；岩样的宏观断裂呈现沿晶和穿晶断裂的组合特征，弹性变形阶段裂纹的萌生和扩展以沿晶断裂为主，而穿晶断裂的扩展和贯通往往导致局部失稳，表现为非弹性变形阶段。