基于可视化试验的风积沙宏细观力学特性分析

薛嘉麟， 杜 强

(1.内蒙古工业大学 土木工程学院， 内蒙古 呼和浩特 010051；2.沙旱区地质灾害与岩土工程防御自治区高等学校重点实验室， 内蒙古 呼和浩特 010051)

内蒙古自治区存在规模较大的沙漠，如库布齐沙漠、毛乌素沙漠、腾格里沙漠等。沙漠地区矿藏资源丰富，却因为交通不便，影响其开发和利用。近年来，关于沙漠地区风积沙的研究是个热点问题。张俊[1]通过对风积沙颗粒分析和含水率等方面研究，验证了公路换填风积沙垫层的可行性。张宏等[2]对风积沙进行了矿物成分、化学成分、酸碱度、含盐量与级配特征等方面的研究。王丽英等[3]研究了风积沙改性土基坑填料力学性质和结构受盐-冻作用的影响规律。阮波等[4]通过对低温养护条件下水泥改良风积沙无侧限抗压强度试验进行研究，提出低温条件下，不同掺量水泥风积沙的刚度损失率。张飞阳[5]对低含水率下风积沙地基的抗拔承载特性做出研究，并得出风积沙地基的上拔承载力随着基础底板尺寸及含水率的增加而增加，风积沙地基加载前中期开始是受到张拉破坏，最后的破坏是由于张拉破坏和剪切破坏联合作用导致的。吴克雄等[6]对西北地区高性能风积沙固化剂的研制及工程应用做出研究。董瑞鑫等[7]对干湿循环与风沙吹蚀作用下风积沙混凝土的抗硫酸盐耐久性进行了研究。上述文献大多对风积沙的物理化学性质和风积沙土性改良等方面进行研究，对风积沙的力学特性研究较少。此外，风积沙级配特殊，如果从细观角度进行研究可能有助于进一步解释其力学机理。

为了探明风积沙的细观力学特性，研究受力过程中土颗粒的运动情况，本论文对现有的直剪仪、压缩仪进行可视化改进，并结合数字照相变形量测技术，对各密实度下直剪试验和压缩试验中试样的变形进行分析，得出土样位移场和应变场云图。结合离散元数值模拟程序，建立风积沙的细观数值模型，对风积沙的室内直剪试验、压缩试验和三轴试验进行数值模拟，研究风积沙的细观力学特性。

1 风积沙的颗粒分析试验

试验所用风积沙取自库布齐沙漠，风积沙呈浅红棕色，颗粒粒径均匀，内部掺杂黑色或白色细小颗粒，见图1。

图1 风积沙样本Fig.1 Aeoliansand sample

对风积沙进行颗粒分析试验，得到颗粒级配曲线，并与标准砂进行对比分析，见图2。由风积沙级配曲线可知风积沙粒径分布集中在0.075～0.25 mm，以颗粒较大的砂粒为主，其不均匀系数Cu=2.30，曲率系数Cc=0.709，级配不良。

图2 级配曲线Fig.2 Classification curve

2 可视化试验

为了更好地观察剪切和压缩试验过程中土颗粒的位移及试样的变形情况，对直剪仪和压缩仪进行了可视化改进，改进后的试样盒见图3。利用对比试验的方法，验证新型可视化试样盒的可靠性，并为后续分析提供图像和基本力学参数。

图3 改进后的试验器材Fig.3 Improved test equipment

2.1 直剪试验

利用改进后的直剪仪对风积沙进行直剪试验，并与常规直剪试验进行对比，试验曲线见图4。由试验可得：改进后的可视化直剪盒所得结果与常规直剪试验结果基本一致，风积沙内摩擦角大约为26°，黏聚力在1～3 kPa范围内，风积沙的抗剪强度主要靠其内摩擦力构成，黏聚力可以忽略不计。

图4 不同垂直压力下的剪应力-位移曲线Fig.4 Stress-strain curve under different vertical pressures

将库布齐沙漠风积沙的抗剪强度指标与文献[8]中级配较好的砂土进行对比，见表1。可知砂土的内摩擦角不仅受不均匀系数影响，还与曲率系数有关，风积沙的不均匀系数较低，曲率系数小于1,与级配良好的砂土相比，内摩擦角小于级配良好的砂土。但随着土体密实度的增加，其抗剪强度有所提高。

表1 风积沙与级配良好砂土的对比

2.2 压缩试验

利用改进后的压缩仪对风积沙进行压缩试验，并与常规压缩试验进行对比，试验曲线见图5。由试验结果可知，改进后的可视化压缩盒所得结果与常规压缩试验结果基本一致。本次试验所使用的风积沙初始孔隙比在0.615左右， 压缩系数为0.25 MPa-1。

图5 改进前后压缩试验对比Fig.5 Comparison of improved front and rear compression tests

对比文献[9]中的结果可知，颗粒粗细越悬殊，试样压缩性越高,而细粒含量越多，压缩性越低。粒径的增大、粒径变化范围的拓宽、颗粒形态的不规则均会导致压缩性的升高。库布齐风积沙的颗粒粒径均匀，颗粒较细，因此压缩性较低。

通过上述试验结果可知，改进后的可视化直剪仪与压缩仪，所得结果与常规土工试验仪器所得结果基本一致，证明了改进后仪器的可行性。

3 可视化数字图像分析研究

在可视化直剪与压缩试验过程中，采用高清数码相机对试验过程进行拍摄，并利用PhotoInfor软件对所采集照片进行数字图像分析。得到各级垂直荷载作用下风积沙剪切和压缩过程中研究区域内土体的位移和应变变化过程。图6中红色方框为PhotoInfor分析区域。

图6 可视化分析区域Fig.6 Visual analysis area

3.1 直剪试验

由于在施加不同法向应力的情况下得到的剪切变形规律大致相同，此处以法向压力为 100 kPa和400 kPa的试验结果为例进行分析，见图7。

图7 直剪试验分析区域位移场云图Fig.7 Cloud map of regional displacement field analyzed by direct shear test

通过对最大剪位移云图分析发现，风积沙的剪切变形不均匀，表现为剪切接触处位移较大，上部和下部位移较小，且越远离接触面位移越小。在剪切接触面两边位移较小，中央位移较大。说明直剪破坏由剪切面中心开始，沿剪切面延申，最终形成一条贯穿的剪切带。试验还表明在不同法向荷载下试样剪切变形的方向变化和剪切带的形成过程大致相同，这进一步验证了可视化仪器结果的可靠性。

3.2 压缩试验

图8为不同压力下压缩试验中的位移场云图。由图8可知分析区域内压缩变形由上透水石中间附近位置开始。

图8 不同压力下压缩试验分析区域位移场云图Fig.8 Regional displacement field cloud map of the compression test analysis at different pressures

中间蓝色区域位移量为0是由于上下盒之间的亚克力板反光造成，为不可避免误差，在任意荷载等级下此位置位移均为0，因此可忽略此影响。由于风积沙颗粒均匀细小且无粘聚力，故土体颗粒均匀压缩，无论在何种荷载作用下，各个水平面的土体颗粒位移都大致相同。在任意一级荷载作用下，都是试样上部的位移变化最大，且中间位移相比边缘较大。随着荷载等级的增大，风积沙颗粒间的变形带逐渐形成，试样上部位移场云图出现明显的弧形分布。且随着荷载的增加，位移变化也逐渐增大且变化范围逐渐向周围进行扩伸，整体位移呈规律性的增大。到了加荷后期，可以看出变形范围已经基本保持稳定，位移最大增量仍主要集中在试样整体上部的中间区域，向外拓展延伸的趋势减小，此时风积沙基本已压缩稳定。

4 细观数值模拟研究

为从细观方面进一步分析风积沙的力学特性，基于可视化力学试验结果，采用数值模拟软件Particle Follow Code，简称PFC，对直剪试验、压缩试验和三轴试验进行数值模拟。

由颗分试验结果可知，风积沙实际粒径范围为0.075～0.25 mm，故数值模拟试验中拟生成风积沙颗粒粒径同样在0.08～0.25 mm范围内，赋予风积沙颗粒为黄色，赋予墙体为黑色，赋予坐标轴为棕色。给墙体赋予108N/m的强度和刚度，设定风积沙的重度为2 650 N/m3，重力加速度为9.8 m/s2，通过不断运行模型得出相应试验结果，使之与实际可视化直接剪切试验的结果进行对比，经过多次试算校核，得出较为准确的风积沙细观参数，最后确定风积沙强度、刚度均为2×107N/m，摩擦系数为1.2。

4.1 直剪试验

建立七面墙作为约束的试验范围，模拟上下剪切盒生成两个封闭的区域。以应变控制式直剪仪的实际尺寸确定墙体长度，即整个剪切盒长度为6.18 cm，高度为4 cm，模拟二维平面内的直剪试验。上下两个剪切盒内各生成3 600个风积沙颗粒。直剪试验数值模拟模型如见图9。

图9 直剪试验数值模拟模型图Fig. 9 Numerical simulation model diagram of direct shear test

由于PFC2D中无法通过命令流直接给墙体施加荷载，所以在模拟过程中，首先对顶部墙体施加竖直向下1×10-7m/s的速度，使风积沙试样发生竖向压缩，从而土颗粒对墙体产生反力，经公式转化后的法向应力分别达到100 kPa、200 kPa、300 kPa、400 kPa时停止压缩，与之相对应的压缩量依次为1.68 mm、2.32 mm、2.93 mm、3.50 mm。

控制剪切位移最大为6 mm，记录垂直压力在100 kPa、200 kPa、300 kPa、400 kPa下风积沙颗粒间的接触力变化情况以及剪切应力与剪切位移的数值关系曲线，共进行四次模拟试验，得到如图10所示的试验数据。

图10 不同垂直压力下竖向反力和切应力变化曲线Fig.10 Change curves of vertical reverse force and shear stress under different vertical pressure

根据以上数据可以得出风积沙在进行数值模拟时，当垂直压力为100 kPa、200 kPa、300 kPa、400 kPa时对应的抗剪强度分别为48 kPa、105 kPa、142 kPa、196 kPa，由此便可以绘制风积沙PFC2D数值模拟的剪切包线，见图11。

图11 室内可视化直接剪切试验与PFC2D数值模拟抗剪强度包线对比Fig. 11 Indoor visualization direct shear testbeing compared with PFC2D numerical simulated shear strength envelope

由抗剪强度的库仑定律可以得出风积沙的黏聚力为2.5 kPa，内摩擦角为25°。由PFC2D数值模拟得出的试验结果与室内可视化直接剪切试验的试验结果很接近。与室内直接剪切试验相比，在数值模拟的过程中，为避免颗粒散落，在剪切面上增设了水平方向的两个墙体，可能是导致数值模拟试验结果中内摩擦角偏小、黏聚力偏大的主要原因。

图12是直剪试验颗粒间的接触力链。由图12可知，剪切试验刚开始进行时，土样中颗粒分布分散，只有部分颗粒产生一定的位移。继续施加水平剪切力，上下盒发生错动开始产生剪切带，剪切带位于上下盒交界处。

图12 直剪试验的接触力链Fig.12 Contact force chain for the direct shear test

另外，剪切带随着接触面上的风积沙不断的平移、错动，最终稳定、贯通。在剪切过程中剪切带以上土颗粒向上竖直运动。推测土颗粒的移动是由于在剪切过程中，土颗粒之间相互挤压造成，表现为剪胀性。直剪试验完成后可以观察到接触力链呈对角分布，角度大致为45°。颗粒细观运动结果与可视化直剪试验得到的宏观力学现象相一致。

4.2 压缩试验

建立四面墙作为约束的试验范围，以压缩仪的实际尺寸确定墙体长度，即长度为6.18 cm，高度为2 cm，模拟二维平面内的压缩试验。由于不能直接对边界墙施加荷载，故对顶面墙体施加y轴负方向1×10-7m/s的速度，将速度等效替代垂直压力，便可以模拟压缩试验，并监控测量墙体随试验进行在竖直方向的反力，即压缩试验的荷载，控制最大压缩量为4 mm。在试样范围内设置三个测量圆，测定压缩过程中土体孔隙率的变化。得出压缩荷载、孔隙率随压缩进程的变化曲线，经过公式转换，达到孔隙比与荷载之间的关系曲线，即e-p曲线，得出风积沙的压缩特性。压缩试验数值模拟模型见图13。

图13 压缩试验颗粒模型图Fig.13 Model of compression test

对三个代表颗粒Ball 1(x=0.01，y=-0.01)、Ball 2(x=0.02，y=-0.01)、Ball 3(x=0.03，y=-0.01)的水平速度、竖直速度进行监控。限于篇幅，此处只列出Ball 1图，见图14。

图14 Ball 1水平方向和竖直方向速度曲线Fig.14 Velocity curves of in horizontal and vertical direction of Ball 1

在模拟压缩过程中发现颗粒体水平方向的速度均比较小，少数颗粒会突变，而竖直方向速度基本等同于压缩速度。见图15。

图15 压缩后颗粒接触力链Fig.15 Grain contact force chain after compression

同一高度上的颗粒体水平方向的位置变化很小，而高度越靠上的颗粒体竖直方向位置下降情况越明显。颗粒细观运动结果反映了宏观变化的内在机理。压缩结束后模型应力链呈倾斜的散乱分布，倾角与水平面小于45°，中部较为密集，两边稍疏。

由于PFC2D是二维平面范围内的数值模拟，其生成的所有颗粒形心全部处于同一平面内。相反，PFC3D模型生成的颗粒情况并不是这样，而是更接近于室内可视化试验实测的孔隙率，因此PFC3D模型中生成的颗粒间孔隙率要远大于PFC2D模型生成的颗粒间孔隙率。研究表明，用相同的球体填充PFC墙体范围时，在球与球之间接触十分紧密的情况下，PFC3D生成的颗粒间孔隙率是0.259 5，PFC2D模型生成的颗粒间孔隙率是0.093 1。即室内可视化试验的实测孔隙率是PFC2D模型生成的颗粒间孔隙率的2.787倍，本论文中采用该系数进行孔隙率转换，再进行两者的对比分析，e-p曲线见图16。

图16 现场试验与模拟试验e-p曲线对比Fig.16 Comparison of e-p curve of field test and simulation test

从整体上讲，数值模拟得出的试验结果与室内可视化压缩试验的试验结果比较接近。室内试验的e-p曲线变化幅度较小，而数值模拟试验的e-p曲线变化幅度较大，造成这种现象的主要原因是风积沙的初始孔隙率不同以及压缩速率不同。

4.3 三轴试验

为更好地了解风积沙在剪切过程中的细观变化，使用PFC3D进行了三轴试验的离散元数值模拟研究。三轴试样的高度和横截面的直径的长度比为2∶1，因此模拟三轴试样的高度也应该是横截面直径的二倍。又因为要考虑进行模拟验算时计算机所能承担的最大计算量，最终假定的三轴试样尺寸为直径15 mm，高30 mm，见图17。

图17 三轴试验数值模拟建模图Fig.17 Numerical simulation and modeling diagram of the triaxial test

PFC3D在运行结束之后会获得有关土颗粒的多组图像，如颗粒的速度矢量图、位移矢量图、接触力分布图等。这些图像对试验结果处理有着很重要的作用。图18～20为围压100 kPa下各应力的变化曲线。

图18 土体应变变化Fig.18 change of soil strain

图19 围压变化Fig.19 Surrounding pressure changes

本试验以围压为50 kPa，100 kPa，200 kPa下的轴向应力应变情况的观测为主，发现其对应的最大剪应力分别为142 kPa，252 kPa，366 kPa，且各围压下风积沙轴向应力应变情况趋势相同，只是大小有所不同。

图20 轴向应力变化Fig.20 Changes in axial stress

5 结 论

1) 在误差允许范围内，经过改进的方形可视仪器具有可行性，且在试验观测上具有其独特的优势。

2) 由宏细观试验可知，库布齐沙漠的风积沙存在着结构松散、级配不良、孔隙率大、透水性强、保水性较差、水稳性好、粘聚力小甚至无粘聚力、抗剪强度低的特点。

3) 在剪切过程中，风积沙颗粒初期分布散乱，大粒径颗粒先发生运动，逐步形成稳定、贯通的剪切带，且具有剪胀性。

4) 压缩过程中，颗粒体水平方向的速度均比较小，少数颗粒会突变，而竖直方向速度基本等同于压缩速度。同一高度上的颗粒水平方向的位置变化很小，土体的压缩主要发生在土样的上部。

由风积沙细观力学试验可知，风积沙的抗剪强度会随着土体密实度的提高而提高。风积沙的压缩性较低，且压缩变形主要集中在土体上部。因此在工程实践中，可以通过压实浅部土体的方法，来提高强度，减小变形，满足工程对强度和变形的要求。