DOI: 10.11817/j.issn.1672-7207.2016.05.032

堆石体在真三轴应力状态下的非共轴性与剪胀特性

马刚^{1, 2}，刘嘉英^{1, 2}，常晓林^{1, 2}，周伟^{1, 2}

 (1. 武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉，430072；

2. 武汉大学 水工岩石力学教育部重点实验室，湖北 武汉，430072)

摘要：基于连续-离散耦合分析方法，发展随机颗粒不连续变形分析方法。采用不规则多面体模拟实际堆石颗粒，制备初始各向同性的堆石料数值试样，进行等平均静水压力p、等中主应力系数b应力路径的真三轴数值试验。分析堆石料在真三轴应力状态下的剪胀特性，以及应变增量与应力增量的非共轴性。研究结果表明：刚开始加载时，应变增量与应力增量共轴，随着加载的进行，应变增量的方向发生偏转，偏转量与中主应力系数b有关，在三轴压缩应力路径(b=0)和三轴拉伸应力路径(b=1)时，几乎不发生偏转；对比Rowe，Roscoe，改进的Roscoe以及Lagioia剪胀模型，改进的Roscoe剪胀模型预测的剪胀曲线与数值试验结果较吻合，在改进的Roscoe剪胀模型上引入1个角隅函数，以反映中主应力对剪胀特性的影响。

关键词：堆石料；真三轴应力状态；连续-离散耦合分析；非共轴性；剪胀特性

中图分类号：TV641          文献标志码：A         文章编号：1672-7207(2016)05-1697-11

Non-coaxiality and dilatancy of rockfill materials under true triaxial stress condition

MA Gang^{1, 2}, LIU Jiaying^{1, 2}, CHANG Xiaolin^{1, 2}, ZHOU Wei^{1, 2}

 (1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University,

Wuhan 430072, China;

2. Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of Ministry of Education,

Wuhan University, Wuhan 430072, China)

Abstract: Based on the combined finite-discrete element method, the stochastic granular discontinuous deformation method was developed. Numerical samples composed of irregular polyhedras were compressed to the isotropical state and then subjected to true triaxial loading with the constant mean stress and the constant intermediate principal stress ratio. The dilatancy of rockfill materials was analyzed, as well as the non-coaxiality of the stress increment and the strain increment. The results show that the strain increment and stress increment are coaxial at first, and then the direction of the strain increment has some deflection which is related to the intermediate principal stress ratio. No deflection is shown in the triaxial compression and extension path. The comparison of the Rowe’s, Roscoe’s, modified Roscoe’s and Lagioia’s models indicates that the dilatancy curve predicted by modified Roscoe’s model fits well with the simulated results. A ridge function is introduced to describe the influence of the intermediate principal stress ratio on the dilatancy.

Key words: rockfill materials; true triaxial stress condition; combined FEM/DEM; non-coaxiality; dilatancy

随着我国社会经济的快速发展以及西部水电开发进程的加快，西南地区正在或即将建设一批调节性能好的高堆石坝，坝高大多在200 m以上，有些甚至超过300 m，属超高堆石坝工程。超高堆石坝的建设对堆石料力学特性的研究提出了更高的要求。目前，常规三轴试验仍是研究堆石料力学特性的主要手段，现有的本构模型也大多基于堆石料的常规三轴试验结果建立。常规三轴试验只能模拟轴对称应力状态，无法考虑中主应力的影响，而实际工程中坝体的各个部位都处于三向不等的应力状态，因此，有必要进行堆石料的真三轴试验，研究堆石料在真三轴应力状态下的变形和强度特性。国内外很多学者对砂土进行了真三轴试验^[1]，但是堆石料粒径比砂土的粒径要大得多，颗粒力学性质也与砂土存在较大差异，因此不能简单地将砂土的真三轴试验结果外推到堆石料中。由于缺乏大尺寸的真三轴试验设备，目前堆石料的真三轴试验成果很少。施维成^[2]采用中型真三轴仪TSW-40对粗粒土进行了一系列的真三轴试验，研究了中主应力对粗粒土强度和变形特性的影响， 并在此基础上提出了一个粗粒土的三维破坏准则及双屈服面三维弹塑性模型。施维成等^[3]采用小型真三轴仪ZSY-1对砾石料进行了等、等b试验，研究了中主应力对砾石料变形和强度的影响。XIAO等^[4]整理了施维成^[2]的粗粒土真三轴试验结果，研究了围压、中主应力系数、试样密度对粗粒土剪胀特性的影响，采用最小二乘法拟合试验数据，提出了一个新的剪胀模型。除了试验研究外，一些学者采用离散单元法进行了颗粒集合体在真三轴加载情况下的细观数值模拟^[5-9]，研究了颗粒集合体的细观组构和接触力在加载过程中的演化特性，从颗粒层面解释了复杂力学特性的细观机理。基于连续-离散耦合分析方法(combined finite-discrete element method，FDEM)^[10]，周伟等^[11]发展了随机颗粒不连续变形分析方法(stochastic granular discontinuous deformation method，SGDD)。考虑实际堆石料一般由棱角状和亚棱角状的颗粒组成，提出了不规则多面体颗粒的随机生成算法。马刚等^[12]采用SGDD方法进行了堆石料常规三轴剪切试验的细观数值模拟，再现了堆石料的非线性、压硬性、剪胀和剪缩等主要力学特性。通过选择合适的细观参数，基于SGDD的细观数值试验具有较好的预测能力，为再现堆石料的宏观力学特性、揭示其细观力学机理提供了一条新的途径。周伟等^[13]采用SGDD进行了等、等b应力路径的真三轴数值试验，从宏细观2个层面研究了中主应力对堆石体变形和强度特性的影响。本文作者采用SGDD方法进行了堆石料在等p、等b应力路径的真三轴数值试验，采用伺服控制加载保持试样的应力状态按照预定的应力路径变化。为了避免试样初始各向异性的影响，各向等速压缩松散颗粒集合体制备各向同性的数值试样，颗粒形状为不规则的凸多面体以接近真实的堆石颗粒形状。着重分析了堆石料在真三轴应力状态下的剪胀特性，以及应变增量与应力增量的非共轴性。

1  连续-离散耦合分析方法

连续-离散耦合分析方法结合了有限单元法和离散单元法，将基于有限单元法的连续介质力学分析与基于离散单元法的接触检索、接触力计算和显式动力学求解融合在一起。连续-离散耦合分析方法的提出是为了解决由大量可变形颗粒组成的颗粒集合体的瞬时动力学问题，颗粒间的接触满足互不侵入条件并传递法向和切向接触力，在大多数情况下颗粒会发生断裂、破碎和磨损，导致颗粒形状的变化和颗粒数量的增多。

在连续-离散耦合分析中，通过变分形式简化颗粒间接触的理论假定，认为法向接触力是法向侵入量的函数，而切向接触力是法向接触力和接触状态的函数。将接触的边界问题转化为在接触边界域构造泛函及其变分形式，在接触边界域，接触颗粒的位移场满足：

                 

接触问题的变分形式需要在接触边界域上构造一个泛函，通过寻找泛函的驻值来满足不可贯入条件。为了接触边界域上满足接触约束条件，定义泛函为

  

式中：p为罚参数。

由于：

≥0            

若泛函在接触边界域上为最小值，则罚参数p必须为正值。通过求解式(2)中修正泛函的极小值，近似满足接触约束条件。罚参数p越大，接触约束条件的满足程度越好，当罚参数无穷大时，接触约束条件能够精确满足。在静态或隐式动力学问题中，通过迭代求解的方法来精确满足不可贯入条件。而在瞬时动力学问题中，放弃完全不可贯入条件，而采用足够大的罚函数，使接触的侵入量相对于颗粒尺寸来说可以忽略不计。

采用Munjiza-NBS算法检索颗粒间的接触，基于颗粒的有限元网格离散并结合接触势的概念进行接触力分析。由于每个颗粒都被离散为单独的有限元网格，因此在接触力分析中，可以方便地使用有限元节点的几何坐标来描述接触颗粒的几何形状，并且接触面上接触力的分布更加真实。更重要的是，大大改善了接触边界附近的局部应变场的数值畸变性，当考虑颗粒材料的脆性断裂和破碎时，这一点尤为重要。

采用二阶四面体单元离散颗粒，如图1所示。单元的中任意1点的位置矢量x和位移矢量u用形函数表示为：

     

     

式中：x_i，y_i和z_i为第i个节点在全局坐标系下的坐标；u_i，v _i和w_i为第i个节点在全局坐标系下的位移；N_i为第i个节点的形函数。

定义单元中任意1点在全局坐标系中的初始位置矢量x_i为

             

式中：x_in为单元节点的初始位置矩阵；N为形函数   矩阵。

单元中任意1点在全局坐标系中的当前位置矢量x_c为

              

式中：x_cn为单元节点的当前位置矩阵。

此时变形梯度张量F可表示为

           

定义变形梯度张量F的行列式为Jacobian行列式：

                

由变形梯度张量F定义右置的Cauchy-Creen张量为

                         

采用Neo-Hookean模型来计算单元的应力，Neo-Hookean模型是各向同性线弹性模型的扩展，适用于可压缩的Neo-Hookean材料在大变形情况下的应力计算。此时，Cauchy应力矢量T可表示为

         

式中：和为拉梅常数。

            

式中：E为弹性模量；为泊松比。

图1  二阶四面体单元

Fig. 1  Quadratic tetrahedron element

2  真三轴数值试验

2.1  数值试样

采用随机颗粒生成程序RPG生成形状随机的不规则多面体颗粒，如图2所示。定义与颗粒等体积球的直径为其等效直径，式中V为三维凸多面体颗粒的体积。定义试样的体积分数为等效粒径小于某值的颗粒体积与颗粒总体积之比， (D^*＜d^*)，在本文的数值试验中，试样的体积分数在最大等效粒径 mm与最小等效粒径 mm之间线性变化。试样的相对密度为0.9，对应的孔隙比e=0.518，孔隙率n=0.342。

制备数值试样时，首先采用随机颗粒生成程序RPG在1个较大的立方体空间内生成不接触的松散颗粒集合体。将颗粒集合体信息导入堆石料细观数值模拟软件SGDD中，为了避免由制样产生的初始各向异性，在试样的各个方向采用位移控制等速地压缩试样直至目标大小，如图3所示。在此过程中颗粒间的滑动摩擦角和重力加速度都设为0，且颗粒不发生损伤和破碎。最终生成的数值试样如图4(a)所示，试样形状为立方体，其长×宽×高为300 mm×300 mm× 300 mm，共包含8 927个不规则凸多面体颗粒，试样中颗粒的等效粒径分布如图4(b)所示，采用二阶四面体单元离散为142 277个单元，401 590个节点。由于采用各向等速地压缩制备试样，其组构和接触力的空间分布没有表现出明显的各向异性。采用傅里叶函数来拟合颗粒间接触法向、粒间法向接触力和切向接触力与角度的关系，其数学表达式为：

      

      

式中：为细观组构量与剪切方向的夹角；f₀为所有颗粒法向接触力的平均值；和分别为接触法向和法向接触力各向异性的主方向；a_c和a_n为傅里叶系数，其数值分别反映接触法向和法向接触力的各向异性程度。

图2  不规则多面体颗粒示意图

Fig. 2  Sketch map of irregular polyhedral particle

图3  各向等压制备试样

Fig. 3  Sample preparation by isokinetic triaxial compressing loose particle assembly

图4  数值试样及试样中颗粒的等效粒径分布

Fig. 4  Numerical sample and equivalent particle size distribution

接触法向分布在3个平面上的各向异性系数a_c分别为0.022，0.044和0.018，法向接触力分布在3个平面上的各向异性系数a_n分别为0.012，0.011和0.008，因此可以认为数值试样是初始各向同性的。

2.2  细观参数

数值试验所需要的细观参数较多，其中部分参数可以通过常规物理力学试验直接确定，比如颗粒密度、颗粒弹性模量和泊松比等，还有部分参数可参照一般岩石的取值范围来确定，比如颗粒母岩的内摩擦角、单轴抗压与抗拉强度之比。

除此之外还有颗粒间摩擦角、法向接触刚度K_n、切向接触刚度K_s、颗粒单轴抗压强度f_c和损伤阀值R_n等。本文的真三轴数值试验不针对具体工程的堆石料，只是在一般堆石料的参数取值范围内选取了一套参数进行数值模拟，所用参数见表1。

表1  堆石料真三轴数值试验所用的细观参数

Table 1  Micro-parameters of rockfill in true triaxial numerical test

2.3  应力应变不变量

描述堆石料真三轴应力状态的应力不变量有广义剪应力q、平均静水压力p和应力罗德角：

                

   

，           

式中：中主应力系数b是反映中主应力与大主应力、小主应力相对大小的参数，其取值范围为    0≤b≤1.00。当b=0时，=，应力罗德角=0，试样处于三轴压缩状态；当b=1.00时，，应力罗德角，试样处于三轴拉伸状态。

2.4  加载路径

首先给试样施加三向等压应力直至达到预定的围压值，然后再进行等p、等b应力路径的剪切试验。剪切时，在试样轴向进行位移控制加载，在试样的2个侧面施加应力控制边界条件，在此过程中保持静水压力p和中主应力系数b不变，直至试验结束。在加载过程中，试样的轴向为大主应力方向，大主应力为

                  

式中：F₁为作用在试样轴向加载板上的集中力；，为试样轴向的截面积。

联立以下2式：

             

                 

可得施加在试样2个侧面的中主应力和小主应力分别为：

       

            

在加载过程中，通过伺服控制机制动态调整作用在2个侧面上的集中力荷载F₂和F₃：

         

        

分别进行了平均静水压力p为2.4 MPa，中主应力系数b为0，0.25，0.50，0.75和1.00的等p、等b应力路径的真三轴数值试验，三维应力空间下的应力路径见图5。其中，b=0对应三轴压缩试验(TC)，b= 0.50对应简单剪切试验(SS)，b=1.00对应三轴拉伸试验(TE)。

图5  真三轴数值试验应力路径

Fig. 5  Stress paths of true triaxial tests

3  宏观应力应变关系

图6所示为广义剪应力q、体积应变和偏应力随大主应变的变化曲线。图中体积应变以剪缩为正，剪胀为负。

随着中主应力系数b的增大，颗粒集合体的剪应力与大主应变曲线的初始段斜率越来越大，峰值剪应力及其对应的大主应变却逐渐减小。最大的峰值剪应力出现在b=0时，此时为三轴压缩应力路径，而b=1.00时，颗粒集合体的峰值剪应力最小，此时为三轴拉伸应力路径，这个规律与已有的室内试验和离散元模拟结果相似。不同中主应力系数b时，试样的体积响应不同，在经过短暂而微小的压缩变形后，试样进入剪胀状态，进入剪胀时对应的大主应变和剪胀变形都与中主应力系数b有关。具体来说，随着中主应力系数b的增大，颗粒集合体更快的进入剪胀状态，且剪胀体变更大。偏应力随大主应变的演化曲线与剪应力比较相似，区别在于b=0和b=0.25时偏应力差别不大，而不像剪应力应变曲线那样呈现出明显的单调变化趋势，比如在BARRETO等^[8]的离散元模拟中，b=0.25时的偏应力就大于b=0时的偏应力。

图6  b不同时数值试验结果与大主应变的关系曲线

Fig. 6  Simulated behaviors plotted as functions of major principal strain for different b values

本文的真三轴数值试验以及THORNTON等^{[5, 7-8]}进行的离散元数值模拟试验，颗粒集合体的剪应力均表现出微小的峰后软化。产生这种现象的原因有2个：一是上述研究中所用数值试样的宽高比均为1.0，LADE^[1]研究了试样的宽高比对应变局部化或剪切带的影响，他发现当试样的高宽比为1.0时，剪切带会与试样顶部的加载板和底部的基座相交，相交会约束剪切带的发展并使试样的应变更加均匀；另一个可能是，本文数值试验采用的是不规则的多面体颗粒，颗粒之间咬合作用较强，为颗粒集合体提供了一个更加稳定的细观结构。

图7所示为主应变之间的关系曲线，中主应变、小主应变与大主应变之间均为非线性关系。不同中主应力系数b时，小主应变始终为负值，表明在加载过程中试样在小主应力方向始终发生膨胀，膨胀变形量随中主应力系数b的增大而增大。试样在中主应力方向的变形方向与中主应力系数b有关，在平面应变条件下，试样在中主应力方向既不膨胀也不收缩，此时对应的中主应力系数为b_ps。当b＜b_ps时，中主应变对应膨胀变形；当b＞b_ps时，中主应变对应收缩变形。

图7  b不同时主应变之间的关系曲线

Fig. 7  Relationship between principal strains for different b values

图8所示为偏主应变之间的关系曲线。与主应变之间的非线性关系不同，偏主应变之间近似为线性关系。偏主应变之间的关系曲线可以用下式拟合：

                         

式(25)预测的结果与数值模拟结果拟合得非常好。

图8  b不同时偏主应变之间的关系曲线

Fig. 8  Relationship between principal deviatoric strains for different b values

4  非共轴性

大量的工程实测结果和计算分析均表明：堆石坝内堆石料在填筑期的应力路径近似为等应力比的路 径^[14]。蓄水期上游堆石料内的小主应力方向接近于大坝上游面法向方向，水库蓄水过程中，由于水荷载的作用方向与竣工期坝体内小主应力方向大体一致，随着水荷载的增加，使坝轴线上游侧小主应力增大，而偏应力减小，大小主应力比发生明显的变化，主应力方向也将发生明显的旋转，这将导致主应力方向和主应变率方向的不一致，即非共轴性。目前对堆石料所开展的试验研究几乎都是常规三轴剪切试验，属于等比例加载条件，故不能反映堆石料的非共轴性。

定义加载过程中试样的应力罗德角和应变罗德角为(如图9所示)：

        

         

图10所示为静水压力2.4 MPa时，不同中主应力系数b的应力罗德角和应变罗德角随大主应变的变化曲线。由于采用等p、等b加载应力路径，所以，数值试验中应力罗德角不发生变化。应变罗德角刚开始时与应力罗德角重合，随着加载的进行，应变罗德角开始偏移应力罗德角，开始出现非共轴性。本文的数值试验没有发生明显的应变软化，所以应变罗德角没有明显减小。在加载过程中，应变罗德角对应力罗德角的偏移量可以量化为，不同中主应力系数的见图10(f)。由图10(f)可以看出：在三轴压缩和三轴拉  伸应力路径下，也即b=0和b=1.00时，偏转量小于2°，可以认为试样没有出现非共轴性。

图 9  应力罗德角和应变罗德角的定义

Fig. 9  Definitions of stress lode angle and strain lode angle

图10  应力罗德角和应变罗德角随大主应变的变化

Fig. 10  Evolutions of stress lode angle and strain lode angle with major principal strain

5  剪胀特性

在土力学中，常用应力比q/p和增量塑性应变率来描述堆石料的剪胀特性。基于摩擦性材料的能量耗散特点，ROSCOE等^[15]提出经典剪胀模型：

              

式中：和分别为体积应变增量和剪应变增量的塑性部分；q/p，为应力比；M为对应零剪胀时的应力比。

ROWE^[16]分析了颗粒材料的剪胀特性，将剪胀因子d表示为

          

在Roscoe剪胀模型中，剪胀因子d随应力比线性变化，导致预测的剪胀特性与试验结果差别较大。为了改进Roscoe剪胀模型，JEFFERIES等^[17-20]在式(28)中引入了1个常数：

           

此外，LAGIOIA等^[21]提出了1个适应性更好的剪胀模型：

        

式中：和为模型参数，当趋近0，d趋向于无穷，表明在各向同性加载情况下出现单纯的体积变形。

图11所示为以上4个剪胀模型预测的剪胀曲线与数值试验结果的对比。Roscoe剪胀模型和Rowe剪胀模型的预测能力较差，与数值试验结果相差较远。而Roscoe剪胀模型的2个改进形式，式(30)和式(31)的预测曲线与数值试验结果较吻合。考虑到式(30)的简洁性以及较好的预测能力，在接下来的研究中将采用这个剪胀模型来描述堆石料的剪胀特性。

图12所示为静水压力为2.4 MPa时，不同中主应力系数b的应力比q/p与增量应变比的关系曲线。不同中主应力系数b时，应力比与增量应变比之间近似为线性关系。采用式(30)拟合各子图中的数值试验结果，拟合曲线的斜率为，与轴的交点为特征应力比M。可见：不同中主应力系数b时，拟合剪胀曲线的特征应力比和斜率不同，表现出了中主应力相关性。为了反映中主应力的影响，将特征应力比M和斜率表示为中主应力系数b或者应力罗德角的函数：

               

式中：和分别为在三轴压缩情况下，即b=0时的特征应力比和曲线斜率。式(32)中的为角隅函数，在b=0时，有=1。

图11  不同剪胀模型预测剪胀特性与数值试验结果对比

Fig. 11  Comparison of simulated dilatancy behavior with predications by different dilatancy models

图12  不同中主应力系数的应力剪胀曲线

Fig. 12  Stress-dilatancy relations along different stress paths

图13  剪胀模型的预测值与数值试验结果对比

Fig. 13  Comparisons of dilatancy parameters between predictions and simulated results

采用一个统一的角隅函数来描述特征应力比和曲线斜率与中主应力系数的关系：

        

式中：a为拟合参数；k为三轴拉伸情况下的特征应力比与三轴压缩情况下的比值。由图13可以看出：本文采用的角隅函数拟合效果较好。

6  结论

1) 随着中主应力系数b的增大，颗粒集合体的剪应力与大主应变关系曲线的初始段斜率越来越大，峰值剪应力及其对应的大主应变却逐渐减小。最大的峰值剪应力出现在b=0时，此时为三轴压缩应力路径，而b=1.00时，颗粒集合体的峰值剪应力最小，此时为三轴拉伸应力路径。

2) 在经过短暂而微小的压缩变形后，试样进入剪胀状态，进入剪胀时对应的大主应变和剪胀变形的大小都与中主应力系数b有关。具体来说，随着中主应力系数b的增大，颗粒集合体更快地进入剪胀状态，且剪胀体变更大。

3) 应变罗德角刚开始时与应力罗德角重合，随着加载的进行，应变罗德角开始偏移应力罗德角，开始出现非共轴性。偏转量与中主应力系数b有关，在三轴压缩应力路径(b=0)和三轴拉伸应力路径(b=1.00)时，几乎不发生偏转。

4) 对比了Rowe，Roscoe，改进的Roscoe以及Lagioia剪胀模型，改进的Roscoe剪胀模型预测的剪胀曲线与数值试验结果较吻合，在改进的Roscoe剪胀模型上引入1个角隅函数，以反映真三轴应力状态下的剪胀特性。

参考文献：

[1] LADE P V. Assessment of test data for selection of 3-D failure criterion for sand[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(4): 307-333.

[2] 施维成. 粗粒土真三轴试验与本构模型研究[D]. 南京: 河海大学岩土工程科学研究所, 2008: 1-2.

SHI Weicheng. True triaxial tests on coarse-grained soils and study on constitutive model[D]. Nanjing: Hohai University. Geotechnical Research Institute, 2008: 1-2.

[3] 施维成, 朱俊高, 刘汉龙. 中主应力对砾石料变形和强度的影响[J]. 岩土工程学报, 2008, 30(10): 1449-1453.

SHI Weicheng, ZHU Jungao, LIU Hanlong. Influence of intermediate principal stress on deformation and strength of gravel[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(10): 1449-1453.

[4] XIAO Y, LIU H L, ZHU J G, et al. Dilatancy equation of rockfill material under the true triaxial stress condition[J]. Science China Technological Sciences, 2011, 54(Suppl 1): 175-184

[5] THORNTON C. Numerical simulations of deviatoric shear deformation of granular media[J]. Geotechnique, 2000, 50(1): 43-53

[6] NG T T. Macro-and micro-behaviors of granular materials under different sample preparation methods and stress paths[J]. International Journal of Solids and Structures, 2004, 41(21): 5871-5884

[7] MAHMUD SAZZAD M, SUZUKI K, MODARESSI- FARAHMAND-RAZAVI A. Macro-micro responses of granular materials under different b values using DEM[J]. International Journal of Geomechanics, 2012, 12(3): 220-228

[8] BARRETO D, O’SULLIVAN C. The influence of inter-particle friction and the intermediate stress ratio on soil response under generalised stress conditions[J]. Granular Matter, 2012, 14(4): 505-521.

[9] 周伟, 谢婷蜓, 马刚, 等. 基于颗粒流程序的真三轴应力状态下堆石体的变形和强度特性研究[J]. 岩土力学, 2012, 33(10): 3006-3012.

ZHOU Wei, XIE Tingting, MA Gang, et al. Stress and deformation analysis of rockfill in true triaxial stress conditions based on PFC[J]. Rock and Soil Mechanics, 2012, 33(10): 3006-3012.

[10] MUNJIZA A. The combined finite–discrete element method[M]. New York: John Wiley & Sons Ltd, 2004: 1-2.

[11] 周伟, 常晓林, 周创兵, 等. 堆石体应力变形细观模拟的随机散粒体不连续变形模型及其应用[J]. 岩石力学与工程学报, 2009, 28(3): 491-499.

ZHOU Wei, CHANG Xiaolin, ZHOU Chuangbing, et al. Stochastic granule discontinuous deformation model of rockfill and its application[J]. Chinese Journal of Rock Mechanics and Engineering 2009, 28(3): 491-4995.

[12] 马刚, 周伟, 常晓林, 等. 堆石体三轴剪切试验的三维细观数值模拟[J]. 岩土工程学报, 2011, 33(5): 746-753.

MA Gang, ZHOU Wei, CHANG Xiaolin, et al. 3D mesoscopic numerical simulation of triaxial shear tests for rockfill[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(5): 746-753.

[13] 周伟, 刘东, 马刚, 等. 基于随机散粒体模型的堆石体真三轴数值试验研究[J]. 岩土工程学报, 2012, 34(4): 748-755.

ZHOU Wei, LIU Dong, MA Gang, et al. Numerical simulation of true triaxial tests on mechanical behaviors of rockfill based on stochastic granule model[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(4): 748-755.

[14] 柏树田, 周晓光, 晁华怡. 应力路径对堆石变形特性的影响[J]. 水力发电学报, 1999(4): 76-80.

BAI Shutian, ZHOU Xiaoguang, CHAO Huayi. Effects of stress path on the deformation of rockfill materials[J]. Journal of Hydroelectric Engineering, 1999(4): 76-80.

[15] ROSCOE K H, SCHOFIELD A N, THURAIRAJAH A. Yielding of clays in states wetter than critical[J]. Geotechnique, 1963, 13(3): 211-240

[16] 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 A: Mathematical, Physical and Engineering Sciences. London: The Royal Society, 1962: 500-527.

[17] JEFFERIES M G. Nor-Sand: a simple critical state model for sand[J]. Geotechnique, 1993, 43(1): 91-103

[18] GAJO A, WOOD M. Severn-trent sand: a kinematic-hardening constitutive model: the q–p formulation[J]. Geotechnique, 1999, 49(5): 595-614

[19] LI X S, DAFALIAS Y F, WANG Z L. State-dependant dilatancy in critical-state constitutive modelling of sand[J]. Canadian Geotechnical Journal, 1999, 36(4): 599-611

[20] YANG Y, MURALEETHARAN K K. The middle surface concept and its application to the elasto-plastic behaviour of saturated sands[J]. Geotechnique, 2003, 53(4): 421-431

[21] LAGIOIA R, PUZRIN A M, POTTS D M. A new versatile expression for yield and plastic potential surfaces[J]. Computers and Geotechnics, 1996, 19(3): 171-191.

(编辑  杨幼平)

收稿日期：2015-06-26；修回日期：2015-08-25

基金项目(Foundation item)：国家自然科学基金资助项目(51379161，51509190)；中央高校基本科研业务费专项资金资助项目(2042015kf0022)；博士后科学基金面上资助项目(2015M572195) (Projects(51379161, 51509190) supported by the National Natural Science Foundation of China; Project(2042015kf0022) supported by the Fundamental Research Funds for the Central Universities; Project(2015M572195) supported by China Postdoctoral Science Foundation)

通信作者：马刚，博士(后)，讲师，从事颗粒材料宏细观多尺度力学特性研究；E-mail: magang630@whu.edu.cn