土石混合体随机结构模型生成与直剪强度数值试验研究

罗伟，王优，张帅浩，程肖，黄姗，荣耀

土石混合体随机结构模型生成与直剪强度数值试验研究

罗伟1, 2，王优3，张帅浩3，程肖3，黄姗3，荣耀2

(1. 华东交通大学 土木建筑学院，江西 南昌 330013；2. 江西省交通科学研究院，江西 南昌 330200；3. 中南大学 土木工程学院，湖南 长沙 410075)

为获取土石混合体力学数值试验的样本，提出基于AutoCAD二次开发的土石混合体随机结构模型生成方法。利用数字图像处理技术提取块石轮廓和尺寸数字信息，并将其存储为可供AutoCAD读取的脚本文件格式，得到可供缩放的AutoCAD图形交互文件。利用AutoCAD强大的块处理功能进行VBA编程建块，以及AutoCAD软件平台具有的相交和内含快速判别功能，实现了土石混合体中块石的高效率投放和土石混合体数值模型建立。针对某一土石混合体级配曲线进行直剪数值试验，验证了本文方法的有效性。试验结果表明：土石混合体材料在初始阶段存在明显的“绕石效应”，在塑性剪切应变逐步扩展至贯通过程中，土石混合体逐步出现塑性变形而发生剪切破坏。

土石混合体；数字图像处理；块石投放；随机结构模型；直接剪切强度；数值试验

岩土体细观结构的构建与力学性质分析一直是岩土工程领域研究的热点，国内外学者分别从多个层面开展了一系列研究[1−4]。而随着计算机技术，尤其是数字图像处理技术[5]的发展，为学者在岩土体细观结构定量分析方面提供了一种新的研究思路。Tovey等[6]基于数字图像处理技术，开发了一种土体微结构分析系统，根据不同类型土体的扫描图像，对土体的微观结构特性进行定量分析，并得到了一些有价值的结论。在土石混合体研究方面，徐文杰等[7−9]采用数字图像处理技术对土石混合体的细观结构特征与力学性质之间的关系进行了定量分析。数字图像处理技术在岩土体细观结构的定量研究中得到了越来越广泛的应用。大尺度的直剪试验研究结果表明[8]，土石混合体的内部块石含量及块体特征，在很大程度上控制了其变形和破坏模式。获取随机块石图形交互文件，是一个涉及多学科、极为复杂的过程，主要有2种思路。一种是通过数字图像处理获得单颗粒块石信息录入块石库，然后通过数值投放获得图形交互文件。基于此思想，油新华[10]采用规则几何体圆形随机生成了土石混合体的二维结构，Graziani等[11]使用基于规则形状岩石块(如：圆形，三角形，矩形等)对二维土石混合体随机模型进行了一系列的数值模拟，这虽使工作量减小，但与实际情况存在较大差异；WANG等[12]提出“提取与放置”(take and place)的方法，即提取一个多边形，通过判断轮廓上的点是否位于其他多边形内以达到判别块体间是否出现相交和内含的目的。但由于轮廓上的已知点数有限，所以此法不能保证将所有相交情况囊括在内，且所取点越多，结果越精确，这无疑会大大增加计算量，因此，此方法有待进一步改进。第2种思路是通过对土石混合体断面进行图像处理直接获得相应的图形交互文件。徐文杰等[13]采取此法虽使得块石的形状由规则可以变得不规则，但由于其相交判断法则的局限性使得最后投放的块体只能为凸多边形，尚未准确描述实际块石的形状。基于此，本文采取上述第1种思路，编制MATLAB程序对拍摄的单个颗粒块石图片进行二值化、降噪处理，提取像素点坐标，并将其转化为块石轮廓的离散点坐标。随后求解块石中心点坐标和对应的“极径”，并将块石数据点坐标写成AutoCAD可以读取的形式。然后，基于AutoCAD中强大的块处理功能，在二次开发环境VBA中编制相关程序[14]来实现建块过程，并对块石进行缩放。最后根据MATLAB中计算出的每类块石的投放数目，在AutoCAD中实现投放。由于在AutoCAD里可以方便准确地判断块体是否相交和内含，解决了投放时遇到的最大问题，因而极大提高了投放的效率。投放结束后生成图形交互文件，方便土石混合体数值模型的建立。另外，本文针对某一土石混合体级配曲线进行试验，得到的投放结果应用于直剪数值实验，验证了此法的有效性。

1 基于数字图像处理的二维块石 构建

1.1 块石图像预处理

数字图像，即二维图像用有限数值的像素来表示。在灰度图像中每个像素点对应一个整数值，称为灰度值。常见的有二值化和256色，其对应的灰度值分别为0~1和0~255。灰度值的不同表征了图像所包含的信息不同，这些离散的像素点对应的灰度值构成了数字图像处理的基础。本文以对单颗粒块石拍摄数字图像的方式获得块石可视图像，通过编制MATLAB程序，将可视图像转化为二值图像并对其进行降噪处理，然后提取像素点坐标，得到块石轮廓的离散点坐标，将块石轮廓线信息写入Excel表格，将其导入VBA中的blockbuilding程序，便可得到一个整体的块石图块。

在块石图像预处理过程中，首先在彩色图像上进行处理，通过去掉背景来去除噪声和部分阴影的影响。然后对图像进行高斯模糊，高斯模糊的主要功能是把某一高斯曲线附近的像素值进行统计加权得到该曲线的像素值，便于得到物体轮廓。根据多次调试，模糊半径取15个像素时获得的块石轮廓最为清晰。MATLAB程序中的二值化主要是通过自动选取图像分割阈值的方法来实现，最终可将图像分割为目标和背景。最后对二值化后图像中的目标进行颗粒分析，提取其轮廓线和尺寸信息，并将其存储到可供AutoCAD读取的脚本文件格式，读入AutoCAD得到图形交互文件，如图1所示。

1.2 像素与实际尺寸转化比例

为获取块石实际尺寸，在拍照过程中加入“标定块”(本文中取1 cm×1 cm)，其像素面积为2 383，实际面积为1 cm2，据此可推算其他块体面积，计算结果如表1所示。

(a) 原始块石图像；(b) 去背景处理后块石图像；(c) 高斯模糊处理后块石图像；(d) 二值化处理后块石图像；(e) 块石轮廓线；(f) 块石图形交互文件(AutoCAD)

表1 块石像素面积与实际面积

1.3 批量创建图块

AutoCAD中的图块是一组图形实体的总称，在应用的过程中，图块是一个独立的对象。可以根据需要通过缩放图块和调整图块角度将其插入到区域中的任一位置。

将块石的实体写入图块，能够将每个块石视为一个独立完整的对象，而不仅仅是多条多段线或多个三角面。并且块石创建图块后，能够便捷地实现块石的投放(将图块插入到任一指定位置)，辅以相应的相交内含判断准则便能够实现土石混合体模型的构建。而块石的大小和主轴的方向均可在图块插入时定量控制。

2 块石随机投放

在对多边形块石形态分析之后，多边形块石的所有顶点坐标，将从局部坐标转化为全局坐标，进行块石的投放。该过程的实现在已有文献[12]大多采用“提取与放置”(take and place)的方法，即提取一个块石，并将其随机放置于虚拟试件(全局坐标系)中的某一位置。这一过程的实现往往伴随着较长的计算时间，投放效率较低，且从理论上不能完全避免块石的重叠。

基于此，本文块石的投放是通过在AutoCAD二次开发环境VBA中编制相关程序来实现，Auto CAD具有强大的块处理功能，从理论上保证了块石不发生重叠。具体步骤如下：

1) 根据块石的粒径确定块石种类，然后由含石量、块石中心点的位置范围以及每种块石的级配系数等参数得到每一种块石所需投放的数目及面积；

2) 使用上述数字图像处理的方法得到一批符合要求的石头；

3) 将块石轮廓上散点的坐标写入Excel文件并导入VBA中，完成块的创建；

4) 将每种块石的信息写入取石报告中，运行VBA程序便可以实现块石的随机投放。

由于较大粒径的块石更容易和已经存在的其他块石发生重叠，因此，在实际操作中，总是根据块石级配曲线，优先投放较大粒径的块石。块石投放实现流程图如图2所示。

2.1 块石相交判定

图2 块石投放流程图

图3 块石相交情形图示

在块石相交判断的基础上进行块石内含判断，块石内含情形如图4所示。运用AutoCAD内置的GetBoundingBox函数得到块石边框的最大、最小点坐标值，如图5所示。若2块石没有相交，且max2min1&min2>min1，则2块石为内含关系。反之，则不出现内含。

图4 块石内含情形图示

图5 块石内含判定图示

由于本文真实块石来源为道砟，因此所取颗粒的粒径集中分布于一个较窄的区间内，这与自然界中的实际情况不符，所以需要对块石进行缩放以符合级配要求。

图6 不同级配的块石粒径曲线

2.3 块石缩放

自然情况下，石头粒径分布范围很广。本文通过拍照获得的块体粒径往往相差较小，投放时会与真实情况不符，因此在根据一条级配曲线进行数值投放时需要对块体进行缩放以让投放结果尽可能接近实际情况。对某一条具体的级配曲线，可以先将粒径等间距地分为若干段，然后计算每一段的平均值，接着根据前述已经处理的块石粒径的均值，以两者之比作为缩放比例。通常，缩放的思路如下：将平面块体划分为以原点为公共顶点的三角形，然后求出每个三角形的形心坐标和面积，最后便可得出整个块体的形心位置。为检验投放效果，本文选取了3条级配曲线如图6所示。图7为对应的投放结果，3个模型的含石量均为40%。

(a) 级配1；(b) 级配2；(c) 级配3

3 数值模型的建立

在得到块石投放的交互文件后，我们便可根据需要进行相应的数值实验。在块石边缘轮廓线提取过程中，为保证识别精度，块石边缘由大量像素点连接构成，在网格划分时每一个像素点将对应网格模型的一个节点，这势必将导致网格数量急剧增加，使得数值试验效率降低，甚至失败。即便在经过前述的高斯模糊处理后每个块石上仍有1 000多个点。

图8(a)为边缘未简化情况下对块石直接进行网格划分，在划分过程中可以发现该块石轮廓线由1 233条线段组成，划分完成后发现该块石单元数量为35 855，节点数目为36 469，单个块石像素点数目如此庞大，当模型块石含量较多时，很难开展相应的数值试验。所以对块石边缘像素点进行适当的简化显得极为必要。本文通过编制相应的MATLAB程序，有效减少了块石边缘轮廓线像素点。主要思想是利用插值，即计算出各轮廓点相对于形心的坐标然后在0到2π内以一定步长进行极径的插值，这样就能在保持块体原有特征的前提下，缩减轮廓点个数。对边缘轮廓线像素点进行简化，虽然可能会使得块石识别精度有所降低，却大大提高了建模效率，同时也使得包含大量块石的数值模型试验成为可能。在有效保证识别精度的前提下，经过多次尝试将块石轮廓线上的像素点数目降至36个，以保证建模及数值试验效率。

添加边缘简化程序后块石网格划分如图8(b)所示，该块石单元数量为126个，节点数目为75个。图9为对边缘像素点进行简化后的块石在100 cm×50 cm的区域内某次的投放效果图。

(a) 未简化；(b) 简化

3.2 土石混合体数值模型直剪试验

为了模拟土石混合体这类特殊地质体在外荷载作用下的变形破坏特征，本文将运用有限差分软件FLAC3D对上述建立的土石混合体模型进行直接剪切试验[15]。

图9 边缘简化块石某次投放结果矢量图

3.2.1 几何模型

图10所示为截取出的土石混合体剪切计算模型及网格模型，网格在划分时采用四节点四边形单元，划分完成后再延纵向(垂直纸面方向)拉伸一个单元长度，形成三维网格模型(假三维)，该模型对应的单元数量为95 460，节点总数为49 917。

(a) 直剪试验计算模型；(b) 土石混合体网格模型

3.2.2 边界条件

计算中模型边界条件为：模型底部为纵向、横向及法向全约束；模型下半部分两侧进行横向约束；在模型上半部分左侧施加了一定的位移荷载；在模型的顶部施加相应的法向荷载，如图10(a)所示。土石混合体试样与剪切盒之间用接触面模拟。按照作用机理，FLAC3D中包含3种不同类型的接触面，分别是黏性接触面、滑动性黏结接触面和库伦滑移接触面，按照工程实际中土石接触情况，我们选取库伦滑移接触面来进行研究。库伦滑动接触面是一种仅有完好和破坏2种对立状态的黏性接触面。发生破坏时，接触面单元的力学行为由内摩擦角、黏聚力和接触刚度决定，黏结强度在没有设定时其值默认为0，破坏后的黏结单元不能承受拉 应力。

库伦剪切强度的表达式如下：

式中：为接触面表面的黏聚力；为接触面表面内摩擦角；为接触面面积；F为接触面的法向力；为孔隙水压力。本文中不考虑水的作用，即设置0。界面强度参数如表2所示。

表2 界面强度参数取值表

3.2.3 材料参数

本文基于已有文献[16−17]，确定了土体和块石参数取值，如表3所示。

表3 土石混合体计算参数

3.2.4 结果分析

在剪切试验过程中土体及块石均采用摩尔−库伦本构模型。为便于对比分析，对均质土体及土石混合体同时进行直剪试验。

从图11和图12可以看出，均质土体的剪切带是一条规则的直线，土石混合体的剪切带则较为宽厚和曲折，且模型内部塑性区具有明显的“绕石效应”，即塑性区绕开石头延伸。

(a) 级配1；(b) 级配2；(c) 级配3

图12 均质土体剪切试验结果

通过在接触面上设置监测点来获取剪切过程中的位移和应力。以模型的左下角为原点，监测点的坐标分别为(20,25)，(40,25)，(60,25)，(80,25)(单位：cm)，最终的位移和应力值为4个监测点的平均值。实验所得剪切应力和剪切位移关系曲线见图13。

在剪切的初始阶段，土石混合体和均质土体的剪切曲线几乎是直线，表明两者皆处于弹性阶段，但土石混合体的直线斜率较大，说明土石混合体的抗剪强度较大。土石混合体和均质土体的剪切曲线在经历直线段后变化缓和，进入塑性屈服阶段，并逐渐达到峰值强度，随后开始下降，并逐渐达到相应的残余强度。在整个剪切过程中，土石混合体的剪切曲线整体在均质土体之上，说明土石混合体的抗剪强度由于块石的存在而得到了提升。在相同的含石量下，不同级配对土石混合体的剪切特性均有较大影响。

图13 剪切应力随剪切位移的变化曲线

4 讨论

1) 本文没有考虑块石与块石之间直接接触。虽然块石分离与块石接触的投放判断过程有较大差异，且块石接触有点点接触、点线接触、线线接触，情况更为复杂，但按照本文研究思路可以构建考虑块石接触的数值模型。块石与块石接触对土石混合体的力学性质产生的影响将在后续工作中开展。

2) 采用本文思路可实现真三维数值分析模型的构建，其中如何实现三维块石轮廓重构及三维块石相交判定是研究的难点。目前，基于多视几何方法，通过对三维块石多个面拍照，实现照片对块石实体全覆盖。对照片进行特征点提取及配准，重建出多个三维空间点，从而可以构建出块石的三维网格。另外，在块石相交判定方面，虽然判定原理相似，由于构成三维块石的基本元素相比二维块石更为庞大，因此其实现过程也更为复杂。

5 结论

1) 通过数字图像处理技术可以提取到块石轮廓上的点，再通过MATLAB编程将坐标写成AutoCAD可以读取的格式，从而实现建块及投放过程。

2) 利用AutoCAD的内置函数编制的相交及内含判定程序原理较为完善，因此可以获得较理想的投放结果。

3) 非均质的土石混合体以及均质土体的直剪数值实验表明，均质土体的剪切带是一条规则的直线，土石混合体的剪切带则较为宽厚和曲折，且样本内部塑性区具有明显的“绕石效应”。

4) 在含石量相同的情况下，级配对土石混合体的塑性区扩展及剪切特性均有较大的影响。

[1] GONG Jian, LIU Jun. Analysis on the mechanical behaviors of soil-rock mixtures using discrete element method[J]. Procedia Engineering, 2015, 102: 1783−1792.

[2] XU Wenjie, XU Qiang, HU Ruilin. Study on the shear strength of soil–rock mixture by large scale direct shear test[J]. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(8): 1235−1247.

[3] Vallejo L E, ZHOU Y. The mechanical properties of simulated soil–rock mixtures[C]// Proceedings of the 13th International Conference on Soil Mechanics and Foundation Engineering. New Delhi, India; [s.n.], 1994: 365–368.

[4] WANG Y, LI X, WU YF, et al. Experimental study on meso-damage cracking characteristics of RSA by CT test[J]. Environmental Earth Sciences, 2015, 74(6): 5023−5037.

[5] 涂望明, 魏友国, 施少敏. MATLAB 在数字图像处理中的应用[J]. 微计算机信息, 2007, 23(6): 299−300.TU Wangming, WEI Youguo, SHI Shaomin. Application of MATLAB to digital image processing[J]. Microcomputer Information, 2007, 23(6): 299−300.

[6] Tovey N K, Krinsley D H, Pye K, et al. Backscattered scanning electron microscopy and image analysis of sediments and sedimentary rocks[J]. Journal of Paleolimnology, 2000, 23(2): 223−224.

[7] 徐文杰, 岳中琦, 胡瑞林. 基于数字图像的土、岩和混凝土内部结构定量分析和力学数值计算的研究进展[J]. 工程地质学报, 2007, 15(3): 289−313. XU Wenjie, YUE Zhongqi, HU Ruilin. Current status of digital image based quantitative analysis of internal structures of soils, rocks and concretes and associated numerical simulation[J]. Journal of Engineering Geology, 2007, 15(3): 289−313.

[8] 徐文杰, 胡瑞林, 岳中琦, 等. 基于数字图像分析及大型直剪试验的土石混合体块石含量与抗剪强度关系研究[J]. 岩石力学与工程学报, 2008, 27(5): 996−1007. XU Wenjie, HU Ruilin, YUE Zhongqi, et al. Research on relationship between rock block proportion and shear strength of soil-rock mixtures based on digital image analysis and large direct shear test[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(5): 996− 1007.

[9] 徐文杰, 胡瑞林, 岳中琦, 等. 土石混合体细观结构及力学特性数值模拟研究[J]. 岩石力学与工程学报, 2007, 26(2): 300−311. XU Wenjie, HU Ruilin, YUE Zhongqi, et al. Mesostructural character and numerical simulation of mechanical properties of soil-rock mixtures[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(2): 300−311.

[10] 油新华. 土石混合体的随机结构模型及其应用研究[D]. 北京: 北京交通大学, 2001. YOU Xinhua. Stochastic structural model of the earth- rock aggregate and its application[D]. Beijing: Beijing Jiaotong University, 2001.

[11] Graziani A, Rossini C, Rotonda T. Characterization and DEM modeling of shear zones atalarge dam foundation[J]. International Journal of Geomechanics, 2012, 12(6): 648– 664.

[12] WANG Z M, Kwan A K, Chan H C. Mesoscopic study of concrete: Generation of random aggregate structure and finite element mesh[J]. Computers and Structures, 1999, 70(5): 533−544.

[13] 徐文杰, 胡瑞林, 岳中崎. 土–石混合体随机细观结构生成系统的研发及其细观结构力学数值试验研究[J]. 岩石力学与工程学报, 2009, 28(8): 1652−1665. XU Wenjie, HU Ruilin, YUE Zhongqi. Development of random mesostructure generating system of soil-rock mixture and study of its mesostructural mechanics based on numerical test[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(8): 1652−1665.

[14] 张帆, 郑立楷, 王华杰. AutoCAD VBA开发精彩实例教程[M]. 北京: 清华大学出版社, 2004. ZHANG Fan, ZHENG Likai, WANG Huajie. Development of excellent example tutorials based on AutoCAD VBA[M]. Beijing: Tsinghua University Press, 2004.

[15] HU Q J, SHI R D, ZHENG L N, et al. Progressive failure mechanism of a large bedding slope with a strain- softening interface[J]. Bulletin of Engineering Geology and the Environment, 2018, 77(1): 69−85.

[16] XU W J, WANG S, ZHANG H Y, et al. Discrete element modelling of a soil-rock mixture used in an embankment dam[J]. International Journal of Rock Mechanics & Mining Sciences, 2016, 86: 141−156.

[17] GU T F, ZHANG M S, WANG J D, et al. The effect of irrigation on slope stability in the Heifangtai platfom, Gansu province, China[J]. Engineering Geology, 2018, 248(8): 346−356.

A method for generating random structure model of soil-rock mixtures and study of its direct shear strength based on numerical test

LUO Wei1, 2, WANG You3, ZHANG Shuaihao3, CHENG Xiao3, HUANG Shan3, RONG Yao2

(1. School of Civil Engineering and Architecture, East China Jiaotong University, Nanchang 330013, China; 2. Jiangxi Transportation Institute, Nanchang 330200, China; 3. School of Civil Engineering, Central South University, Changsha 410075, China)

In order to obtain samples of the numerical mechanics tests of soil-rock mixtures, a method for generating random structural model of soil-rock mixtures was proposed based on the secondary development of AutoCAD. The contour and size information of multi rock blocks were extracted using the digital image processing technology. It was stored as a script file format for AutoCAD read, and a scalable AutoCAD graphic interactive file was obtained. Blocks were built by VBA programming based on the powerful block processing function of AutoCAD. The high efficiency distribution of rock blocks and the establishment of the numerical model of the soil-rock mixture were achieved by using the intersecting fast distinguishing function of the AutoCAD software platform. Direct shear numerical tests on a gradation curve verify the effectiveness of the proposed method. Shear failure tests of the soil-rock mixtures show that there is an obvious “rock-around effect” in the initial stage of the soil-rock mixtures. During the gradual expansion of the plastic shear strain, the soil-rock mixtures appear plastic deformation and eventually failures.

soil-rock mixtures; digital image processing; distribution of block stones; random structural model; direct shear strength; numerical test

TU457

A

1672 − 7029(2019)07− 1681 − 09

10.19713/j.cnki.43−1423/u.2019.07.11

2018−10−10

国家自然科学基金资助项目(51878668)；贵州省交通运输厅科技项目(2017123033)；江西省交通运输厅科技项目(2016C0058，2016C0007)

罗伟(1986−)，男，江西东乡人，博士，从事道路与铁道工程、岩土构筑物稳定性分析；E−mail：luoweicsu@126.com

(编辑 涂鹏)