压-剪复合应力下非球形颗粒材料空心圆柱剪切试验的离散元模拟

史旦达，王 威，薛剑峰，邵 伟

（1.上海海事大学 海洋科学与工程学院，上海 201306；2.新南威尔士大学 工程与信息学院，堪培拉 澳大利亚 2612；3.四川大学 水力学与山区河流开发保护国家重点实验室，四川 成都 610065）

1 研究背景

在实际岩土工程问题中，建筑物地基或者基坑土体往往处于压-剪复合应力状态，即其土单元体会受到正应力和剪应力的耦合作用［1］。空心圆柱扭剪仪（Hollow Cylinder Apparatus，HCA）是室内研究复杂应力条件下土体静、动力特性的重要试验工具，针对砂土材料，Yang等［2］利用HCA研究了Toyoura砂在应力主轴旋转条件下的不排水剪切强度，重点讨论了中主应力和组构各向异性对土体强度规律的影响。Cai等［3］采用HCA探讨了压-剪复合应力条件下砂土的非共轴变形特性，其室内试验重点研究了加载路径对非共轴特性的影响。Chen等［4］在空心圆柱试验研究基础上，分析了多种不同主应力轴旋转路径下饱和粉土的动强度及其不排水各向异性特性。然而，空心圆柱室内试验一般只能得到土体宏观的力学响应，很难分析土体内部细微观结构的变化，因此难以建立其宏-细（微）观之间的内在关联［5］。离散单元法（Discrete Element Method，DEM）是目前开展宏细（微）观土力学研究的重要手段，近年来，基于DEM的空心圆柱数值模拟方法得到了发展，例如，Li等［6-8］针对散粒材料，采用DEM方法开展了空心圆柱剪切模拟，通过分析数值试样内部孔隙率和剪应变率的分布探讨了剪切带的演化规律，并与玻璃微珠HCA试验结果进行了对比。Farhang等［9］采用离散元TRUBAL程序建立了空心圆柱试样，通过DEM模拟分析了主应力方向变化对试样抗剪强度的影响。史旦达等［10］基于PFC3D（Particle Flow Code in 3-Dimension）离散元模拟平台构建了空心圆柱型数值模型，通过引入“分离式簇墙”技术近似模拟了HCA室内试验橡皮膜的柔性边界效应，并通过动态更新顶部扭矩层颗粒的旋转速度较好的控制了剪应力的施加，提升了空心圆柱DEM模拟的效果。然而，目前针对砂土等散粒材料的空心圆柱DEM模拟均仅局限于纯球形颗粒单元，没有考虑颗粒形状对试样宏细观力学响应的影响。

鉴于以上分析，本文在前期研究基础上，引入PFC3D“团聚颗粒”方法（CLUPM方法）构建三维椭球形颗粒单元。通过对空心圆柱数值试样同时施加竖向正应力和环向剪应力实现压-剪复合加载条件，并实施大主应力方向角α=30°恒定不变的单调剪切模拟。通过DEM模拟，重点分析颗粒形状和围压对数值试样应力-应变关系、抗剪强度、剪切带形态和组构各向异性演化规律的影响，并与已有的实际砂土室内试验结果展开对比。

2 数值试样建模

2.1 空心圆柱试样的受力状态压-剪单调剪切条件下，空心圆柱试样的典型受力状态如图1（a）所示。图1（a）中，W为竖向轴力，T为环向扭矩，pi、po为内、外围压。单元体的应力状态如图1（b）所示，σz、σr、σθ分别为竖向、径向和环向正应力，τzθ为剪应力。由各应力分量可进一步得到各主应力分布，如图1（c）所示，σ1、σ2、σ3分别为大、中、小主应力，α为大主应力方向角（即大主应力方向与竖直方向的夹角）。

图1 空心圆柱试样的应力状态

由pi、po、W、T并结合加载过程中测得的试样内径d和外径D，可计算得到各应力分量值，计算公式为：

由式（1）—（4），可进一步计算得到各主应力大小及大主应力方向角：

试样峰值内摩擦角φp的计算公式为：

式中（σ1/σ3）p为峰值应力比。

类似主应力的计算方法，可以得到大主应变ε1和小主应变ε3，进一步定义偏应变ε13为：

2.2 试样制备及参数设置数值试样为空心圆柱型，内径、外径、高度分别为60、100和200 mm。内、外边界均采用20面簇墙形式来近似模拟室内试验橡皮膜的柔性边界效应，每面簇墙的高度均为10 mm，上、下边界均为刚性墙体。颗粒单元为椭球形，在PFC3D中采用“团聚颗粒”方法（即“CLUMP”方法）构建，从纯球形颗粒至椭球形颗粒之间的转换满足“体积等效”和“质量等效”原则［11］，采用椭球形颗粒的长轴和短轴之比Se来反映颗粒形状效应，文中Se取1.0、1.2、1.5三种情况，Se=1.0即退化为纯球形颗粒，如图2所示。颗粒粒径范围0.53～3.3 mm，级配服从均匀分布，平均粒径d50为1.92 mm。考虑到计算效率，与实际砂土相比，DEM颗粒粒径做了适当放大，关于颗粒尺寸效应的讨论详见文献［10］。

颗粒之间的接触模型选择接触刚度模型，法向和切向的接触刚度均为5×104N·m-1。根据Procter和Barton［12］的试验研究，石英砂颗粒之间的摩擦系数约为0.49，因此本文DEM模型中颗粒之间的摩擦系数设为0.5。此外，考虑到实际砂土颗粒形状各异且颗粒之间存在较为显著的咬合作用，在DEM建模中引入线性旋转阻抗模型（linear rolling resistance model）从而一定程度上弥补颗粒过度旋转问题，关于旋转阻抗模型的详细介绍可以参考相关文献［13-14］，本文中旋转阻抗系数取0.1。PFC3D中提供了二种常用的阻尼机制，分别为局部阻尼和黏性阻尼机制，粘性阻尼适用于强夯、爆炸等高速率冲击问题的模拟，而局部阻尼主要用于静力或拟静力问题的分析。本文模型属于静力加载问题，因此采用局部阻尼，阻尼系数取0.7。

图2 颗粒形状

图3 数值试样（Se=1.5）

数值试样采用重力沉积制样，不同颗粒形状试样的初始孔隙率n0均控制为0.388，生成的数值试样效果如图3所示。重力沉积完成后，Se=1.0、1.2、1.5三种试样的初始平均接触数分别为5.8、6.8、7.1，关于颗粒集合体平均接触数的定义和计算方法详见文献［15］。试样生成后，通过内、外簇墙和上、下墙体施加各向均等围压σ3c，文中，σ3c取100、200和400 kPa三种情况。围压施加后，开始对试样实施大主应力方向角α=30°恒定不变的定向单调剪切模拟，竖向应力增量通过上、下墙体相对运动的方式予以施加，控制上、下墙体的运动速度均为ν=1 mm/s。剪应力通过扭矩层颗粒（如图4（a）所示）绕中心轴的旋转予以施加，通过目标剪应力值与实测剪应力值的偏差来动态调整扭矩层颗粒的旋转速度，扭矩层厚度取为颗粒平均粒径d50的2倍，即3.84 mm。当偏应变ε13达到10%，单调剪切终止。三种颗粒形状下，大主应力方向角α的实际控制效果如图4（b）所示。由图4（b）可见，对于球形和非球形颗粒的模拟，本文采用的数值建模技术均能较好的控制定向剪切的模拟效果。

图4 扭矩层设置及大主应力方向角控制效果

3 宏观剪切强度规律分析

3.1 各应力分量-偏应变关系图5给出了围压σ3c=100、200和400 kPa条件下数值试样的各应力分量（σz、σr、σq、τzθ）-偏应变ε13关系曲线。由图5可知，在定向剪切过程中，环向正应力σq和径向正应力σr均维持初始围压不变，而竖向正应力σz和剪应力τzθ逐渐增加，且随着偏应变的增加，两者基本保持相同的变化趋势。围压大小对σz和τzθ增长曲线的形态有一定影响，低围压下（100 kPa和200 kPa），σz和τzθ曲线呈现出较为明显的应变软化特征，在本文采用的围压水平下，σz和τzθ应力峰值出现在4.4%～7.8%偏应变范围内。高围压下（400 kPa），σz和τzθ曲线逐渐趋于应变硬化特征，σz和τzθ曲线没有明显的应力峰值点，但在加载后期，σz和τzθ的增长幅度明显变缓。颗粒形状对σz和τzθ曲线走势的影响不大，但颗粒形状会影响σz和τzθ的峰值大小；以σ3c=100 kPa为例，当Se从1.0增大到1.5时，σz峰值从168 kPa增加到180 kPa（增幅为7%），τzθ峰值从58 kPa增加到77 kPa（增幅为33%）。

图5 各应力分量随偏应变变化曲线

3.2 偏应力-偏应变关系图6给出了颗粒形状Se=1.0、1.2、1.5下数值试样的偏应力-偏应变（q-ε13）关系曲线。从图6可知，围压大小对偏应力-偏应变曲线的影响较为显著；随着围压的增加，偏应力-偏应变曲线的初始斜率也即试样的初始变形模量显著增加。Janbu［16］最早提出了土体初始变形模量E0与围压σ3c之间的定量关系式：

式中：K、n分别为初始变形模量系数和初始变形模量指数；pr为参考压力，可取100 kPa。

图6 偏应力-偏应变曲线

表1给出了各工况下数值试样初始变形模量E0的数值大小。利用Janbu公式对E0进行进一步分析，将数值试样的E0/pr～σ3c/pr关系绘于双对数坐标上并进行线性拟合，可得初始变形模量系数K和指数n的数值。Se=1.0、1.2、1.5试样的初始变形模量指数n分别为0.51、0.45、0.43。根据文献［1］的试验研究，福建标准砂的初始变形模量指数n在0.4～0.5之间，因此，本文数值试样的n值符合福建标准砂n值范围。需要特别说明的是，颗粒形状对剪切强度的影响较为显著，但对n值的影响并不十分显著，即使采用纯球形颗粒模拟其n值也基本位于实际砂土n值范围内，文献［9］的数值模拟也得到了与本文相同的规律。

表1 数值试样力学特性参数汇总

3.3 抗剪强度变化规律图7为围压σ3c=100、200和400 kPa条件下应力比-偏应变（σ1/σ3～ε13）关系曲线，大、小主应力的数值可由前文公式（5）计算得到。分析图7可知，与图5各应力分量的变化规律一致，在低围压下（100 kPa和200 kPa），三种颗粒形状下的应力比曲线均呈现出较为明显的应变软化特征；而在高围压条件下（400 kPa），相应的应力比曲线转变为应变硬化型。在图7中提取各工况下的峰值应力比（σ3c=400kPa工况下取10%偏应变水平所对应的应力比为峰值应力比）代入前文公式（8），可以计算得到数值试样的峰值内摩擦角φp，结果汇总于表1。分析表1可知，无论纯球形（Se=1.0）还是非球形颗粒模拟，随着围压的增加，数值试样的峰值内摩擦角均逐渐减小，这一规律符合实际砂土排水剪切时抗剪强度随围压变化的一般规律［17-18］。颗粒形状显著影响试样的抗剪强度，以σ3c=100 kPa为例，当椭球颗粒长短轴比Se从1.0增加到1.5时，数值试样的φp可增加7.9°。文献［11］曾对颗粒形状影响试样抗剪强度的细观机理进行过研究，其结果表明，试样抗剪强度的高低与颗粒形状差异引起的试样初始平均接触数多少密切相关，初始平均接触数越多，试样后继的抗剪强度就越高。由2.2节可知，在初始孔隙率n0相同的条件下，Se=1.5试样的初始平均接触数要明显高于Se=1.0试样，因此其表现出的宏观抗剪强度也越高。

图7 应力比—偏应变曲线

4 细观力学响应分析

4.1 剪切带形态及演化规律剪切带是颗粒材料试样应变局部化的重要表征，也是构建土体渐进破坏理论的重要基础［5］。本节采用颗粒位移场和对应的二维孔隙率云图来分析数值试样剪切带的形态及其演化规律。试样的内部孔隙率通过在数值试样内部布置测量球测得，从上至下共布置10层测量球，每层12个，量测球直径为20 mm，如图8所示。

图8 测量球布置

限于篇幅，以Se=1.5、σ3c=100 kPa工况为例来分析剪切带的演化过程，图9给出不同偏应变水平下数值试样的颗粒位移场矢量图及其对应的内部孔隙率云图，内部孔隙率云图可通过将轴对称空心圆柱试样沿平面展开后得到，具体展开方式详见文献［6］。分析图9可知，在荷载施加前（ε13=0），由于各向均等围压的作用，颗粒的位移矢量均指向试样内部，且围压施加后试样的局部孔隙率要略小于初始孔隙率（n0=0.388）。当试样达到峰值强度时（ε13=5%），从位移矢量图中可以看出试样已出现清晰可辨的剪切带区域，且从内部孔隙率云图中可以发现孔隙率数值较大的局域集中出现在试样斜向部位，说明剪切带已开始形成。当偏应变ε13达到为7.5%时（图9（c）），从位移矢量图中可以看出剪切带发展更为明显且形态更为稳定，从孔隙率云图中可以发现斜向剪切带已经贯通，试样带内的局部孔隙率要明显高于带外。当加载至ε13=10%（图9（d）），即达到残余状态时，试样的剪切带形态已非常稳定，剪切带与水平方向的夹角（下文中也称剪切带倾角θc）约为24.4°；由孔隙率分布云图可知，剪切带内最大局部孔隙率可达0.47以上，明显大于试样的初始孔隙率，反映了带内颗粒的重排列和剪胀效应。文献［19］曾借助X射线显微CT技术分析了Leighton Buzzard干砂（粒径为0.6～1.18 mm）在三轴剪切时剪切带的微观演化过程，试验研究表明，砂土试样的剪切带在应力峰值阶段开始形成，且达到残余状态时（轴向应变约15%），剪切带已经非常明显。因此，本文DEM数值试样呈现出的剪切带演化规律与上述实际砂土的规律保持一致。

图9 数值试样的剪切带演变规律（Se=1.5、σ3c=100kPa）

下文讨论颗粒形状和围压对剪切带倾角的影响，各工况下数值试样的剪切带倾角汇总于表1，表1中小括号里的数值为Mohr-Coulomb理论值。由图10可知，在主应力轴偏转条件下，由Mohr-Coulomb理论得到的剪切带倾角θc的理论值为：

图10 基于Mohr-Coulomb理论的剪切带倾角判定

分析表1中的数值，可得两点规律：（1）颗粒形状影响数值试样的剪切带倾角，椭球颗粒长短轴比Se越大，剪切带倾角越大，且这一规律与所施加的围压水平关系不大；（2）数值试样实测的剪切带倾角均小于Mohr-Coulomb理论值，其原因为Mohr-Coulomb理论值为剪切带倾角的上限值，没有考虑试样的初始各向异性。文献［20］曾以玻璃微珠为材料，进行了大主应力方向α=30°条件下的空心圆柱定向剪切试验，实测的剪切带倾角为22°，而本文数值试样的剪切带倾角范围为16°～25°，与文献［20］室内试验所得的数值非常接近。

4.2 组构各向异性演化分析在描述砂土细观组构的各参量中，颗粒之间的接触法向（contact normal）是影响其宏观强度和变形特性的主要细观参数［11］。Rothenburg和Bathurst［21］提出可采用二阶傅里叶函数来定量表述接触法向的各向异性演化规律，其公式为：

式中：a为各向异性系数，其数值大小反映各向异性程度；θa为接触法向各向异性主方向角。

对接触法向进行统计时，首先将三维接触矢量向竖向平面进行投影，再在平面内按10°一个区间进行分布密度的统计并绘制成玫瑰图，最后采用式（12）对统计结果进行拟合，得到各向异性分布参数a和θa的数值。其中，θa为接触法向各向异性主方向与水平轴正向的夹角。

图11给出了围压200 kPa条件下，三种颗粒形状数值试样在剪切加载过程中接触法向演化玫瑰图，图11中虚线表示式（12）的傅里叶拟合结果。对于每种颗粒形状的试样，应变水平均取加荷前、峰值强度前、峰值强度点、残余强度点四种状态进行分析。由图11可知，在加载开始前，纯球形颗粒试样（Se=1.0）的接触法向各向异性程度较小（a=0.074），且其各向异性主方向偏于水平方向（θa=182.8°）；而随着颗粒长短轴比Se的增加（图11（b1）和图11（c1）），数值试样接触法向分布的初始各向异性程度明显增强，且各向异性主方向开始偏于竖直方向。随着荷载的施加，无论采用纯球形颗粒还是非球形颗粒，接触法向各向异性的主方向均发生了明显偏转，且均逐渐靠近加载时的大主应力方向（α=30°定向剪切时，大主应力方向与水平轴正向的夹角为120°），体现了数值试样细观组构的应力诱发各向异性演化规律。

图11 接触法向各向异性演变规律（σ3c=200kPa）

图12 接触法向各向异性主方向角随偏应变变化规律

进一步将所有模拟工况下，接触法向各向异性主方向角θa随偏应变ε13的变化规律绘制于图12。分析图12（a）可知，当采用纯球形颗粒（Se=1.0）模拟时，围压大小对θa的影响较为显著，与高围压（400 kPa）情况相比，在低围压（100 kPa和200 kPa）下，当加载至残余状态时，数值试样的接触法向各向异性主方向越容易接近大主应力方向，体现出更强的应力诱发各向异性。出现这一现象的原因主要为低围压条件下数值试样的强度曲线呈应变软化型，数值试样呈现出明显的剪胀特性，剪切过程中颗粒之间出现较为剧烈的接触法向的重构，因此各向异性主方向更易发生大的偏转，从而更易接近大主应力方向；而在高围压条件下，数值试样的强度曲线呈应变硬化型，数值试样以剪缩为主，剪切过程中颗粒之间未发生剧烈的接触法向重组，因此θa的变化程度不如低围压条件下显著。进一步分析图12（b）和图12（c）可知，当引入非球形颗粒模拟时，围压对θa的影响开始变小，但与纯球形颗粒情况类似，在低围压下数值试样的接触法向各向异性主方向越靠近大主应力方向。

5 结论

本文针对空心圆柱试样的特点，基于PFC3D离散元平台开展了固定大主应力方向角α=30°不变的单调剪切模拟，研究了压-剪复合应力条件下数值试样的宏细观力学响应，重点讨论了颗粒形状和围压的影响。主要结论有：（1）数值试样的初始变形模量能够反映实际砂土的压硬性规律，初始变形模量随着围压的增大而增大，且初始变形模量指数其数值与福建标准砂试验结果较为接近。（2）数值试样的偏应力-偏应变曲线在低围压下表现为应变软化型，而在高围压下逐渐趋于应变硬化；其峰值内摩擦角随着颗粒长短轴比Se的增大而增大，随着围压的增大而减小；在细观机理上，颗粒形状对峰值强度的影响与试样的初始平均接触数密切相关。（3）数值试样的剪切带形成于应力峰值阶段，达到残余状态时其形态趋于稳定；在本文模拟工况下，其剪切带倾角在16°～25°范围内，与已有文献中玻璃微珠空心圆柱试验的实测结果吻合较好，且随着Se的增大，剪切带倾角越大。（4）无论接触法向初始各向异性的主方向分布如何，随着剪切荷载的施加，接触法向各向异性的主方向θa均逐渐接近大主应力方向，也即试样发现明显的应力诱发各向异性；与高围压情况相比，低围压下数值试样的θa越靠近大主应力方向；当采用非球形颗粒时，与纯球形颗粒相比，围压对θa变化规律的影响程度减弱。