颗粒材料定向剪切特性的空心圆柱离散元模拟

史旦达，杨晨捷，薛剑峰，王 威

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

1 研究背景

在基坑开挖、路堤填筑和预制桩贯入等岩土工程实际问题中，地基土体常常处于复杂应力状态，其土单元体的水平面上会受到正应力和剪应力的共同作用，此时主应力方向不再处于竖直或水平方向，即存在主应力方向的偏转[1］。由于土体的原生各向异性，在不同主应力方向加载条件下，土体的剪切力学响应会产生一定差异[2］。空心圆柱仪（Hollow Cylinder Apparatus，HCA）是目前研究复杂应力条件下土体各向异性剪切特性的常用试验工具，针对砂土等散体颗粒材料，Yang等[3］、冷艺等[4］、许成顺等[5］、于艺林等[6］和Cai等[7］均采用空心圆柱仪研究了不同主应力方向定向剪切条件下砂土的强度和变形特性，试验研究总结了很多有益的规律。然而，常规的HCA试验只能分析试样宏观的力学响应，尚不能从细观层面分析试样内部剪切带、接触力链及细观组构的演化规律，这不利于土体各向异性宏细观规律的深入研究。鉴于此，近年来针对空心圆柱仪的构造和加荷特点，采用离散单元法（Discrete Element Method，DEM）对空心圆柱试验展开宏细观模拟的研究受到了广泛关注。例如，Li等[8-10］基于离散元PFC3D（Particle Flow Code in 3-Dimension）数值模拟平台构建空心圆柱数值试样，通过捕获颗粒试样内部孔隙率和剪应变率的分布探讨了剪切带的演化规律，并且分析了主应力方向和中主应力对非共轴特性的影响。Farhang等[11］采用离散元TRUBAL程序建立了空心圆柱数值试样，数值模拟分析了主应力方向变化对试样宏观剪切强度的影响，但其研究中没有考虑拉伸状态下试样的力学响应。除上述空心圆柱模拟外，Li等[12］也曾通过对三维多面体离散元模型进行非比例加载来实现主应力方向的偏转，数值模拟揭示了主应力轴循环旋转下颗粒试样的宏观变形特性和组构各向异性演化规律。然而，与三轴、直剪等常规的单元试验细观模拟不同，由于空心圆柱试验比较复杂，关于空心圆柱的离散元数值建模研究国内外均刚刚起步，在侧向柔性边界条件、扭矩施加与控制方式、空心圆柱试样三维组构分析等方面的研究均有待深入。

本文在借鉴已有的空心圆柱离散元建模经验基础上，引入“动态更新法”对扭矩的施加效果进行改进，并在此基础上，进行5种不同大主应力方向角条件下颗粒材料数值试样的定向剪切模拟，重点分析大主应力方向角对试样剪切强度的影响，并与实际砂土空心圆柱试验结果进行对比。最后，从细观角度探讨试样剪切带的演化规律，并通过分析试样内部局部孔隙比和配位数的变化揭示剪切带的演化机理。

2 空心圆柱离散元建模

2.1 空心圆柱试样单元体的应力状态空心圆柱扭剪试验时，其单元体的典型受力状态如图1所示。图1（a）中，pi、po、W和T分别为施加在空心圆柱试样上的内围压、外围压、竖向荷载和环向扭矩；图1（b）中，r、θ、z坐标分别表示径向、环向和竖向，σz、σθ、σr和 τzθ分别为作用在单元体上的竖向、环向、径向正应力和剪应力；图1（c）中，σ1、σ2、σ3和α分别表示大、中、小主应力和大主应力方向角（也即大主应力与竖直方向的夹角）。

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

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

由式（1）—式（4），可进一步计算得到各主应力大小及大主应力方向角：试样峰值内摩擦角φp的计算公式为：

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

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

2.2 数值试样制备及参数设置DEM数值试样亦为空心圆柱型，初始内径d0、外径D0和高度H0分别为60、100和200 mm。为了最大程度上消除空心圆柱试样应力应变分布的不均匀性，Sayao等[13］曾提出了如下3点建模准则：（1）D0-d0=40～120 mm；（2）0.65≤d0/D0≤0.82；（3）1.8≤H0/D0≤2.2。对照以上3点，本文数值试样的尺寸符合准则（1）和（3），而内外径之比（d0/D0=0.6）略小于准则（2），但数值模型的d0/D0与文献[4，7，14］所采用的实际空心圆柱砂样的尺寸是一致的。

考虑到计算效率，与实际砂土的级配相比，颗粒粒径做了适当的放大，数值试样的粒径范围为1.4～2.4 mm，平均粒径d50=1.9 mm。针对土力学单元试验的DEM模拟，一般可通过控制试样外观尺寸与颗粒平均粒径的比值（也称特征长度比值）大于某一定值来消除粒径放大可能带来的颗粒尺寸效应[15］。孙其诚等[16］研究表明，对于常规三轴试验的模拟，当特征长度比值大于50时，可以有效消除颗粒尺寸效应对试样抗剪强度的影响。在本文中，数值试样初始外径与颗粒平均粒径的比值D0/d50=52.6＞50，借鉴文献[16］三轴模拟总结的经验，在本文中可以近似忽略颗粒尺寸效应的影响。

数值试样的上、下边界为刚性墙体，为了能够有效模拟室内空心圆柱试验中橡皮膜的柔性边界效应，内、外边界各采用20面簇墙来模拟，每面簇墙能独立施加伺服围压，如图2所示。Zhao等[17］最早在三轴试验离散元模拟中采用分离式簇墙方法来反映橡皮膜效应，取得了较好的模拟效果。Li等[8］在空心圆柱模拟中也采用了类似的建模方法，并讨论了簇墙数量对剪切带形成的影响，其结果表明，簇墙数量越多模拟效果越好，但计算效率会显著降低；当簇墙高度Hs与颗粒平均粒径d50的比值小于10时，数值模拟已经能够较好的捕捉试样内部剪切带的发生。在本文中，内外簇墙均为20面，每面簇墙的高度Hs=10 mm，簇墙高度与颗粒平均粒径的比值Hs/d50=5.26＜10，可认为能够近似反映橡皮膜的柔性边界效应。

图2 空心圆柱数值试样的边界墙体建模

数值试样制备采用重力沉积法，初始孔隙比e0=0.563，颗粒数量20 227个，颗粒形状为三维纯球形颗粒，颗粒密度为2650 kg·m-3。颗粒与颗粒之间、颗粒与墙体之间均采用线性接触刚度模型，且法向与切向接触刚度均为2×106N·m-1。颗粒与颗粒之间的摩擦系数为0.8，根据Procter等[18］的室内试验研究，石英砂颗粒之间的平均粒间摩擦系数为0.49，由于本文数值试样采用的纯球形颗粒不能反映实际不规则形状砂土颗粒之间的旋转阻抗效应，因而粒间接触摩擦系数的取值比实际石英砂稍大，以近似体现颗粒旋转阻抗的影响。为了反映光滑墙面的加载效果，颗粒与墙体之间的摩擦系数均为0。离散元PFC3D中提供了二种阻尼机制来模拟土体的阻尼特性，一种为黏性阻尼机制，较适用于冲击荷载、强夯等动力问题模拟；另一种为局部阻尼机制，可用于静力和拟静力问题分析。本文采用局部阻尼机制，该方法通过在颗粒接触力合力的反方向施加一定大小的静态阻尼力以使颗粒系统尽快达到静力平衡状态，阻尼力的大小与设定的局部阻尼系数α有关，本文中α=0.7[19］。图2（c）给出了由DEM生成的空心圆柱数值试样。

2.3 荷载的施加方式与室内HCA试验相同，在数值模拟中，可以通过施加内围压、外围压、竖向荷载和环向扭矩来实现对数值试样的定向剪切。内、外围压的施加可通过PFC3D的围压伺服机制实现[20］。竖向荷载的施加可通过对上、下刚性墙体赋予一定的运动速率实现，本文中采用双向加载方式，上、下墙体以相同的速率相向运动，运动速率均为1 mm/s。环向扭矩的施加是实现主应力方向偏转的关键，Li等[8-10］和Farhang等[11］均采用在试样顶部设置一定厚度的颗粒加载层（下文也称“扭矩层”），并对加载层内的颗粒施加一定的角速度来实现环向扭矩的施加，但他们的模拟均采用了恒定角速度的方法（下文也称“定值法”），该方法在模拟效果上尚存在一定的不足。为了取得更好的模拟效果，本文对角速度的施加方式进行了适当改进。

如图3（a）所示，数值试样顶部的颗粒扭矩层厚度取为颗粒平均粒径d50的2倍，即为3.8 mm。在加载过程中，随着轴向应变的增加，初始扭矩层厚度内的颗粒会溢出初始扭矩层从而导致扭矩施加的不均匀。因此，为了保证整个加载过程中扭矩施加的均匀性，每加载100个时步由自编的程序自动动态更新一次扭矩层。通过对扭矩层内的颗粒施加绕着试样中心的旋转速度可以实施剪应力的施加，图3（b）给出了数值试样扭矩层颗粒的速度矢量图。

图3 数值试样顶部扭矩层的设置

根据上文所述，已有的数值模拟均对扭矩层内的颗粒施加定值角速度（例如，Li等[9］对扭矩层内的颗粒施加0.2 rad/s的恒定角速度），定值法的优势是实施方便、模拟效率较高，但缺陷是大主应力角的维持效果较差。针对这一缺陷，本文通过自编的扭矩控制函数来实现对扭矩层颗粒角速度的动态调整，在加载过程中，扭矩控制函数每隔100个时步计算一次当前的大主应力方向角α值，再由当前计算值与设定的目标值之间的差值，赋予扭矩层颗粒新的角速度，以保证在整个加载过程中，α值能持续逼近设定的目标值。以α=30°模拟工况为例，图4给出了定值法和本文动态更新法的效果对比。由图4可见，定值法在加载后期α值会明显偏离设定的目标值，而动态更新法可以在整个加载过程中较好地控制大主应力方向角。

数值模拟控制α=0°、30°、45°、60°、90°共5种工况。剪切开始前，通过围压伺服机制对数值试样施加各向均等围压，初始围压大小为100 kPa。

图4 扭矩施加方式的改进效果对比

3 宏观模拟结果分析

3.1 各应力分量-偏应变关系图5给出了5种大主应力角α条件下，数值试样各应力分量（σz、σq、σr、τzθ） 随偏应变ε13的变化规律。

图5 各应力分量随偏应变的变化规律

（1）由图5（a）可知，当α=0°时，试样处于纯压缩状态（也即常规三轴压缩状态），施加于试样上的剪应力τzθ接近于0；竖向应力σz在加载初期增长较快，试样表现出较大的初始刚度，当ε13=1%左右时，σz的增长幅度明显变缓，当ε13=8%左右时，σz达到峰值状态，最大值为505 kPa；在竖向加载过程中，内、外围压通过围压伺服机制保持恒定不变，以此来保证环向应力σθ和径向应力σr恒定不变（也即σθ=σr=100 kPa）。（2）由图5（b）可知，当α=30°时，试样处于扭剪压缩状态，通过 2.3节介绍的扭矩施加方式对数值试样施加剪应力τzθ，且竖向应力σz保持与剪应力τzθ相同的增长规律，τzθ和σz的峰值均出现在ε13=8.5%左右，峰值大小分别为81和247 kPa；在加载过程中，与α=0°工况一致，σθ和σr仍保持100 kPa恒定不变。（3）由图5（c）可知，当α=45°时，试样处于纯扭剪状态，试样上仅施加剪应力增量，而各正应力分量（σz、σθ、σr）均维持初始各向均等围压值不变；在整个加载过程中，剪应力的增长幅度较为平缓，剪应力τzθ的峰值出现在ε13=8.3%左右，峰值大小为60 kPa。（4）由图5（d）可知，当α=60°时，试样处于扭剪拉伸（伸长）状态，试样受到侧向正应力增量和环向剪应力增量的共同作用，通过调整围压伺服机制实现对内、外侧墙施加相同的压力增量，以保证在侧向增压过程中σθ=σr，且竖向应力σz维持初始值不变；剪应力τzθ和环向应力σθ的峰值均出现在ε13=8.7%左右，峰值大小分别为67和213 kPa。（5）由图5（e）可知，当α=90°时，试样处于纯拉伸状态（也即常规三轴拉伸状态），试样上的剪应力τzθ接近于0，通过侧向增压且维持竖向应力σz不变来实现竖向拉伸，σθ和σr的最大值出现在ε13=7.5%左右，峰值大小均为385 kPa。

3.2 大主应力方向角-偏应变关系由图5并结合式（7），可以计算得到数值模拟过程中实际的大主应力角α，以分析定向剪切的控制效果。图6（a）绘出了DEM模拟中大主应力方向角随偏应变的变化曲线，为了便于与实际砂土HCA试验结果进行对比，图6（b）同时给出了文献[21］福建标准砂HCA定向剪切试验得到的实测剪切应力路径（大主应力方向角α-加载时间t关系曲线）。分析图6可知，DEM模拟得到的剪切应力路径与砂土室内试验实测的应力路径十分相似，说明本文采用的空心圆柱离散元建模技术能够较好的用于砂土等颗粒材料定向剪切的宏细观模拟。

需要说明的是，图6中，与室内试验的实测曲线相比，30°和60°工况下DEM模拟曲线存在小幅波动，其原因与数值模拟和室内试验剪应力的不同施加方式有关。室内HCA试验时，剪应力通过顶部（或底部）加载板的旋转予以施加；而DEM模拟不能通过墙体的旋转直接施加剪应力，而是通过设置顶部颗粒扭矩层的方式施加剪应力，此时，剪应力的大小会受到扭矩层颗粒与试样颗粒之间接触力和配位数的波动影响，从而不可避免会引起一定的小幅波动，但这一小幅波动不会影响试样总体的宏细观力学响应。

图6 数值模拟中大主应力方向角的控制效果及与试验对比

3.3 试样剪切强度变化规律由图5所示的各应力分量变化规律并结合式（5），可计算得到大、小主应力（σ1、σ3）的大小，进一步得到应力比σ1/σ3随偏应变ε13的变化规律，如图7（a）所示。分析图7（a）可知，对比文献[7］所列的密实Portaway砂（相对密度Dr≈90%）定向剪切强度曲线规律（文献[7］中图4（a）），本文数值试样表现出近似实际密砂的应力-应变响应；α=0°工况下的强度曲线明显高于其他4种模拟工况，且试样的残余强度也明显高于其他4种工况；对于α=30°、45°、60°工况，由于剪应力的施加以使试样处于扭剪应力状态，其强度曲线的波动要明显大于纯压缩（α=0°）和纯拉伸（α=90°）状态。由图 7（a）可知，对于α=0°、30°、45°、60°、90°共 5种工况，峰值应力比（σ1/σ3）p的大小依次为5.042、4.405、4.001、3.545和3.854。

图7 数值试样的剪切强度变化规律

由图7（a）得到的上述应力比峰值，代入式（8），可计算得到数值试样的峰值内摩擦角φp。

图7（b）给出了峰值内摩擦角φp随大主应力方向角α的变化规律，为了与实际砂土试验结果展开对比，图 7（b）中还给出了已有文献中 Toyoura砂[14］、南京云母砂[6］、福建标准砂[4］的定向剪切试验结果，其中，Toyoura砂和南京云母砂为密实状态，福建标准砂为松散状态。分析图7（b）可知，由于DEM数值试样在试样密度、颗粒形状等方面与实际砂土试样存在一定的差异，因此从定量上分析，DEM试样的峰值内摩擦角数值与实际砂土尚存在差异；但从定性规律上看，DEM试样表现出与实际砂土一致的定向剪切强度变化规律：当α从0°增至60°时，φp随着α的增加而逐渐减小；而当α从60°增至90°时，φp值则又略有增加，也即α=60°工况下试样的抗剪强度为最低。Cai等[7］曾对砂土的空心圆柱剪切特性进行了总结，结果表明，由于砂土的初始各向异性特征，砂样的排水剪切强度最低值一般出现在α=60°～75°范围内。因此，本文DEM数值模拟较好的反映了实际砂土原生各向异性对后继剪切强度的影响规律。

4 细观力学响应分析

4.1 颗粒位移场及剪切带演化在本文采用的5种模拟工况中，α=30°工况，也即试样处于扭剪压缩状态时，试样剪切带的演化最为明显。图8（a）—（d）给出了不同应变水平下α=30°试样的颗粒位移场矢量图，从颗粒位移场的演化规律中可以清晰的看出内部剪切带的生成和演化过程，在加载初期，当应变水平较低时（图8（a）），试样以压缩变形为主且位移场分布较为均匀；随着应变水平的增加（图8（b）），试样内部逐渐产生变形的局部化，试样内部逐渐生成斜向的剪切带；当达到峰值剪切强度时（图8（c）），剪切带的形态已经清晰可辨且基本定型；当进入峰后阶段时（图8（d）），剪切带的厚度略微减小而倾角基本保持不变。李博等[22］曾在GDS空心圆柱仪上进行了平均粒径0.2 mm的玻璃微珠密实试样的定向排水剪切试验，图8（e）给出了该试验得到的α=30°工况下试样的剪切破坏照片。对比分析图8（d）和8（e）可知，数值试样整体的破坏特征以及剪切带的形态（厚度、倾角等）均与文献[22］玻璃微珠室内试验结果吻合较好，以上结果也从细观规律上反映了离散元模拟的有效性。

4.2 局部孔隙比和配位数演化通过在数值试样内部设置量测球可以分析试样内部局部孔隙比和配位数的演化规律。在距离试样顶面20、60和100 mm处设置了3层共12个量测球，每层4个量测球呈对称布置，量测球的直径均为16 mm，量测球的布置如图9所示。考虑到本文采用双向加载模式，沿着斜向的剪切带，试样下部与上部的变形规律基本对称，因此，考虑到监测和分析效率，量测球仅设置在试样上半部分。仍以α=30°工况为例，分析图8（d）所示的剪切带形态及图9中量测球的位置可知，7号、10号和12号量测球落入试样剪切带内，其他量测球位于剪切带外。

图8 数值试样剪切带演化规律及与室内试验对比（α=30°）

图9 数值试样内部量测球的布置形式

图10分别给出了各量测球测得的局部孔隙比和配位数随偏应变的变化曲线。分析图10，可以发现：试样内部局部孔隙比和配位数的变化规律可以从细观上较好的揭示剪切带的生成和演化过程，对于位于剪切带内的7号、10号和12号量测球，由其测得的局部孔隙比在2%～4%偏应变水平之间发生了激增，而对应的配位数则发生明显的递减。这一应变水平正是试样剪切带触发和逐渐生长的阶段；而位于剪切带外的其他量测球，由于没有发生明显的变形局部化，在整个加载过程中，局部孔隙比均缓慢增加，对应的配位数则均缓慢减少。需要指出的是，图10（a）所示的各量测球局部初始孔隙比均略小于试样的总体初始孔隙比（e0=0.563），其原因与试样施加100 kPa均等围压后，试样的体积发生一定的收缩有关（类似于室内砂样固结排水产生的体积收缩）。

5 结论

本文基于PFC3D离散元模拟平台，通过构建空心圆柱型数值试样，研究了定轴剪切条件下颗粒材料数值试样的宏细观力学响应，重点讨论了大主应力方向角α变化对数值试样剪切强度的影响，主要结论有：（1）与以往采用定值法相比，本文引入的扭矩层颗粒动态更新法能够较好的实现剪应力的施加，从而保证剪切过程中主应力方向的定轴模拟效果；（2）数值模拟能够较好的再现室内空心圆柱砂样定向剪切试验的应力路径和应力-应变响应，数值试样反映了实际砂土的初始各向异性特性，在α=60°工况下数值试样的抗剪强度为最低，数值试样呈现出与室内砂样一致的定向剪切强度变化规律；（3）数值模拟能够较好的捕获试样剪切带的发生与演化规律，且数值试样呈现出的剪切带形态与室内玻璃微珠试样较为一致，在细观上，数值试样的剪切带演化可由内部局部孔隙比和配位数的变化来表征。

需要说明的是，由于空心圆柱模拟的复杂性，本文所得的模拟结果还较为初步，后续将重点开展颗粒形状、主应力轴循环旋转等复杂条件下空心圆柱模拟的相关研究。