0.   引言

堆石坝在大坝填筑和蓄水运行过程中，堆石料颗粒破碎会引起流变变形，且流变变形会持续数年甚至数十年之久。目前，一般选用元件模型^[1-2]或经验流变模型^[3]来反映堆石料的长期变形特性，其中又以经验流变模型为主。沈珠江^[3]最早开展了堆石料流变试验的研究，基于西北口灰岩堆石料试验结果提出了3参数经验流变模型，该模型将体积流变量和剪切流变量分开考虑，模型中假定体积流变量仅与围压有关，剪切流变量仅与应力水平有关，同时采用指数函数描述流变速率的衰减过程。后续研究者对该模型进行了各种形式的修正，如郭兴文等^[4]依据水布垭面板堆石坝灰岩堆石料在围压较高时体积流变量呈非线性增长的特性，提出了修正形式的3参数经验流变模型；李国英等^[5]、米占宽^[6]针对糯扎渡心墙堆石坝筑坝所用砂岩和泥质砂岩的堆石料流变试验结果，考虑了围压和剪应力对体积流变的综合影响，提出了包含剪应力相关项的6参数和7参数流变模型；李海芳等^[7]，于浩等^[8]基于九甸峡主堆石料的流变试验结果提出了幂函数形式的7参数流变模型；程展林等^[9]、左永振等^[10]，周伟等^[11]，Zhou等^[12-13]基于水布垭面板堆石坝灰岩流变试验结果，提出了9参数幂函数流变模型；王海俊等^[14-15]提出了基于双曲线型流变速率衰减的4参数流变模型；Fu等^[16-17]考虑了体积流变量与平均主应力和应力比的关系，提出了8参数流变模型；黄耀英等^[18]认为单指数形式的经验流变模型难以准确描述流变速率的衰减过程，提出了组合形式的3参数流变模型。鉴于经验流变模型表达形式简单，在采用有限元法进行数值分析时可采用初应变法，易于编程实现且可以与现有的瞬变模型结合^{[16-17, 19-26]}，因此在土石坝、高填方路基等工程的流变分析中得到了广泛应用。

堆石体的流变在宏观上表现为高接触应力-颗粒破碎和颗粒重新排列-应力释放、调整和转移的循环过程，期间伴随着能量耗散，经验流变模型只是从量值上反映了堆石料的流变变形，而没有阐释堆石料发生流变的机理。对红石岩堰塞料开展了大型三轴压缩和流变试验，分析了流变过程中的能量耗散，提出了基于颗粒破碎耗能的堆石料流变模型，与现有的经验流变模型相比，本文提出的流变模型不仅反映了流变变形的机理，而且具有参数少、物理意义明确等优点。

1.   堆石料大型三轴压缩试验与流变试验方案

1.1   大型三轴压缩试验

试验所用堆石料取自红石岩堰塞料，岩性为白云质灰岩，现场堰塞料级配较宽且不均匀，孤石最大粒径可达数米到数十米，剔除掉大尺寸孤石后统计的级配如图 1所示，最大粒径为0.4 m。试样尺寸Ф300×700 mm，采取平均线级配缩尺后进行室内试验，制样干密度为2.20 g/m³，对应的相对密度为0.95，孔隙率为22%。试验围压分别为300，600，1000，1500 kPa。试验按照《土工试验方法标准》GB/T50123—2019进行。

图  1  堰塞料现场和试验级配曲线

Figure  1.  Field and test grain-size distribution curves

下载: 全尺寸图片 幻灯片

1.2   流变试验

根据三轴压缩试验得到的破坏时的偏应力${({\sigma _1} - {\sigma _3})_{\rm f} }$，每级围压下分别开展应力水平SL为0.0，0.2，0.4，0.8的流变试验（SL=$({\sigma _1} - {\sigma _3})/{({\sigma _1} - {\sigma _3})_{\rm f} }$），以24 h轴向应变增量小于总流变量的1‰作为停机标准。

2.   堆石料大型三轴压缩与流变试验结果

2.1   大型三轴压缩试验

试验得到的应力应变曲线如图 2所示，应力应变曲线基本呈现出应变硬化型，仅在300 kPa低围压下出现轻微的软化现象，中高围压下均为应变硬化型。300，600 kPa低围压下出现剪胀现象，高围压下则始终处于剪缩状态。根据试验结果整理得到的破坏应力比M_f见表 1，M_f=(q/p)_f，其中q，p分别为偏应力和球应力。M_f随围压升高逐渐降低，也即峰值摩擦角随围压逐渐降低。可见，堰塞料应力应变特性与其它堆石料类似。

图  2  三轴压缩试验应力应变关系曲线

Figure  2.  Stress-strain curves of test materials

下载: 全尺寸图片 幻灯片

表  1  流变试验结果

Table  1.  Creep strains of tests

σ3/kPa	Mf	SL	η=q/p	t/d	εvf/%	εsf/%
300	1.87	0.0	0.000	4.09	0.341	0.000
		0.2	0.801	5.89	0.401	0.184
		0.4	1.178	6.76	0.458	0.409
		0.8	1.705	8.27	0.494	1.593
600	1.79	0.0	0.000	4.76	0.429	0.000
1000	1.73	0.0	0.000	5.12	0.532	0.000
		0.2	0.667	6.50	0.614	0.188
		0.4	1.112	8.08	0.674	0.546
		0.8	1.596	10.40	0.747	1.843
1500	1.69	0.0	0.000	6.25	0.635	0.000
		0.2	0.668	7.19	0.722	0.218
		0.4	1.035	8.30	0.788	0.552
		0.8	1.530	11.88	0.861	1.964

下载: 导出CSV 

| 显示表格

2.2   流变试验

流变试验结果见表 1，表中列出了每一级围压σ₃及相应应力水平SL对应的流变时间t、最终体积流变量${\varepsilon _{{\text{vf}}}}({\varepsilon _{{\text{vf}}}} = {\varepsilon _{{\text{1f}}}} + 2{\varepsilon _{{\text{3f}}}})$和最终剪切流变量${\varepsilon _{{\text{sf}}}}({\varepsilon _{{\text{sf}}}} = {\varepsilon _{{\text{1f}}}} -$ ${\varepsilon _{{\text{vf}}}}{\text{/}}3)$，其中围压为600 kPa时仅开展了SL=0的流变试验。从表 1中可以看出，围压相同时，体积流变量和剪切流变量随着应力水平增大逐渐增大，应力水平相同时，体积流变量和剪切流变量随着围压增大逐渐增大；图 3给出了围压1000 kPa下流变时间过程，从图 3可以看出，随着流变时间增加，流变速率逐渐降低，流变变形最终趋于稳定，且应力水平越高流变变形持续的时间也越长。将围压300，1000 kPa下三轴压缩试验曲线与流变试验曲线绘制在同一张图中，如图 4所示。从图 4中可以看出：流变试验过程中应变方向（图中蓝色箭头）与三轴压缩试验过程中应变方向（图中黑色箭头）明显不同，体积流变方向出现转折，与静力加载时低围压下会出现先剪缩后剪胀的趋势不同，不同围压和应力水平下流变试验过程中体积流变均处于持续剪缩状态，这一结论与Fu等^[16]、Zhang等^[25]的研究结论一致。

图  3  轴向流变和体积应变与时间关系曲线(1000 kPa)

Figure  3.  Axial creep strains and volumetric creep strains against time at confining pressure of 1000 kPa

下载: 全尺寸图片 幻灯片

图  4  三轴压缩与流变试验结果对比(300，1000 kPa)

Figure  4.  Comparison of triaxial compression test and creep test results at confining pressures of 300 and 1000 kPa

下载: 全尺寸图片 幻灯片

3.   堆石料流变过程中的应变能分析

目前，在堆石料的流变研究中，颗粒破碎被认为是产生流变的主要原因^{[5-6, 11, 27]}，但是目前常用的经验流变模型，由于直接拟合试验数据，不能反映流变产生的机理，因此本节从能量的角度，首先分析流变过程中的输入能与应力条件的关系，然后结合流变过程中颗粒破碎耗能，建立一个新的流变模型。

3.1   流变应变能及其与围压和归一化应力比的关系

将流变试验过程中应力与应变的乘积定义为输入能，并将流变全过程累积的输入能称为流变应变能W_creep：

由于流变过程中平均主应力p、广义剪应力q保持不变，且流变应变表现为剪缩，因此式(1)可表达为

从式（2）可知，可以从能量分析的角度，分析流变应变能的变化规律。采用式（2）计算的红石岩堰塞料流变过程中的应变能W_creep与围压${\sigma _3}$和归一化应力比η/M_f的关系如图 5（a）所示，其中η为应力比，表达为η=q/p。可以看出，当η/M_f不变时，流变应变能W_creep随着围压${\sigma _3}$增加逐渐增加；当围压不变时，流变应变能W_creep随着η/M_f增加而增加。为验证其普适性，图 5中尚分别给出了某抽水蓄能水电站花岗岩堆石料^[16]、茨哈峡砂岩堆石料^[28]、新疆大石峡砂砾石面板坝砂砾料^[27]流变应变能W_creep与围压${\sigma _3}$和归一化应力比η/M_f的关系图，虽然堆石料母岩岩性、级配和孔隙率不同，但W_creep与${\sigma _3}$和η/M_f的关系是一致的。堰塞料、花岗岩堆石料、砂岩堆石料、砂砾料的W_creep与${\sigma _3}$, η/M_f的关系均落到类似的平滑曲面上，曲面随着${\sigma _3}$, η/M_f的增大而不断抬升。因此在后续内容中，根据图 5中试验数据建立起W_creep与${\sigma _3}$, η/M_f的关系式。

图  5  流变应变能与围压和归一化应力比之间的关系

Figure  5.  Relationship among creep strain energy and confining pressure and normalized stress ratio

下载: 全尺寸图片 幻灯片

3.2   流变应变能的分解

分析式（2）中流变应变能W_creep的计算公式，可以将其分为两部分：①平均主应力引起的应变能W_p= pε_vf；②广义剪应力引起的应变能W_q=qε_sf。当应力水平SL=0时，p=${\sigma _3}$，定义此时的流变应变能为固结流变应变能${W_{{\sigma _3}}}$。固结流变应变能${W_{{\sigma _3}}}$与围压${\sigma _3}$的关系如图 6所示，二者之间关系可采用幂函数的形式描述：

图  6  固结流变应变能与围压之间的关系

Figure  6.  Relationship between consolidation creep strain energy and confining pressure

下载: 全尺寸图片 幻灯片

式中：c₁, n₁为拟合参数，分别为拟合曲线的截距和斜率。

随着应力水平的增大，平均主应力p和广义剪应力q随之增大，此时的流变应变能W_creep与${W_{{\sigma _3}}}$的比值表示了应力水平的提高对流变应变能的影响。流变应变能与固结流变应变能的比值W_creep/${W_{{\sigma _3}}}$与归一化应力比η/M_f的关系如图 7所示，堰塞料、花岗岩堆石料、砂岩堆石料和砂砾料等4种材料的W_creep/${W_{{\sigma _3}}}$与η/M_f的规律一致，随着η/M_f的增加，W_creep/${W_{{\sigma _3}}}$从1.0迅速增加，可以用双曲线公式描述二者间的关系：

图  7  流变应变能和固结流变应变能的比值与应力比之间的关系

Figure  7.  Relationship between ratio of creep strain energy to consolidation creep strain energy and stress ratio

下载: 全尺寸图片 幻灯片

式中：c₂为拟合参数，表示W_creep/W_σ3增加的速率。

将式（3）代入式（4），可得流变应变能W_creep的完整表达式：

式中：${c_1}$, ${n_1}$, ${c_2}$为参数。式（5）反映了流变应变能与围压和归一化应力比的相关关系。

4.   堆石料的流变剪胀方程和流变模型

4.1   流变剪胀比

定义流变过程中体应变增量与轴应变增量比值为切线体积比μ_t，

切线体积比μ_t可反映流变过程中的应变路径，下面结合图 4的流变试验进行说明，把图 4中流变体变-轴变曲线$({\varepsilon _{\rm v}}{\text{ - }}{\varepsilon _1})$绘制到同一张图中，如图 8所示。

图  8  围压为300，1000 kPa时流变应变路径

Figure  8.  Creep strain paths at confining pressures of 300 and 1000 kPa

下载: 全尺寸图片 幻灯片

从图 8中可看出，μ_t表示了体变-轴变曲线的斜率，流变过程中的应变路径较为复杂，流变路径不是一成不变的，即切线体积比μ_t是变化的。在应力水平SL=0时，试样处于等向固结状态，此时无剪应变产生，切线体积比恒定为3.0；当围压相同时，随着应力水平的增大，切线体积比逐渐减小，曲线上抬。

从应力水平为0.2，0.4，0.8的流变试验结果来看，此时的切线体积比并不为恒定值，其可能的原因是流变试验过程中试样端部的横截面面积是不断变化的，但受目前试验水平的限制，流变试验采用的是应力控制，也即试验过程中轴向荷载并不会随试样端部横截面面积的变化而变化，已不处于严格的应力恒定状态^[29]。此外，由于颗粒尺寸效应，内部组构的演化规律等都会对流变路径产生影响，以目前的研究水平还难以准确描述这一过程^{[29, 30-31]}。因此在目前的流变研究中，一般假定流变路径是不变的，即认为不同应力水平下切线体积比是恒定值。傅中志等^[29]、王占军等^[32]以流变结束时体应变值与轴应变值的比值作为切线体积比${\mu _{\rm t}}$的近似值，式（6）可化简为

采用式（7）计算的流变应变路径如图 8中黑色虚线所示，假定流变过程中流变应变方向不变是可以接受的。

假定流变过程中应变不可恢复，流变应变全部为塑性应变。切线体积比μ_t与剪胀比d_g具有如下的关系：

分析式（8）可知，当SL=0，此时试样处于等向固结流变状态，剪应变为零，d_g=∞；随着应力水平SL逐渐增大，应力路径逐渐向破坏线靠近，剪应变的增量远大于体积应变，此时流变剪胀比d_g逐渐趋于0。整理图 4中堰塞料在不同围压和应力水平的剪胀比d_g与应力比η的关系如图 9所示（应力水平SL=0时的剪胀比为∞，故在图 9中舍弃不予绘制）。从图 9中可知，当应力比较小时流变剪胀比远大于三轴压缩试验的剪胀比，随着应力比增大，流变剪胀比逐渐减小，此外与三轴压缩试验的剪胀比不同，不同应力水平下的流变剪胀比均大于0，也即流变过程中始终处于剪缩状态，这一规律与Zhang等^[25]、傅中志等^[29]、王占军等^[32]和Lade等^[33]的研究结论一致。

图  9  三轴压缩和流变试验的剪胀比与应力比的关系对比

Figure  9.  Comparison of relationship between dilatancy ratio and stress ratio in triaxial compression tests and creep tests

下载: 全尺寸图片 幻灯片

4.2   基于破碎耗能的流变剪胀方程

堆石料在静力加载过程和流变过程均产生了颗粒破碎。因此堆石料在流变过程中的剪胀方程，仍可采用静力加载过程中剪胀方程的推导方式^[34-35]。

考虑颗粒破碎耗能的ROWE^[35]剪胀方程为

式中：M为临界状态应力比；dE_b为颗粒破碎耗能增量。

式（9）中等号左边可视为流变过程中的输入能，等号右边可视为流变过程中能量耗散，由于dE_b的表达式未知，由式（9）不能直接求得剪胀方程。整理式（9）得到颗粒破碎耗能增量dE_b：

采用式（10）计算图 8中围压1000 kPa下流变试验过程中颗粒破碎耗能增量，得到颗粒破碎耗能E_b与轴向应变${\varepsilon _1}$的关系如图 10所示。

图  10  颗粒破碎耗能与轴向应变之间的关系

Figure  10.  Relationship between particle breakage energy consumption and axial strain

下载: 全尺寸图片 幻灯片

定义图 10中E_b-ε₁曲线的斜率为颗粒破碎耗能因子A^*，A^*=dE_b/dε₁，A^*表示流变过程中，产生单位轴向应变所需的颗粒破碎能耗。从图 11可知，随着应力水平的提高，A^*逐渐降低，产生这一现象的主要原因在于高应力水平下试样已经接近破坏状态，此时试样单位轴向应变下变形所需的能量降低。由于流变过程中颗粒破碎耗能E_b与轴应变${\varepsilon _1}$相关系数较高，可采用线性关系表达。在三轴试验条件下，dE_b和dε₁可分别表达为

图  11  A^*/(3+M)/p与η的关系图

Figure  11.  Relationship between A^*/(3+M)/p and η

下载: 全尺寸图片 幻灯片

式（11）代入式（9），可以得到颗粒破碎耗能因子A^*及剪胀方程d_g表达式：

在同一围压和应力水平下的流变试验，平均主应力p、应力比η和剪胀比d_g均保持不变，因此在流变过程中A^*亦为定值。也即流变过程中颗粒破碎耗能E_b与轴应变${\varepsilon _1}$呈线性关系。

由于应力水平SL=0时，η=0，d_g=∞，可以得出此时A^*值：

整理图 5中4种堆石料的A^*与应力比η的关系如图 11所示，在图 11中对A^*进行了归一化处理。

由图 11看出，堰塞料、花岗岩堆石料、砂岩堆石料、砂砾料等4种材料归一化后的A^*/(3+M)/p随着应力比η的升高而逐渐降低，呈现出非线性下降的趋势。在应力比η=0时，A^*/(3+M)/p=1.0；随着应力比η逐渐增大，A^*/(3+M)/p的值逐渐降低，趋近于0。因此通过以下的双曲线公式来描述A^*/(3+M)/p与η的关系：

式中：c₃，n₃为拟合参数。c具有明显的围压依赖性，随着围压升高，c值逐渐降低，c的大小控制着拟合曲线的弯曲程度。

采用式（14）对图 11中A^*/(3+M)/p拟合，拟合结果如图 11中虚线所示。式（14）可以反映出A^*/(3+M)/p随着η增加而降低的趋势。

利用式（14）可计算得到颗粒破碎耗能因子A^*，代入剪胀方程式（12）计算得到的剪胀比与试验数据对比如图 12所示。从图 12可看出，式（12）计算出的剪胀方程可以体现出剪胀比随应力比升高而由无穷大逐渐降低的趋势；同时，由于引入式（12），（14）计算出的剪胀方程可以反映围压的影响。

图  12  不同围压下流变试验的剪胀比计算结果与试验结果对比

Figure  12.  Comparison between calculated and test results of dilatancy ratio under different confining pressures

下载: 全尺寸图片 幻灯片

4.3   流变模型及验证

联立式（2），（5），（8）可以求得最终的体积流变量${\varepsilon _{\rm vf}}$与剪切流变量${\varepsilon _{\rm sf}}$表达式为

由式（12），（14），本文所建立的流变模型包含c₁，n₁，c₂，c₃，n₃等5个参数，与传统的堆石料经验流变模型相比，模型参数均具有明确的物理意义。其中c₁，n₁用于反映围压${\sigma _3}$对固结流变应变能W_σ3的影响；c₂用于反映流变应变能W_creep与固结流变应变能W_σ3二者之间的比例关系；${c_3}$，${n_3}$用于反映剪胀比d_g与应力p和应力比η的关系。整理图 5中4种堆石料流变模型参数如表 2所示。采用表 2中参数，利用式（16）对4种堆石料体积流变量${\varepsilon _{\rm vf}}$与剪切流变量${\varepsilon _{\rm sf}}$进行计算，计算结果如图 13，14所示。

表  2  流变模型参数汇总

Table  2.  Parameters of creep model

堆石料	c1	n1	c2	c3	n3	M
堰塞料	1.3864	3.0899	1.4560	0.1960	-0.022	1.52
花岗岩	1.4199	2.6800	1.5390	0.3014	-0.072	1.60
砂岩	1.4307	2.1930	1.6283	0.6676	-0.425	1.51
砂砾石	1.3438	2.1841	0.8968	0.7880	-0.296	1.57

下载: 导出CSV 

| 显示表格

图  13  体积流变计算结果与试验结果对比

Figure  13.  Comparison between calculated and test results of volumetric creep strain

下载: 全尺寸图片 幻灯片

图  14  剪切流变计算结果与试验结果对比

Figure  14.  Comparison between calculated and test results of deviatoric creep strain

下载: 全尺寸图片 幻灯片

由图 13，14，本文所提流变模型可较好的拟合流变试验结果，可以反映出体积流变量${\varepsilon _{{\text{vf}}}}$与剪切流变量${\varepsilon _{{\text{sf}}}}$对围压和应力依赖性，且可以反映${\varepsilon _{{\text{vf}}}}$与${\varepsilon _{{\text{sf}}}}$随着围压和应力比升高而升高的现象。从拟合效果来看，虽然对剪切流变的拟合效果良好，但是，对于低围压下较高应力水平时体积流变的拟合效果较差。可能的原因有：①流变试验过程中试样体积变化，试样内部不处于严格的恒定应力状态；②堆石料流变的体变测量一般采用排水法测量，流变体变量值较小，且流变持续时间较长导致体变测量结果精度较差。但是，本文提出的流变模型避免了以往经验流变模型的缺点；从能量耗散的角度解释了流变发生的机理；参数较少，且均具有明确的物理意义；拟合结果总体上与试验结果规律一致，为流变模型的研究开辟了新的研究思路。

4.4   流变数值计算方法

与经验流变模型中^[3-6]流变过程中流变应变与流变时间的关系类似，流变应变能与时间可以建立如下的衰减关系式：

式中：α为流变速率参数，控制着应变能发展速率的快慢。

对式（16）求导，结合指数函数性质，得到时间t时刻的应变能速率：

对于t时刻应变能增量表达式为

式（18）代入流变剪胀方程d_g，得到t时刻的体积应变增量和剪切应变增量表达式：

式（19）仅具有理论意义，在数值计算时需要划分具体的时间步开展计算。取时间t时刻后，增量时间步为Δt内，流变引起的体积应变增量和剪切应变增量分别是Δε_vt，Δε_st：

采用Prandtl-Reuss流动法则，在Δt内流变应变增量以张量形式表示为

式中：I=[1 1 1 0 0 0]^T, s为偏应力张量，q为广义剪应力。

由式（16）~（21）的分析可知，本文所提的流变模型在数值编程时与常规的经验模型最主要的区别在于，不需要分别编程流变体积应变和流变剪切应变与时间的关系函数（即不需要分别建立式（15）中应变与时间的关系函数），而是从流变应变能公式和剪胀方程入手，在具体计算时，增量体积应变和增量剪切应变通过增量形式的流变应变能，并结合剪胀方程分配得到，最后结合Prandtl-Reuss流动法则开展具体的计算。

5.   结论

针对红石岩堰塞坝堆石料开展了三轴压缩试验和流变变形试验，同时分析了花岗岩堆石料、砂岩堆石料、砂砾料等不同类型堆石料的流变试验成果，得到3点结论。

（1）分析了流变应变能与围压和应力比的相关关系，发现固结流变应变能与围压之间关系可采用幂函数的形式表达，以及流变应变能和固结流变应变能的比值与应力比呈非线性相关关系，据此给出了流变应变能的表达式。

（2）分析了流变过程中体积应变与轴向应变的相关关系，发现不同围压和不同应力水平下流变过程中体积流变均表现为剪缩，且切线体积比保持不变，且应力水平为0时剪胀比为无穷大，随着应力水平的提高剪胀比逐渐降低。

（3）在上述研究的基础上，基于考虑颗粒破碎耗能的ROWE剪胀方程，推导得到了最终体积流变和剪切流变的表达式，与现有的经验流变模型相比，本文提出的流变模型反映了堆石料的流变机理，具有较大的理论价值，且在进行数值分析时仍可采用初应变法，易于编程实现，具有较大的推广价值。