砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究

1、第 32 卷 第 10 期岩土工程学报Vol.32 No.102010 年10 月Chinese Journal of Geotechnical EngineeringOct. 2010砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究史旦达 1, 3，周健 2, 3，刘文白 1，邓益兵 1（1. 上海海事大学海洋环境与工程学院，上海 200135；2. 同济大学地下建筑与工程系，上海 200092； 3. 同济大学岩土及地下工程教育部重点实验室，上海 200092）摘 要：基于 PFC2D 非圆颗粒单元的二次开发，对砂土直剪力学过程进行了非圆颗粒仿真模拟，分析了数值试样的应 力剪胀关系并与实际砂

2、土进行对比，探讨了颗粒位移与颗粒旋转特征及其与剪切带演化的内在关联，研究了主应力 与主应变增量的非共轴效应，揭示了细观组构各向异性的演化规律及其与宏观剪切强度之间的宏细观关联。研究结果 表明，数值试验能够较好的模拟实际砂土的应力剪胀关系和剪切过程主应力与主应变增量的非共轴效应；剪切带的 演化与颗粒位移和颗粒旋转密切相关，颗粒形状影响剪切带的厚度；试样宏观的剪切强度主要受控于粒间法向接触力 的分布及其各向异性演化；整个加荷过程中，剪切带内大主应力的偏转方向与法向接触力各向异性的主方向保持了良 好的一致性。关键词：直剪试验；非圆颗粒；剪切带；非共轴效应；组构各向异性中图分类号：TU441.5文献标

3、识码：A文章编号：10004548(2010)10155709 作者简介：史旦达(1979 )，男，浙江舟山人，博士，副教授，主要从事细观岩土力学、海洋岩土力学、港口软基处 理等方面的研究工作。E-mail: 。Exploring macro- and micro-scale responses of sand in direct shear tests by numerical simulations using non-circular particlesSHI Dan-da1, 3, ZHOU Jian2, 3, LIU Wen-bai1, DENG Yi-bing1(1. Colleg

4、e of Ocean Environment and Engineering, Shanghai Maritime University, Shanghai 200135, China; 2. Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China; 3. Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai 2000

5、92, China)Abstract: The main objective of this study is to explore the macro- and micro-scale responses of sand in the direct shear tests by discrete element simulation. Non-circular particle element is self-developed in particle flow code in 2-dimension (PFC2D). The stress-dilatancy relationship of

6、 numerical samples is compared with that of real sand. The evolution of shear band is studied in relation with the variation of particle displacement and rotation. The effect of non-coaxiality between the directions of principal stress and principal strain increment is examined. Special attention is

7、 focused on the evolution of stress-induced fabric anisotropy during the shear loading. It is found that the discrete element simulation has a good ability to reproduce the stress-dilatancy relationship and non-coaxiality effect of real sand. The formation of shear band is strongly dependent on the

8、distributions of particle displacement and rotation. A strong correlation between the shear strength and the evolution of normal contact force is evident in the simulation results. Coincidence between the orientation of normal contact force anisotropy andthe direction of the major principal stress h

9、as been observed during all stages of shear loading.Key words: direct shear test; non-circular particle; shear band; non-coaxiality effect; fabric anisotropy0引言直剪试验是室内研究土体剪切力学性状的重要试 验方法，然而，常规的直剪仪只能从宏观上测得土体 的抗剪强度指标，不能分析剪切过程试样内部应力应 变的变化。最近，刘文白等1开发了半膜式细观结构 观测直剪仪，利用数字照相变形量测技术，研究了直 剪过程中砂土试样应力应变场的分布规律，但该试

10、验没有分析剪切面处颗粒的滑移、旋转等相对运动规律。 数值仿真模拟作为对室内试验的有力补充，目前已被 广泛用于包括直剪试验在内的土体力学性质的仿真分 基金项目：国家自然科学基金项目（50909057，90815008）；上海高校 选拔培养优秀青年教师专项基金项目（shs08004）；同济大学岩土及地 下工程教育部重点实验室开放基金项目（KLE-TJGE-0904）；上海市教 委科研创新项目（11YZ132）收稿日期：201001111558岩 土 工 程 学 报2010 年析。Potts 等2、Tejchman 和 Bauer3分别采用各向同 性弹塑性模型和亚塑性模型，利用有限元法分析了砂 土的

11、直剪力学性状，研究了剪切带内的应变局部化问 题。有限元法可以较好的分析剪切过程的应力应变关 系，但无法从颗粒层面揭示剪切强度和变形的细观力 学机理，也无法反映颗粒形状、孔隙率等试样要素变 化对模拟结果的影响。离散元法在分析岩土体大变形 破坏问题及其宏细观机理研究方面具有独特的优势。 刘斯宏和徐永福4在铝棒堆积体直剪试验的基础上， 利用离散元分析了粒状材料的直剪力学行为。Masson 和 Martinez5基于 PVC 材料细观力学特性，对松、密 试样的直剪特性进行了离散元对比分析。 Cui 和 OSullivan6针对均一粒径钢球集合体，利用三维离散 元分析了模拟颗粒材料的直剪宏细观力学响应。

12、Wang 等7则针对均一粒径玻璃圆球集合体，利用离散连 续耦合分析方法分析了剪切带的演化规律。总结以上 文献可知，目前关于直剪试验的仿真模拟大多针对铝 棒、钢球等理想代用材料，数值模拟也以纯圆形颗粒 居多，而对于实际的砂土，其材料属性不同于理想材 料，且具有明显的角粒特征，针对实际砂土的直剪行 为开展非圆形颗粒的仿真模拟研究，具有重要的研究 意义和应用价值。本文利用颗粒离散元 PFC2D 程序，开发了二维椭 圆形颗粒单元，并将其用于砂土直剪力学性状的仿真 模拟。数值模拟着重研究以下几个关键问题：宏观 应力剪胀关系及与实际砂土的对比；颗粒位移、 颗粒旋转与剪切带演化；主应力与主应变增量的非 共轴

13、效应；组构各向异性的演化规律；宏细观 力学量之间的关联。1数值试样制备1.1椭圆形颗粒单元关于 PFC2D 中椭圆形非圆颗粒和角粒形非圆颗粒 的开发工作，文献8，9已进行了详细的阐述，限于 篇幅，本文在此不再重复赘述。非圆颗粒单元的开发1.2 重力沉积法制样数值试样制备仍采用文献8中的重力沉积法，具 体步骤为：首先生成四周墙体，在高 30 mm、宽 60 mm 的二维空间内，按初始孔隙率 0.4 生成最大粒径 0.8 mm、最小粒径 0.4 mm，级配服从均匀分布的初始 纯圆颗粒，如图 2（a）所示；删去纯圆颗粒，在同 一位置用等效的椭圆形颗粒代替，颗粒长轴的初始方 向在 0180范围内随机定

14、向，如图 2（b）所示；施加重力，椭圆形颗粒在重力作用下沉积，并达到 初始平衡状态，如图 2（c）所示；在高度 20 mm 的 位置施加顶部墙体，同时删去顶墙以上的颗粒，如图 2（d）所示；卸除重力，施加上覆固结压力，等待 施加剪切荷载。图 2 数值试样制备（Se=1.4）Fig. 2 Sample preparation for numerical simulations (Se=1.4)直剪试验二维数值试样的外观尺寸与常规室内直 剪仪的尺寸相仿（直剪仪高 20 mm、直径 60 mm）。 为了减少生成颗粒的数量以提高计算效率，数值试样 平均粒径（d50=0.6 mm）略大于福建标准砂平均粒

15、径（d50=0.34 mm），但级配的分布与实际标准砂保持一致10。不同 Se 试样的初始试样信息汇总于表 1。由表 1 可知，不同 Se 试样的初始孔隙率 n0 基本保持在 0.175 左右，但非圆颗粒试样的初始平均接触数(Cn)ini 明显 大于纯圆颗粒试样，随着 Se 的增加，试样的(Cn)ini 逐 渐增加。表 1 数值试样初始信息Table 1 Initial parameters of numerical samples需要遵循“面积等效原则”和“质量等效原则”，对于长短轴比 Se1.01.21.41.6椭圆形颗粒，定义颗粒的长短轴之比为形状指标，用颗粒数量340734063402

16、3413Se 来表示。全文中取 Se=1.0（纯圆），1.2，1.4，1.6 四固结后高度/mm19.9619.9019.8619.91种情况进行对比分析。图 1 给出了椭圆形颗粒单元的示意图。初始孔隙率 n0初始平均接触0.1750.1760.1770.173数(Cn)ini3.7074.1244.2394.302图 1 椭圆形颗粒Fig. 1 Elliptical particles1.3细观参数设置数值模拟中，颗粒与颗粒之间、颗粒与墙体之间 的接触特性均采用线性接触模型来表述，需要设定的 细观参数有颗粒/颗粒法向接触刚度和切线接触刚度第 10 期史旦达，等. 砂土直剪力学性状的非圆颗粒模

17、拟与宏细观机理研究1559(kn)p，(ks)p，颗粒/墙体法向接触刚度和切向接触刚度 (kn)w，(ks)w，颗粒/颗粒和颗粒/墙体摩擦系数 fp，fw， 颗粒密度rs。本文作者在文献10中已完成了福建标准 砂二维模拟的细观参数标定工作，本文数值模拟仍然 采用文献10的结果。考虑到直剪仪墙体可视为刚性 体，所以数值模拟中颗粒/墙体接触刚度取为颗粒/颗 粒接触刚度的 2 倍。细观参数的具体数值汇总于表 2。表 2 细观参数取值Table 2 Values of micro-mechanical parameters2宏观力学响应分析2.1宏观力学特性曲线数值试验可以得到剪应力比水平位移(t/s

18、Nu) 和竖向位移水平位移(hu)关系曲线，如图 4 所示。 图 4（b）中竖向位移 h 为正表示体积剪胀。图 4 中实 心点标示处表示剪应力比峰值点和对应的体积剪胀峰 值点。细观(kn)p/参数(Nm1)数值(ks)p/fp(Nm1)(kn)w/ (Nm1)(ks)w/fw(Nm1)210811080.5410821080.5大小3细观rs/(kgm )参数Se=1.0Se=1.2Se=1.4Se=1.6数值2643170117381805大小1.4加荷条件文献4研究表明，由于剪切盒内壁摩擦的影响， 室内直剪试验测得的剪切强度在土体剪胀时，比土体 实际的强度要高；而在土体剪缩时，则比土体实际

19、的 强度要低。所以，计算剪切面上的剪应力比t/sN 时应 该恰当考虑侧壁摩擦效应的影响。图 3 给出了数值模 拟剪切过程示意图，本文不针对固结应力水平展开重 点讨论，所有试样施加同一上覆压力 p=200 kPa。剪 切过程中，围成下盒的、号墙体以 0.01 mm/min 的剪切速率向右移动，围成上盒的、号墙体保持固定，顶部号墙体由伺服机制控制保 持上覆压力不变，号墙体发生的竖向位移用 h 表示； 当水平剪切位移 u 达到 10 mm 时，加荷终止。图 3 数值模拟加荷示意图Fig. 3 Schematic diagram of shear loading in numerical simula

20、tions数值试验可以全程监测各个墙体的位移量和墙体上的作用力，剪切面上的剪应力比t/sN 可由下式计算：图 4 数值试样宏观力学特性曲线Fig. 4 Macro-mechanical response of numerical samples分析图 4，可以得到以下几点规律：二维数值 试样在 0.175 左右的初始孔隙率条件下，表现出典型 密实砂样的剪切性状，剪切强度曲线的峰后软化和剪 切过程的体积剪胀现象均表现的十分明显；颗粒形 状影响剪切强度，随着 Se 的增大，试样的峰值强度和 残余强度（u=10 mm）均明显增加，颗粒形状对初始 剪切刚度的影响并不大，不同 Se 试样强度曲线的初始

21、斜率基本相同；对于不同 Se 试样，峰值强度大概出 现在水平位移 23 mm 之间，且随着 Se 的增大，出 现峰值强度时的位移水平逐渐增加；对于不同 Se 试 样，除加荷初期发生少量体积剪缩外，剪切过程试样 均发生明显的体积剪胀，对应峰值强度点时，不同 Se 试样的体积剪胀量较为接近，但当加荷进行至残余强 度时，纯圆试样（Se=1.0）的体积剪胀已基本趋于稳 定，而非圆试样仍有较明显的剪胀发生。2.2应力剪胀关系（1）Rowe-Davis 应力剪胀公式Rowe 于 1962 年提出了平面应变条件下颗粒集合 t = T = N4 - N6 - T5，(1)体的应力剪胀关系11：s NNN5 +

22、 T4 + T6sin f+ siny cr 式中，N4、T4，N5、T5，N6、T6 分别表示、号墙体上所受的法向作用力和切向作用力，见图 3。sinfps =1 + sinfcrsiny，(2)式中，fps 表示平面应变(plain strain)条件下的内摩擦1560岩 土 工 程 学 报2010 年角，fcr 表示临界状态(critical state)内摩擦角，y表示剪 胀角。Davis12最早提出了直剪试验内摩擦角与平面应 变试验内摩擦角之间的换算关系，即力剪胀公式。值得注意的是，在图 5 和表 3 中，无 论是对于实际的砂土还是数值试样，在实测值中，横 坐标(tan(fds)r+

23、sinyp)/cosyp 的数值均略大于纵坐标 tan(fds)p，对于非圆数值试样，这一数值上的差异表现tan fds= cosy sinfps 1 - siny sin fps，(3)的更为明显，其原因与试样到达残余强度时，其剪胀角并非完全为零有关，尤其对于非圆数值试样，由 2.1式中，fds 表示直剪试验（direct shear）的内摩擦角。 将式（2）代入式（3），消去fps，得到fds，fcr，y三者之间的关系节可知，残余强度时试样仍有较明显的体积剪胀发生。tanfds= sinfcr + siny cosy。(4)当试样处于临界状态时，通常认为剪胀角y=0， 由式（4）可得rta

24、n (fds )= sin fcr，(5)式中，(fds)r 表示直剪试验的残余内摩擦角。 将式（5）代入式（4），得到tan (fds )+ siny图 5 应力剪胀关系tanfds=r cosy。(6)Fig. 5 Stress-dilatancy relationship式（6）即为描述砂土直剪应力剪胀关系的Rowe-Davis 公式。（2）数值试样的应力剪胀关系引用上述 Rowe-Davis 公式来分析数值试样的应 力剪胀关系。图 4 中，数值试验内摩擦角fds 和剪胀 角y分别可由下式计算： t 2.3剪切强度与初始平均接触数之间的关系在文献8中，本文作者通过砂土双轴试验数值模 拟初

25、步探讨了抗剪强度与初始平均接触数之间的关 系，本文在此进一步分析直剪试验过程中剪切强度与 初始平均接触数之间的关联。图 6 给出了数值试样峰 值内摩擦角(fds)p、残余内摩擦角(fds)r 与试样初始平均 接触数(Cn)ini 之间的关系。stan fds =N，(7)tany = dhdu。(8)计算图 4 中不同 Se 试样的峰值内摩擦角(fds)p、峰 值剪胀角yp 和残余内摩擦角(fds)r，并将数值汇总于表 3，表 3 中 up 表示峰值强度对应的水平位移大小。表 3 中还列出了文献13，14中关于实际砂土的分析结果。图 5 为公式（6）理论计算值与数值试验实测值之 间的对比关系。

26、由图 5 可知，对于不同 Se 试样，数值 试验数据点与理论预测值均非常接近，数值试样表现 出的应力剪胀关系符合实际砂土的 Rowe-Davis 应图 6 剪切强度与初始平均接触数关系Fig. 6 Shear strength versus initial averaged contact number表 3 应力剪胀关系实测值Table 3 Measured values of stress-dilatancy relationship试样(fds)p(fds)rypup/mmtan(fds)p(tan(fds)r+sinyp)/cosypSe=1.021.315.96.52.230.390

27、.40Se=1.224.220.37.12.500.450.50Se=1.428.824.27.72.750.550.59Se=1.631.027.58.13.020.600.6640.529.714.70.850.85oura sand1440.733.011.8 0.860.87Leight Buzzard sand13 (e0=0.565, p=170 kPa)Toy (e0=0.65, p=49 kPa)第 10 期史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究1561分析图 6 可知，试样的峰值内摩擦角和残余内摩 擦角与初始平均接触数之间成良好的线性拟合关系， 本文结果

28、和文献8的研究结果进一步说明，无论是平 面应变状态还是直剪状态，试样所能发挥的宏观剪切 强度与试样初始获取的粒间平均接触数的多少密切相 关。3细观力学响应分析3.1颗粒位移与剪切带剪切带的发生和演化与试样内部应变局部化现象 密切相关。数值试验的优势在于可以直观的分析颗粒 的水平位移和竖向位移，研究内部剪应变和体应变的 局部化现象，揭示剪切带的发生机制与演化过程。对 于纯圆试样，文献4，5均已作过研究。本文侧重于 讨论颗粒形状的影响，分析颗粒形状变化对剪切带厚 度的影响。以 Se=1.4 试样为例，图 7 给出了 u=10 mm 时，不 同竖向位置颗粒发生的水平位移和竖向位移分布。分 析图 7（

29、a）可以清晰的看出，变形的局部化发生在狭 长的剪切带内，且剪切带厚度的分布是不对称的，位 于下盒的剪切带厚度要明显大于上盒，对于 Se=1.4 试 样，剪切带的厚度约为 7.6 mm。分析图 7（b）可以 发现，较大的竖向位移集中在剪切带内，说明试样宏 观的体积剪胀主要由剪切带内颗粒的运动控制。值得 注意的是，竖向位移存在负值，表示有些颗粒发生了 向下运动，这说明颗粒的运动轨迹是不规则的，但向 下运动的最大位移量要比向上运动的最大位移量小很 多。图 7 颗粒位移分布图Fig. 7 Distribution of particle displacement目前，关于剪切带厚度的研究大多关注于颗粒

30、平 均粒径的影响，而关于颗粒形状影响的研究相对缺乏。 对于实际砂土，每个颗粒形状各异，很难对形状指标 变化对剪切带的影响展开研究，而在数值试验中可以 很快捷的生成单一颗粒形状的仿真试样，为这一领域 的深入研究提供了条件。在图 7（a）中，取竖向间隔 为 1 mm，计算每 1 mm 试样高度范围内颗粒水平位移的平均值，并将不同 Se 试样的结果绘于同一幅图中， 如图 8 所示。分析图 8 可知，随着 Se 的增大，试样剪 切带厚度逐渐增加，Se=1.0，1.2，1.4，1.6 时，相应 的剪切带厚度 D 分别约为 6.0，6.5，7.6，8.0 mm，这 一结果表明，当颗粒的角粒特征增强时，不仅

31、表现出 更大的剪切强度和体积剪胀特征，而且剪切过程中试 样内部应变局部化的范围也会扩大。图 8 颗粒平均水平位移曲线Fig. 8 Curves of averaged particle horizontal displacement3.2颗粒旋转除颗粒位移外，颗粒旋转也是表征剪切带演化过 程的重要细观力学量。同样以 Se=1.4 试样为例，图 9（a）给出了 u=10 mm 时，不同竖向位置颗粒的颗粒旋转量分布。图 9（a）中，正值表示颗粒逆时针旋转， 负值表示颗粒顺时针旋转。由图 9（a）可知，较大的 颗粒旋转主要集中在剪切带区域内，而剪切带之外， 特别是接近上下端面位置，颗粒旋转量逐渐减小

32、。综 合分析不同 Se 试样情况，在图 9（a）中，取竖向间隔 为 1 mm，计算每 1 mm 试样高度范围内颗粒旋转量（正、负均取绝对值）的平均值，并将不同 Se 试样的 结果绘于同一图中，如图 9（b）所示。分析图 9（b） 可得以下几点规律：纯圆试样的颗粒旋转要明显大 于非圆试样；颗粒旋转剪切带内最大，并向剪切带 外两侧逐渐递减；无论是纯圆试样还是非圆试样， 剪切带内的颗粒旋转峰值均出现在下盒，其原因可能 与 3.1 节分析中，下盒的剪切带厚度较大有关。图 9 颗粒旋转分布图Fig. 9 Distribution of particle rotation3.3主应力与主应变增量的非共轴效

33、应（1）平均应力与应变率1562岩 土 工 程 学 报2010 年k k 在 PFC2D 中，可以通过设置在试样内部的量测圈 直接测得试样平均应力和应变率。平均应力计算公式 为a = 1212s xyarctan()o xx - s yy2de xy，(12)sij 1 N=fi l j(i, j=x, y)，(9)b =arctan(2de)xx - de yy。(13)V k =1式中，V 表示量测圈的体积（面积），Nc 表示量测圈内图 11 给出了加荷过程中a，b的变化过程。分析图 11 可得几点规律：固结状态时，竖向应力大于水的颗粒接触数总数，fik 表示接触点 k 处接触力矢量在平应

34、力，加荷开始后，大主应力方向从90左右开始i 方向上的分量，ljk 表示接触点 k 处颗粒支矢量（连接两个接触颗粒圆心的矢量）在 j 方向上的分量。 应变率计算公式为偏转，峰值强度后，大主应力的方向逐渐趋于稳定；加荷开始后，应变增量主轴的偏转要比应力主轴的 偏转快的多，在很小的剪切位移水平下，大主应变增eij= dn id(i, j=x, y)，(10)量方向即发生了很大的偏转，并趋于稳定，这一变化x j式中，dvi 表示量测圈内两相邻颗粒 t 时刻在 i 方向上 的速度矢量差，dxj 表示两相邻颗粒 t 时刻颗粒支矢量 在 j 方向上的分量。由式（10）可计算得到 dt 时步内的应变增量为规

35、律与文献15，16的试验结果是一致的，说明数值 试验能够较好的模拟非共轴效应；非共轴效应在峰 值强度前表现的较为明显，峰值强度后，随着塑性变 形的增加，非共轴效应逐渐减小；颗粒形状对大主 应变增量方向角b的影响要小于对大主应力方向角adeij = eij dt。(11)的影响，对于不同 Se试样，峰值强度后，a角大致在由式（9）和（11）可以计算得到各应力分量sxx，syy，sxy 和应变增量分量 dexx，deyy，dexy 的大小，进 一步可计算得到主应力和主应变增量。本文量测圈的布置共分 L1，L2，L3 三层，每层布 置 8 个量测圈，量测圈的直径为 6 mm，如图 10 所示。图 1

36、0 量测圈的布置Fig. 10 Layout of measurement circles（2）非共轴效应分析Stroud15、Tatsuoka 等16均通过砂土的单剪试验 (simple shear test)研究证实，伴随剪应力的施加，应变 增量主轴的旋转并不与应力主轴的旋转同步，应力主 轴的旋转滞后于应变增量主轴，主应力与主应变增量 之间存在非共轴效应。数值试验能否有效模拟非共轴 效应是本节研究重点。Jewell17试验研究表明，直剪试验中，剪切带内 砂土的变形近似于单剪试验中砂土的变形模式。本文 数值试验中，布置在 L2 层的第 9 至第 16 号量测圈基 本覆盖了剪切带区域（如图 1

37、0 所示）。取第 9 至第 16 号量测圈各应力分量和应变增量分量的平均值来分析 加荷过程中应力主轴和应变增量主轴的旋转情况。4050区间范围内，b角在 4045区间范围 内。图 11 大主应力与大主应变增量方向角Fig. 11 Inclination of major principal stress and major principal strain increment along horizontal plane3.4组构各向异性演化规律对于散粒体材料，加荷过程试样表现出的宏观强 度和变形特性与其内部细观组构的演化规律密切相 关。Rothenburg 和 Bathurst18提出可用傅

38、里叶函数来 近似描述单调剪切荷载作用下粒间接触法向、粒间法 向接触力和切向接触力各向异性的演化规律，其数学 表达式分别为用a，b分别表示大主应力方向和大主应变增量方向与水平方向的夹角，则a，b分别可由下式计算E (q ) = 1 1 + a cos 2 (q - q2a)，(14)第 10 期史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究1563接触力系（a）加荷前（b）峰值强度时（u=2.75 mm）（c）残余强度时（u=10 mm）图 12 组构各向异性演化（Se=1.4）Fig. 12 Evolution of fabric anisotropy（Se=1.4）fn (q

39、) = f0 1 + an cos 2 (q - qn )，(15)接触力、切向接触力的各向异性分布玫瑰图。玫瑰图ft (q ) =- f0 at sin 2 (q - qt )。(16)的绘制每 10分一个区间，统计接触法向时，M 表示式中 f0 是相对于所有接触的平均法向接触力；qa，qn，qt 分别是接触法向各向异性、法向接触力各向异 性和切向接触力各向异性的主方向；a，an，at 是傅立 叶系数，它们的数值大小反映对应细观组构参量的各 向异性程度。以 Se=1.4 试样为例，分析加荷过程中细观组构各 向异性的演化规律。组构各向异性的分析着重针对剪 切带区域，数值试验加荷过程中，对 L2

40、 层内（见图 10）颗粒之间的接触法向、法向接触力和切向接触力 的演化进行统计分析。图 12（a）（c）分别给出了 加荷前、峰值强度和残余强度时粒间接触法向、法向接触点总数、DM(q)表示接触法向落入该角度区间的接触点个数；统计接触力时，fn(q)，ft(q)分别取接触法 向落入该角度区间内所有接触点处法向接触力和切向 接触力的平均值。图 12 中虚线为公式（14）（16） 傅里叶函数的拟合结果，图中箭头方向表示拟合得到 的组构参量各向异性的主方向。分析图 12 可知，接触法向的拟合曲线呈“椭圆状” 图形，在整个剪切过程中，接触法向的主方向都偏于 竖直方向，但接触法向分布于各个方向的各向异性程

41、 度并不大。法向接触力和切向接触力的拟合曲线分别 呈“花生状”和“花瓣状”图形，加荷前，数值试样 由重力沉积生成并施加了竖向固结压力，由图 12（a）1564岩 土 工 程 学 报2010 年可知，试样内的主接触力系分布于竖直方向，粒间法 向接触力和切向接触力的主方向都偏于竖直方向，初 始接触力的分布并非各向同性；剪切荷载施加后，随 着主接触力系逐渐向斜 45方向倾斜，法向接触力和 切向接触力的主方向随之倾斜，且各向异性的程度明 显增强；剪切至残余强度时，法向接触力和切向接触 力的各向异性程度减小，且粒间平均法向接触力 f0 减 小。图 13 给出了组构各向异性参数的变化规律。图13(b)中实

42、线所示为 3.3 节计算得到的大主应力方向角a。分析图 13 可得几点宏细观力学量之间的关联：与 an 相比，整个加荷过程中，a 和 at 的数值和变化 量均较小，说明试样宏观的剪切强度主要受控于粒间 法向接触力的分布及其各向异性演化；an，at 出现 峰值时的位移水平与试样呈现峰值强度时的位移水平 十分接近，说明试样宏观强度的变化与内部接触力各 向异性的演化保持同步；对比分析qn 和a的变化规 律，剪切带内大主应力的偏转方向与法向接触力各向 异性的主方向保持了良好的一致性，进一步说明宏观 应力与粒间接触力之间存在内在关联。图 13 组构各向异性参数的变化（Se=1.4）Fig. 13 Var

43、iation of fabric anisotropy parameters（Se=1.4）4 结 论本文基于 PFC2D 非圆颗粒单元的二次开发，对砂 土直剪力学过程进行了非圆颗粒模拟，得到的主要结 论有：（1）数值试样表现出密实砂样的直剪力学性状， 强度的峰后软化和体积剪胀均十分明显，加荷过程表 现出的应力剪胀关系符合实际砂土的 Rowe-Davis 应力剪胀关系。颗粒形状影响峰值强度和残余强度， 颗粒形状的影响在细观机理上与试样的初始平均接触数有关。（2）颗粒的水平位移分布清晰的表明，变形的局 部化发生在狭长的剪切带内，较大的竖向位移和颗粒 旋转均集中在剪切带内，试样宏观的体积剪胀主要由

44、 剪切带内颗粒的运动控制。颗粒形状影响剪切带厚度， Se=1.0，1.2，1.4，1.6 时，剪切带厚度 D 分别约为 6.0， 6.5，7.6，8.0 mm，说明颗粒的角粒特征增强时，试 样内部应变局部化的范围将扩大。（3）数值试验能够较好的模拟主应力与主应变增 量的非共轴效应，非共轴效应在峰值强度前表现的较 为明显，峰值强度后逐渐减小；非共轴效应中，应力 主轴的旋转滞后于应变增量主轴。颗粒形状对大主应 变增量方向角b的影响要小于对大主应力方向角a的 影响。（4）剪切过程中，接触法向的主方向都偏于竖直 方向，但各向异性的程度并不大。加荷前，粒间接触 力的主方向偏于竖直方向；剪切荷载施加后，伴

45、随主 接触力系逐渐向斜 45方向倾斜，法向接触力和切向 接触力的主方向随之倾斜，且各向异性的程度增强； 剪切至残余强度时，法向接触力和切向接触力的各向 异性程度减小，且粒间平均法向接触力 f0 减小。（5）各向异性参数的定量研究表明，试样宏观的剪切强度主要受控于粒间法向接触力的分布及其各向 异性演化，且宏观强度的变化与内部接触力各向异性 的演化保持同步；整个加荷过程中，剪切带内大主应 力的偏转方向与法向接触力各向异性的主方向保持了 良好的一致性。参考文献：1 刘文白, 张 辉, 邓一兵. 基于 DPDM 技术的砂土直剪试 验剪切过程的应力场分析J. 中国水运, 2008, 8(7): 2352

46、39. (LIU Wen-bai, ZHANG Hui, DENG Yi-bing. Analysisof sand stress field during direct shearing by DPDM techniqueJ. China Water Transport, 2008, 8(7): 235239. (in Chinese)2 POTTS D M, DOUNIAS G T, VAUGHAN P R. Finite elementanalysis of the direct shear box testJ. Gotechnique, 1987,37(1): 1123.3 TEJCH

47、MAN J, BAUER E. FE-simulations of a direct and a true simple shear test within a polar hypoplasticityJ. Computers and Geotechnics, 2005, 32(1): 116.4 刘斯宏, 徐永福. 粒状体直剪试验的数值模拟与微观考察J. 岩石力学与工程学报, 2001, 20(3): 288292. (LIUSi-hong, XU Yong-fu. Numerical simulation for a direct box shear test on granular ma

48、terial and microscopic第 10 期史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究1565considerationJ. Chinese Journal of Rock and Mechanics and Engineering, 2001, 20(3): 288292. (in Chinese)(ZHOU Jian, SHI Dan-da, JIA Min-cai, et al. Numerical simulation of mechanical response on sand under monotonic5 MASSON S, MARTINEZ J.

49、 Micromechanical analysis of theshear behavior of a granular materialJ. Journal of Engineering Mechanics, 2001, 127(10): 10071016.11loading by particle flow codeJ. Journal of Tongji University, 2007, 35(10): 12991304. (in Chinese)ROWE P W. The stress dilatancy relation for static6 CUI L, OSULLIVAN C. Exploring the macro- andequilibrium of an assembly of particles in contactC/micro-scale response of an idealized granular material in the direct shear apparatusJ. Gotechnique, 2006, 56(7): 455Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1962: