基于数字图像相关技术的砂土全场变形测量及其DEM数值模拟

王鹏鹏 ，郭晓霞 ，桑 勇 ，邵龙潭 ，陈之祥 ，赵博雅

(1.大连理工大学工业装备结构分析国家重点实验室，辽宁，大连 116024；2.大连理工大学工程力学系，辽宁，大连 116024；3.大连理工大学精密与特种加工教育部重点实验室，辽宁，大连 116024；4.大连理工大学机械工程学院，辽宁，大连 116024)

土的局部变形如剪切带的形成规律和渐进破坏特性是土体区别于常规连续介质的重要性质[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

表2 试验条件及试验测量结果Table 2 Experimental conditions and measurement

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

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)应用于平面应变试验，可以得到试样宏观的应力-应变关系曲线。并且可获得试样的体积变形曲线及局部破坏形态。相较常规试验，此方法可以量化更多的宏观力学特性，以此获得的细观参数更可靠。