0.   引言

近年来,中国交通基础设施建设不断向山区延伸。在此区域内,工程性质优良的土质填料较为缺乏,而爆破开采所得粗粒填料料源充足。为了就地取材、缩短工期和节省投资,普遍采用粗粒料（CGMs）填筑山区交通基础设施路基,而山区高填方路基工后变形及其所致路基病害问题已成为当前岩土工程研究热点课题。目前,学者们常用三轴蠕变仪研究粗粒料（包括砂砾、碎砾石、土石混合料和堆石料等）的蠕变行为,并通过非线性拟合构建蠕变经验模型。国内外学者^[1-4]开展了粗粒料流（蠕）变试验,并给出了各类流（蠕）变模型。目前,常用经验模型包括双曲线^[5-6]、指数^[7-8]和对数模型^[9-10]。另外,有些学者^[11-13]建立了考虑环境因素影响的粗粒料流（蠕）模型。在上述经验模型中,最终蠕变常被表达为围压和应力水平的函数^[7-8,14],这只适用于硬化材料,而不适于材料软化阶段和临界状态下蠕变行为描述。

当前计入粗粒料流（蠕）变的本构模型主要包括元件组合模型^[15-16]和弹黏塑性模型^[17-21]两类。前者以流变力学为基础,只能描述某些流（蠕）变现象,例如应力松弛、弹性后效与流动特性等,而不能较好地反映流变本质,不太适合于粗粒料流（蠕）变行为描述。后者以弹塑性理论为基础,耦合流（蠕）变经验表达式,可以较为全面和合理地描述粗粒料的弹性和（黏）塑性力学行为,比较符合粗粒料实际流变特性且便于工程应用。

有鉴于此,笔者首先依据土的蠕变性质分析给出了蠕变、蠕变速率与最终蠕变的表达式；其次,将上述蠕变表达式嵌入弹塑性理论构建粗粒料弹黏塑性本构模型,并采用软岩料三轴蠕变试验结果予以初步验证；最后,采用上述本构模型描述粗粒料本构特性,针对一山区高填方路基实体工程进行数值模拟,分析了从路基填筑至竣工后15 a这个过程中的路基应力、沉降和蠕变变化规律。研究成果可为深入了解粗粒料高填方体工后长期变形规律提供有益参考。

为统一表达,文中所有应力皆为有效应力。

1.   蠕变特性描述

土的蠕变是指土在某一级荷载作用下发生的与时间相关的变形。除去土的自身因素之外,土的蠕变与应力水平、应力路径和加载历史等外部因素相关。在多级加载工况下,相邻两级荷载下土的蠕变与瞬时变形互相隔断而不发生连续累积。对于单级加载,土的蠕变即为总量；对于多级加载,土的蠕变即为增量。

1.1   蠕变函数

目前,通常采用三轴蠕变仪研究土的蠕变特性。已有试验资料^[2,8]表明：在蠕变试验开始后约1 h内,材料发生瞬时变形,之后才有蠕变产生。因此,试验开始后的总时间t′是瞬时变形时间t₀与实际蠕变时间t_c之和。为便于数学描述,定义无量纲时间t为

式中,t₀=1 h。

鉴于指数函数求导较为方便,笔者选用负指数模型描述粗粒料的蠕变规律：

而蠕变关于时间的变化率（蠕变速率）为

式中,${\epsilon }^{\text{c}}$,${\dot{\epsilon }}^{\text{c}}$分别为蠕变与蠕变速率,a,b为材料参数。

依据初始条件和蠕变变量的物理意义,可由式（2）,（3）获得如下推论：

（1）蠕变函数

当$t\to 0$,${\epsilon }^{\text{c}}\to {\epsilon }^{\text{ci}}=\text{0}$；

当$t\to \infty$,${\epsilon }^{\text{c}}\to a={\epsilon }^{\text{cf}}$。

式中,${\epsilon }^{\text{ci}}$,${\epsilon }^{\text{cf}}$分别为初始与最终蠕变应变。

（2）蠕变速率

当$t\to 0$,${\dot{\epsilon }}^{\text{c}}\to {\dot{\epsilon }}^{\text{ci}}=ab=b{\epsilon }^{\text{cf}}$；

当$t\to \infty$,${\dot{\epsilon }}^{\text{c}}\to {\dot{\epsilon }}^{\text{cf}}\to 0$。

式中,${\dot{\epsilon }}^{\text{ci}}$,${\dot{\epsilon }}^{\text{cf}}$分别为初始与最终蠕变速率。

重新整理式（2）,（3）可得

可见,${\dot{\epsilon }}^{\text{ci}}$,${\epsilon }^{\text{cf}}$为描述土的蠕变特性的表征参数。

1.2   最终蠕变

针对土的三轴蠕变试验,目前学者们^[2-4,11]只研究了应变硬化（阶段）的蠕变特性,而未涉及应变软化（阶段）和临界状态,已有测试结果表明：在每级荷载下,最终体积蠕变${\epsilon }_{\text{v}}^{\text{cf}}$和偏蠕变${\epsilon }_{\text{d}}^{\text{cf}}$均随围压${\sigma }_3$与应力水平S_L（S_L=$({\sigma }_1-{\sigma }_3)/{({\sigma }_1-{\sigma }_3)}_{\text{f}}$）增加而增加。其中,${\sigma }_1-{\sigma }_3$为偏应力,${({\sigma }_1-{\sigma }_3)}_{\text{f}}$为峰值偏应力。

然而,对于应变软化土,即使软化阶段应力水平与硬化阶段相同,由于内部状态（可用状态参数ψ=e-e_c^[22]描述,e,e_c分别为土的当前与临界状态孔隙比）不同,土的蠕变必然不同（S_L1=S_L2,ψ₁≠ψ₂,ε^c1≠ε^c2）。

依据临界状态土力学理论,在临界状态下土的体积应变${\epsilon }_{\text{v}}^{}$保持恒定（${\epsilon }_{\text{v}}^{}$=const且$\text{d}{\epsilon }_{\text{v}}^{}$=0）,而偏应变${\epsilon }_{\text{d}}^{}$则会持续发展（${\epsilon }_{\text{d}}^{}$→∞且$\text{d}{\epsilon }_{\text{d}}^{}$→∞）。因此,可以认为：在临界状态下,体积蠕变${\epsilon }_{\text{v}}^{\text{c}}$为0（${\epsilon }_{\text{v}}^{\text{c}}$=0）,而剪切（偏）蠕变${\epsilon }_{\text{d}}^{\text{c}}$趋于无穷（${\epsilon }_{\text{d}}^{\text{c}}$→∞）。

此外,在整个三轴剪切过程中,不同应力水平下的最终蠕变应是连续变化的（即最终蠕变函数连续可导）。由此推知：土的最终蠕变与围压、应力水平和状态参数相关；在某一围压下,最终体积蠕变随应力水平和状态函数的变化规律为先增加而后减小,至临界状态时变为0。

为使最终蠕变${\epsilon }^{\text{cf}}$表达更具普适性（适用于各类应力路径,如各向等压、三轴剪切与纯剪切等）,采用归一化平均主应力$p/p_{\text{a}}$代替${\sigma }_3/p_{\text{a}}$,应力水平$\eta$/M_f代替S_L,则最终蠕变${\epsilon }^{\text{cf}}$的通用表达式为

式中,$\eta$为应力比($\eta$=q/p),p为平均主应力,q为偏应力,M_f为峰值应力比,$p_{\text{a}}$为大气压力。

依据已有三轴蠕变测试结果可知,体积蠕变终值${\epsilon }_{\text{v}}^{\text{cf}}$主要受围压或平均主应力影响较大,而偏蠕变终值${\epsilon }_{\text{d}}^{\text{cf}}$则主要受应力水平影响较大。因此,土在某一级荷载下的最终蠕变可表达为

式中,ψ₀为状态参数初值,k_v,k_d为材料常数。

2.   弹黏塑性本构模型

针对粗粒料变形的时间效应,本次研究只考虑蠕变特性。依据各类变形成因机理,粗粒料总应变被划分为弹性${\epsilon }^{\text{e}}$、塑性${\epsilon }^{\text{p}}$和蠕变${\epsilon }^{\text{c}}$3个组成部分：

定义非弹性应变${\epsilon }^{\text{in}}$为塑性与蠕变之和：

为耦合蠕变与弹塑性应变,笔者基于CSSM（critical sate soil mechanics）临界状态土力学与CBFD（coupling breakage and fiction dissipation）弹塑性本构模型框架^[23],借鉴屈服面流变模型表达方法,在硬化规则中引入时间变量t,使之拓展为适合描述粗粒料蠕变行为的弹黏塑性本构模型。

2.1   弹性

在该模型中,粗粒料弹性特性依然服从广义Hooke定律。在弹性矩阵中,剪切模量G与体积模量K采用Richart等^[24]的表达式描述：

式中,G₀为一材料常数,$\nu$为泊松比（常取0.30）。

2.2   非弹性

在CBFD模型框架^[21]中,塑性势函数g为

而非关联屈服函数f (即f≠g)为

该模型中的相关变量说明如下：

（1）$\beta$为一状态函数($\beta$=exp(m₁ψ)),其中临界状态孔隙比e_c=(e_c,min+(e_c,max-e_c,min)exp(-m(p/$p_{\text{a}}$)^k)。

（2）M_c为临界状态应力比(M_c =nA(p/$p_{\text{a}}$)^n-1)。

（3）H为一新的软化/硬化参数(H= (1-η/M_p)∫(dW^p/ p)),其中,dW^p为塑性功增量,M_p为一可描述软化/硬化特性的虚拟应力比（M_p =${\beta }^{\prime }$M_c）且${\beta }^{\prime }$为另一状态函数(${\beta }^{\prime }$=$\text{exp}(m_2\langle -\psi \rangle )$),符号$\langle \rangle$定义的运算规则：当x≥0时$\langle x\rangle$=x,而x<0时$\langle x\rangle$=0。

在上述表达式中,参数分别是与剪胀/剪缩、软化/硬化、临界状态强度和临界状态孔隙比等剪切特性关联的材料常数（m₁和$\alpha$；m₂,$\lambda$和$\kappa$；n和A；e_c,min,e_c,max,m和k）。

为描述粗粒料的非弹性特性（inelastic behaviors）重新定义塑性功增量为

式中,$\text{d}{\epsilon }_{\text{v}}^{\text{in}}=\text{d}{\epsilon }_{\text{v}}^{\text{p}}+\text{d}{\epsilon }_{\text{v}}^{\text{c}}$且$\text{d}{\epsilon }_{\text{d}}^{\text{in}}=\text{d}{\epsilon }_{\text{d}}^{\text{p}}+\text{d}{\epsilon }_{\text{d}}^{\text{c}}$。

因此,软化/硬化参数H是应力、塑性应变和蠕变的函数,即$H=H(\sigma ,{\epsilon }^{\text{in}})=H(\sigma ,{\epsilon }^{\text{p}},{\epsilon }^{\text{c}})$,其中,$\sigma$为应力张量。

$\text{d}{\epsilon }_{\text{v}}^{\text{p}}$和$\text{d}{\epsilon }_{\text{d}}^{\text{p}}$分别为塑性体积应变与塑性偏应变增量,可由塑性流动法则计算得到：

其中,$\Lambda$为塑性乘子。而$\text{d}{\epsilon }_{\text{v}}^{\text{c}}$和$\text{d}{\epsilon }_{\text{d}}^{\text{c}}$分别为体积蠕变与偏蠕变增量,

将上述系列方程代入式（13）,即可获得考虑粗粒料蠕变效应的黏塑性硬化规则。相关蠕变导致的屈服面硬化规律参见Zhang等^[3]的论文。

式（11）～（13）,（15）和（16）一起构成了粗粒料弹黏塑性本构方程。

3.   三轴蠕变试验验证

为了验证上述弹黏塑性本构模型,针对一山区高填方体所用填料开展了不同围压和不同应力水平条件下的大型三轴蠕变试验研究。该填料为采石场山体爆破开采所得粗粒料,母岩岩性为弱风化石英砂岩。该填料的颗粒级配如图1所示,其级配特征参数包括：d₁₀为1.94 mm,d₃₀为9.97 mm,d₅₀为21.36 mm,d₆₀为28.02 mm,C_u为14.44,C_c为1.83。其他基本物理力学性质：土粒相对密度G_s为2.72,相对密实度D_r为96%,制样干密度为2.10 g/cm³,制样孔隙比e₀为0.295,母岩抗压强度（饱和）为86.9 MPa,压缩模量E_s4,8为93.5 MPa,压缩指数C_c为1.13×10^-1,再压缩指数C_s为1.48×10^-3,渗透系数为0.227 cm/s。

图  1  粗粒料级配曲线

Figure  1.  Grain-size distribution curve of coarse-grained materials (CGMs)

下载: 全尺寸图片 幻灯片

试验设备为国产大型三轴流变仪,最大围压3 MPa。圆柱形试样尺寸为直径×高度=300 mm×600 mm。依据最大填方高度80 m,确定测试围压为300,600,900,1200 kPa,每级围压下采用的应力水平S_L为0.2,0.4,0.6和0.8。

在测试过程中,先施加一定的围压使试样固结；之后保持围压不变,施加偏应力进行试样剪切；至预设应力水平时,保持围压和偏应力不变,对试样进行该级荷载下的蠕变试验,至变形稳定时结束试验,每个试样试验耗时约10 d。

采用自编本构程序对上述三轴蠕变试验进行了数值模拟,粗粒填料弹黏塑性模型参数见表1。

表  1  粗粒料弹黏塑性模型参数

Table  1.  Parameters of elasto-visco-plasticity model for CGMs

参数	弹性常数	CSL常数				CSE常数		剪胀常数		软化/硬化常数			蠕变常数
符号	G0	ec,min	ec,max	mc	k	n	A	m1	α	m2	λ	κ	b	εv0	kv	εd0	kd
取值	153.5	0.10	0.39	0.148	0.63	0.81	2.89	5.63	1.76	2.48	1.56	29.8	0.48	0.062	0.81	0.019	0.38

下载: 导出CSV 

| 显示表格

数值预测与试验测试结果对比分析如下所示。

（1）应力–应变–体变曲线

图2为填料的应力–应变–体变曲线,由图2可知：在本次试验围压范围内,填料呈现较为微弱的软化与剪胀特性；除在较高围压时预测曲线与实测结果有所偏差之外,模型预测结果与三轴试验结果较为吻合。可见,该模型能够较为合理地描述三轴剪切过程中粗粒填料的软化/硬化和剪胀/剪缩特性。

图  2  粗粒料应力–应变–体变预测结果与测试结果对比

Figure  2.  Comparison between model predictions and test results of stress-strain-volume change of CGMs

下载: 全尺寸图片 幻灯片

（2）总体变–时间曲线

图3为填料的总体变–时间曲线,限于篇幅只给出上、下限围压对应结果。由图3可知：在达到预定应力水平后1 h内,填料体变陡然增加,而后其体变随时间缓慢增长并最终趋于稳定；在较高围压时模型预测结果稍大于试验测试结果,而在较低围压时预测结果与测试结果较为吻合。因此,该模型可较好地预测粗粒填料在瞬时加载与蠕变阶段的变形发展规律。

图  3  粗粒料体积蠕变-时间曲线预测结果与测试结果对比

Figure  3.  Comparison between model predictions and test results of volumetric creep of CGMs

下载: 全尺寸图片 幻灯片

4.   山区高填方体工后变形数值模拟

4.1   计算模型

针对一山区高填方实体工程,采用ABAQUS软件建立了最大填方高度处横剖面数值模型（如图4所示）。该路基底部以下为风化基岩,最大填方高度80 m,顶面宽度25 m,底部宽度195 m。为保持边坡稳定,两侧按照上陡下缓放坡,不同高度处边坡坡比如图4所示。为了解高填方体工后沉降规律,在图4所示的20,40,60 m的3个填方高度处设置系列沉降监测点,监测自2017年2月工程竣工时开始,截止目前已持续监测近3 a之久。

图  4  粗粒料高填方体有限单元网格

Figure  4.  Finite element mesh of high embankment filled with CGMs

下载: 全尺寸图片 幻灯片

如图4所示,代表性单元508位于20 m填方高度左侧边坡附近（简称左侧单元）,单元423位于40 m填方高度右侧边坡附近（简称右侧单元）,单元291位于60 m高度填方体中线下（简称轴线单元）。

在数值模拟中,针对不同材料采用不同本构模型描述。其中,风化基岩采用线弹性模型,弹性模量E=5.32 GPa,泊松比$\nu$=0.12；粗粒填料采用ABAQUS二次开发的材料子程序UMAT,模型参数见表2。

4.2   施工进度

图5为施工进度曲线,水平分层（分期）填筑方案如下：填方体分8级填筑,每级填筑高度10 m,施工期为30 d。工期按前期慢而后期快的原则予以安排。每2级填筑及其后施工间歇期为1个施工期,8级填筑分4个施工期完成,例如第①期包括第1和第2级填筑及其施工间歇期。第①、②、③和④期的施工间歇期分别为90,60,30,0 d。因此,第①、②、③和④期结束时的时间为150,270,360,420 d。第④期后进入路面铺筑和运营阶段。

图  5  粗粒料高填方体填筑施工进度曲线

Figure  5.  Construction history of high embankment filled with CGMs

下载: 全尺寸图片 幻灯片

4.3   结果分析

（1）应力

在竣工时和工后15 a,各应力分量空间分布规律相似；在填方体中,填筑高度越大处S₁₁和S₂₂越小,最大值（S₁₁为0.61,0.85 MPa和S₂₂为1.44,1.39 MPa）均位于中心线下填方体底面；S₁₂随着至坡面距离的增加而减小,绝对值的最大值（S₁₂为0.29,0.46 MPa）位于路基边坡坡脚。

由图6可知：在填筑阶段,代表性单元的各个应力分量随填土荷载的施加而持续增加；在工程竣工后,各个应力分量略有波动,在工后5 a后各个应力分量保持稳定。N508,N423和N291的S₁₁,S₂₂和S₁₂分别为0.35,0.57,0.21 MPa；0.18,0.44,0.12 MPa；0.10,0.47,0.08 MPa；N291位于填方体中心线下,因此剪应力S₁₂→0。模型预测所得应力的时空变化规律与高填方体现场实测变化趋势基本吻合。

图  6  粗粒料高填方体代表性单元应力–时间全过程曲线

Figure  6.  Stress history of representative elements in high embankment filled with CGMs

下载: 全尺寸图片 幻灯片

（2）沉降

在竣工时和工后15 a,高填方体总沉降分布类似；最大总沉降(U₂)_max位于路基中心线下约3/4填方体高度处,离最大总沉降位置越远,总沉降量越小。竣工时(U₂)_max为1.35 m,而工后15 a的(U₂)_max为1.81 m,（蠕变所致）工后沉降为0.46m,约占总沉降量的1/4,占填方体高度的1/200以上。

由图7可知：在填筑阶段,代表性单元的总沉降随填土荷载的施加而持续增加；在工程竣工后,总沉降亦继续增长,但变化速率远小于施工阶段,其中N291最大,N423居中,而N508最小,并于工后约1a趋于稳定；至于蠕变所致工后沉降,N508为0.05 m,约占总沉降的9.1%,占填方体高度的0.06%；N423为0.26 m,约占总沉降的21.3%,占填方体高度的0.33%；N291为0.46 m,约占总沉降的25.4%,占填方体高度的0.58%。这表明粗粒料高填方体工后沉降较为显著,与总沉降处于同一数量级。

图  7  粗粒料高填方体代表性单元沉降–时间全过程曲线

Figure  7.  Settlement history of representative elements in high embankment filled with CGMs

下载: 全尺寸图片 幻灯片

由图8可知：在高填方体横断面上,越靠近中心线位置,监测点的工后沉降越大,反之越小；在竣工后近3 a内,对于代表性单元,工后沉降模型预测结果与现场监测结果较为吻合,由此表明该模型预测粗粒料高填方体工后（蠕变）变形具有较高可靠性。

图  8  粗粒料高填方体监测点与代表性单元工后沉降

Figure  8.  Post-construction settlement of monitoring points and representative elements in high embankment filled with CGMs

下载: 全尺寸图片 幻灯片

（3）蠕变

在高填方体内部,竣工时竖向蠕变应变C_E22较小,最大值为0.43%,且从底面至顶面由大到小分布；工后15 a,C_E22增长1个数量级,最大值为2.32%,与最大沉降位置（路基中心线下约3/4路基高度处）一致,其分布围绕最大蠕变点向周边衰减；竣工后15 a内,C_E22增长约1.89%,约占总量的4/5以上。

由图9可知：在工程竣工后,各个代表性单元的竖向蠕变C_E22与其沉降U₂的变化规律相似；N508工后蠕变C_E22增长约0.38%,约占总量的64.4%；N423工后蠕变C_E22增长约1.02%,约占总量的89.5%；N291工后蠕变C_E22增长约2.23%,约占总量的97.8%。以上分析表明：粗粒料蠕变效应较为显著,是相关高填方体工后变形（工后沉降）的内因。

图  9  粗粒料高填方体代表性单元竖向蠕变–时间全过程曲线

Figure  9.  Vertical creep strain history of representative elements in high embankment filled with CGMs

下载: 全尺寸图片 幻灯片

5.   结论

（1）分析土的蠕变特性给出了蠕变与时间的指数型经验模型,同时分析不同应力水平下软化/硬化材料的最终蠕变变化规律并给出了状态相关数学表达式。

（2）借鉴屈服面流变模型表达方法,将蠕变表达式嵌入CBFD弹塑性模型的硬化规则之中,从而构建了适合于描述粗粒料蠕变行为的弹黏塑性本构模型。

（3）对比分析模型模拟与三轴测试结果可知：该模型能够较好地模拟剪切过程中粗粒料的软化/硬化和剪胀/剪缩特性,并预测蠕变阶段变形发展规律。

（4）依据一山区粗粒料高填方体ABAQUS模拟结果分析可知：模型预测所得应力的时空变化规律与高填方体现场实测变化趋势基本吻合；在竣工后近3 a时间内,高填方体沉降的模型预测结果与现场监测结果较好吻合；粗粒料蠕变效应较为显著,是相关高填方体工后长期变形的重要诱因。这再次表明该模型预测工后（蠕变）变形具有较高可靠性。

有待进一步开展粗粒料高填方体工后沉降的现场长期监测工作,以便全面揭示粗粒料高填方体工后长期变形时空变化规律与粗粒料蠕变的本质联系。