较好地改善了Gudehus-Bauer亚塑性模型模拟的堆石料三轴压缩体积应变曲线。利用改进的Gudehus-Bauer亚塑性本构模型对面板堆石坝进行了非线性有限元应力变形分析，将获取的数值结果与双屈服面模型和E-B模型对比分析，结果表明：改进的Gudehus-Bauer亚塑性模型计算的堆石体沉降量、水平位移分布、主应力分布与沈珠江双屈服面模型较为接近，同时克服了蓄水期邓肯E-B模型面板底部顺坡向拉应力偏大的缺点。摘 要： 介绍了Gudehus-Bauer亚塑性本构模型及其参数，并建立了拟合指数α和β与围压的关系。用Gudehus-Bauer亚塑性模型模拟了堆石料三轴压缩试验曲线，并与堆石料试验值进行比较。根据模拟结果，在模型中增加了体积应变控制项较好地改善了Gudehus-Bauer亚塑性模型模拟的堆石料三轴压缩体积应变曲线。利用改进的Gudehus-Bauer亚塑性本构模型对面板堆石坝进行了非线性有限元应力变形分析，将获取的数值结果与双屈服面模型和E-B模型对比分析，结果表明：改进的Gudehus-Bauer亚塑性模型计算的堆石体沉降量、水平位移分布、主应力分布与沈珠江双屈服面模型较为接近，同时克服了蓄水期邓肯E-B模型面板底部顺坡向拉应力偏大的缺点。
近年来，一种新型的亚塑性本构模型被用于描述堆石料应力-应变关系。该模型以非线性连续介质力学为理论基础，以张量分析为计算工具，直接建立了颗粒体材料的应力率与应变率的关系。同时，该模型还减少了一些人为的假定，较为客观地反映了材料复杂的力学特性。
目前国内外研究人员已建立了较多的亚塑性本构模型，该模型最早由Kolymbas[1]针对砂料提出，并被广泛推广。Gudehus[2]、Bauer[3]针对砂料提出了Gudehus-Bauer亚塑性本构模型，在模型中引入了密度因子和刚度因子，考虑了砂粒料的孔隙比变化特征和临界状态，较好地描述了砂料的主要力学特性。Herle等[4]研究了亚塑性模型，并对其参数确定方法作了详细的介绍。Fuentes等[5]基于亚塑性理论提出颗粒间应变的空间屈服函数，反映了应变历史状态。Wang等[6]提出一种改进的边界面亚塑性本构模型，描述了松砂或者密砂不同阶段的循环加载的力学特性，考虑了材料状态相关剪胀性理论和累积塑性应变对塑性模量的影响。Lin等[7]基于微亚塑性模型理论，采用离散元法和连续统法模拟颗粒料的双轴试验和简单剪切试验，获得了较好的成果。岑威钧等[8-9]在亚塑性模型中增加体变控制项，且取-0.04，并用于描述堆石料应力变形特性，改善了模型模拟的堆石料体积应变曲线。吴长彬等[10]、岑威钧[11]建立了拟合参数与平均压力的关系。Xiang等[12]考虑了不同应力路径的影响，提出了增量非线性亚塑性本构模型，获得了试验验证。明华军等[13]引入考虑颗粒破碎堆石料特征孔隙比的表达式，改进的模型较好地描述堆石料的力学特性。陈泽钦等[14]在Wang-Wu亚塑性模型的基础上，提出了七参数模型，较好描述了堆石料应力变形特性。王仁超等[15]研究了亚塑性模型自适应隐式和显式两种不同的积分算法，并在ABAQUS平台开发UMAT-VUMAT接口，有效地实现了亚塑性模型数值分析应用。唐扬[16]推导了主应力空间的亚塑性模型稳定面，并应用于滑坡稳定性数值计算，较好地分析了潜在不稳定区域。Osinov等[17]基于亚塑性模型进行隧道的动力数值分析，研究了衬砌和土体的应力和变形。Hleibieh等[18]通过与不同边界条件和初始应力状态下获取的砂土试验数据对比，认为亚塑性模型可以有效地描述不同动力响应条件下砂土的应力-应变特性。本文利用改进的Gudenhus-Bauer亚塑性本构模型对面板堆石坝进行了非线性有限元应力分析。
Gudehus-Bauer亚塑性模型由Gudehus等[2-3]提出，其模型为：
(1)
式中：
为应变率张量；
为Jaumann应力率张量；
为规范化应力张量，且
为规范化应力偏张量，
为应变率张量的Euclidean范数；tr(·)为张量求迹；a1为无量纲标量因子；fd为密度因子；fs为刚度因子。
其中Gudehus-Bauer模型中的刚度因子fs和密度因子fd表达式如式(2)和式(4)所示：
(2)
(3)
(4)
式中：α和β为模型的拟合参数；ps为平均压力，ps=-tr(σ)/3；ei0、ec0和ed0分别为零压力下ei、ec、ed对应的界限孔隙比；hs为颗粒硬度，kPa；n为拟合指数。
下界孔隙比ed、临界孔隙比ec和上界孔隙比ei随平均压力ps增加而减小[2]，各孔隙比的变化规律如式(5)所示：
(5)
在无量纲因子a1中引入了Lode角θ，则a1为偏应力方向的函数，如式(6)所示。在常规三轴试验时，无量纲因子a1为常数[19]，即a1=1/c1，
(6)
式中：
为摩擦角；
因此，Gudehus-Bauer亚塑性模型共有8个参数，分别为α、β、φc、hs、n、ei0、ed0、ec0。
根据文献[4, 20]介绍的方法计算了坝体填筑料的Gudehus-Bauer亚塑性本构模型参数中的6个参数，即颗粒硬度hs、拟合指数n、不同零压力下界限孔隙比ei0、ed0、ec0和摩擦角φc等参数，如表1所示。
对于拟合指数α的确定，考虑常规三轴压缩试验峰值状态时
并代入式(1)和式(4)即可求得拟合指数α。依据此结果，吴长彬等[10]、岑威钧[11]建立了指数α与平均应力p的关系。本文在求取拟合指数α后，进行适当的调整，并假定拟合指数α与ln(σ3/pa)具有式(7)的关系：
α=k1ln(σ3/pa)+b1。
(7)
式中：k1为斜率，b1为截距，均为拟合参数。
运用侧限压缩试验求取拟合指数β，边界条件为
引入应力比K=σ11/σ33且K=1-sin φc，增量刚度模量
代入式(1)和式(2)中即可得到拟合指数β，如式(8)所示：
β=ln A/ln B。
(8)
式中：B=ei/e；
。
由式(8)可知，β与ps具有较为复杂的关联，为简化这一关系，假定β与ln(σ3/pa)具有式(9)的关系：
β=k2ln(σ3/pa)+b2。
(9)
式中：k2为斜率；b2为截距，均为拟合参数。
由以上方法即可求取Gudehus-Bauer亚塑性模型的8个参数结果，如表1所示。
表1 模型参数求取结果
Table 1 The result of the model parameters
料层ec0ei0ed0nφc /(°)hs/kPab1k1b2k2λ次堆料0.6680.8350.2720.99245.140 0320.003 3-0.005 10.439 60.092 1-0.02主堆料0.6580.8230.2561.03645.438 929-0.000 5-0.001 20.506 00.071 7-0.03过渡料0.6370.7970.2201.13541.132 6000.055 6-0.016 20.250 20.340 7-0.05垫层料0.5950.7440.1840.96641.846 6980.051 8-0.013 10.049 20.261 40.10
根据表1的参数结果，用Gudehus-Bauer亚塑性模型模拟了主堆石料三轴压缩试验，如图1所示。
图1 主堆石料的三轴试验值与模型模拟结果
Figure 1 The test values and the model simulated
values of major rockfill materials in triaxial test
由图1可知，Gudehus-Bauer亚塑性模型描述的堆石料偏应力曲线与试验值吻合较好；描述的体积应变曲线与试验值的吻合程度较差，该模型不能较好地反映堆石料的应变状态，需进一步改进。
在模型中引入了体积应变项
改进后的模型方程为：
(10)
式中：λ为体积应变调整参数，为拟合参数。
通过对比岑威钧等[8-9]对Gudehus-Bauer亚塑性模型的研究成果，即在模型方程中引入
且b取-0.04，可得b=λσ/tr(σ)，此时b为一变量，且与堆石料的应力状态有关。
用改进的Gudehus-Bauer亚塑性模型模拟了堆石料三轴压缩试验，体积应变曲线与试验值更为吻合，有效改善了原模型这一缺点。
某抽水蓄能电站大坝坝型为面板堆石坝，坝顶高程299.90 m，正常蓄水位294.00 m，最大坝高75.1 m。面板厚度采用渐变的方式，顶部厚度0.3 m，底部厚度0.5 m，坝坡及坝体材料分区如图2所示。
图2 大坝河床断面(m)
Figure 2 The riverbed section of the dam(m)
坝体分为41个横断面，坝体有限元计算网格划分如图3所示，单元总数为7 929，结点总数为8 420。在有限元计算中，采用28级进行加载以模拟坝体施工过程，前18级模拟堆石体施工过程，19级和20级分别模拟面板施工和坝顶施工，第28级大坝蓄水至正常蓄水位为294 m。
图3 堆石坝三维网格
Figure 3 FEM mesh of dam
混凝土面板和坝体堆石料采用三维实体单元模拟,面板垂直缝和周边缝采用无厚度Goodman连接单元模拟，单元的参数取值参考天生桥一级面板堆石坝的力-位移模型试验成果[21]。在材料性质相差较大的接触面设置无厚度Goodman接触面单元模拟，如混凝土与垫层料接触面，趾板与堆石体接触面，单元参数从文献[21]取值。
将改进的Gudehus-Bauer亚塑性模型的分析成果与邓肯E-B模型和沈珠江双屈服面模型数值成果进行对比。其中，采用中点增量法求解非线性方程组，双屈服面模型的屈服面中广义剪应力q采用八面体剪应力，邓肯E-B模型采用平均主应力p和广义剪应力q等效代替σ3和(σ1-σ3)[21]。
竣工期坝体仅受自重时，坝体沉降量和水平位移的最大值及其变形规律与面板堆石坝常见的变形规律较为吻合，如图4和图5所示。
图4 竣工期沉降量(cm)
Figure 4 The vertical displacement in
construction period(cm)
图5 竣工期水平位移(cm)
Figure 5 The horizontal displacement in
construction period(cm)
由图4可知，改进的Gudehus-Bauer亚塑性模型和双屈服面模型计算的沉降量相近，最大值分别为22.9 cm和24.0 cm，E-B模型的最大沉降量为32.1 cm。3种模型计算的坝体水平位移基本呈对称分布。
由图5可知，改进的Gudehus-Bauer亚塑性模型和双屈服面模型计算的水平位移分布规律相近。改进的Gudehus-Bauer亚塑性模型、邓肯E-B模型和双屈服面模型计算的最大上游向水平位移分别为3.43、2.84、5.97 cm，最大下游向水平位移分别为4.45、3.53、6.93 cm。
蓄水期的坝体沉降量和水平位移如图6和图7所示。
图6 蓄水期沉降量(cm)
Figure 6 The vertical displacement in storage
period(cm)
图7 蓄水期水平位移(cm)
Figure 7 The horizontal displacement
in storage period(cm)
由图6可知，蓄水期坝体的沉降量明显有所增大。改进的Gudehus-Bauer亚塑性模型、双屈服面模型和邓肯E-B模型的最大沉降量分别为26.4、26.5、34.8 cm，比竣工期增加3.5、2.5、2.7 cm，约占坝高的0.35%、0.35%和0.46%。
由图7可知，改进的Gudehus-Bauer亚塑性模型和双屈服面模型计算的蓄水期坝体水平位移相近，最大上游向位移分别为1.41 cm、2.98 cm，最大下游向位移为5.31 cm、7.92 cm。E-B模型计算的水平位移发生较大变化，蓄水期产生上游的位移基本消失。
蓄水期面板挠度分布规律，如图8所示。
图8 蓄水期面板挠度
Figure 8 Deflection of face slabs
由图8可知，蓄水期面板在水压力作用下产生明显的下游向位移，分布规律为：面板的最大挠度出现在河谷段面板下部约1/3坝高处，且挠度值向四周逐渐减小；改进的Gudehus-Bauer亚塑性模型与双屈服面模型计算的最大挠度较为接近，分别为0.094 m和0.078 m；而邓肯E-B模型的计算值偏大，最大为0.138 m。
蓄水期面板坝轴线向位移，如图9所示。由如图9可知，面板的总体趋势是岸坡段面板向河床位移，3种模型计算的面板最大坝轴线向位移较为接近，改进的Gudehus-Bauer亚塑性模型、邓肯E-B模型和双屈服面模型计算的最大左岸向位移分别为0.017、0.015、0.011 cm，最大右岸向位移分别为1.6、1.4、1.0 cm。
由图10可知，3种模型计算的蓄水期河床段附近的面板中部顺坡向主要承受压应力。改进的Gudehus-Bauer亚塑性模型和双屈服面模型计算的面板底部和顶部顺坡向处于受拉状态，最大值分别为0.1、0.2 MPa，远小于邓肯E-B模型面板底部顺坡向最大拉应力。邓肯E-B模型计算的面板底部顺坡向基本上都处于受拉状态，最大拉应力为1.70 MPa，出现在河床段面板底部。
图9 面板轴线向位移
Figure 9 Axial displacement of face slabs
图10 面板顺坡向应力
Figure 10 Stress along slop direction of
face slabs
由如图11可知，蓄水期河床段面板中下部承担了最大的轴线向压应力，改进的Gudehus-Bauer亚塑性模型、邓肯E-B模型和双屈服面模型计算的最大坝轴线向压应力分别为2.24、3.60、1.68 MPa，面板坝轴线向基本上处于受拉状态，但拉应力不大，分别为0.49、0.69、0.40 MPa。
图11 面板轴线向应力
Figure 11 Axial stress of face slab
根据改进的堆石料Gudehus-Bauer亚塑性本构模型，将其应用于面板堆石坝的应力变形分析，得到如下结论：
(1)改进的Gudehus-Bauer亚塑性模型计算的堆石体沉降量、水平位移以及坝体的主应力等主要结果与沈珠江双屈服面模型的计算结果较为接近。
(2)改进的Gudehus-Bauer亚塑性模型计算的蓄水期面板应力和变形分布规律与沈珠江双屈服面模型和邓肯E-B模型基本一致，且克服了蓄水期邓肯E-B模型面板底部顺坡向拉应力偏大的缺点。因此，改进的Gudehus-Bauer亚塑性模型为堆石坝数值计算提供了参考。
[1] KOLYMBAS D. A rate-dependent constitutive equation for soils[J]. Mechanics research communications, 1977, 4(6):367-372.
[2] GUDEHUS G.A comprehensive constitutive equation for granular materials[J].Soils and foundations,1996,36(1):1-12.
[3] BAUER E.Calibration of a comprehensive hypoplastic model for granular materials[J].Soils and foundations,1996,36(1):13-26.
[4] HERLE I,GUDEHUS G.Determination of parameters of a hypoplastic constitutive model from properties of grain assemblies[J].Mechanics of cohesive-frictional materials,1999,4(5):461-486.
[5] FUENTES W, TRIANTAFYLLIDIS T. ISA model:a constitutive model for soils with yield surface in the intergranular strain space[J]. International journal for numerical and analytical methods in geomechanics, 2015, 39(11):1235-1254.
[6] WANG G, XIE Y. Modified bounding surface hypoplasticity model for sands under cyclic loading[J]. Journal of engineering mechanics, 2014, 140(1):91-101.
[7] LIN J, WU W. A comparative study between DEM and micropolar hypoplasticity[J]. Powder technology, 2016, 293:121-129.
[8] 岑威钧, 王修信. 堆石料本构建模新途径[J]. 河海大学学报(自然科学版), 2008, 36(1):102-105.
[9] 岑威钧, 王修信, BAUER E,等. 堆石料的亚塑性本构建模及其应用研究[J]. 岩石力学与工程学报, 2007, 26(2):312-322.
[10] 吴长彬, 张岩, 许小东. 应用于堆石坝的弹塑性模型与亚塑性模型的理论对比分析[J]. 灾害与防治工程, 2009(2):8-14.
[11] 岑威钧. 堆石料亚塑性本构模型及面板堆石坝数值分析[D]. 南京：河海大学, 2005.
[12] XIANG B, ZHANG Z L, CHI S C. An improved hypoplastic constitutive model of rockfill considering effect of stress path[J]. Journal of central south university of technology, 2009, 16(6):1006-1013.
[13] 明华军, 孙开畅, 徐小峰,等. 考虑颗粒破碎对特征孔隙比影响的堆石体亚塑性本构模型[J]. 岩土力学, 2016, 37(1):33-40.
[14] 陈泽钦, 刘国明, 林伯江. 堆石料Wang-Wu亚塑性本构模型参数求取的改进[J]. 福州大学学报(自然科学版), 2017, 45(2):199-205.
[15] 王仁超, 曹婷婷, 刘严如. 亚塑性模型不同积分算法的数值实现[J]. 岩土力学, 2017, 38(5):1510-1516.
[16] 唐扬. 亚塑性本构模型及三峡库区塘角滑坡稳定性判据研究[D]. 北京：中国地质大学, 2019.
[17] OSINOV V A, CHRISOPOULOS S, TRIANTAFYLLIDIS T. Numerical analysis of the tunnel-soil interaction caused by an explosion in the tunnel[J]. Soil dynamics and earthquake engineering, 2019，122:318-326.
[18] HLEIBIEH J,HERLE I.The performance of a hypoplastic constitutive model in predictions of centrifuge experiments under earthquake conditions[J].Soil dynamics and earthquake engineering,2019,122:310-317.
[19] Von WOLFFERSDORFF A. A hypoplastic relation for granular materials with a predefined limit state surface[J]. Mechanics of cohesive-frictional materials, 1996, 1(3):251-271.
[20] 傅华, 米占宽,朱秩,等. 福建仙游抽水蓄能电站面板堆石坝坝体填筑料静力特性试验研究[R]. 南京:南京水利科学院, 2004.
[21] 顾淦臣, 沈长松,岑威钧. 土石坝地震工程[M]. 北京: 中国水利水电出版社, 2009.