土的局部变形如剪切带的形成规律和渐进破坏特性是土体区别于常规连续介质的重要性质[1-2]。许多学者在多孔介质的局部变形理论、应力应变试验方法和数值模拟方面做了一系列的研究工作[3-5]。研究土体的试验和模拟手段，对于完善计算土力学理论和防范地基工程灾害具有重要的理论意义和实用价值。
土体是典型的多相多孔材料，其宏观力学响应与土体复杂的微观结构密切相关。岩土力学问题研究中一般是依据室内试验确定土体参数，进而采用数值计算或理论模型预演工程实际[5-6]。室内试验是进行土体结构参数确定的必要手段，数值计算是预演土体结构受力变形行为的实用途径[6]。为了准确确定数值计算中的模型参数，实现室内试验中碎散、多相土体应力变形的整体-局部测量，许多学者将各种试验方法应用在土体的局部变形测量中。例如，局部位移传感器[7-8]，CT扫描[9-11]，数字图像测量[12-15]等。相较于接触式和外辐射式测量，图像测量技术具有无接触、无辐射、布置灵活、适用性高等优点，因而成为测量变形场的主流方法[15-16]。图像测量技术的后处理是影响测量结果和获取土体参数的重要环节，常用的材料变形场测量技术有数字图像相关技术(DIC)和粒子图像测速法(PIV)[16-20]，两者原理相同，可以在测试区域的数字图像中得到大量的可追踪点，通过图像中区域的灰度值相关算法进而可准确地反映出材料的变形场[16-17]，其中粒子图像测速法(PIV)[18]常被应用于具有颗粒流动特征的试验中。基于数字图像相关技术(DIC)的测量方法也被一些学者应用在岩石与岩土研究中。如，Tang等[19]利用三维数字图像相关技术(3D-DIC)研究岩石在不同围压下的渐进破坏过程。王学滨等[20]利用DIC在无侧限压缩试验中观测砂土试样在不同含水率下剪切带的发展过程。Yang等[21]利用3D-DIC技术研究混凝土中集料尺寸与形状对混凝土不均匀收缩的影响。赵燕茹等[22]利用DIC技术观测三点弯试验中切口混凝土梁全场变形并分析混凝土梁在冻融条件下断裂破坏过程中水平位移和应变的变化规律。同时随着计算机技术的发展，基于离散元理论的数值模型逐步被用于研究岩土、岩石等颗粒材料的破坏机制研究中[23-25]。数值模型中颗粒与颗粒之间采用相应的接触模型粘结起来，颗粒之间的相互作用受到接触模型细观参数的约束，颗粒的整体运动轨迹反映材料的宏观变化。其中，唐欣薇等[26]利用数字图像技术获得花岗岩细观结构的表征，结合离散元技术，建立反映岩石细观非均质组构特性的颗粒离散元模型。
为了更真实地反映平面应变状态下土的受力变形特性，合理确定数值计算中所需的细观参数。研制了一种新型的平面应变加载设备，该设备通过对试样的侧向(围压方向上)施加柔性荷载来降低常规平面应变试验中刚性加载所造成的边界约束影响。同时，搭建了能够得到表面变形识别的数字图像采集系统。在此基础上，利用研制的平面应变设备结合二维数字图像相关技术(2D-DIC)，给出了基于图像测量的福建标准砂的宏细观参数确定方法，从而搭建了室内试验-数据处理-参数确定-离散元模拟的计算系统。基于以上系统组成，对福建标准砂的平面应变特性进行了室内试验和离散元模拟。
1 平面应变试验系统
1.1 试验设备
试验仪器是由大连理工大学工业装备结构分析国家重点实验室主导，自主研发的平面应变试验设备，仪器的结构部分如图1(a)~图1(b)所示。试样形状为长方体，在高度、宽度与厚度方向上的尺寸为120 mm×60 mm×100 mm，试样高宽比为2∶1，试验尺寸及其在仪器中的放置状态如图1(c)所示。试样由黑色的乳胶膜包裹，试样前后由刚性板约束来限制试样ε2方向的变形。试样前表面安装了一块透明的钢化玻璃，数字相机可通过玻璃视窗记录试验的变化过程。
与传统的平面应变试验机中各个方向都采用刚性加载的方式不同，本试验机在σ3方向的加载方式采用柔性加载。土样的两个侧面与两个乳胶膜制成的方形柔性水囊接触，水囊内部充满无气水，并由精密水压力控制器对柔性水囊内部施加压力，为试样提供围压σ3。通过试验测试，柔性水囊可以承受400 kPa范围内的水压力。采用柔性水囊加载存在以下优点：一是可以减小加载过程中试样变形时试样与加载装置的摩擦；二是柔性水囊提供的受力分布接近实际土层中的受力状态。在试验机的轴向(σ1方向)方向上，在压力室的上部与底部分别布置2个称重传感器，用2个传感器相互校准减小测量中摩擦力的影响。
图1 平面应变试验装置
Fig.1 Plane strain experimental apparatus
1.2 试样制备及试验条件
试验用砂土采用福建标准砂(中国)，福建标准砂的颗粒分布曲线如图2所示，其基本物理指标如表1所示。采用激振法分四层制平面应变试验试样，制好的试样由试样膜包裹，安置在压力室内，为了减小乳胶膜和试样前后面板的摩擦力，在试样乳胶膜前后涂抹透明的硅脂。试验采用位移控制的方式加载，加载速率为0.2 mm/min。试验的编号从T-1至T-3，试验围压分别为50 kPa、100 kPa、150 kPa。试验状态为干砂剪切试验，如表2所示。
图2 福建标准砂颗粒分布曲线
Fig.2 Grain size distribution curves for Fujian standard sand
表1 试样材料的基本特性
Table 1 Properties of experimental sands
初始密度ρd/(g/cm3)平均粒径D50/mm最大空隙比emax 最小空隙比emin 1.57 0.37 1.04 0.60相对密度Dr/(%)比重Gs不均匀系数Cu曲率系数Cc 87 2.65 5.49 0.68
表2 试验条件及试验测量结果
Table 2 Experimental conditions and measurement
试验编号围压/Kpa试验状态摩擦角φ′/(°)剪切带倾角β/(°) 剪切带宽度T/D50 T-150 干砂41 58 18 T-2100干砂41 58 17 T-3150干砂41 55 14
1.3 数字图像测量系统
采用数字相机通过玻璃视窗记录试样的变形过程，如图1所示，在试样黑色乳胶膜表面喷上白色的散斑图案，便于后续使用数字图像相关法程序的分析。用于采集的相机采用型号为：Pike F-100B/C，像素分辨率为1000×1000，每像素约代表0.11 mm (根据每次试验标定的情况而定)，每5 s中保存一张图像。为了保证采集的照片光照均匀度更好，采用了两个LED灯，从两侧为试样照明。
2 砂土的平面应变试验全场变形测量
2.1 测量区域的选取
数字图像相关法是一种基于图像灰度特征来实现追踪的图像测量方法。该方法的优点是计算局域内每一个像素的灰度信息都会被使用到，是一种适合变形场的稠密点跟踪算法[16, 27-28]。本试验中，采用的数字图像相关法的精度可以达到0.01像素。数字图像相关法是以参考图像为基准，分析变形后的图像相对于参考图像的变化。本试验中均选用初始状态试验图像作为参考图像，试样变形场的计算区域大于整个试样表面的70%，试样数据分析区域如图3所示。
图3 试验数据分析的相关区域
Fig.3 The analyzed area of tests
图3中试样形状为长方体，在厚度方向上试验的变形被刚性板约束在100 mm。在此平面应变试验中，假设厚度方向上的变形与前表面一致，只考虑试样高度(竖向方向)和宽度(水平方向)的变形，利用二维的数字图像相关法可以确定图3中所示变形场计算区域的竖向应变(εy)、竖向位移(Vy)、水平应变(εx)、水平位移(Vx)。采用图3所示的变形场计算区域的面积变形来表示试样的体积变形。试样应力应变曲线中的轴向应变用图3中所示的轴向应变计算区域竖向变形的平均值来表示，该区间选择在试样的中间部分，以消除试样的端部约束。
2.2 试样变形场的测量结果
图4给出了T-3试验中试样在拉格朗日坐标系下的竖向位移场与应变场和渐进破坏过程，图4(a)中标注出了测量区域的范围(单位/mm)。可看出，该试样的破坏是从一个小区域开始，慢慢地贯穿整个试样。变形场以试验开始时初始图像为基准，以压缩方向为正方向。在利用数字图像相关法分析时，采用的子区大小为21像素×21像素，用于跟踪变形区域的目标点超过40000个。图4中的区域是图3中所示的变形场计算区域，图4(a)~图4(e)所示的位移场与应变场对应的是图4(f)中应力-应变曲线上红色标记点的位置，每幅图下方的应变数值是试样整体轴向应变值，取试样中间部分的轴向变形平均值，在图4(a)~图4(e)中左侧图的是位移场(竖向位移用V表示，单位为mm)，最大值用Vmax表示；右侧的图是应变场(竖向应变用ε表示，单位为%)，最大值用εmax表示。由图4可以得出试样的整体变形与局部变形数值上相差很大，特别是在剪切后期，如在图4(e)中，试样的整体变形为7.98%时，局部的最大应变已经为21.04%。从变形场中也可以看出，在剪切前期，变形量较小，试样整体变形表现相对均匀。随着荷载的增大，试样中一些区域开始出现应变集中，渐渐地应变集中区扩展成一条剪切带贯穿整个试样，试样被分成了两个部分。剪切带形成后试样的变形主要集中在剪切带中，剪切带外的变形在达到4%左右时，基本不再发生变化，这个数值要比应力-应变曲线峰值对应轴向变形4.88%要小。基于分叉理论可知，在峰值前试样内部就出现了破坏的区域，即剪切带开始出现；剪切带则完全形成与峰值之后。
2.3 试样的宏观力学特性
在图像测量得到的应变场的基础上，平面应变试验中试样的体应变(εv)由计算区域面积的变化确定，即由式(1)计算得出：
假设同一平面上受到的应力是相等的，试样的整体应力由式(2)~式(5)计算确定。在试样上取了靠近试样中间的6个平面，如图3所示的分布在试样中间部分的6个截面。
式中：σi表示第i个截面的偏应力(i=1~6)；F表示上置轴向力传感器和下置轴向传感器的平均值；F1表示上置传感器的力；F2表示下置传感器的力；Ai表示第i个截面加载过程中的面积；A0是截面的初始面积；εxi表示第i个截面的水平应变；σm表示试样的整体应力。试样在不同围压下应力-应变曲线与体积变形曲线，如图5所示。
图4 试样在竖直方向上不同时刻的位移场与应变场
Fig.4 Vertical displacement contours and strain distribution diagrams through the course of test
图5 应力-应变曲线与体积变形曲线
Fig.5 Stress-strain curves and volumetric deformation
由位移传感器(LVDT)测量得到的应力-应变曲线与数字图像相关法(DIC)测量得到的应力-应变曲线对比可知，采用图像测量方法的峰值出现的要比传统的测量方法早(图5(a))，因此采用传统方法得出的结果有可能低估了试样内的局部变形。基于数字图像相关技术不仅可以得到试样宏观的应力-应变曲线而且也可以得到试样的其他参数，如摩擦角φ′，剪切带的倾角β与剪切带的宽度T，结果如表2中所示。其中：摩擦角φ′利用应力应变曲线数据画出莫尔圆根据摩尔-库伦准则可以得到；剪切带的倾角β可以通过应变场的数据直接得出；剪切带的宽度T获取方法如图6(a)所示。取与剪切角方向垂直的5条线段，根据数字图像相关法提取数据，可以得到这5条线段上的应变值。把5组应变加和平均，以应变曲线上斜率最大的点为分界点。左右两个分界点之间的区域为剪切带区，外侧的为剪切过渡区，如图6(b)所示。
图6 剪切带的计算宽度示意图
Fig.6 Principle of obtaining the shear band thickness
2.4 福建标准砂的宏观指标选取
根据2.3节的得到的结果可知，基于数字图像相关技术得到了福建标准砂的宏观应力-应变曲线，结果与常规方法得到的结果有所差异。由此在数字图像测量技术下测得的福建标准砂宏观力学反映可以得到一些宏观力学参数包括抗压强度、弹性模量、体积变形。对于体积变形，常规试验只能得到试样饱和状态下排水剪切试验的体积变形。而基于数字图像测量技术不但可以得到各种试验条件下试样的体积变形，还可得到砂土渐进破坏的过程，这是常规试验方法无法实现的。基于测量获得的全场变形结果，可直接得到福建标准砂试样在平面应变条件下的剪切带宽度与倾角。在此基础上，利用离散元技术对试验结果进行数值模拟，既需要符合福建标准砂的宏观应力应变关系、体积变形，又需要与试验中剪切带的形态保持一致，使得细观参数的获取更加可靠。
3 平面应变试验离散元分析
3.1 颗粒级配与试样尺寸
离散元技术可以预演岩土材料颗粒之间的相互作用。为了使模拟接近真实的岩土材料，首先选择对应的理论模型，然后依据试验结果对模型参数进行标定。在福建标准砂的平面应变试验中，试样沿着剪切带发生错动，在剪切带内部，颗粒一定伴随着滚动与滑动，因此数值模型采用抗滚动摩擦模型(Rolling Resistance Linear Model)[29]，基于这种离散元模型，需要标定的参数为颗粒的法向刚度kn，颗粒的切向刚度ks，颗粒的摩擦系数f及颗粒的滚动摩擦系数fr。数值模拟的试样尺寸为120 mm×60 mm的二维矩形，与实际试样的高度宽度一致。在实际试样中尺寸过小的颗粒对试样力学性能影响很小，为了提高计算效率，去掉了尺寸相对较小的颗粒。在本次离散元模拟中采用颗粒形状为圆形颗粒，根据颗分曲线，颗粒半径在0.158 mm~0.75 mm之间均匀分布，即模拟中颗粒的最大半径Rmax为0.75 mm，颗粒最小半径Rmin为0.158 mm。由于本试验采用的新型平面应变加载设备为试样提供的围压是柔性荷载，为了使数值模拟与试验的边界条件接近，本模拟中采用Wang等[30]与蒋明镜等[31]提出的基于接触黏结联结的微小颗粒来模拟试验中的侧向的柔性膜，如图7所示。
图7 离散元模拟模型
Fig.7 Simulation model of discrete element method
岩土材料的宏观力学特性其实是由颗粒之间的细观特性决定的，而颗粒的细观特性并不能直接测量得到，因此为了得到土颗粒的细观特性，常用的方法是由试验宏观力学行为标定出土颗粒的细观参数，然后用标定后细观参数代入数值模型模拟多种试验条件以验证微观参数的合理性。在本文中抗滚动摩擦模型需要确定的4个参数，采用试错法，经过多次的模拟并且与宏观试验结果相对比，确定出基于数值模拟的细观参数如表3所示。
表3 福建标准砂细观参数表
Table 3 Microscopic parameters of Fujian sand
注：表中kn表示颗粒的法向刚度；ks表示颗粒的切向刚度；f表示颗粒的摩擦系数；fr表示颗粒的滚动摩擦系数；Rmax表示颗粒的最大半径；Rmin表示颗粒最小半径；n表示颗粒的孔隙率。
kn/(N/m) ks/(N/m) f fr Rmax/mm Rmin/mmn 7×107 4.7×107 0.42 0.050.75 0.158 0.15
3.2 离散元(DEM)模拟结果分析
利用数字图像相关法(DIC)可以观测到试样的整体变形过程和局部变形过程，但是由于试验条件的限制，试样被包裹在试样膜中，观测试样的变形是通过试样膜的变形来表现的。当剪切带形成以后，试样被分为上下两个部分并沿着剪切带滑动，因而无法直接观测到颗粒的运动形态。借助PFC2D软件，基于试验数据对离散元参数标定后，可以观测到土颗粒的运动形态。
利用fish语言生成试样颗粒，对生成的模型固结后开始加载试验，上下加载板同时向试样中心缓慢加载。把表3中的细观参数放入模型中进行离散元的数值模拟。数值模拟中试样的应力、应变、体积变形利用9个测量圆进行测量，测量圆的范围与图3中变形场测量区域接近，整体应力与应变由所有的测量圆加权平均所得，模拟得到的各个围压下的应力-应变曲线与试验的应力-应变曲线对比如图5(b)所示，试验得到的试样体积变形与数值模拟得到的体积变形随轴向应变的曲线如图5(c)与图5(d)所示。100 kPa侧向压力下试样数值模拟的轴向位移如图8所示。
图8 离散元模拟结果
Fig.8 Simulation result of discrete element method
通过对3个围压条件下的平面应变试验与模拟结果对比表明，利用标定的细观参数模拟得到的应力-应变曲线和体积变形与试验结果吻合。试验中由于采用了数字测量技术，因此能给利用更多的信息来标定细观参数，除了对比应力-应变曲线和体积变形曲线也可以对比试验与模拟的变形场。图8(a)为数值模拟获取的100 kPa侧向压力下试样的轴向位移图，由于采用上下加载板同时向试样中心加载的方式，因此试样位移值有正值和负值。图8(b)是试验获取的100 kPa侧向压力下剪切带形成后欧拉坐标系下的轴向位移图，试验加载时自上而下加载，位移值都以竖向向下为正。从图8(a)与图8(b)中可以看出，两者的剪切带倾角基本相同，位移趋势较为一致，数值模拟的柔性边界变形与试验也较为相似。因此，这也表明所确定的细观参数是准确的。不可否认，由于制样精度和土料均匀性并非完全理想，模拟中得到的X型剪切带与试验中出现的是单条剪切带并不能够完全吻合，这也是下一步改进计算理论、提升制样精度的研究方向。
4 结论
研制了一种新型的平面应变加载设备，搭建了能够得到表面变形识别的数字图像采集系统。在此基础上，利用研制的平面应变设备结合二维数字图像相关技术(2D-DIC)，给出了基于图像测量的福建标准砂的宏细观参数确定方法，基于以上研究得到主要结论如下：
(1) 基于轴向传感器的位移测量方法只能得到试样的整体变形，不能描述岩土材料的非均匀变形且低估了岩土材料内部的局部变形量。数字图像测量技术可以量化加载过程中材料的非均匀变形。
(2) 加载过程中剪切带是由应变集中区域慢慢扩展形成的。剪切带产生后试样的变形主要集中在剪切区域，剪切带外侧区域变形不大。
(3) 通过将数字图像相关技术(DIC)应用于平面应变试验，可以得到试样宏观的应力-应变关系曲线。并且可获得试样的体积变形曲线及局部破坏形态。相较常规试验，此方法可以量化更多的宏观力学特性，以此获得的细观参数更可靠。
[1]Thomas T.Plastic flow and fracture in solids [J].Journal of Mathematics and Mechanics, 1958, 7: 291-322.
[2]Shao L T, Liu G, Zeng F T, et al.Recognition of the stress-strain curve based on the local deformation measurement of soil specimens in the triaxial test [J].Geotechnical Testing Journal, 2016, 39(4): 658-672.
[3]Zhang X, Li L, Chen G, et al.A photogrammetry-based method to measure total and local volume changes of unsaturated soils during triaxial testing [J].Acta Geotechnica, 2014, 10(1): 55-82.
[4]杨智勇, 曹子君, 李典庆, 等.颗粒接触摩擦系数空间变异性对颗粒流双轴数值试验的影响[J].工程力学,2017, 34(5): 235-246.Yang Zhiyong, Cao Zijun, Li Dianqing, et al.Effect of spatially variable friction coefficient of granular materials on its Macro-mechanical behaviors using biaxial compression numerical simulation [J].Engineering Mechanics, 2017, 34(5): 235-246.(in Chinese)
[5]姜浩, 徐明.碎石料应力路径大型三轴试验的离散元模拟研究[J].工程力学, 2014, 31(10): 151-157, 180.Jiang Hao, Xu Ming.Study of stress-pathdependent behavior of rockfills using discrete element method [J].Engineering Mechanics, 2014, 31(10): 151-157, 180.(in Chinese)
[6]Medina-Cetina Z, Rechenmacher A.Influence of boundary conditions, specimen geometry and material heterogeneity on model calibration from triaxial tests [J].International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34(6): 627-643.
[7]Scholey G K, Frost J D, Presti D L, et al.A review of instrumentation for measuring small strains during triaxial testing of soil specimens [J].Geotechnical Testing Journal, 1995, 18(2): 137-156.
[8]陈超斌, 叶冠林.基于LVDT的小应变三轴仪研制及其软土试验应用[J].岩土力学, 2018, 39(6): 2304-2310.Chen Chaobin, Ye Guanlin.Development of small-strain triaxial apparatus using LVDT sensors and its application to soft clay test [J].Rock and Soil Mechanics, 2018,39(6): 2304-2310.(in Chinese)
[9]Alikarami R, Andò E, Gkiousas-Kapnisis M, et al.Strain localisation and grain breakage in sand under shearing at high mean stress: insights from in situ X-ray tomography[J].Acta Geotechnica, 2014, 10(1): 15-30.
[10]姚志华, 陈正汉, 李加贵, 等.基于CT技术的原状黄土细观结构动态演化特征[J].农业工程学报, 2017,33(13): 134-142.Yao Zhihua, Chen Zhenghan, Li Jiangui, et al.Meso-structure dynamic evolution characteristic of undisturbed loess based on CT technology [J].Transactions of the Chinese Society of Agricultural Engineering, 2017, 33(13): 134-142.(in Chinese)
[11]Cheng Z, Wang J.Quantification of the strain field of sands based on X-ray micro-tomography: A comparison between a grid-based method and a mesh-based method [J].Powder Technology, 2019, 344:314-334.
[12]周健, 史旦达, 吴峰, 等.基于数字图像技术的砂土液化可视化动三轴试验研究[J].岩土工程学报, 2011,33(1): 81-87.Zhou Jian, Shi Danda, Wu Feng, et al.Visualized cyclic triaxial tests on sand liquefaction using digital imaging technique [J].Chinese Journal of Geotechnical Engineering, 2011, 33(1): 81-87.(in Chinese)
[13]Mousa A.Revisiting the Calibration Philosophy of Constitutive Models in Geomechanics [J].International Journal of Geomechanics, 2017, 17(8): 06017002.
[14]Sun L, Abolhasannejad V, Gao L, et al.Non-contact optical sensing of asphalt mixture deformation using 3D stereo vision [J].Measurement, 2016, 85: 100-117.
[15]Li L, Zhang X.Factors influencing the accuracy of the photogrammetry-based deformation measurement method[J].Acta Geotechnica, 2019, 14(2): 559-574.
[16]Sutton M, Wolters W, Peters W, et al.Determination of displacements using an improved digital correlation method [J].Image and Vision Computing, 1983, 1(3):133-139.
[17]Pan B, Dafang W, Yong X.Incremental calculation for large deformation measurement using reliability-guided digital image correlation [J].Optics and Lasers in Engineering, 2012, 50(4): 586-592.
[18]Keane R D, Adrian R J.Theory of cross-correlation analysis of PIV images [J].Applied scientific research,1992, 49(3): 191-215.
[19]Tang Y, Okubo S, Xu J, et al.Progressive failure behaviors and crack evolution of rocks under triaxial compression by 3D digital image correlation [J].Engineering Geology, 2019, 249(31): 172-185.
[20]王学滨, 杜亚志, 潘一山, 等.基于数字图像相关方法的等应变率下不同含水率砂样剪切带观测[J].岩土力学, 2015, 36(3): 625-632.Wang Xuebin, Du Yazhi, Pan Yishan, et al.Measurement of shear bands of sand specimens with different water contents under constant strain rate based on digital image correlation method [J].Rock and Soil Mechanics, 2015,36(3): 625-632.(in Chinese)
[21]Yang C, Wei J, Huang H, et al.Application of 3D-DIC to characterize the effect of aggregate size and volume on non-uniform shrinkage strain distribution in concrete [J].Cement and Concrete Composites, 2018, 86: 178-189.
[22]赵燕茹, 王磊, 韩霄峰, 等.冻融条件下玄武岩纤维混凝土断裂韧度研究[J].工程力学, 2017, 34(9): 92-101.Zhao Yanru, Wang Lei, Han Xiaofeng, et al.Fracture toughness of basalt-fiber reinforced concrete subjected to cyclic freezing and thawing [J].Engineering Mechanics,2017, 34(9): 92-101.(in Chinese)
[23]肖军华, 张德, 王延海, 等.基于DEM-FDM耦合的普通铁路碎石道床-土质基床界面接触应力分析 [J].工程力学, 2018, 35(9): 170-179.Xiao Junhua, Zhange De, Wang Yanhai, et al.Study on interface stress between ballast and subgrade for traditional railway based on coupled DEM-FDM [J].Engineering Mechanics, 2018, 35(9): 170-179.(in Chinese)
[24]胡光辉, 徐涛, 陈崇枫, 等.基于离散元法的脆性岩石细观蠕变失稳研究 [J].工程力学, 2018, 35(9): 26-36.Hu Guanghui, Xu Tao, Chen Chongfeng, et al.A microscopic study of creep and fracturing of brittle rocks based on discrete element method [J].Engineering Mechanics, 2018, 35(9): 26-36.(in Chinese)
[25]Zhang F, Li M, Ming P, et al.Three-dimensional DEM modeling of the stress–strain behavior for the gap-graded soils subjected to internal erosion [J].Acta Geotechnica,2019, 14(2): 487-503.
[26]唐欣薇, 黄文敏, 周元德, 等.华南风化花岗岩劈拉断裂行为的试验与细观模拟研究 [J].工程力学, 2017,34(6): 246-256.Tang Xinwei, Huang Wenmin, Zhou Yuande, et al.Experimental and meso-scale numerical modeling of splitting tensile behavior of weathered granites from south China [J].Engineering Mechanics, 2017, 34(6):246-256.(in Chinese)
[27]Pan B, Qian K, Xie H, et al.Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review [J].Measurement Science and Technology, 2009, 20(6): 062001.
[28]Pan B, Xie H, Wang Z.Equivalence of digital image correlation criteria for pattern matching [J].Applied Optics, 2010, 49(28): 5501-5509.
[29]Jiang M, Yu H, Harris D.A novel discrete model for granular material incorporating rolling resistance [J].Computers and Geotechnics, 2005, 32(5): 340-357.
[30]Wang Y, Leung S.A particulate-scale investigation of cemented sand behavior [J].Canadian Geotechnical Journal, 2008, 45(1): 29-44.
[31]蒋明镜, 张望城, 孙渝刚, 等.理想胶结砂土力学特性及剪切带形成的离散元分析[J].岩土工程学报, 2012,34(12): 2162-2169.Jiang Mingjing, Zhang Wangcheng, Sun Yugang, et al.Mechanical behavior and shear band formation in idealized cemented sands by DEM [J].Chinese Journal of Geotechnical Engineering, 2012, 34(12): 2162-2169.(in Chinese)