粗粒土细观组构分析的影响因素研究

, ,2,,,

(1.浙江工业大学 岩土工程研究所,浙江 杭州 310023;2.浙江同济科技职业学院，浙江 杭州 311231)

土的力学性质的复杂性，包括应力路径相关性、各向异性和时间效应等，是与其细观组构特征密切相关的，因此对土的细观组构的形成和演化规律进行研究是认识土的力学行为机理的重要途径.在土的力学特性和细观组构联合研究方面，Oda等[1-3]通过研究发现散粒材料的初始各向异性主要表现在接触应力的各向异性分布、孔隙空间的优势方向和非圆球体颗粒的优势方向3个方面，基于图像分析法提出了使用二阶组构张量来描述诱发各向异性，并发现了颗粒材料的孔隙分布与外力加载的方向具有一致性[4].Zhang等[5]研究了在循环剪切荷载条件下二维光弹颗粒结构中应力发展，利用光弹技术检测到接触网格与接触力，并且研究了组构张量、应力张量以及应力在法向和切向的分布规律.杨仲轩等[6-9]利用三轴试验研究砂土的各向异性并给出了定量方法，发现砂土力学特性受到材料各向异性和应力各向异性的耦合影响；周健等[10-14]利用了离散元程序对颗粒材料的力学变形特性进行了数值模拟研究，并且在颗粒流理论的基础上，还引入了不同颗粒的接触本构模型，探索适用于砂土和黏性土等不同颗粒材料的颗粒模型，这为深入分析土的细观力学机理提供了可行的途径和研究基础.Hu等[15]采用PFC2D模拟了循环荷载作用下砂土的诱发各向异性，研究发现试样的初始不等向固结方式、颗粒材料的法向与切向刚度比及循环荷载峰值等对砂土的诱发各向异性都会产生一定程度的影响.

研究表明：土的组构各向异性是颗粒材料最重要的组构特征，而采用离散元细观数值模拟的方法可以对土的各向异性组构和受力特性进行联合分析，是有效的方法.但是在现有的离散元细观组构研究中，采用的试样颗粒级配、试样尺寸和数值制样方法等都不统一，这对各种组构分析研究结果之间的相互参考性提出了质疑.为此，通过采用离散元数值模拟，对不同的颗粒级配、不同的试样尺寸、不同的制样方法以及不同的加载方式等各种因素下粒状土的细观参数变化及各向异性组构特征等进行分析，比较各种不同条件下土样的颗粒接触分布、平均配位数和组构偏量的变化规律，并且结合试样的应力

应变特征分析了应用离散元方法研究粗粒土各向异性组构特征时需考虑的前期影响因素，可为室内试验和理论研究提供参考.

1 数值模拟试验设计及细观组构参数

1.1 颗粒材料参数及接触模型选择

为了研究制样方法和试样的物理及几何特征对试样的组构及初始接触的影响，采用颗粒流程序PFC2D对数值试样进行了双轴等向压缩和双轴剪切试验.

在数值试验中合理选择颗粒间的接触模型及模型参数对正确分析颗粒材料组构特征是至关重要的.颗粒流程序中提供了线性接触刚度模型、粘结模型以及滑动等模型，其中线性接触刚度模型能较好的反映颗粒材料在接触点处的挤压与摩擦作用，是一种可用于模拟散体颗粒单元间作用的简单且有效的模型[16]，为此，采用线性接触刚度模型来模拟砂土的力学性质.

在线性接触刚度模型中，法向接触力与相对位移的关系为

(1)

式中：Kn为颗粒的法向接触刚度；Un为接触的法向相对位移；ni为接触单位法向向量.切向接触力与相对位移的关系为

(2)

(3)

由于颗粒材料的细观模拟参数很难由试验精确确定，为此，参考文献[15-18]中的取值方法，经过反复试算调整颗粒间法、切向刚度比.选取土样的细观模拟参数如表1所示.

表1 PFC2D数值模拟基本参数Table 1 Parameters used in PFC2D simulations

1.2 试样制备及数值模拟方案

为了研究试样的颗粒级配、制样方法、试样尺寸和加载方式等条件的影响，采用的具体数值模拟方案如表2所示.其中各个影响因素的控制方法如下：

1) 颗粒级配.采用均匀分布和高斯分布(不均匀分布)方法生成了两种颗粒试样，分别表示为均匀试样(U)和不均匀试样(G)，其中前者控制颗粒粒径分布在较小的区间内(0.3～0.399 mm)；后者的颗粒粒径的分布范围较大(0.05～0.92 6mm).

2) 几何尺寸.采用相同形状的2种不同尺寸，分别为80 mm×30 mm(B) ，40 mm×15 mm(L).

3) 制样方法.以孔隙比为标准，采用压实法(P)和振动法(Z)2种不同方式制备试样.前者是通过上边界墙体以一定速度向下运动压实试样，当边界应力达到一定值时停止运动，稳定一定时间后再继续压实，如此循环若干次；后者是控制4个边界以相同的速度向同一个方向运动，待其位移达到一定的值时，改变其速度的方向，达到原位置时继续反向以相同速度压实试样，如此反复若干次.

4) 加载方式.采用双轴等向压缩(C)和双轴剪切(S)两种加荷方式.前者控制试样的上下边界和左右边界以相同的应变率向内移动压缩试样，即压缩的速率与运动方向上的试样尺寸成正比，实现试样的等向压缩；后者是通过控制上下边界以一定的速率相对移动进行压缩，而左右边界同时以一定的速率相对扩张，通过伺服控制以维持均应力稳定，实现试样的纯剪切加载.

表2 试验方案1)Table 2 Test Program

注：1)b×h为试样尺寸，其中b为试样宽，h为试样高；试样编号中的字母为本节中试验研究因素的符号.

1.3 细观组构参数

颗粒流方法提供了一种从细观层面上观察和分析散体颗粒材料在荷载作用下的力学变形规律的研究手段.为了分析砂样的接触及其细观组构的变化规律，采用如下统计方法：

1) 接触力网络又可称颗粒间的接触力链，可以直观地观察试样内部颗粒间接触力的分布以及接触力的优势方向.

2) 接触力方向分布图可以反映颗粒间接触力在不同方向上概率分布情况，也可以反映在不同方向上平均接触力的大小.

3) 平均接触数Z(又称配位数)指散体颗粒间接触的密度，反映了颗粒之间的连通性与结构性，可表示为

(4)

式中：Nc为颗粒间的总接触数；N为颗粒数.

4) 组构张量.组构张量是反映散体颗粒及接触的几何分布形式和颗粒间相互作用的综合特征，可用于描述散体颗粒材料的组构特征[19]，其表达式为

(5)

式中：Nc为颗粒间总接触数；ni为接触法向的单位矢量.

5)组构偏量指组构张量特征值φ1和φ2之差，即(φ1-φ2).对各向同性散体材料时，组构偏量值为0，各向异性越明显，组构偏量值也越大[20].

2 PFC2D数值模拟验证

为了验证离散元方法分析颗粒材料组构及接触特征的合理性，对Majmudar等[[21]进行的等向压缩光弹试验进行数值模拟.采用线性接触模型，并根据光弹试验模型确定PFC2D模型中试样几何尺寸、颗粒粒径及细观参数，具体如表3所示.

表3 PFC2D数值模拟基本参数Table 3 Parameters used in PFC2D simulations

通过统计等向压缩条件下，颗粒法向和切向接触力大小的概率分布如图1所示，横坐标为无量纲化的法向或切向接触力，即法向/切向接触力与平均法向/平均切向接触力的比值(F/Fmean)；纵坐标表示法向或切向接触力的概率密度P(F/Fmean).从图1中可以看出：光弹试验[21]和数值试验中颗粒间的法向接触力大小的概率密度在平均法向接触力附近达到峰值，大于平均法向接触力时，法向接触力的概率密度分布逐渐减小，而切向接触力大小的概率密度分布一直随着切向接触力的增大而减小.

图2给出了试样在等向压缩条件下大于平均接触力的接触方向分布.对比图2(a,b)的结果可以发现：在等向压缩条件下，大于平均接触力的接触，其方向分布规律也具有明显的相似性.

图1 颗粒法向和切向接触力大小的概率密度分布Fig.1 Probability distribution of normal and tangential contact force

图2 大于平均接触力的接触方向分布Fig.2 Direction distribution of contact force

综上分析可知：两类试验得到的接触力大小的概率密度分布规律是一致的，且接触方向的分布规律也是相似的.因此，利用PFC2D模拟砂样的细观组构及结构变化过程是可行且可靠的.

3 数值试验结果分析与讨论

为了分析各种条件对PFC2D方法研究砂土组构各向异性特征的影响，开展了表2中16个试样的数值加载试验，通过统计加载后试样内部的平均接触数和组构偏量的变化等细观参数，明确各种因素对模拟结果的影响.

3.1 颗粒级配的影响

图3给出了2种颗粒级配试样在等向压缩和剪切加载条件下的应力—应变关系，从图3可以看出：等向压缩条件下的应力应变曲线基本相同，在剪切加载条件下的偏应力轴应变关系曲线的规律一致，在后期略有区别，表明不同颗粒级配对试样的应力应变关系的影响不大.

图3 不同颗粒级配试样的加载曲线Fig.3 The stress-strain curves of specimens with different gradations

图4为2种不同颗粒级配试样采用压实法制样后加载前的接触力网络，图4中线段粗细表示接触力的大小，把大于平均接触力接触定义为主要接触，通过对比图4(a,b)中在加载前的接触力网络可以发现不同颗粒级配试样的接触力网络略有不同，前者主要接触力较大，分布不均匀，而后者主要接触力相对较小，分布比较均匀.原因可能是不均匀分布粒径分布范围较大，起骨架作用的颗粒数较少，而均匀分布恰好相反，起到主要骨架作用的颗粒数较多.

图5给出了2种不同颗粒级配试样的接触力方向分布图，图5中黑色线表示接触力在各方向上的平均接触力，虚线表示各个方向上出现接触数的概率.从图5中可以看出：均匀分布试样在不同方向上的接触力分布规律近似圆形，也可以说明均匀分布试样的接触力在不同方向较为均匀，接触力分布的优势方向不明显.表4中均匀分布试样的平均接触数大于不均匀分布，同时从表5中对比2种不同颗粒级配的试样在加载前后组构偏量的差值可以发现，不均匀分布试样的组构偏量变化值大于均匀分布的试样.以上分析尽管从应力-应变分析来看，试样颗粒级配的影响不明显，但是从组构分析的角度看，颗粒级配对试样的初始细观组构及演化有一定的影响.

图5 不同颗粒级配试样法向接触力方向分布Fig.5 Distribution of normal contact force with different gradation

试样大小加载方式平均接触数压实制备加载前加载后差值振动制备加载前加载后差值不均匀大试样等向压缩双轴剪切3．5924．1440．5523．6000．0083．5503．9980．4483．410-0．140不均匀小试样等向压缩双轴剪切3．6514．1000．4493．6520．0013．5234．0640．5413．5430．020均匀大试样等向压缩双轴剪切5．2075．2100．0035．2260．0195．2435．2730．0305．242-0．001均匀小试样等向压缩双轴剪切5．1625．166-0．0045．152-0．0145．3185．315-0．0035．276-0．042

表5 不同试样在加载前后对组构偏量的影响Table 5 Fabric deviator of different samples before and after loading

3.2 试样几何尺寸的影响分析

图6给出了2种不同尺寸的试样在加载前的接触力链，可以发现二者的差异并不明显.从表4中分析2种不同几何尺寸试样在加载前后平均接触数的变化，发现加载前大尺寸试样的平均接触数总是略大于小尺寸试样的平均接触数，但相差不大.加载后，2种不同几何尺寸的试样的平均接触数相差很小.同时，从表5中对比分析不同几何尺寸对组构偏量的影响，可以发现试样几何尺寸的不同对组构偏量的影响不大.从图7可以看出：不同几何尺寸试样在加载过程中的偏应力-轴应变的变化规律虽然具有相似性，但数值上有些区别.由此可知试样尺寸的大小对试样的应力—应变关系分析有一定的影响，但从试样细观组构研究的角度来看，试样几何尺寸的影响不大.

图6 不同尺寸试样接触力链Fig.6 Contact force chain of different sample size

图7 不同尺寸试样的加载曲线Fig.7 The stress-strain curves of specimens with different sizes

3.3 试样制备方法的影响分析

图8给出了压实法与振动法两种制样方法获得试样的接触力链.从图8中可以发现：试样经压实后接触力在各方向上分布较为均匀，经振动后接触力分布在竖轴方向上形成了明显方向优势.从图9中可以看出：2种不同制备方法下试样法向接触力方向分布具有明显区别，压实制备下法向接触力在各个方向上分布近似为圆形，即分布较为均匀，而振动制备下试样在90°和270°附近的平均法向接触力较大，在0°和180°附近的平均法向接触力较小，说明振动制备过程引起了试样颗粒接触力链的重新调整，接触力在竖轴方向上集中发展.

图8 不同制备方法下试样法向接触力链Fig.8 Contact force chain with different preparing methods

对比分析表4中不同制备方法下试样颗粒的平均接触数的变化，可以发现压实制备下与振动制备下试样的平均接触数比较接近.同时从表5中不同制备方法生成的试样在加载前后组构偏量的变化看出，加载前振动制备试样的组构偏量大于压实制备试样的组构偏量，双轴等向压缩加载后，组构偏量都减小，双轴剪切加载后，组构偏量都增大，且最终振动制备的试样组构偏量越大，试样的各向异性越显著.说明不同的制备试样的方法对材料初始组构及结构演化均有较大影响.从图10中可以看出：压实与振动制备试样在双轴等向压缩加载条件下的应力—应变关系曲线的发展规律具有明显的一致性，但数值有明显差别，因此，无论从研究土的应力—应变变化的角度还是从研究土的细观组构演化的角度讲，选择合适的制样方法都是重要的.

图10 不同制备方法下试样的加载曲线Fig.10 Stress-train curves of different samples with different preparing methods

3.4 试样的加载方式的影响分析

图11给出了试样经等向压缩和双轴剪切两种加载方式下的接触力链，可以看出试样在等向压缩后的接触力方向分布更加均匀，而在双轴剪切作用下接触力方向分布在竖轴方向上形成了优势方向.图12(a,c,e)和图12(b,d,f)分别表示压实与振动制备试样在双轴等向压缩和双轴剪切加载前后的法向接触力方向分布，可以看出在等向压缩作用下，试样的法向接触力方向分布越均匀，试样的颗粒的排列更紧密，在双轴剪切作用下，颗粒间法向接触力在竖轴方向上出现明显优势方向.法向接触力的分布规律与试样接触力网络的分布规律是一致的.

图11 不同制备方法下试样法向接触力链Fig.11 Contact force chain with different preparing methods

图12 不同加载方式下试样法向接触力方向分布Fig.12 Distribution of contact force direction with different preparing methods

从表4中试样加载前后接触平均数变化来看，双轴等向压缩加载后试样的平均接触数要比双轴剪切作用后试样的平均接触数大.原因在于等向压缩作用后，试样中起到骨架作用的颗粒数变多.同时，分析表5中不同加载方式作用前后试样的组构偏量的变化可以发现，在双轴等向压缩加载后，试样的组构偏量变小，而在双轴剪切加载后，试样的组构偏量增大.说明加载方式对试样的细观组构有较大的影响.

综上分析可知：不同的加载方式决定了试样的细观组构演化特征，尽管土样的颗粒级配和土样的制备方法不同会对试样的初始组构有一定的影响，但是，在相同的加载路径下，试样的组构演化规律是一致的，而且试样的应力-应变发展特征也是一致的.其中双轴等向压缩作用使得试样的组构偏量减小，颗粒间的接触力分布较为均匀，而双轴剪切作用使得试样组构偏量的增大，颗粒间的接触力在加载方向上出现优势方向，使得试样在加载后各向异性显著.

4 结 论

在验证离散元数值模拟研究粗粒土细观组构合理性的基础上，研究试样颗粒级配、试样尺寸、制样方式和加载方式等条件对进行砂土细观组构分析的影响，结果表明：试样颗粒级配对相同加载下土样的应力应变响应影响不大，但是对土样的细观组构分析有一些影响，主要表现在颗粒的接触数和组构偏量取值有所不同.因此，分析不同土类在特定应力路径下的细观组构特征时，应注意土样的来源、选取有代表性的颗粒级配；试样尺寸对土样的应力-应变响应有一定影响，但是对分析试样细观组构及其演化影响不明显.因此，从砂土细观组构分析的角度，采用离散元分析可以采用较小的试样尺寸来提高数值模拟的效率；试样制备方法对试样的初始组构及加载后的组构演化有较显著的影响.因此当以土样的细观组构及各向异性演化为研究内容时，须考虑土体的生成特征及沉积环境的影响，采取相应的方法制备土样；加载方式对土样的组构特征和各向异性演化有重要影响，因此为了分析土样的组构特征和演化，须考虑土样的应力历史及加载路径的影响；采用等向压缩的方法可以消除土样由于制样不同形成的各向异性的影响，因此研究正常固结土和超固结土在不同加载路径下的组构变化时，可以通过先施加等向压力再施加固结压力的方式对试样进行预压.

参考文献：

[1] ODA M. Initial fabrics and their relations to mechanical properties of granular material[J].Soils and foundations, tokyo,1972，12(1):17-36.

[2] ODA M. The mechanism of fabric changes during compressional deformation of sand[J].Soils and foundations,tokyo,1972,12(2):1-18.

[3] ODA M. Deformation mechanism of sand in triaxial compression tests[J]. Soils and foundations,tokyo,1972,12(4):45-63.

[4] ODA M, KAZAMA H, KONISHI J. Effects of induced anisotropy on the development of shear bands in granular materials[J]. Mechanics of materials,1998,28:103-111.

[5] ZHANG J, MAJMUDAR T S, TORDESILLAS A,et al. Statistical properties of a 2D granular material subjected to cyclic shear[J].Granular matter,2009,12(2):159-172.

[6] 杨仲轩,李相崧,明海燕.砂土各向异性和不排水剪切特性研究[J].深圳大学学报(理工版),2009(2):158-163.

[7] YAMADA Y, ISHIHARA K. Anisotropic deformation characteristics of sand under three-dimentional stress conditions[J].Soils and foundation,1979,19(2):79-94.

[8] 陈春霖,张惠明.饱和砂土三轴试验中的若干问题[J].岩土工程学报,2000(6):659-663.

[9] 刘延志,胡敏云,沈映,等.杭州市粉土的各向异性室内试验研究[J].浙江工业大学学报,2012,40(2):188-192.

[10] 周健,池毓蔚,池永,等.砂土双轴试验的颗粒流模拟[J].岩土工程学报,2000(6):701-704.

[11] 蒋明镜,付昌,刘静德,等.各向异性结构性砂土离散元分析[J].岩土力学,2015(S1):577-584.

[12] 蒋明镜,孙渝刚.人工胶结砂土力学特性的离散元模拟[J].岩土力学,2011(6):1849-1856.

[13] 李学丰,孔亮,袁琪,等.宏细观结合各向异性砂土的变形特性模拟[J].同济大学学报(自然科学版),2016(2):173-179.

[14] 周凤玺,赖远明.冻结砂土力学性质的离散元模拟[J].岩土力学,2010(12):4016-4020.

[15] HU M, O’SULLIVAN C, JARDINE R R, et al. Stress-induced anisotropy in sand under cyclic loading[J]. Ganular matter,2010,12:469-476.

[16] 史旦达.单调与循环加荷条件下砂土力学性质细观模拟[D].上海：同济大学,2007.

[17] 沈一帆,胡敏云,计国贤.基坑水泥搅拌桩桩身水平位移影响因素的颗粒流数值模拟[J].浙江工业大学学报,2011,39(3):273-278.

[18] 孙新坡,何思明,于忆骅.基于离散元法崩塌体动力破碎分析[J].浙江工业大学学报,2015,43(4):464-467.

[19] ODA M. Fabrics and their effects on the deformaiton behaviours of sand[J]. Department of foundation engineering, saitama university,1976，15(4):75.

[20] THORNTON C. Numerical simulations of deviatoric shear deformation of granular media[J]. Geotechnique，2000,47(2):319-329.

[21] MAJMUDAR T S, BEHRINGER R P. Contact force measurements and stress-induced anisotropy in granular materials[J] . Nature,2005，435：1079.