第32卷 第10期 岩 土 工 程 学 报 Vol.32 No.10 2010年10月 Chinese Journal of Geotechnical Engineering Oct. 2010 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究
史旦达1, 3，周 健2, 3，刘文白1，邓益兵1
（1. 上海海事大学海洋环境与工程学院，上海 200135；2. 同济大学地下建筑与工程系，上海 200092；
3. 同济大学岩土及地下工程教育部重点实验室，上海 200092）
摘要：基于PFC2D非圆颗粒单元的二次开发，对砂土直剪力学过程进行了非圆颗粒仿真模拟，分析了数值试样的应力–剪胀关系并与实际砂土进行对比，探讨了颗粒位移与颗粒旋转特征及其与剪切带演化的内在关联，研究了主应力与主应变增量的非共轴效应，揭示了细观组构各向异性的演化规律及其与宏观剪切强度之间的宏细观关联。研究结果表明，数值试验能够较好的模拟实际砂土的应力–剪胀关系和剪切过程主应力与主应变增量的非共轴效应；剪切带的演化与颗粒位移和颗粒旋转密切相关，颗粒形状影响剪切带的厚度；试样宏观的剪切强度主要受控于粒间法向接触力的分布及其各向异性演化；整个加荷过程中，剪切带内大主应力的偏转方向与法向接触力各向异性的主方向保持了良好的一致性。
关键词：直剪试验；非圆颗粒；剪切带；非共轴效应；组构各向异性
中图分类号：TU441.5 文献标识码：A 文章编号：1000–4548(2010)10–1557–09
作者简介：史旦达(1979– )，男，浙江舟山人，博士，副教授，主要从事细观岩土力学、海洋岩土力学、港口软基处理等方面的研究工作。E-mail: shidanda@https://www.doczj.com/doc/306066405.html,。
Exploring macro- and micro-scale responses of sand in direct shear tests by
numerical simulations using non-circular particles
SHI Dan-da1, 3, ZHOU Jian2, 3, LIU Wen-bai1, DENG Yi-bing1
(1. College 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 200092, 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 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 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 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 and the direction of the major principal stress has been observed during all stages of shear loading.
Key words: direct shear test; non-circular particle; shear band; non-coaxiality effect; fabric anisotropy
0 引 言
直剪试验是室内研究土体剪切力学性状的重要试验方法，然而，常规的直剪仪只能从宏观上测得土体的抗剪强度指标，不能分析剪切过程试样内部应力应变的变化。最近，刘文白等[1]开发了半膜式细观结构观测直剪仪，利用数字照相变形量测技术，研究了直剪过程中砂土试样应力应变场的分布规律，但该试验没有分析剪切面处颗粒的滑移、旋转等相对运动规律。数值仿真模拟作为对室内试验的有力补充，目前已被广泛用于包括直剪试验在内的土体力学性质的仿真分───────
基金项目：国家自然科学基金项目（50909057，90815008）；上海高校选拔培养优秀青年教师专项基金项目（shs08004）；同济大学岩土及地下工程教育部重点实验室开放基金项目（KLE-TJGE-0904）；上海市教委科研创新项目（11YZ132）
收稿日期：2010–01–11
1558 岩土工程学报 2010年
析。Potts等[2]、Tejchman和Bauer[3]分别采用各向同性弹塑性模型和亚塑性模型，利用有限元法分析了砂土的直剪力学性状，研究了剪切带内的应变局部化问题。有限元法可以较好的分析剪切过程的应力应变关系，但无法从颗粒层面揭示剪切强度和变形的细观力学机理，也无法反映颗粒形状、孔隙率等试样要素变化对模拟结果的影响。离散元法在分析岩土体大变形破坏问题及其宏细观机理研究方面具有独特的优势。刘斯宏和徐永福[4]在铝棒堆积体直剪试验的基础上，利用离散元分析了粒状材料的直剪力学行为。Masson 和Martinez[5]基于PVC材料细观力学特性，对松、密试样的直剪特性进行了离散元对比分析。Cui和O′Sullivan[6]针对均一粒径钢球集合体，利用三维离散元分析了模拟颗粒材料的直剪宏细观力学响应。Wang 等[7]则针对均一粒径玻璃圆球集合体，利用离散–连续耦合分析方法分析了剪切带的演化规律。总结以上文献可知，目前关于直剪试验的仿真模拟大多针对铝棒、钢球等理想代用材料，数值模拟也以纯圆形颗粒居多，而对于实际的砂土，其材料属性不同于理想材料，且具有明显的角粒特征，针对实际砂土的直剪行为开展非圆形颗粒的仿真模拟研究，具有重要的研究意义和应用价值。
本文利用颗粒离散元PFC2D程序，开发了二维椭圆形颗粒单元，并将其用于砂土直剪力学性状的仿真模拟。数值模拟着重研究以下几个关键问题：①宏观应力–剪胀关系及与实际砂土的对比；②颗粒位移、颗粒旋转与剪切带演化；③主应力与主应变增量的非共轴效应；④组构各向异性的演化规律；⑤宏–细观力学量之间的关联。
1 数值试样制备
1.1 椭圆形颗粒单元
关于PFC2D中椭圆形非圆颗粒和角粒形非圆颗粒的开发工作，文献[8，9]已进行了详细的阐述，限于篇幅，本文在此不再重复赘述。非圆颗粒单元的开发需要遵循“面积等效原则”和“质量等效原则”，对于椭圆形颗粒，定义颗粒的长短轴之比为形状指标，用S e来表示。全文中取S e=1.0（纯圆），1.2，1.4，1.6四种情况进行对比分析。图1给出了椭圆形颗粒单元的示意图。
图1 椭圆形颗粒
Fig. 1 Elliptical particles 1.2 重力沉积法制样
数值试样制备仍采用文献[8]中的重力沉积法，具
体步骤为：①首先生成四周墙体，在高30 mm、宽60 mm的二维空间内，按初始孔隙率0.4生成最大粒径0.8 mm、最小粒径0.4 mm，级配服从均匀分布的初始
纯圆颗粒，如图2（a）所示；②删去纯圆颗粒，在同
一位置用等效的椭圆形颗粒代替，颗粒长轴的初始方
向在0°～180°范围内随机定向，如图2（b）所示；
③施加重力，椭圆形颗粒在重力作用下沉积，并达到
初始平衡状态，如图2（c）所示；④在高度20 mm的
位置施加顶部墙体，同时删去顶墙以上的颗粒，如图
2（d）所示；⑤卸除重力，施加上覆固结压力，等待
施加剪切荷载。
图2 数值试样制备（S e=1.4）
Fig. 2 Sample preparation for numerical simulations (S e=1.4)
直剪试验二维数值试样的外观尺寸与常规室内直
剪仪的尺寸相仿（直剪仪高20 mm、直径60 mm）。
为了减少生成颗粒的数量以提高计算效率，数值试样
平均粒径（d50=0.6 mm）略大于福建标准砂平均粒径（d50=0.34 mm），但级配的分布与实际标准砂保持一
致[10]。不同S e试样的初始试样信息汇总于表1。由表
1可知，不同S e试样的初始孔隙率n0基本保持在0.175
左右，但非圆颗粒试样的初始平均接触数(C n)ini明显
大于纯圆颗粒试样，随着S e的增加，试样的(C n)ini逐
渐增加。
表1 数值试样初始信息
Table 1 Initial parameters of numerical samples 长短轴比S e 1.0 1.2 1.4 1.6 颗粒数量3407 3406 3402 3413 固结后高度/mm19.96 19.90 19.86 19.91 初始孔隙率n00.175 0.176 0.177 0.173 初始平均接触
数(C n)ini
3.707
4.124 4.239 4.302 1.3 细观参数设置
数值模拟中，颗粒与颗粒之间、颗粒与墙体之间
的接触特性均采用线性接触模型来表述，需要设定的
细观参数有颗粒/颗粒法向接触刚度和切线接触刚度
第10期 史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究
1559
(k n )p ，(k s )p ，颗粒/墙体法向接触刚度和切向接触刚度(k n )w ，(k s )w ，颗粒/颗粒和颗粒/墙体摩擦系数f p ，f w ，颗粒密度ρs 。本文作者在文献[10]中已完成了福建标准砂二维模拟的细观参数标定工作，本文数值模拟仍然采用文献[10]的结果。考虑到直剪仪墙体可视为刚性体，所以数值模拟中颗粒/墙体接触刚度取为颗粒/颗粒接触刚度的2倍。细观参数的具体数值汇总于表2。
表2 细观参数取值
Table 2 Values of micro-mechanical parameters
细观参数 (k n )p / (N·m -1) (k s )p / (N·m -1) f p (k n )w / (N·m -1) (k s )w / (N·m -1)f w 数值大小 2×108
1×108
0.5
4×108
2×108
0.5
ρs /(kg·m -3)
细观参数 S e =1.0 S e =1.2
S e =1.4
S e =1.6
数值大小
2643 1701 1738 1805
1.4 加荷条件
文献[4]研究表明，由于剪切盒内壁摩擦的影响，
室内直剪试验测得的剪切强度在土体剪胀时，比土体实际的强度要高；而在土体剪缩时，则比土体实际的强度要低。所以，计算剪切面上的剪应力比τ/σN 时应该恰当考虑侧壁摩擦效应的影响。图3给出了数值模拟剪切过程示意图，本文不针对固结应力水平展开重点讨论，所有试样施加同一上覆压力p =200 kPa 。剪切过程中，围成下盒的①、②、③、⑦号墙体以0.01 mm/min 的剪切速率向右移动，围成上盒的④、⑥、⑧号墙体保持固定，顶部⑤号墙体由伺服机制控制保持上覆压力不变，⑤号墙体发生的竖向位移用h 表示；当水平剪切位移u 达到10 mm 时，加荷终止。
图3 数值模拟加荷示意图
Fig. 3 Schematic diagram of shear loading in numerical
simulations
数值试验可以全程监测各个墙体的位移量和墙体上的作用力，剪切面上的剪应力比τ/σN 可由下式计算：
4
65
N 546N N T T N N T T τσ??==++ ， (1) 式中，N 4、T 4，N 5、T 5，N 6、T 6分别表示④、⑤、⑥号墙体上所受的法向作用力和切向作用力，见图3。
2 宏观力学响应分析
2.1 宏观力学特性曲线
数值试验可以得到剪应力比–水平位移(τ/σN –u )和竖向位移–水平位移(h –u )关系曲线，如图4所示。图4（b ）中竖向位移h 为正表示体积剪胀。图4中实心点标示处表示剪应力比峰值点和对应的体积剪胀峰值点。
图4 数值试样宏观力学特性曲线
Fig. 4 Macro-mechanical response of numerical samples
分析图4，可以得到以下几点规律：①二维数值试样在0.175左右的初始孔隙率条件下，表现出典型密实砂样的剪切性状，剪切强度曲线的峰后软化和剪切过程的体积剪胀现象均表现的十分明显；②颗粒形状影响剪切强度，随着S e 的增大，试样的峰值强度和残余强度（u =10 mm ）均明显增加，颗粒形状对初始剪切刚度的影响并不大，不同S e 试样强度曲线的初始斜率基本相同；③对于不同S e 试样，峰值强度大概出现在水平位移2～3 mm 之间，且随着S e 的增大，出现峰值强度时的位移水平逐渐增加；④对于不同S e 试样，除加荷初期发生少量体积剪缩外，剪切过程试样均发生明显的体积剪胀，对应峰值强度点时，不同S e 试样的体积剪胀量较为接近，但当加荷进行至残余强度时，纯圆试样（S e =1.0）的体积剪胀已基本趋于稳定，而非圆试样仍有较明显的剪胀发生。 2.2 应力–剪胀关系
（1）Rowe-Davis 应力–剪胀公式
Rowe 于1962年提出了平面应变条件下颗粒集合体的应力–剪胀关系[11]：
cr ps cr sin sin sin 1sin sin φψ
φφψ+=+ ， (2)
式中，φps 表示平面应变(plain strain)条件下的内摩擦
1560 岩 土 工 程 学 报 2010年 角，φcr 表示临界状态(critical state)内摩擦角，ψ表示剪胀角。 Davis [12]最早提出了直剪试验内摩擦角与平面应变试验内摩擦角之间的换算关系，即
ps
ds ps cos sin tan 1sin sin ψφφψφ=? ， (3)
式中，φds 表示直剪试验（direct shear ）的内摩擦角。
将式（2）代入式（3），消去φps ，得到φds ，φcr ，ψ三者之间的关系
cr ds sin sin tan cos φψ
φψ+= 。 (4)
当试样处于临界状态时，通常认为剪胀角ψ=0，由式（4）可得
()ds cr r tan sin φφ= ， (5) 式中，(φds )r 表示直剪试验的残余内摩擦角。
将式（5）代入式（4），得到
()ds r ds tan sin tan cos φψ
φψ+= 。 (6)
式（6）即为描述砂土直剪应力–剪胀关系的Rowe-Davis 公式。
（2）数值试样的应力–剪胀关系
引用上述Rowe-Davis 公式来分析数值试样的应力–剪胀关系。图4中，数值试验内摩擦角φds 和剪胀角ψ分别可由下式计算：
ds N
tan τ
φσ= ， (7) d tan d h
u
ψ=
。 (8) 计算图4中不同S e 试样的峰值内摩擦角(φds )p 、峰值剪胀角ψp 和残余内摩擦角(φds )r ，并将数值汇总于表3，表3中u p 表示峰值强度对应的水平位移大小。表3中还列出了文献[13，14]中关于实际砂土的分析结果。
图5为公式（6）理论计算值与数值试验实测值之间的对比关系。由图5可知，对于不同S e 试样，数值试验数据点与理论预测值均非常接近，数值试样表现出的应力–剪胀关系符合实际砂土的Rowe-Davis 应
力–剪胀公式。值得注意的是，在图5和表3中，无论是对于实际的砂土还是数值试样，在实测值中，横
坐标(tan(φds )r +sin ψp )/cos ψp 的数值均略大于纵坐标tan(φds )p ，对于非圆数值试样，这一数值上的差异表现的更为明显，其原因与试样到达残余强度时，其剪胀角并非完全为零有关，尤其对于非圆数值试样，由2.1节可知，残余强度时试样仍有较明显的体积剪胀发生。
图5 应力–剪胀关系 Fig. 5 Stress-dilatancy relationship
2.3 剪切强度与初始平均接触数之间的关系 在文献[8]中，本文作者通过砂土双轴试验数值模拟初步探讨了抗剪强度与初始平均接触数之间的关系，本文在此进一步分析直剪试验过程中剪切强度与初始平均接触数之间的关联。图6给出了数值试样峰值内摩擦角(φds )p 、残余内摩擦角(φds )r 与试样初始平均接触数(C n )ini 之间的关系。
图6 剪切强度与初始平均接触数关系
Fig. 6 Shear strength versus initial averaged contact number
表3 应力–剪胀关系实测值
Table 3 Measured values of stress-dilatancy relationship
试样 (φds )p
(φds )r ψp u p /mm tan(φds )p (tan(φds )r +sin ψp )/cos ψp S e = 1.0 21.3 15.9 6.5 2.23 0.39 0.40 S e = 1.2 24.2 20.3 7.1 2.50 0.45 0.50 S e = 1.4 28.8 24.2 7.7 2.75 0.55 0.59 S e = 1.6 31.0 27.5 8.1 3.02 0.60 0.66
Leight Buzzard sand [13]
(e 0=0.565, p =170 kPa) 40.5 29.7 14.7 0.85 0.85
Toyoura sand [14]
(e 0=0.65, p =49 kPa)
40.7 33.0 11.8 0.86 0.87
第10期史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究 1561
分析图6可知，试样的峰值内摩擦角和残余内摩擦角与初始平均接触数之间成良好的线性拟合关系，本文结果和文献[8]的研究结果进一步说明，无论是平面应变状态还是直剪状态，试样所能发挥的宏观剪切强度与试样初始获取的粒间平均接触数的多少密切相关。
3 细观力学响应分析
3.1 颗粒位移与剪切带
剪切带的发生和演化与试样内部应变局部化现象密切相关。数值试验的优势在于可以直观的分析颗粒的水平位移和竖向位移，研究内部剪应变和体应变的局部化现象，揭示剪切带的发生机制与演化过程。对于纯圆试样，文献[4，5]均已作过研究。本文侧重于讨论颗粒形状的影响，分析颗粒形状变化对剪切带厚度的影响。
以S e=1.4试样为例，图7给出了u=10 mm时，不同竖向位置颗粒发生的水平位移和竖向位移分布。分析图7（a）可以清晰的看出，变形的局部化发生在狭长的剪切带内，且剪切带厚度的分布是不对称的，位于下盒的剪切带厚度要明显大于上盒，对于S e=1.4试样，剪切带的厚度约为7.6 mm。分析图7（b）可以发现，较大的竖向位移集中在剪切带内，说明试样宏观的体积剪胀主要由剪切带内颗粒的运动控制。值得注意的是，竖向位移存在负值，表示有些颗粒发生了向下运动，这说明颗粒的运动轨迹是不规则的，但向下运动的最大位移量要比向上运动的最大位移量小很多。
图7 颗粒位移分布图
Fig. 7 Distribution of particle displacement
目前，关于剪切带厚度的研究大多关注于颗粒平均粒径的影响，而关于颗粒形状影响的研究相对缺乏。对于实际砂土，每个颗粒形状各异，很难对形状指标变化对剪切带的影响展开研究，而在数值试验中可以很快捷的生成单一颗粒形状的仿真试样，为这一领域的深入研究提供了条件。在图7（a）中，取竖向间隔为1 mm，计算每1 mm试样高度范围内颗粒水平位移的平均值，并将不同S e试样的结果绘于同一幅图中，如图8所示。分析图8可知，随着S e的增大，试样剪切带厚度逐渐增加，S e=1.0，1.2，1.4，1.6时，相应的剪切带厚度D分别约为6.0，6.5，7.6，8.0 mm，这一结果表明，当颗粒的角粒特征增强时，不仅表现出更大的剪切强度和体积剪胀特征，而且剪切过程中试样内部应变局部化的范围也会扩大。
图8 颗粒平均水平位移曲线
Fig. 8 Curves of averaged particle horizontal displacement 3.2 颗粒旋转
除颗粒位移外，颗粒旋转也是表征剪切带演化过程的重要细观力学量。同样以S e=1.4试样为例，图9（a）给出了u=10 mm时，不同竖向位置颗粒的颗粒旋转量分布。图9（a）中，正值表示颗粒逆时针旋转，负值表示颗粒顺时针旋转。由图9（a）可知，较大的颗粒旋转主要集中在剪切带区域内，而剪切带之外，特别是接近上下端面位置，颗粒旋转量逐渐减小。综合分析不同S e试样情况，在图9（a）中，取竖向间隔为1 mm，计算每1 mm试样高度范围内颗粒旋转量（正、负均取绝对值）的平均值，并将不同S e试样的结果绘于同一图中，如图9（b）所示。分析图9（b）可得以下几点规律：①纯圆试样的颗粒旋转要明显大于非圆试样；②颗粒旋转剪切带内最大，并向剪切带外两侧逐渐递减；③无论是纯圆试样还是非圆试样，剪切带内的颗粒旋转峰值均出现在下盒，其原因可能与3.1节分析中，下盒的剪切带厚度较大有关。
图9 颗粒旋转分布图
Fig. 9 Distribution of particle rotation
3.3 主应力与主应变增量的非共轴效应
（1）平均应力与应变率
1562 岩 土 工 程 学 报 2010年
在PFC 2D 中，可以通过设置在试样内部的量测圈直接测得试样平均应力和应变率。平均应力计算公式为
c 11N k k
ij i j k f l V σ==∑ (i , j =x , y ) ， (9)
式中，V 表示量测圈的体积（面积），N c 表示量测圈内的颗粒接触数总数，f i k 表示接触点k 处接触力矢量在i 方向上的分量，l j k 表示接触点k 处颗粒支矢量（连接两个接触颗粒圆心的矢量）在j 方向上的分量。
应变率计算公式为
d d i ij j x ν
ε
= (i , j =x , y ) ， (10) 式中，d v i 表示量测圈内两相邻颗粒t 时刻在i 方向上的速度矢量差，d x j 表示两相邻颗粒t 时刻颗粒支矢量在j 方向上的分量。
由式（10）可计算得到d t 时步内的应变增量为
d d ij ij t εε= 。 (11) 由式（9）和（11）可以计算得到各应力分量σxx ，
σyy ，σxy 和应变增量分量d εxx ，d εyy ，d εxy 的大小，进一步可计算得到主应力和主应变增量。
本文量测圈的布置共分L 1，L 2，L 3三层，每层布置8个量测圈，量测圈的直径为6 mm ，如图10所示。
图10 量测圈的布置
Fig. 10 Layout of measurement circles
（2）非共轴效应分析
Stroud [15]、Tatsuoka 等[16]均通过砂土的单剪试验(simple shear test)研究证实，伴随剪应力的施加，应变增量主轴的旋转并不与应力主轴的旋转同步，应力主轴的旋转滞后于应变增量主轴，主应力与主应变增量之间存在非共轴效应。数值试验能否有效模拟非共轴效应是本节研究重点。
Jewell [17]试验研究表明，直剪试验中，剪切带内砂土的变形近似于单剪试验中砂土的变形模式。本文数值试验中，布置在L 2层的第9至第16号量测圈基本覆盖了剪切带区域（如图10所示）。取第9至第16号量测圈各应力分量和应变增量分量的平均值来分析加荷过程中应力主轴和应变增量主轴的旋转情况。
用α，β分别表示大主应力方向和大主应变增量方向与水平方向的夹角，则α，β分别可由下式计算
21arctan(
)2
xy
xx yy
σασσ=? ， (12) 2d 1
arctan(
)2
d d xy
xx yy
εβεε=? 。 (13)
图11给出了加荷过程中α，β的变化过程。分析图11可得几点规律：①固结状态时，竖向应力大于水平应力，加荷开始后，大主应力方向从90°左右开始偏转，峰值强度后，大主应力的方向逐渐趋于稳定；②加荷开始后，应变增量主轴的偏转要比应力主轴的偏转快的多，在很小的剪切位移水平下，大主应变增量方向即发生了很大的偏转，并趋于稳定，这一变化规律与文献[15，16]的试验结果是一致的，说明数值试验能够较好的模拟非共轴效应；③非共轴效应在峰值强度前表现的较为明显，峰值强度后，随着塑性变形的增加，非共轴效应逐渐减小；④颗粒形状对大主应变增量方向角β的影响要小于对大主应力方向角α的影响，对于不同S e 试样，峰值强度后，α角大致在40°～50°区间范围内，β角在40°～45°区间范围内。
图11 大主应力与大主应变增量方向角
Fig. 11 Inclination of major principal stress and major principal
strain increment along horizontal plane
3.4 组构各向异性演化规律
对于散粒体材料，加荷过程试样表现出的宏观强度和变形特性与其内部细观组构的演化规律密切相关。Rothenburg 和Bathurst [18]提出可用傅里叶函数来近似描述单调剪切荷载作用下粒间接触法向、粒间法向接触力和切向接触力各向异性的演化规律，其数学表达式分别为
()()a 11cos 22π
E a θθθ??=+??? ， (14)
第10期 史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究 1563
Y
()()n 0n n 1cos 2f f a θθθ=+????? ， (15) ()()t 0t t sin 2f f a θθθ=?? 。 (16)
式中 f 0是相对于所有接触的平均法向接触力；θa ，θn ，θt 分别是接触法向各向异性、法向接触力各向异性和切向接触力各向异性的主方向；a ，a n ，a t 是傅立叶系数，它们的数值大小反映对应细观组构参量的各向异性程度。
以S e =1.4试样为例，分析加荷过程中细观组构各向异性的演化规律。组构各向异性的分析着重针对剪切带区域，数值试验加荷过程中，对L 2层内（见图10）颗粒之间的接触法向、法向接触力和切向接触力的演化进行统计分析。图12（a ）～（c ）分别给出了加荷前、峰值强度和残余强度时粒间接触法向、法向
接触力、切向接触力的各向异性分布玫瑰图。玫瑰图的绘制每10°分一个区间，统计接触法向时，M 表示接触点总数、?M (θ)表示接触法向落入该角度区间的接触点个数；统计接触力时，f n (θ)，f t (θ)分别取接触法向落入该角度区间内所有接触点处法向接触力和切向接触力的平均值。图12中虚线为公式（14）～（16）傅里叶函数的拟合结果，图中箭头方向表示拟合得到的组构参量各向异性的主方向。
分析图12可知，接触法向的拟合曲线呈“椭圆状”图形，在整个剪切过程中，接触法向的主方向都偏于竖直方向，但接触法向分布于各个方向的各向异性程度并不大。法向接触力和切向接触力的拟合曲线分别呈“花生状”和“花瓣状”图形，加荷前，数值试样由重力沉积生成并施加了竖向固结压力，由图12（a ）
接触力系
1564 岩土工程学报 2010年
可知，试样内的主接触力系分布于竖直方向，粒间法向接触力和切向接触力的主方向都偏于竖直方向，初始接触力的分布并非各向同性；剪切荷载施加后，随着主接触力系逐渐向斜45°方向倾斜，法向接触力和切向接触力的主方向随之倾斜，且各向异性的程度明显增强；剪切至残余强度时，法向接触力和切向接触力的各向异性程度减小，且粒间平均法向接触力f0减小。
图13给出了组构各向异性参数的变化规律。图13(b)中实线所示为3.3节计算得到的大主应力方向角α。分析图13可得几点宏－细观力学量之间的关联：①与a n相比，整个加荷过程中，a和a t的数值和变化量均较小，说明试样宏观的剪切强度主要受控于粒间法向接触力的分布及其各向异性演化；②a n，a t出现峰值时的位移水平与试样呈现峰值强度时的位移水平十分接近，说明试样宏观强度的变化与内部接触力各向异性的演化保持同步；③对比分析θn和α的变化规律，剪切带内大主应力的偏转方向与法向接触力各向异性的主方向保持了良好的一致性，进一步说明宏观应力与粒间接触力之间存在内在关联。
图13 组构各向异性参数的变化（S e=1.4） Fig. 13 Variation of fabric anisotropy parameters（S e=1.4）
4 结 论
本文基于PFC2D非圆颗粒单元的二次开发，对砂土直剪力学过程进行了非圆颗粒模拟，得到的主要结论有：
（1）数值试样表现出密实砂样的直剪力学性状，强度的峰后软化和体积剪胀均十分明显，加荷过程表现出的应力–剪胀关系符合实际砂土的Rowe-Davis 应力–剪胀关系。颗粒形状影响峰值强度和残余强度，颗粒形状的影响在细观机理上与试样的初始平均接触数有关。
（2）颗粒的水平位移分布清晰的表明，变形的局部化发生在狭长的剪切带内，较大的竖向位移和颗粒旋转均集中在剪切带内，试样宏观的体积剪胀主要由剪切带内颗粒的运动控制。颗粒形状影响剪切带厚度，S e=1.0，1.2，1.4，1.6时，剪切带厚度D分别约为6.0，6.5，7.6，8.0 mm，说明颗粒的角粒特征增强时，试样内部应变局部化的范围将扩大。
（3）数值试验能够较好的模拟主应力与主应变增量的非共轴效应，非共轴效应在峰值强度前表现的较为明显，峰值强度后逐渐减小；非共轴效应中，应力主轴的旋转滞后于应变增量主轴。颗粒形状对大主应变增量方向角β的影响要小于对大主应力方向角α的影响。
（4）剪切过程中，接触法向的主方向都偏于竖直方向，但各向异性的程度并不大。加荷前，粒间接触力的主方向偏于竖直方向；剪切荷载施加后，伴随主接触力系逐渐向斜45°方向倾斜，法向接触力和切向接触力的主方向随之倾斜，且各向异性的程度增强；剪切至残余强度时，法向接触力和切向接触力的各向异性程度减小，且粒间平均法向接触力f0减小。
（5）各向异性参数的定量研究表明，试样宏观的剪切强度主要受控于粒间法向接触力的分布及其各向异性演化，且宏观强度的变化与内部接触力各向异性的演化保持同步；整个加荷过程中，剪切带内大主应力的偏转方向与法向接触力各向异性的主方向保持了良好的一致性。
参考文献：
[1] 刘文白, 张辉, 邓一兵. 基于DPDM技术的砂土直剪试
验剪切过程的应力场分析[J]. 中国水运, 2008, 8(7): 235–
239. (LIU Wen-bai, ZHANG Hui, DENG Yi-bing. Analysis
of sand stress field during direct shearing by DPDM technique[J]. China Water Transport, 2008, 8(7): 235–239.
(in Chinese))
[2] POTTS D M, DOUNIAS G T, V AUGHAN P R. Finite element
analysis of the direct shear box test[J]. Géotechnique, 1987,
37(1): 11–23.
[3] TEJCHMAN J, BAUER E. FE-simulations of a direct and a
true simple shear test within a polar hypoplasticity[J].
Computers and Geotechnics, 2005, 32(1): 1–16.
[4] 刘斯宏, 徐永福. 粒状体直剪试验的数值模拟与微观考察
[J]. 岩石力学与工程学报, 2001, 20(3): 288–292. (LIU Si-hong, XU Yong-fu. Numerical simulation for a direct box
shear test on granular material and microscopic
第10期史旦达，等. 砂土直剪力学性状的非圆颗粒模拟与宏细观机理研究1565
consideration[J]. Chinese Journal of Rock and Mechanics and Engineering, 2001, 20(3): 288–292. (in Chinese))
[5] MASSON S, MARTINEZ J. Micromechanical analysis of the
shear behavior of a granular material[J]. Journal of Engineering Mechanics, 2001, 127(10): 1007–1016.
[6] CUI L, O′SULLIV AN C. Exploring the macro- and
micro-scale response of an idealized granular material in the direct shear apparatus[J]. Géotechnique, 2006, 56(7): 455–
468.
[7] WANG J, DOVE J E, GUTIERREZ M S. Discrete-continuum
analysis of shear banding in the direct shear test[J].
Géotechnique, 2007, 57(6): 513–526.
[8] 史旦达, 周健, 刘文白, 等. 砂土单调剪切特性的非圆颗
粒模拟[J]. 岩土工程学报, 2008, 30(9): 1361–1366. (SHI Dan-da, ZHOU Jian, LIU Wen-bai, et al. Numerical simulation for behaviors of sand with non-circular particles under monotonic shear loading[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(9): 1361–1366. (in Chinese))
[9] 史旦达, 周健, 刘文白, 等. 循环荷载作用下砂土液化特
性的非圆颗粒数值模拟[J]. 水利学报, 2008, 39(9): 1074–
1082. (SHI Dan-da, ZHOU Jian, LIU Wen-bai, et al.
Numerical simulation of liquefaction of sands with non-circular particles under cyclic loading[J]. Chinese Journal of Hydraulic Engineering, 2008, 39(9): 1074–1082.
(in Chinese))
[10] 周健, 史旦达, 贾敏才, 等. 砂土单调剪切力学性状的
颗粒流模拟[J]. 同济大学学报, 2007, 35(10): 1299–1304.
(ZHOU Jian, SHI Dan-da, JIA Min-cai, et al. Numerical simulation of mechanical response on sand under monotonic loading by particle flow code[J]. Journal of Tongji University, 2007, 35(10): 1299–1304. (in Chinese))
[11] ROWE P W. The stress dilatancy relation for static
equilibrium of an assembly of particles in contact[C]// Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1962: 500–527. [12] DA VIS E H. Theories of plasticity and the failure of soil
masses[M]. In Soil Mechanics: Selected topics(. Lee K I, ed), London: Butterworth, 1968: 341–380.
[13] LINGS M L, DIETZ M S. An improved direct shear
apparatus for sand[J]. Géotechnique, 2004, 54(4): 245–256. [14] SHIBUYA S, MITACHI T, TAMATE S. Interpretation of
direct shear box testing of sands as quasi-simple shear[J].
Géotechnique, 1997, 47(4): 769–790.
[15] STROUD M A. The behavior of sand at low stress levels in
the simple shear apparatus[D]. Cambridge: University of Cambridge, 1971.
[16] TATSUOKA F, PRADHAM T B S, Horii N. Discussion on
direct shear tests on reinforced sand[J]. Géotechnique, 1988, 38(1): 148–153.
[17] JEWELL R A. Direct shear tests on sand[J]. Géotechnique,
1989, 39(2): 309–322.
[18] ROTHENBURG L, BATHURST R J. Analytical study of
induced anisotropy in idealized granular materials[J].
Géotechnique, 1989, 39(4): 601–614.