堆石料的三维应力分数阶本构模型

吉华　孙逸飞

摘 要：为合理反映粗粒土的状态依赖非关联应力应变特性，提出应力分数阶塑性力学模型。已有模型基于三轴试验结果，无法对堆石料真三轴条件下的应力应变特性进行预测，为解决这一问题，基于特征应力法对已有分数阶塑性力学模型进行完善。进一步选取不同初始状态条件下堆石料的真三轴压缩试验数据对模型进行验证，结果表明，三维化后的分数阶岩土塑性力学模型可以合理地模拟不同初始状态的堆石料在真三轴压缩条件下的应力应变行为。与传统塑性力学模型相比，提出的分数阶塑性模型在描述堆石料非关联流动时不需要额外引入塑性势函数，仅需对已有屈服面求解分数阶导数。此外，模型在特征应力空间中推导完成再映射到原应力空间，可描述土体的三维强度特性，无需额外采用三维强度准则。

关键词：堆石料;三维应力;塑性力学;分数阶微积分;状态参数;本构模型

中图分类号：TU452 文献标志码：A 文章编号：2096-6717（2022）04-0027-08

Three-dimensional stress-fractional constitutive model for rockfill

JI Hua1， 2， SUN Yifei3

（1. Huadong Engineering Corporation Limited， PowerChina， Hangzhou 310014， P. R. China; 2. College of Civil and Tranportation Engineering， Hohai University， Nanjing 210098， P.R. China; 3. Faculty of Civil and Environmental Engineering， Ruhr University Bochum， Bochum 44801， Germany）

Abstract： The stress-fractional plasticity model has been proposed to model the state-dependent stress-strain behavior of granular soils. However， the previous model was established based on triaxial test， which cannot be easily used to predict the three-dimensional （3D） stress-strain characteristics of rockfill. To solve this problem， an attempt is made to modify the fractional plasticity model for rockfill by using a novel concept of characteristic stress to make the model capable of predicting the 3D stress-strain behavior of rockfill. A series of true triaxial test results of rockfills under different initial states are simulated， from which a good agreement between the model predictions and the corresponding testing results is found. Comparing with the classical plasticity model， no additional plastic potential is required; in addition， the model can describe 3D strength characteristics by mapping from the characteristic stress space to the normal stress space， where no specific 3D strength criterion is needed.

Keywords：rockfill; three-dimensional stress; plasticity; fractional calculus; state parameter; constitutive model

堆石料常常被用于水利工程，诸如堆石坝、堤坝等的建造[1]。因此，在此类工程的设计、建造和运营全过程中，堆石料的强度和变形特性显得尤为重要。由于河流水荷载的复杂性，堆石料在其全寿命周期内经常会遭受不同应力路径的外部荷载作用。为了分析不同应力路径下堆石料的强度和变形行为，学者们[2-6]开展了一系列试验和理论研究，并建立了相应的本构模型。Xiao等[5]通过对不同围压和初始密度条件下堆石料的强度和剪胀行为进行研究，提出了状态依赖的堆石料边界面塑性力学本构模型。陈生水等[7]通过研究颗粒破碎对堆石料临界状态特性等的影响，建立了一个可以考虑堆石料颗粒破碎特性的广义塑性力学模型。此外，在实际工程中，堆石料等粗粒土常遭受三维荷载作用，中主应力对其强度和变形行为影响较大[8]。为此，Xiao等[9]提出了基于统一破坏准则的堆石料三维边界面塑性本构模型;姚仰平等[10-11]、Tian等[12]提出了基于转换应力法[13]的软土和粗粒土三维统一硬化模型。为了考虑材料的非关联塑性流动特性，以上這些模型在建立过程中均额外引入了不同于屈服面函数的塑性势面函数。这些塑性势面函数往往基于经验假定。那么，是否可以在不引入塑性势的条件下建立堆石料的三维本构模型呢？为此，Sun等[14]通过分析岩土材料的三轴试验特性，提出了土体的应力分数阶塑性力学方法[14-15]。Lu等[16]基于特征应力法，提出了软土的三维分数阶塑性本构模型。Qu等[17]建立了适用于岩石的分数阶塑性力学模型。李海潮等[18]进一步引入温度效应，拓展建立了软岩的分数阶热塑性模型。这些模型均处于“应力分数阶塑性力学”框架，还有一些诸如“时间分数阶塑性力学”框架下的模型[19]，仍然采用屈服面和塑性势面不一致的假定，这里不再赘述。不同于传统塑性力学，应力分数阶塑性力学基于非局域性的分数阶微分算子，在某一应力点求导所得塑性流动方向不仅与该点的应力状态有关，还与到达该点的加载历史或相对临界状态的距离有关。在实际建模过程中，仅需对屈服面进行一阶和分数阶求导，即可建立堆石料的非关联弹塑性本构模型。然而，目前提出的分数阶塑性力学模型[14， 20]仅对粗粒土的三轴应力应变特性进行了分析，对土体真三轴应力应变特性尚未进行预测和模拟。58794B6A-604A-455F-A538-39F6135C5C96

为解决这一问题并同时完善现有分数阶岩土塑性力学理论体系，有必要建立一个三维应力分数阶塑性力学本构模型。笔者拟通过特征应力法[21]对现有分数阶塑性力学模型进行修正，使其可以合理地预测土体的三维应力应变关系。

1 基本定义

1.1 应力应变的定义

根据临界状态土力学[22]，约定堆石料是均匀各向同性材料;受压为正，受拉为负;模型推导所用到的应力均为有效应力。有效增量应力张量（σ·′）定义为

式中：i、j = 1、2、3;σ·′ij为应力分量;ei、ej为基底矢量。相应的增量变张量（ε·）定义为

式中：ε·ij为应变分量。此外，约定在主应变空间中ε·1、ε·2和ε·3分别代表第一、第二和第三主应变增量;相应的主应力空间中σ·′1、σ·′2和σ·′3分别是第一、第二和第三有效主应力增量。

除此之外，还有平均有效主应力（p′）、偏应力（q）、中主应力系数（b）这几个需要用到的定义。

式中：偏应力张量s=σ′-p′I，I为单位张量。

1.2 分数阶导数

Caputo的分数阶左导数和右导数分别定义为[23]

式中：D =α/xα;α∈（n-1，n]，是求导阶数，其中n是正整数，取值范围为1～2;σ′是加载应力;Γ（x）=∫∞0e-ττx-1dτ。

1.3 特征应力法

根据文献[21]，用特征应力来代替原有强度准则中的常规应力，使得修正后的强度准则可以反映土体的三维强度特性。特征应力（σ′ij）定义为

σ′ij=paσ′ijpaβ （7）

式中：pa = 100 kPa是大气压强。β∈（0，1]是特征应力参数，反映了中主应力比对土体强度的影响。根据Lu等[21]的研究，β的大小不受土体密度和围压等因素影响，对于给定材料，它是一个定参数。据此，修正后的特征有效平均主应力（p′）和特征偏应力（q）分别为

2 三维分数阶应力应变关系

增量主应变可以分解为弹性应变和塑性应变两部分，即

式中：ε·e为弹性应变张量;

ε·p为塑性应变张量，根据非线性K-G弹性理论，可得

式中：G是剪切模量;K是体积模量，分别定义为

式中：G0是弹性常数;ν是泊松比;e是孔隙比。此外，塑性应变张量（ε·p）可以计算为

式中：Kp为硬化模量，其定义将在后面的本构模型部分给出。n和m分别为加载方向和塑性流动方向，定义为[24]

在知道屈服面的情况下，分数阶塑性力学模型无需再额外假定一个塑性势面函数，而是通过分数阶微分直接得到一个非关联的流动法则，从而使模型得到简化。由文献[15]可知，当α越大时，所得到的非关联性越明显。具体参见图1，当α由1.05增加到1.3时，同一应力水平下加载方向与流动方向的夹角变大，即非关联程度增加。

3 三维分数阶塑性力学模型

3.1 屈服面

由传统塑性力学可知，加载方向是屈服面的法向。选取特征应力条件下的修正剑桥模型屈服面函数来描述土体的受荷屈服。

式中：p′0是特征应力空间屈服面与p′轴的交点。M是特征应力空间的临界状态参数，可定义为

式中：φc是常规三轴压缩时的临界状态摩擦角。如图2所示，特征应力法修正后的屈服面在常规主应力空间的屈服面不再是一个圆锥，而是一个在π平面的投影，是曲边三角形的曲面，如图3所示。因此，可以反映土体在真三轴应力条件下的屈服特性。

3.2 加载张量与流动张量

将式（19）分别代入式（17）、式（18），进一步求得df和dg为将式（21）和式（22）代回至式（15）和式（16），即可得到塑性加载张量和塑性流动张量的具体表达式

比較式（23）和式（24）可知，当α=1时，模型退化为传统关联流动法则。由图4可知，对屈服面进行分数阶求导后，所得子午平面上塑性流动方向不再垂直于屈服面，而是与加载方向存在一个夹角。

此外，为了可以反映受荷变形的状态依赖特性，Sun等[14-15]建议采用状态依赖的分数阶数

式中：Δ>0，是材料常数。式（25）是经验公式，来源于基于试验现象的考虑：土体所处状态会影响其塑性流动特性，而分数阶数决定了塑性流动方向。因此，运用式（25）直观地考虑了加载过程中土体所处状态对其塑性流动方向的影响。ψ是状态参数，可以定义为

式中：ec是临界状态孔隙比，定义为[25]

式中：eΓ、λ和ξ是e-p′空间的临界状态参数。

3.3 硬化模量

硬化模量描述了土体的硬化特性。Sun等[14]采用Li等[25]提出的常规三轴应力条件下的硬化模量。为了使建立的模型可以反映真三维应力条件下的硬化特性，基于特征应力法修正了Li等[25]的硬化模量，见式（28）。

式中：特征峰值应力比（Mp）定义为

式中：h1和h2是材料常数，描述了硬化参数对孔隙比的依赖;k是材料常数，描述了特征峰值应力比对土体物质状态的依赖。由式（28）可知，硬化模量的大小受多重因素影响（e、e0、p′、ψ、G），且在剪切尚未开始时（η=0），Kp=+∞，在临界状态和峰值应力点时，Kp=+∞。因此，所提三维硬化模量合理。为进一步验证所提出的模型，选择某堆石料的真三维试验数据对模型进行验证。

4 模型参数识别与模型验证

4.1 模型参数识别

模型含有11个参数，分别为：φc、λ、eΓ、ξ、Δ、β、k、h1、h2、G0以及ν。这些模型参数均可以通过三轴试验获得。

参数φc、λ、eΓ和ξ为临界状态参数。φc可以通过拟合p′-q平面临界状态线的斜率得到;λ、eΓ、ξ可以通过拟合e-p′空间的临界状态数据得到[5]。图5所示为临界状态参数对模型预测结果的影响，可见，随着φc和eΓ数值的增加，模型预测得到的同一应变水平下的应力比逐渐增加;随着λ的增加，模型预测的同一应变水平下的应力比逐渐减小。58794B6A-604A-455F-A538-39F6135C5C96

参数Δ决定了土体的塑性流动特性，可以由试样状态转换时的剪胀关系得到，即d=0。将式（25）代入式（22），得

式中：ψd和ηd分别是ψ和η在状态转换点的数值。参数k可以通过试样在峰值应力状态点的应力应变关系得到，即Kp=0时

式中：ψp和ηp分别是ψ和η在峰值应力点的数值。

参数Δ和k对模型预测的影响见图6。分析发现：随着参数Δ的增加，模型预测得到的同一應变水平下的应力比逐渐增加;随着参数k的增加，模型预测得到的同一应变水平下的应力比逐渐降低。

特征应力参数β可由土体相关三轴压缩（φc）和拉伸（φe）时的临界状态摩擦角关系式获得[21]。β≤0.1时，模型在π平面的屈服轨迹非常接近于Matsuoka-Nakai准则[26]。因此，β近似取0.1。

参数h1和h2可以通过拟合土体受荷时ε1-q关系曲线获得。具体参见文献[15， 25]，这里不再赘述。图7所示为参数h1和h2对模型预测结果的影响，分析可知：当参数h1增加时，同一应变水平下模型预测的应力比增加;当参数h2增加时，同一应变水平下模型预测的应力比减小。

两个弹性参数（G0、ν）可以通过拟合堆石料初始加载阶段的应力应变曲线获得，具体参见文献[27]。由于受荷后堆石料弹性变形相对塑性变形较小，因此，不再探讨弹性参数对模型预测的影响。

4.2 模型验证

选取不同初始状态的堆石料[28]在真三轴压缩条件下的试验数据进行拟合和对比，以验证所提出三维分数阶岩土塑性力学模型的合理性。

所选堆石料为两河口心墙土石坝筑坝料，试验初始孔隙比为0.68。试验时固定第三主应力和中主应力比不变，第二和第一主应力按照设定中主应力比值的大小比例增加。其他有关材料和试样准备资料见文献[28]，不再赘述。根据前述方法，确定具体模型参数为：G0 = 90、ν=0.25、φc=46°、λ=0.11、eΓ=0.404、ξ=0.1、Δ=0.2、β = 0.1、k = 0.1、h1 = 1.2、h2 = 0.3。

图8～图10为所提出分数阶塑性力学模型对真三轴试验数据的模拟。对比分析可以发现：模型可以合理地模拟堆石料在不同围压和中主应力比条件下的真三维应力应变特性。具体来说，随着第一主应变的增加，试样的应力比不断增加;第二主应变只有在中主应力比等于0时，即常规三轴试验条件下，数值为负，其他中主应力比条件下均为正;第三主应变在任何中主应力比条件下均为负值。这些应力应变特性均可被较好地模拟。

此外，需要指出的是，中主应力系数较小时，图8～图10中试样在同一偏应力水平所产生的第一主应变大于第二、第三主应变，但随着中主应力系数的增加，第一主应变小于第二、第三主应变;这是因为随着中主应力系数从0变化到1，原先的第二主应力大小慢慢接近于第一主应力，对侧向应变的主导作用逐渐变大。从而在第一主应力方向产生的应变会逐渐小于同等偏应力水平下其他两个主应力方向的应变。

5 结论

基于特征应力法，修正了已提出的分数阶应力塑性力学模型，使其可以对堆石料在一般三维应力条件下的应力应变进行预测。并选择堆石料的真三轴压缩试验数据，对模型进了验证。主要结论如下：

1）基于特征应力修正的屈服面，其在主应力空间不再是一个圆锥，而是一个可以反映土体三维屈服特性的曲面，该面在π平面的投影为曲边三角形。

2）无需额外假定塑性势面，只需要对已有屈服面进行分数阶微分求导，便可以得到非关联的流动法则。进一步得出剪胀方程受分数阶求导阶数、特征应力参数以及中主应力系数的影响。

3）基于特征应力法修正了已有常规三轴应力条件下的硬化模量，使其可以反映三维应力条件下堆石料的硬化特性。

4）对比试验数据，验证了所提模型的合理性;模型可以对不同初始状态试样在不同中主应力比条件下的真三轴应力应变特性进行较好的模拟。

5）建立模型仅考虑了完全饱和堆石料的变形特性，对于干燥或半干燥堆石料遇水饱和后存在的湿化或者软化性能还有待进一步研究。此外，模型也未考虑堆石料在填筑或者运维过程中的颗粒破碎，对于考虑颗粒破碎影响的堆石料三维分数阶塑性模型也有待进一步研究。

参考文献：

[1] 刘汉龙， 孙逸飞， 杨贵， 等. 粗粒料颗粒破碎特性研究述评[J]. 河海大学学报（自然科学版）， 2012， 40（4）： 361-369.

LIU H L， SUN Y F， YANG G， et al. A review of particle breakage characteristics of coarse aggregates [J]. Journal of Hohai University （Natural Sciences）， 2012， 40（4）： 361-369. （in Chinese）

[2] 刘汉龙， 秦红玉， 高玉峰， 等. 堆石粗粒料颗粒破碎试验研究[J]. 岩土力学， 2005， 26（4）： 562-566.

LIU H L， QIN H Y， GAO Y F， et al. Experimental study on particle breakage of rockfill and coarse aggregates [J]. Rock and Soil Mechanics， 2005， 26（4）： 562-566. （in Chinese）

[3] 贾宇峰， 王丙申， 迟世春. 堆石料剪切过程中的颗粒破碎研究[J]. 岩土工程学报， 2015， 37（9）： 1692-1697.

JIA Y F， WANG B S， CHI S C. Particle breakage of rockfill during triaxial tests [J]. Chinese Journal of Geotechnical Engineering， 2015， 37（9）： 1692-1697. （in Chinese）58794B6A-604A-455F-A538-39F6135C5C96

[4] XIAO Y， LIU H L， CHEN Y M， et al. Strength and deformation of rockfill material based on large-scale triaxial compression tests. I： Influences of density and pressure [J]. Journal of Geotechnical and Geoenvironmental Engineering， 2014， 140（12）： 04014070.

[5] XIAO Y， LIU H L， CHEN Y M， et al. Bounding surface model for rockfill materials dependent on density and pressure under triaxial stress conditions [J]. Journal of Engineering Mechanics， 2014， 140（4）： 04014002.

[6] DAFALIAS Y F， TAIEBAT M. SANISAND-Z： Zero elastic range sand plasticity model [J]. Géotechnique， 2016， 66（12）： 999-1013.

[7] 陳生水， 傅中志， 韩华强， 等. 一个考虑颗粒破碎的堆石料弹塑性本构模型[J]. 岩土工程学报， 2011， 33（10）： 1489-1495.

CHEN S S， FU Z Z， HAN H Q， et al. An elastoplastic model for rockfill materials considering particle breakage [J]. Chinese Journal of Geotechnical Engineering， 2011， 33（10）： 1489-1495. （in Chinese）

[8] XIAO Y， LIU H L， CHEN Y M， et al. Influence of intermediate principal stress on the strength and dilatancy behavior of rockfill material [J]. Journal of Geotechnical and Geoenvironmental Engineering， 2014， 140（11）： 04014064.

[9] XIAO Y， LIU H L， ZHU J G， et al. Modeling and behaviours of rockfill materials in three-dimensional stress space [J]. Science China Technological Sciences， 2012， 55（10）： 2877-2892.

[10] YAO Y P， YAMAMOTO H， WANG N D. Constitutive model considering sand crushing [J]. Soils and Foundations， 2008， 48（4）： 603-608.

[11] 姚仰平， 刘林， 罗汀. 砂土的UH模型[J]. 岩土工程学报， 2016， 38（12）： 2147-2153.

YAO Y P， LIU L， LUO T. UH model for sands [J]. Chinese Journal of Geotechnical Engineering， 2016， 38（12）： 2147-2153. （in Chinese）

[12] TIAN Y， YAO Y P. Modelling the non-coaxiality of soils from the view of cross-anisotropy [J]. Computers and Geotechnics， 2017， 86： 219-229.

[13] YAO Y P， WANG N D. Transformed stress method for generalizing soil constitutive models [J]. Journal of Engineering Mechanics， 2014， 140（3）： 614-629.

[14] SUN Y F， XIAO Y. Fractional order plasticity model for granular soils subjected to monotonic triaxial compression [J]. International Journal of Solids and Structures， 2017， 118/119： 224-234.

[15] SUN Y F， SHEN Y. Constitutive model of granular soils using fractional-order plastic-flow rule [J]. International Journal of Geomechanics， 2017， 17（8）： 04017025.

[16] LU D C， LIANG J Y， DU X L， et al. Fractional elastoplastic constitutive model for soils based on a novel 3D fractional plastic flow rule [J]. Computers and Geotechnics， 2019， 105： 277-290.58794B6A-604A-455F-A538-39F6135C5C96

[17] QU P F， ZHU Q Z， SUN Y F. Elastoplastic modelling of mechanical behavior of rocks with fractional-order plastic flow [J]. International Journal of Mechanical Sciences， 2019， 163： 105102.

[18] 李海潮， 马博， 张升， 等. 基于分数阶热弹塑性理论的软岩力学特性描述[J]. 岩石力学与工程学报， 2020， 39（7）： 1311-1320.

LI H C， MA B， ZHANG S， et al. Mechanical behaviors of soft rocks based on the fractional thermal elastic-plastic theory [J]. Chinese Journal of Rock Mechanics and Engineering， 2020， 39（7）： 1311-1320. （in Chinese）

[19] SUN Y， XIAO Y. Fractional order model for granular soils under drained cyclic loading [J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2017， 41（4）： 555-577.

[20] SUN Y， GAO Y， SHEN Y. Mathematical aspect of the state-dependent stress-dilatancy of granular soil under triaxial loading [J]. Géotechnique， 2019， 69（2）： 158-165.

[21] LU D C， MA C， DU X L， et al. Development of a new nonlinear unified strength theory for geomaterials based on the characteristic stress concept [J]. International Journal of Geomechanics， 2017， 17（2）： 04016058.

[22] SCHOFIELD A， WROTH P. Critical state soil mechanics [M]. New York， USA： McGraw-Hill London， 1968.

[23] PODLUBNY I. Fractional differential equations： An introduction to fractional derivatives， fractional differential equations， to methods of their solution and some of their applications [M]. San Diego， California： Academic press， 1998.

[24] SUN Y F， NIMBALKAR S. Stress-fractional soil model with reduced elastic region [J]. Soils and Foundations， 2019， 59（6）： 2007-2023.

[25] LI X S， DAFALIAS Y F. Dilatancy for cohesionless soils [J]. Géotechnique， 2000， 50（4）： 449-460.

[26] MATSUOKA H， NAKAI T R. Stress-deformation and strength characteristics of soil under three different principal stresses [J]. Proceedings of the Japan Society of Civil Engineers， 1974， 232： 59-70.

[27] SUN Y F， WICHTMANN T， SUMELKA W， et al. Karlsruhe fine sand under monotonic and cyclic loads： Modelling and validation [J]. Soil Dynamics and Earthquake Engineering， 2020， 133： 106119.

[28] XIAO Y， LIU H L， ZHU J G， et al. A 3D bounding surface model for rockfill materials [J]. Science China Technological Sciences， 2011， 54（11）： 2904-2915.

（編辑 王秀玲）

收稿日期：2020-07-23

基金项目：中央高校基本科研业务费（2017B05214）;中国博士后科学基金（2017M621607）

作者简介：吉华（1987- ），博士（后），主要从事岩土材料力学特性与本构理论研究，E-mail：ji_h@ecidi.com。

孙逸飞（通信作者），博士，E-mail：yifei.sun@rub.de。

Received：2020-07-23

Foundation items：Fundamental Research Funds for the Central Universities （No. 2017B05214）; China Postdoctoral Science Foundation （No. 2017M621607）

Author brief：JI Hua （1987- ）， postdoctoral fellow， main research interests： mechanical properties and constitutive theory of geomaterials， E-mail： ji_h@ecidi.com.

SUN Yifei （corresponding author）， PhD， E-mail： yifei.sun@rub.de.58794B6A-604A-455F-A538-39F6135C5C96