非饱和土蠕变力学特性试验及分段模拟

肖华杰，胡 鹏，曾彩云，贾 逸，魏建柄

(1.南宁市城市建设投资发展有限责任公司,南宁 530031；2.核工业西南勘察设计研究院有限公司,成都 610061；3.西南交通大学地球科学与环境工程学院,成都 610031；4.中国地质调查局地质灾害防治技术中心,成都 611734；5.陕西铁道工程勘察有限公司,西安 710043)

大量工程实践证明,基坑支护位移、路基沉降和边坡变形等工程问题不仅与岩土体材料性质和外界荷载有关,还很大程度上受荷载作用时间影响[1-3]。外界荷载的作用时间是决定土体变形累积量的重要因素之一,而非饱和土是一种包含固、液、气三相体系的土壤,相比固、液两相体系饱和土,非饱和土的蠕变力学特性更为复杂。在库岸边坡中,由于降雨入渗及库水位的变化,土体在饱和及非饱和状态之间转化,土体具备非饱和特性,其蠕变变形逐渐累积,对库岸边坡的长期稳定性造成潜在威胁[4-5]。现对于非饱和土蠕变特性的研究已有一定的进展,郑俊等[6]以三峡库区某滑坡非饱和滑带土为研究对象,基于伯格模型研究其非饱和蠕变特性；李梦姿[7]采用重力式恒定分级加载法,研究非饱和土在不同基质吸力下的应变特征,分析在不同吸力下蠕变特性的变化规律；祝艳波等[8]开展非饱和碎屑土的三轴压缩蠕变试验,分析基质吸力对其强度与蠕变特性的影响,探索吸力对非饱和碎屑土强度损失率的影响；刘虎虎等[9]开展不同法向应力、含水率下的滑带土试样的剪切蠕变试验,研究不同含水率下蠕变模量、残余抗剪强度指标的变化趋势。

以某库岸边坡中非饱和粉质黏土为研究对象,进行考虑基质吸力控制条件下的三轴压缩固结排水蠕变试验,分析非饱和土的蠕变曲线特征,研究等时偏应力-应变关系,基于此提出蠕变特性分段模拟的思路,分别构建一个考虑基质吸力的弹性体和分数阶黏滞体来描述非饱和土的瞬时、蠕变应变,将其串联得到新的非饱和土非线性蠕变本构模型,分段求解模型参数,进行试验数据与预测曲线的对比,证明该模型的可行性和合理性。再引入相关文献中非饱和粉质黏土相关数据,证明该模型反映非饱和粉质黏土蠕变力学行为的适用性。

1 蠕变试验及成果

1.1 试验材料

非饱和粉质黏土试样取自某水电工程近坝库岸边坡,土样的基本物理力学参数如表1所示。将土样运回实验室后,将其自然风干后碾散,再过2 mm筛。为了使试样更易成形,用蒸馏水将其配制成含水率为20%的土样,再用保鲜膜包裹,静置1 d以使水分扩散均匀,最后通过削土器制备重塑圆柱样,规格为ø60 mm×120 mm。

表1 基本物理力学参数Table 1 Basic physical and mechanical parameters

1.2 试验步骤

通过FSR-6型非饱和土三轴蠕变仪,进行基质吸力控制条件下的非饱和三轴压缩蠕变试验。为了确定蠕变试验加载方案,首先进行三轴排水剪切试验确定排水剪切强度τf。设置围压σ3=100 kPa,基质吸力s分别为100、200、300、400 kPa,试验采用逐级增量加载方式,每级应力水平历时200 h以上。应力水平D从0.55开始,对应偏应力差为0.55τf,每一级递增0.05直至破坏,破坏偏应力为(σ1-σ3)f,其中,σ1和σ3分别为轴向应力和围压,σ1-σ3为偏应力。三轴排水剪切试验结果及蠕变试验加载方案如表2所示。

表2 蠕变试验加载方案Table 2 Creep test loading plan

1.3 成果分析

如图1所示为s=400 kPa时的蠕变试验结果。图1中基于玻尔兹曼叠加原理[10]处理得到分别加载蠕变曲线如图2所示,土样应变数据统计如表3所示。择取图2中1、26、51、76、101、126、151、176、201 h共9个时间节点的偏应力和应变数据,将其绘制成等时偏应力-应变曲线,如图3所示。如图4所示为不同基质吸力下的等时偏应力-应变曲线,由于曲线密集,为便于观察,仅保留时间节点为1 h和201 h的曲线。

图1 分级加载蠕变曲线Fig.1 Graded loading creep curves

图2 分别加载蠕变曲线Fig.2 Separate loading creep curves

表3 应变数据统计Table 3 Strain datas statistics

图3 等时偏应力-应变曲线Fig.3 Isochronous eccentric stress-strain curves

图4 不同基质吸力的等时偏应力-应变曲线Fig.4 Isochronous eccentric stress-strain curves of different matrix suction

由图1可看出,土样在加载瞬间首先表现出一定量的弹性瞬时应变,接着进入衰减蠕变阶段,该阶段蠕变变形不断累积,然后进入稳定蠕变阶段。在下一级应力水平下,继续重复该过程。结合表2可知,土样在s为100、200 kPa时,历经5级加载,s为300、400 kPa时,加载等级为4级。

由图2可看出,非饱和土在同一基质吸力下的每一级蠕变曲线在形态上较为相似,而实际当s为100、200、300、400 kPa的蠕变曲线形态和发展趋势都大致相同。

由表3可知,随着应力水平的提升,土样的瞬时、蠕变应变都呈递增趋势,每一级加载下,瞬时应变都大于蠕变应变,且瞬时应变的增长幅度大于蠕变应变,这可能是由于试验设置的偏应力荷载增长幅度较大,土样在外界荷载作用增加的情况下,瞬时变形的力学响应更为灵敏。

由图3可发现,时间节点为1 h的曲线为线性相关,与26 h及其以后的节点在横轴上差距较大,26～201 h共8个节点的曲线性态和变化规律基本一致,形成曲线簇。该曲线簇有逐渐偏于应变轴的趋势,表现出非线性特征。结合图4可看出,在不同基质吸力条件下,时间节点为1 h的曲线均为线性相关,实际上26～201 h共8个节点的曲线簇均表现出较为明显的非线性特征。

2 基于分段模拟的蠕变模型

2.1 模型构建

采取蠕变曲线分段模拟的思路,分别对非饱和土的瞬时应变和蠕变应变进行模拟。瞬时应变由弹性变形产生,服从Hooke定律,由此可用弹簧体对瞬时性态进行描述。而图2中非饱和蠕变应变包含衰减蠕变阶段和稳定蠕变阶段,一般采用牛顿体描述其蠕变性态。由于非饱和粉质黏土蠕变的非线性特征,传统元件模型难以辨识其蠕变曲线,故在本文试验结果的基础上,考虑瞬时性态和蠕变性态的非定常性,对流变参数进行非定常化,从而建立非线性蠕变模型。为了模型简洁便于计算,仅考虑Maxwell模型(串联)或Kelvin模型(并联)作为基础模型,由于Kelvin模型中H体(弹簧体)和N体(牛顿体)并联,其中H体和N体的应变相等,与本文分段模拟瞬时应变和蠕变应变的思路不符,而Maxwell模型中总应变等于H体和N体的应变之和,总应力与H体应力及N体应力相等,满足本文分段模拟思路,故选用Maxwell模型作为基础模型。由此分别构建分数阶黏滞体和考虑基质吸力的弹性体,参考Maxwell模型的模型结构将其串联后得到本文分段模拟的蠕变本构模型,模型示意图如图5所示。

图5 模型示意图Fig.5 Model schematic diagram

图5中的模型结构为串联形式,其状态方程为

(1)

式(1)中:σve和εve分别为分数阶黏滞体的应力和蠕变应变；σe和εe分别为弹性体的应力和瞬时应变；σ为初始应力。

2.2 基于分数阶微积分的蠕变性态模拟

2.2.1 Riemann-Liouville型分数阶微积分

从时间的角度而言,传统整数阶微积分所表征的是一个物理或力学过程某时刻的变化或某种性质,而分数阶微积分的核心在于其阶数为有理分数、无理数甚至复数,可表征与某一现象的整个发展历史有关的性质[11]。分数阶微积分具有非局部性、非线性,能描述黏弹性及多孔材料的“记忆性”和非线性动力系统的“遗传性”,且物理意义清晰,形式简单[12]。

鉴于分数阶微积分在非线性动力系统中的优越性,故将其引入到非饱和土蠕变性态的辨识。由于时域的不同,产生不同的分数阶微积分定义,常用的有: Grunwald-Letnikov型、Riemann-Liouville型和Caputo 型[12]。其中Riemann-Liouville型因其简练和效果优良而广泛应用于数理分析,依该理论,函数f(t)在可积区间[0,t]的α阶Riemann-Liouville积分定义为

g)α-1f(g)dg

(2)

函数f(t)的α阶微分相应地定义为

(3)

式中:t为时间；f(t)为在可积区间[0,t]的某一函数；α为大于0的分数阶数；g为用于拉普拉斯变换的某一自变量；Γ(α)为伽马函数；n=[α]，n为大于α的最小整数。

Γ(α)定义为

(4)

(5)

2.2.2 基于分数阶微积分的软体元件

(6)

式(6)中:σ为应力；η为黏滞系数。

非饱和土是一种包含固液气三相体系的复杂土体,分数阶软体元件可理解为介于理想固体和理想流体之间的材料,Scott-Blair[13]提出描述该状态的软体元件,其本构方程为

2)从集合预报系统效能评价来看,据降水量调整方案扰动能量最大,表明其对预报场影响最大,最有可能改进预报场效果。各成员之间差异较大,包含真实大气状态的可能性也较大,且集合平均预报误差较小,效果较理想。该方案降水预报对强降水区域的模拟较好,ETS评分在24 h累积降水量≥50 mm等级预报上优于其他方案,相对控制预报来讲,提高精度约为15%,预报指示意义更明确。

(7)

式(7)中:ηα为该软体元件中的黏滞系数。

当应力σ恒定时,采用Riemann-Liouville分数阶微积分算子理论,对式(7)进行分数阶积分得

(8)

式(8)即为基于分数阶微积分的软体元件。

取σ=500 kPa,ηα=15×105kPa·h时,根据式(8)可绘制α的不同取值所对应的蠕变曲线如图6所示。

图6 分数阶软体元件蠕变曲线Fig.6 Creep curves of fractional soft element

由图6可看出,α∈[0,1],当α值从0增大至1的过程中,蠕变曲线逐渐远离横轴。当α在区间[0,1)内逐渐增大时,曲线斜率也逐渐递增至某一固定值,蠕变曲线表征出较为明显的非线性特征；当α增至1时,蠕变曲线表现出完全线性关系,此时该软件元件为牛顿黏壶。

2.2.3 基于分数阶微积分的非定常黏滞体

对于非线性蠕变模型的研究,常用的方法是考虑流变参数的非定常性,将流变参数进行非定常化处理。由于蠕变应变的累积与时间和应力关联密切,故本文将黏滞系数改为与时间和应力的函数,文献[14]提出一种描述黏弹性应变的SN元件,其本构方程为

(9)

式(9)中:η0、λ、c为与蠕变特性相关的黏滞参数。

通过SN元件对分数阶软体元件进行改进,将式(9)代入式(8)可得：

(10)

式(10)即为改进后基于分数阶微积分的变参数黏滞体的蠕变方程。

2.3 考虑基质吸力的弹性体

非饱和土在外界荷载下瞬间响应,产生弹性瞬时应变,弹性变形过程满足Hooke定律,故用弹簧体进行描述,其本构方程为

(11)

式(11)中:E0为初始弹性模量。

由式(11)可看出,弹性应变εe仅与弹簧体的应力σe和初始弹性模量E0相关,由于土样蠕变试验是逐渐增量加载的方式,每一级加载下应力保持恒定。故而在同一应力水平下,E0作为唯一因素影响非饱和土对瞬时加载的弹性响应,探索E0与基质吸力s的关系便可建立瞬时应变εe与s之间的联系。这样既使蠕变参数E0非定常化,又使非饱和土蠕变本构模型可以定量化反映含水率作用,从而将基质吸力作为独立变量体现到模型中,由此在弹簧体的基础上建立应力-基质吸力-瞬时应变关系模型。Janbu[15]研究发现初始切线模型E0和围压σ3在双对数坐标中线性相关,即初始切线模量E0是σ3的幂函数。故此,E0与s之间可采用同样函数

(12)

式(12)中:pa为大气压,pa=101.33 kPa；F和k为材料常数。

将式(12)代入式(11)可得：

(13)

2.4 非线性蠕变模型的建立

将式(10)和式(13)代入式(1)可得：

(14)

式(14)即为本文考虑基质吸力的非线性蠕变本构模型。

3 模型参数求解及验证

3.1 参数F、n的确定

本文建立的考虑基质吸力的非线性蠕变本构模型包含F、n、η0、λ、c、α共6个参数,将这6个参数分为两类进行求解,其中,F和n为瞬时应变参数,η0、λ、c、α为蠕变应变参数。式(13)中,σe/εe与参数F、n之间表现出幂次关系,σe对应蠕变试验中的偏应力σ1-σ3,对表3中s=400 kPa下4级应力水平下的(σ1-σ3)/εe求取平均值,这样可得到s=400 kPa的σe/εe值。同样的方法求取s=100～300 kPa的σe/εe值,与s/pa进行拟合,如图7所示。

由图7可知,曲线幂次拟合效果较好,R2达到0.975 1,取参数F=521.76,n=0.107 3。E0随着基质吸力s的增大而递增。随着s的减小,库岸边坡内非饱和土初始切线模量不断减小,土体变软,说明随着库岸边坡内含水率的增加,土体蠕变变形更加显著。

3.2 参数η0、λ、c、α的确定及模型验证

参数η0、λ、c、α采用一般的非线性最小二乘拟合便可得到,利用数学优化软件1stOpt,基于Levenberg-Marquardt 算法,对本文试验数据进行辨识,限于篇幅,仅给出在s为300 kPa和400 kPa下的参数η0、λ、c、α取值,如表4所示。引用文献[7]中的元件模型与本文模型进行对比验证,以s为300 kPa和400 kPa的试验数据为例,验证结果如图7所示。

图7 s与σe /εe关系曲线Fig.7 Relationship between s and σe /εe

表4 参数η0、λ、c、αTable 4 Parameters of η0、λ、c、α

由图8可看出,基于分段模拟所建蠕变本构模型辨识能力较强,平均R2=0.987 2,文献[7]中元件模型辨识能力稍弱,平均R2=0.930 5。文献[7]模型预测值在衰减蠕变阶段和稳定蠕变阶段结束部分与试验数据差异较大。

图8 试验值与理论值对比曲线Fig.8 Comparison curve between experimental and theoretical values

3.3 模型适用性验证

蠕变试验对象为非饱和粉质黏土,结合试验成果,基于分段模拟的思路建立了非线性蠕变本构模型,为了验证该模型预测非饱和粉质黏土蠕变力学行为的适用性,引用文献[16-17]中非饱和粉质黏土蠕变试验数据,通过本文和文献[7]模型进行拟合验证,拟合结果如图9所示。

图9 试验值与理论值对比曲线Fig.9 Comparison curve between experimental and theoretical values

由图9可看出,文献[7]模型对于非饱和粉质黏土的衰减蠕变阶段和稳定蠕变阶段结束部分辨识能力较差,辨识文献[16-17]数据的平均R2分别为0.941 8和0.937 2,本文所建非线性蠕变模型的平均R2分别为0.990 6和0.988 4。

综合图8、图9分析可得,本文所建非线性蠕变模型对非饱和粉质黏土的辨识能力较强,能够较为准确地描述非饱和粉质黏土的蠕变力学行为。

4 结论

(1)进行了非饱和粉质黏土固结排水三轴压缩蠕变试验,在不同应力条件下该土样蠕变曲线形态较相似,通过对其瞬时、蠕变应变特征的分析,基于分段模拟思路构建考虑基质吸力的弹性体和基于分数阶微积分的黏滞体,由此建立非饱和粉质黏土蠕变本构模型。

(2)在相同加载荷载下,基质吸力的减小会导致非饱和土蠕变变形更为显著。库岸边坡工程中含水率是影响非饱和土蠕变变形累积的重要因素,含水率增大会加剧时效变形。

(3)通过所建模型对本文和相关文献中非饱和粉质黏土蠕变数据进行辨识,对比分析试验曲线和预测曲线,证明本文模型反映非饱和粉质黏土蠕变力学行为的合理性和适用性。

本文试验背景为某水电工程近坝库岸边坡,取边坡粉质黏土研究其非饱和蠕变特性,后续研究中还应对非饱和土的蠕变变形机制及力学模型在工程实践中的应用进行进一步的探讨。