黄河宁蒙河段黄土蠕变特性与分数阶本构模型研究

2023-07-21许旭兵张帆侯佼建张敏杨浩明董成会

华北水利水电大学学报(自然科学版) 2023年3期

许旭兵, 张帆, 侯佼建, 张敏, 杨浩明, 董成会

(1.黄河水利委员会黄河水利科学研究院,河南郑州 450003; 2.同济大学土木工程学院,上海 200092; 3.河南新华五岳抽水蓄能发电有限公司,河南信阳 465450)

黄河宁蒙河段地区人工黄土边坡的破坏和失稳、地基和路基缓慢的竖向变形、渠道水工建筑物地基和土石坝的变形和失稳等工程问题会影响基础设施的正常使用,甚至威胁人民生命和财产安全。这些工程问题具有时间累积效应,在长期附加荷载、干湿循环、冻融、加卸载等外界作用下,土体极易发生蠕变破坏。因此,有必要对宁蒙河段黄土蠕变特性进行研究,为黄土地区地基沉降预测和防控提供参考。

近年来,国内外学者进行了大量的室内土体蠕变试验研究,从无围压的单向控制蠕变试验到可控围压的三轴蠕变试验[1-8]。葛苗苗等[2]使用高压固结仪对不同压实度黄土进行了一维蠕变试验,研究发现:压实度越小,蠕变占总变形的比例越大;随着应力水平的增大,蠕变占总变形的比例减小。罗汀等[3]对重塑黄土开展了一系列的一维蠕变试验,研究表明:蠕变特性曲线分为减速变形和稳态蠕变两个阶段,各级应力条件下的稳态蠕变段曲线平行,且次固结系数基本保持不变。王鹏程等[4]研究了90%压实度下重塑黄土的三轴蠕变特性,研究表明:黄土总体呈现非线性衰减蠕变,围压和含水率均对蠕变特性具有显著影响。CETIN H和GOKOG L A[8]研究了不同排水条件下土体三轴蠕变特性,发现:不排水蠕变试验中,土样会在高应变水平下发生破坏;而排水蠕变试验土样在较低应变水平下就会发生破坏。马林等[9]为研究加筋土的力学特性,对加筋黄土进行了固结不排水三轴试验,分析了不同筋材、不同围压和不同加筋方式等影响下土体的力学特性。

学者们主要从经验蠕变模型、屈服面蠕变模型和元件蠕变模型3个方面,对土体蠕变本构模型开展研究。经验模型具有形式简单、模型参数易获取等特点。常用的经验模型有幂函数型、双曲线型、指数型及多项式型等。如MEI G X和YIN J H[10]采用对数函数建立了固结耦合模型,该模型对土体长期固结试验结果的拟合效果较好。余湘娟等[11]建立了双曲线形式的经验模型,通过实例计算表明,该模型适用于不同荷载作用下的次固结沉降预测。屈服面模型主要研究屈服面、关联准则及硬化定律三要素随时间的变化规律。YIN J H和GRAHAM J[12]通过剑桥模型修正和黏塑性理论耦合建立了各向同性一维弹黏塑性模型,该模型可用于描述土体加速蠕变、加卸载和松弛等情况。MARANHA J等[13]基于下屈服面理论,建立了黏塑性模型,其可用于描述黏土蠕变特性,该模型能够描述屈服面内部产生的黏塑性应变的问题。元件模型通常是由虎克弹性体、牛顿黏性体、圣维南塑性体等元件通过串并联方式组合而成。付艳斌和杨志银[14]对软黏土进行三轴蠕变试验,并采用五元件黏弹性模型描述软黏土在弹性范围内的蠕变特性。XIAO H B等[15]建立了多元件组合的黏弹塑性蠕变模型来描述膨胀土与红黏土的蠕变特性,研究表明:使用的串并联元件越多,试验结果的拟合效果越好。

传统元件蠕变模型虽能通过元件累加提高精度,但是其对土体蠕变试验曲线中的加速蠕变阶段不能全面、精确地描述[16]。研究发现,在传统元件模型中引入分数阶导数元件模型能够很好地解决上述问题[17]。殷德顺等[18]认为,土体介于理想固体和流体之间,将分数阶微积分算子引入到元件模型,建立了分数阶蠕变模型,其与多元件模型进行对比发现,分数阶模型使用较少的参数就能达到很好的模拟效果。刘林超等[19]对经典三元件模型和分数阶三元件模型进行对比发现,分数导数流变模型具有参数少,能在较宽范围内描述软土的流变行为。何利军等[20]针对传统元件模型具有参数多、物理意义不明确和适用性差的缺点,将分数阶偏微分理论引入到Burger模型中,有效地描述了湛江软黏土的蠕变特性。张金良等[21]考虑龄期、压实度以及掺量的影响,对新型固化剂改良的黄土进行系列抗冲刷试验研究,结果表明,新型固化剂改良的黄土抗冲刷性能随压实度的增大而增强,且呈非线性关系。黄土非线性蠕变特性对研究黄土地基变形至关重要,但是目前使用分数阶蠕变本构模型对黄土蠕变特性的研究较少。

使用计算机仿真技术研究岩土力学特性已成为当下的热点问题。陈建永[22]将改进后的Drucker-Prager模型编写成了UMAT子程序,通过对上海某车站基坑开挖数值模拟,验证了该模型的适用性。XU X B和CUI Z D[17]建立分数阶蠕变模型模拟土体三轴蠕变特性,并利用FORTRAN语言对该模型进行编译,通过计算机仿真技术成功模拟了土体三轴蠕变试验。

本文结合黄河宁蒙河段渠道水工建筑物黄土地基和渠堤路基沉降的实际情况,对宁蒙河段黄土开展了三轴蠕变试验,着重分析了偏应力水平和围压对黄土蠕变特性的影响,构建了黄土分数阶蠕变本构模型,并将其运用到计算机仿真技术以模拟黄土蠕变。研究成果可为黄土地区地基长期变形量的确定提供参考。

1 黄土土样性质和三轴试验方案

1.1 黄土基本物理力学参数

试验土样取自宁夏回族自治区吴忠市利通区秦渠渠堤,土体为全新世Q4坡积黄土。将取回的土体去除树枝、树叶、植物根系等杂质,风干、碾碎并过2.0 mm筛后,进行土体基本物理力学参数测定,结果见表1。

表1 黄土基本物理力学参数

图1为土体颗粒分析试验结果。从图1中可知,试验所用黄土的颗粒主要以细粒为主,其中粒径为0.005 mm至0.050 mm之间的颗粒占较大比重且含量超过60%,而小于0.005 mm的黏粒含量小于12%,说明该区黄土主要为粉粒,并含有一定量的粉砂粒。

图1 黄土颗粒分布曲线

1.2 三轴试验方案

采用长春机械科学研究院研制的三轴试验设备对黄土进行三轴蠕变试验。土体三轴流变试验机主要由压力室、轴向加压设备、围压施加系统、体积变化和孔隙水压力量测系统组成。与传统的三轴仪不同,土体三轴流变试验机增加了计算机控制与分析系统。

本试验采用分级加载,对同一个土样施加不同的应力,每一级荷载进行线性加载,加载时间为5 min,当达到预定时间、产生预定变形或达到稳定状态(每小时轴向变形小于0.005 mm,视为达到稳定状态[21])时,开始施加下一级荷载。试验方案见表2。

表2 黄土三轴蠕变试验方案

2 试验结果分析及数据处理

2.1 试验结果分析

在不同围压下(SS01、SS02、SS03),黄土分级加载排水蠕变的应变时程曲线如图2所示。由图2可知:围压越大,同一级荷载作用下产生的变形越小;在同一围压下,随着偏应力的增加,每一级的应变增量和瞬时应变也随之增大。其中:围压为75 kPa(SS01)时,在第1级偏应力作用下,产生的瞬间变形为0.302 mm;在第2级偏应力作用下,产生的瞬间变形为1.185 mm;在第3级偏应力作用下,产生的瞬间变形为1.253 mm。不同围压下黄土的蠕变规律基本相似,蠕变曲线具有非线性特性。此外,在偏应力施加初期,轴向应变变化较快,随后慢慢趋于稳定;当偏应力增加时,轴向应变随之增大,稳定所需要的时间也随之增加。因此,该黄土具有典型的非线性蠕变特性。

图2 不同围压下分级加载试验数据

施加到一定荷载后,土体产生较大瞬间位移,三轴流变仪停止工作,默认此时土体发生破坏,破坏时土体呈腰鼓状,如图3所示。在围压75 kPa作用下,附加荷载增加到第4级时,土体发生破坏;在围压150 kPa作用下,附加荷载增加到第6级时,土体发生破坏。

图3 土体破坏时试样形态

2.2 试验数据处理

蠕变试验中,分级加载的荷载曲线转化为相似曲线簇的处理方法有Boltzmann线性叠加法和陈氏法。Boltzmann线性叠加法假定试样为线性流变体,通过坐标平移法进行叠加处理[23-24];陈氏法考虑非线性黏弹性材料在流变过程中产生的非线性蠕变累积后的影响,通过适当的作图法进行叠加处理[25]。本研究使用陈氏法对分级加载试验数据进行处理分析,处理结果如图4—6所示,其中,SS01-1、SS01-2、SS01-3表示加载的第1阶段、第2阶段、第3阶段。

图4 围压75 kPa陈氏法叠加后的分级加载曲线簇

图5 围压150 kPa陈氏法叠加后的分级加载曲线簇

图6 围压225 kPa陈氏法叠加后的分级加载曲线簇

在低应力水平作用下,黄土的蠕变应变随时间发展逐渐趋于稳定,曲线平缓。在变形前,表现为应变变形速率衰减型的蠕变特征;在变形后期,表现为明显的稳态蠕变特征,达到稳定蠕变阶段所需的时间短。在高应力水平作用下,黄土的蠕变变形在变形前期同样表现出变形速率衰减型的蠕变特性,在变形后期表现出稳态蠕变特征,黄土达到稳态蠕变阶段所需的时间越来越长。

图7为SS02-3处理后数据的轴向应变曲线和曲线变化率随时间的变化情况,试样在施加荷载2个小时内达到稳态蠕变状态。

图7 SS02-3轴向应变曲线和曲线变化率

3 黄土分数阶蠕变模型和参数确定

3.1 分数阶本构模型

Merchant模型如图8(a)所示,该模型为:

图8 模型示意图

(1)

式中:E1和E2为压缩模量;η为黏滞系数;σ为有效应力;t为时间。

分数阶Merchant模型是用非线性Abel黏壶元件代替Merchant模型中的Newton黏壶建立含非线性元件的蠕变模型,如图8(b)所示。分数阶蠕变模型为:

(2)

式中:k为常数;α为分数阶次。

为了验证分数阶Merchant模型的优越性,同时选用经典Merchant模型和Burgers模型对同一组试验结果进行拟合,Burgers模型的本构关系为:

(3)

式中η′为Burgers模型中黏壶的黏滞系数。

使用3种模型分别对同一组试验数据(SS01-1)进行拟合,拟合结果见表3,拟合效果如图9所示。Merchant模型和Burgers模型的拟合优度R2分别为0.841 7和0.939 7,说明Burgers模型的拟合度优于Merchant模型。但与传统模型相比,分数阶Merchant模型具有更高的拟合精度和更好的预测效果。传统的元件模型不能很好地反映黄土的蠕变特性,在以往研究中,将多个元件串并联起来建立的模型的拟合精度有所提高。但是,复杂模型中存在很多参数,这些参数没有明确的物理意义,不利于分析和归纳土体的流变特性。本研究选用分数阶Merchant模型对黄土蠕变特性进行分析研究。

图9 分数阶模型与传统模型对同一试验数据(SS01-1)拟合程度对比

表3 SS01-1不同模型参数

3.2 分数阶蠕变模型的参数确定

含有4个参数的传统分数阶Merchant模型准确描述了黄土的瞬时弹性变形、衰减蠕变和稳态蠕变3个阶段。对处理后的蠕变试验数据进行拟合,得到分数阶Merchant模型参数,见表4。蠕变试验数据和分数阶Merchant模型模拟结果如图10—12所示。从表4中可知:拟合优度R2均大于0.980 0,说明分数阶Merchant模型能很好地模拟黄土的蠕变特性;分数阶次α虽然接近于0,但较整数阶蠕变模型,其很大程度上提高了试验数据的拟合精度。文献[17]使用分数阶Merchant模型对软黏土三轴蠕变试验进行拟合,得到分数阶次α为0.620 0～0.630 0,其他模型参数也有很大差异,说明不同类型土体的模型参数存在一定差异。

图10 分数阶Merchant模型对围压75 kPa下分级加载蠕变拟合结果

图11 分数阶Merchant模型对围压150 kPa下分级加载蠕变拟合结果

图12 分数阶Merchant模型对围压225 kPa下分级加载蠕变拟合结果

表4 不同围压不同偏应力下排水蠕变分数阶Merchant模型参数

4 黄土分数阶蠕变模型的应用

合理的数值模拟可弥补试验的不足,为理论分析提供重要依据。文献[17]基于分数阶Merchant模型实现了本构模型在ABAQUS中的二次开发,并验证该本构模型的适用性。本研究为了验证黄土分数阶Merchant模型的适用性,建立直径39 mm、高度78 mm的圆柱体有限元模型,土样有限元模型网格划分如图13(a)所示。计算分析过程分为3步,即:con分析、load分析和creep分析。模拟分为两步,即:①模拟等压固结过程;②保持围压不变,模拟增加100 kPa主应力的加载过程,荷载的施加速率和作用时间均与室内试验保持一致。模型参数使用表4中拟合所得参数。围压为75 kPa的分级加载蠕变试验计算结果的应变云图如图13(b)所示。

图13 围压75 kPa(SS01)试验网格划分和计算结果云图

图14为75 kPa围压下蠕变试验和数值模拟结果的对比情况。从图14中可以看出,ABAQUS计算结果与蠕变数据基本一致。因此,分数阶Merchant模型的UMAT子程序可以有效地描述黄土的瞬时弹性变形、衰减蠕变和稳态蠕变阶段。

图14 围压为75 kPa(SS01)蠕变试验数值模拟计算结果与试验结果对比

黄河宁蒙河段渠堤往往作为交通便道使用,使得渠堤承受更大附加荷载,从而造成渠堤地基不均匀沉降,引发渠堤土体边坡不稳,影响渠道正常运行。本研究选取一段渠道路堤,采用有限元软件ABAQUS进行建模计算,渠道断面尺寸如图15所示。模型顶面为自由应力边界,底面为固定边界,侧面水平位移受限,顶面施加100 kPa的均布荷载。网格划分如图15所示,计算时调用文献[17]中开发的UMAT子程序,选用SS02的试验参数,即E1=63.42 MPa、E2=156.67 MPa、η=160.29 MPa·h、α=0.001 20。

图15 模型尺寸和网格划分

计算分析过程分为3步,即:geo分析、load分析和creep分析。第1步的类型为Geostatic,在此分析步中对模型y方向施加-10 N/kg的坚向力,进行地应力平衡;第2步类型为load,在模型顶面施加100 kPa的均布荷载;第3步类型为creep,在此步中调用分数阶Merchant模型对土体蠕变变形进行计算。加载120 d后,位移分布的计算结果云图如图16所示。