卵石土蠕变本构模型及其工程应用

刘 健，唐 扬，易 龙，彭元诚，周跃峰

(1.湖北白洋长江公路大桥有限公司，武汉 430000； 2.长江科学院 水利部岩土力学与工程重点实验室，武汉 430010； 3.中交第二公路勘察设计研究院有限公司，武汉 430000)

1 研究背景

第四系卵石层广泛分布于我国西南片区，区内卵石土的物理力学性质研究，已成为对工程活动最富实际意义的岩土工程研究方向之一[1-3]。针对以卵石土为地基的构筑物来说，根据流变学中“万物皆流动”的原则，在一个时长范围内，时间有可能成为影响卵石土的变形或强度特性的有效因素，甚至成为地基沉降破坏的控制因素，故有必要对砂卵石土的流变特性进行系统的研究，为判断时间因素对工程的实际影响程度提供理论依据。

粗粒土流变的研究主要集中在面板坝中堆石料的粗粒土，通过对高应力状态下堆石料的流变原理进行分析，建立起高应力状态下粗粒土流变本构模型。陈晓斌[4]针对红砂岩粗粒土在不同应力状态下的剪胀性和流变特性进行研究，采用大型三轴试验仪，进行了不同应力状态下的三轴流变试验；周墨臻等[5]针对堆石料流变特性，提出了一种多级应力加载和流变的三轴试验方法，可有效避免流变的起算时间问题；李海芳等[6]对两河口水电站混合料进行了三轴流变试验研究，探讨和分析了流变机理和模型;朱晟等[7]结合某工程堆石料的三轴流变试验资料，提出了适合于土石坝筑坝材料流变特性分析的参数流变模型。此外，采用反分析工程实际观测数据及理论方法来建立流变模型，也是较为常见的方法[8-9]。

目前，在卵石土流变领域的研究还相对较少。鉴于此，本文以湖北省内某大桥锚碇工程下覆卵石土为研究对象，采用室内试验手段，研究其在不同应力状态下的流变特性；同时，根据室内试验结果，建立卵石土的流变本构模型。最终，将该本构模型用于锚碇下覆卵石土地基的长期沉降变形计算。

2 卵石土的蠕变试验

卵石土长期变形特性的研究是基于三轴流变试验开展的。流变是指物体受力变形中存在的与时间有关的变形特性。土体具有流变性，常见的流变现象主要包括蠕变、松弛、流动、应变率效应和长期强度效应等，可采用三轴蠕变试验模型来研究卵石土的蠕变现象，揭示其蠕变规律，建立相应理论及计算方法。

以湖北省内某大桥锚碇工程下覆卵石土为母料配制室内试验需要的试样。由于原始级配的最大粒径较大，室内试验由于仪器尺寸的限制，需要对原始级配进行级配缩尺处理。依据《土工试验方法标准》(GB/T 50123—2019)[10]的粗颗粒土缩尺方法，试验级配采用等量替代法进行级配缩尺，用60～5 mm粒组等量替代>60 mm粒组。原始级配与试验级配如图1所示。

图1 蠕变试验级配曲线Fig.1 Particle gradation for creep tests

三轴蠕变试验采用图2所示的大型高压三轴压缩试验仪。该仪器配有专门稳压的高压蓄能罐，可以提供流变试验需要的稳定轴向压力。试样尺寸为Ф300×H600 mm，围压σ3=0.1、0.2、0.3、0.4 MPa共4级，应力水平SL=0.2、0.4、0.6、0.8共4级，采用应力式加载。试样干密度为1.87 g/cm3，处于饱和状态。试验条件如表1所示。

图2 大型静力应力控制式三轴试验仪Fig.2 Large-scale triaxial test equipment

表1 蠕变试验条件Table 1 Conditions of creep tests

蠕变试验中需要监测轴向变形和体积变形，根据《土工试验方法标准》(GB/T 50123—2019)[10]中粗颗粒土三轴流变试验要求，取稳定标准为每24 h内轴向应变变化量<0.05%。根据大型三轴剪切试验得到的最大偏应力值，按围压计算预施加的各级应力水平下的偏应力竖向荷载，并分级施加在试样上。在已知应力条件下，稳定应力状态若干时间，记录不同时刻试样变形，当变形趋于稳定后施加下一级荷载。试验过程中，为了避免温度变化对试验成果的影响，实验室采用空调进行温度控制，试验过程中温度控制在(20±1)℃左右。

根据以往对土体蠕变特性的试验结果，认为指数型衰减曲线能较为准确地反映土料的蠕变特性。因此，将本次蠕变试验得到的分级加载的轴向应变的剩余应变、体积应变的剩余应变与时间的关系整理在双对数坐标系下(图3和图4)。同时，为区分这两部分应变，对蠕变过程中的轴向应变与时间关系曲线直线段的起始时间进行统计，直线段起始时间的平均值约1 h。为了统一，在整理蠕变试验成果时，两部分应变时间以1 h为界，1 h以前的应变为初始弹塑性应变。

图3 蠕变试验轴向应变的剩余应变与时间的关系Fig.3 Residual strain of axial strain versus time in creep test

图4 蠕变试验体积应变的剩余应变与时间的关系Fig.4 Residual strain of volumetric strain versus time in creep test

从图3和图4可知，轴向应变的剩余应变、体积应变的剩余应变与时间关系在双对数坐标下呈现较好的线性关系。且该卵石土的蠕变是一个逐渐趋于停止的衰减曲线，荷载作用后前期的变形较大，几十个小时就能完成变形的80%。对比分级加载的流变曲线，随着应力水平的增加，蠕变趋于停止的时间逐渐增加。

3 卵石土的蠕变本构模型

3.1 蠕变本构模型

按照滞后变形理论，总应变可以分为瞬时产生的弹塑性应变εep和滞后产生的蠕变εL(t)两部分[11-13]，即

ε=εep+εL(t) 。

(1)

式中εep的计算可以采用任何一种现成的弹塑性模型。

对于滞后黏滞变形，根据上述进行的卵石土蠕变试验，发现蠕变量与时间曲线在双对数坐标系下呈很好的线性关系，“剩余蠕变应变”(εf-εL)与时间曲线在双对数座标系下也呈很好的线性关系，可以采用幂函数表达卵石土的蠕变量εL与时间t的关系,即

εL=εf(1-t-λ) 。

(2)

式中εf相当于t→∞时的最终轴向蠕变量，最终轴向蠕变量εf与应力状态有关。根据式(1)，可将式(2)进行如下变换，即

(εf+εep)-ε=εft-λ。

(3)

根据εep以及不同时间t的应变ε，拟合得到εf、λ。图5为卵石土料不同应力水平的εf与围压的关系曲线，εf与围压有很好的线性关系，且εf与围压成正比,其关系式可表示为

εf=α+βσ3。

(4)

图5 蠕变试验轴向蠕变量εf与围压的关系Fig.5 Axial creep εf versus confining pressure in creep test

图6给出了卵石土系数α、β与应力水平SL之间的相互关系。由图6可知，α与应力水平SL之间可以采用幂函数(式(5))表达，β与应力水平SL之间可以采用双曲线函数(式(6))表达，拟合曲线和试验成果有较好的一致性。

图6 蠕变试验α、β与应力水平SL的关系Fig.6 Relations of α and β versus stress level SL in creep test

(5)

(6)

式中a、b、c、d均为待定系数。

将式(5)、式(6)代入到式(4)，可得到εf与应力状态之间的表达式为

(7)

图7是指数λ与围压的拟合曲线，可以看出λ与围压之间服从幂函数关系，即

图7 蠕变试验λ与围压之间的关系Fig.7 Relations of λ versus confining pressure in creep test

(8)

从体变蠕变余量与时间关系曲线(图4)中可以看出，体积蠕变量εLV的时间曲线可以用幂函数表达为

εLV=εfV(1-tλV) 。

(9)

根据不同时间t的体积蠕变量，拟合得到εfV，且εfV为应力状态的函数。图8为指数λV与围压的关系曲线，我们假定λV为常数。

图8 蠕变试验λV与围压之间的关系Fig.8 Relations of λV versus confining pressure in creep test

图9为不同应力水平下的εfV与围压之间的关系曲线，εfV与围压有很好的线性关系，可以采用线性函数拟合为

图9 蠕变试验εfV与围压之间的关系Fig.9 Relations of εfV versus confining pressure in creep test

εfV=αV+βVσ3。

(10)

图10为αV、βV与应力水平SL之间的相互关系，两者可采用幂函数表达为：

图10 蠕变试验αV、βV与应力水平SL之间的关系Fig.10 Relations of αV and βV versus stress level in creep test

(11)

(12)

将式(11)、式(12)代入式(10)，可得应力状态函数εfV表达式为

(13)

综上所述，参数a、b、c、d、η、m、cα、dα、cβ、dβ、λV可以表达卵石土体轴向应变和体积应变的蠕变特征。该模型可以通过二次开发接口(UMAT)嵌入到ABAQUS软件中进行有限元计算。

3.2 蠕变本构模型在ABAQUS软件中的实现

ABAQUS是一款功能强大的通用有限元软件，由于岩土工程问题涉及到较多的非线性问题，而该软件在这方面有较为突出的贡献，故该软件非常适合用于分析岩土工程问题。尽管ABAQUS中的材料模型非常丰富，但计算岩土体蠕变问题时缺乏较好的本构关系，从而限制了ABAQUS解决该类问题。因此，ABAQUS提供了方便的二次开发接口(UMAT)，可以让用户灵活地创建自定义材料模型。

每一个增量步下，ABAQUS提供应变、应力、应变增量和时间增量信息。采用上述公式，求得每一个增量步下的蠕变应变增量和应力增量，分别更新应力、刚度矩阵及状态变量。本研究使用FORTRAN语言，编写蠕变本构模型的数值实现代码，由于篇幅有限，不再给出程序代码。本文采用提出的11参数本构模型，结合有限元计算方法，用于湖北省内大桥锚碇下覆卵石土地基的长期沉降变形计算。

4 某大桥锚碇地基长期沉降变形计算

4.1 有限元计算模型及参数

某工程大桥锚碇长101.9 m，宽71.5 m。锚碇场区内主要分布有粉质黏土及卵石土，下覆基岩为泥质砂岩。施工预将上部8 m厚的粉质黏土全部开挖，将锚碇建在卵石土层上。在充分考虑基坑施工过程对周边地表、以及基坑施工对持力层压缩可能产生的影响，必须对模型3个纵深方向进行适当的延伸。基于地下基岩的埋深位置，模型深度取50 m。锚碇计算模型网格主要采用八节点的六面体等参元，锚碇共计12 290个单元，节点总数14 779个，土体共计75 456个单元，节点总数82 008个。锚碇及土体的三维有限元网格剖分如图11所示。

图11 锚碇和模型整体的三维网格剖分Fig.11 Three-dimensional mesh division of anchorage and the whole foundation

根据室内试验的结果，提取计算所涉及地基的地层分布及材料参数。卵石土的本构关系采用上述推导的11参数卵石土蠕变模型，相应的参数如表2所示。锚碇部分采用弹性模型计算，参数选取参考《公路钢筋混凝土及预应力混凝土桥涵设计规范》(JTG 3362—2018)[14]规定的计算参数。弹性模量为3.25×104MPa，泊松比为0.2。卵石土和混凝土的密度分别为1.87、2.50 g/cm3。

表2 卵石土蠕变本构模型计算参数Table 2 Calculation parameters of the creep constitutive model for cobbly soil

4.2 计算过程

整个计算包括静力计算和蠕变计算。静力计算大体上分为初始应力平衡、基坑开挖、锚碇基础施工、锚碇主体部分施工、施加锚力等过程。随后针对锚碇进行长期沉降变形计算，即在上述计算过程结束之后，采用蠕变本构模型，计算地基沉降随时间的变化情况。下面选取部分要点进行说明。

(1)开挖：开挖过程在计算过程中为一次性开挖，根据现场施工的实际情况，设置开挖的边坡坡比为1∶1.25，开挖深度为8 m，即将上部粉质黏土层全部挖走，将卵石土层作为地基。

(2)接触面设置：计算过程中涉及到锚碇与地基土体的接触作用，在以下地方设置接触面，如锚碇基础底部与卵石土地基、锚碇基础四周与土体。

(3)运行期过程中会施加拉索拉力，根据相关资料，恒载为302 727 kN，施加于锚碇内部。

4.3 计算结果分析

由于本研究中主要讨论卵石土的蠕变特性，故计算结果主要反映蠕变过程中的变形情况。受锚碇结构影响(前轻后重)，锚碇后趾的变形较大，因此，选取模型后趾的变形与相应位置的监测数据进行对比。从整个运行期变形情况(图12、图13)来看，地基在开始阶段沉降变形速率较大，之后会逐渐趋于平缓，蠕变模型很好地反映了地基土体在蠕变初期和蠕变平稳阶段的变形过程。其中的回弹过程主要是施加拉索拉力之后，导致地基土体回弹。

图12 土体竖向变形与时间的关系Fig.12 Relationship between vertical deformation and time

图13 不同施工节点下竖向位移云图Fig.13 Vertical displacement contours of anchorage foundation in different construction periods

表3为计算结果与监测结果的对比情况。相比监测结果，蠕变初期变形阶段的沉降计算结果偏大，蠕变稳定阶段的沉降计算结果偏小。具体反映在施工结束运行305 d时蠕变位移的增量为5.87 cm，而监测结果仅增加了5.1 cm，计算结果偏大。而之后施加恒载运行30 d时蠕变位移的增量为0.03 cm，而监测结果增加了0.2 cm，计算结果偏小。最后，采用该方法预测了锚碇施加恒载运行3 a以后的沉降值为19.77 cm。

表3 地基沉降监测结果与有限元计算结果对比Table 3 Comparison of foundation settlement between monitoring and finite element calculation

综上所述，该卵石土的蠕变本构模型可以很好地反映卵石土的蠕变特性，且在实际工程中，计算的土体长期变形与监测结果基本一致，该模型可以很好地进行工程应用。

5 结 论

(1)针对卵石土的蠕变特性进行了试验研究，得到了流变量与时间的关系，试验发现蠕变初始阶段的变形较大，变化规律符合幂函数关系式，基于试验结果提出了针对卵石土蠕变特性的11参数流变本构模型。

(2)采用提出的11参数蠕变本构模型，用于某大桥锚碇下覆卵石土地基的长期沉降变形计算。结果显示，锚碇施工结束之后305 d地基土体沉降变形增加了5.87 cm，总沉降量为19.87 cm；恒载运行30 d内的地基土体沉降变形增加了0.03 cm，总沉降量为19.07 cm；恒载运行3 a内的地基土体沉降变形增加了0.73 cm，总沉降量为19.77 cm。与地基沉降监测结果相比，蠕变加速变形阶段的沉降计算结果偏大，蠕变稳定阶段的沉降计算结果偏小。