姚兆明, 蹇膨远, 郭梦圆

(1.安徽理工大学土木建筑学院, 淮南 232001;2.矿山地下工程教育部工程研究中心,淮南 232001)

人工冻结法被广泛应用于城市地下交通、地下采矿巷道以及深基坑等地下工程中。随着城市化进程的加快,城市用地日益紧张,地下空间开始逐渐被开发利用,但地下工程开挖深度的增加,冻结法施工的冻结管、冻结壁等设施被破坏的风险也增加。为确保工程、设施和人员的安全,分析冻土的力学特性显得尤为重要。

人工冻土的蠕变行为是冻土最重要的力学性能之一,而蠕变模型是用来表征冻土蠕变特性的主要途径。因此,学者通过研究提出大量的蠕变模型来模拟冻土的蠕变行为,并将其蠕变模型大致分为三类[1]:①经验模型;②元件模型;③黏弹塑性模型。Yao等[2-3]进行了不同温度下的冻土蠕变试验,建立了以温度为自变量的冻土一维蠕变模型,并在三轴蠕变试验的基础上提出具有剪切强度衰减的冻土蠕变模型;刘萌心等[4]在该模型基础上引入应力历史并修正,得到的模型能较好地反映不同温度及压力下的蠕变过程;李昂等[5]建立了基于Burgers模型的单因素、双因素元件蠕变模型,并通过对实验数据拟合得到了考虑时间、含水率和应力状态的多因素经验蠕变模型;周小棚等[6]对西原模型进行改进克服了其无法描述岩土加速蠕变阶段的缺陷,并推导出了蠕变破坏时间;姚兆明等[7]采用温度和加载系数作为内变量,通过开展冻土单轴压缩与分级加载蠕变试验,提出了冻结黏土内变量蠕变模型,并验证了模型的合理性。

上述模型都是整数阶模型,对于蠕变中的非线性渐进过程无法准确描述,而分数阶本构模型可以描述蠕变中应变的非线性渐进过程,并且分数阶模型具有参数少、形式简单等优点。近年来,分数阶微积分理论被越来越多地用于岩土工程领域来描述不同类型岩土材料的力学特性。Yin等[8]将经典西原模型中的阻尼器换为分数阶导数Abel阻尼器,建立了分数阶本构模型并推导出了土的分数阶模型在蠕变、应力松弛、加载和卸载等不同条件下的解析公式;殷德顺等[9]把Hollomon提出的金属塑性拉伸变形方程式[10]引入到岩土中,在分数阶导数的基础上提出了岩土应变硬化指数理论,获得了反映岩土应变硬化能力的参数;肖华杰等[11]通过分段模拟来描述瞬时和蠕变应变,并构建了考虑基质吸力的弹性体和分数阶黏滞体,建立了可考虑基质吸力的非饱和粉质黏土蠕变本构模型;Wang等[12]基于流变学和分数阶微积分建立了分数阶弹性和黏性单元,建立了物理意义明确、形式简单的基于分数阶弹性和黏性单元的黏土蠕变本构模型。

然而,大多数模型参数较多且无法描述冻土蠕变完整的三个阶段,尤其是非线性的加速蠕变阶段[13]。随着损伤力学不断发展与应用,许多学者将损伤效应引入到各蠕变本构模型中。Li等[14]基于Riemann-Liouville型积分函数,建立了考虑温度-损伤-应力耦合的冻结砂岩非线性蠕变本构方程,并基于分数阶理论建立了非线性蠕变损伤方程;李德建等[15]建立了与弛豫时间相关的分数阶变阶函数,依此构造变阶分数阶损伤蠕变模型,且进一步拓展到三轴状态下,使能够合理地描述砂岩蠕变的三个阶段。

分析山西某矿井井筒检查孔黏土不同冻结温度下的单轴蠕变试验曲线,得到温度对冻结黏土蠕变特性的影响规律。在Singh-Mitchell[16]模型的基础上,引入分数阶导数理论,建立分数阶冻土蠕变模型。通过分析蠕变与时间取对数的拟合曲线发现两者具有线性关系,进而得出只与温度有关的模型参数。鉴于建立的分数阶冻土蠕变模型不能反映冻土蠕变加速阶段,基于Weibull概率分布假设,将损伤因子引入建立的分数阶蠕变模型,进而建立人工冻土分数阶损伤蠕变模型。

1 人工冻结黏土单轴抗压、蠕变试验分析

1.1 人工冻土单轴抗压强度试验

试验土样取自山西省某煤矿黏土,取样深度为281.20～296.80 m,黏土含水率为20%。从模具中取出原状土样,将其切碎、烘干、过筛,根据重塑黏土的制样方法严格按照原状土样的含水率等指标进行配比,并按照规范要求加工成直径50 mm,高度100 mm的圆柱体试样。将制备好的试样放入冰箱中,分别在-5、-10、-15 ℃三个温度下恒温养护24 h,如图1(a)、图1(b)所示。在安徽理工大学WDT-100型冻土试验机上进行人工冻土单轴抗压强度试验[17],试验分别在-5、-10、-15 ℃三个温度下(每组温度下三个试样)进行,应变速率设定为1%/min,实验满足以下三个条件之一时,仪器便会自动停止:①应力值下降20%;②应变超过15%;③力峰值后应变增加3%。微机每15 s自动采集数据并显示相应的应力-应变曲线,试验结束后取出已经破坏的土样,如图1(c)所示。

图1 冻结黏土试样Fig.1 Frozen clay test sample

1.2 试验结果与分析

根据单轴抗压强度试验,得到冻结黏土在三个温度水平下的单轴抗压强度如表1所示。

冻结黏土从加载到破坏大致经历三个阶段:线性增长,塑性屈服和破坏阶段。随冻结温度T的降低,冻结黏土的单轴抗压强度逐渐增大。初始阶段,应力-应变曲线基本成线性关系,到达屈服点后,强度继续增大,表现为硬化特点。冻结黏土达到最大应力后,曲线开始下降,此时土体已发生塑性破坏,如图2所示。

对冻结抗压强度与温度进行拟合,可得到人工冻结黏土的抗压强度与温度在一定条件下呈线性关系,如图3所示。温度与冻结抗压强度两者之间满足关系式

σs=2.3-0.052T,R2=0.97

(1)

式(1)中:σs为单轴抗压强度,MPa;T为冻结温度,℃。

由式(1)可知,冻结温度越低,单轴抗压强度则越大,温度每降低1 ℃,单轴抗压强度则会增加约0.052 MPa。

表1 单轴抗压强度试验值Table 1 Uniaxial compressive strength test value

1.3 人工冻土单轴蠕变试验

单轴蠕变试验在安徽理工大学WDT-100型冻土试验机上进行。将加工好的试样分别在-5、-10、-15 ℃三个温度下进行蠕变试验。以冻结黏土的单抗抗压强度试验结果为依据,采用分级加载(σ=0.3σs、0.5σs、0.7σs)的方式进行加载,其中σs为单轴强度抗压强度,即每组温度下三个试样的单轴抗压强度平均值。微机每60 s自动采集数据并显示相应的应变-时间曲线。实验满足以下三个条件之一时,仪器便会自动停止:①应力值下降20%;②应变超过15%;③时间超过10 h。

1.4 单轴蠕变试验结果与分析

根据各温度下单轴抗压强度平均值,得到冻结黏土在不同温度及加载等级下的荷载,如表2所示。

图4所示为不同温度、应力水平下的冻结黏土蠕变试验曲线。从图中可以看出,冻结黏土在施加荷载时有一定的瞬时蠕变。当应力水平为0.3σs和0.5σs时,整个蠕变过程中应变值变化范围随冻结温度的降低而逐渐减小,冻结黏土主要以衰减蠕变和稳定蠕变为主;当应力水平为0.7σs时,蠕变曲线开始出现加速阶段,蠕变损伤值较大,此时土体内部已出现裂缝并随着荷载的持续作用,裂缝不断扩大从而导致土体结构破坏。在高应力水平下,温度也影响着土体的破坏过程。如冻结温度为-5 ℃时,土体蠕变从衰减蠕变经历短暂的稳定蠕变进入加速阶段,这是由于冻结黏土在温度相对较高时抗压强度较低,因此在较高应力水平下呈现一种非稳定状态,土体自身强度抵抗不了施加的恒定荷载,故在较短时间内急速破坏。

图2 各温度下单轴抗压强度关系曲线Fig.2 Uniaxial compressive strength relationship curve at various temperatures

图3 温度与单轴抗压强度关系曲线Fig.3 Relation curve between temperature and uniaxial compressive strength

表2 不同温度和加载等级下的荷载Table 2 Loads at different temperatures and loading levels

图4 各温度和各加载等级下的冻结黏土蠕变曲线Fig.4 Creep curve of frozen clay at various temperatures and loading levels

2 S-M分数阶蠕变模型

2.1 分数阶微积分定义

参看文献[9],其中Riemann-Liouville型分数阶微积分对函数f(t)的β阶积分定义为

(2)

分数阶微分定义为

(3)

式中:D为分数阶微积分算子;t为自变量;τ为分数阶积分后的自变量。β>0,且n-1<β≤n(n为正整数);Γ(*)为Gamma函数,其定义为

(4)

式(4)中:s为复数自变量;Re(s)为复数s的实部。

对于函数f(x)=cx,其中c为常数,在0≤β≤1时,其分数阶微分为

(5)

当取t=ν0τ时,其中ν0为常数,则

(6)

(7)

2.2 分数阶S-M蠕变模型的建立

Singh-Mitchell (S-M)模型是常用于描述岩土蠕变特性的一种经验模型,表示为

(8)

对式(8)积分可得

(9)

冻结黏土的应力应变特性是介于理想固体和理想流体之间某个关系,式(9)未能反映这种特性,而分数阶导数可以很好地描述这种关系。当m=1、t1=1时,对式(9)进行分数阶微分可得

(10)

式中:B=A/(1-m)。

2.3 模型参数确定

试验的土样含水率均为20%,为使S-M模型考虑温度的影响,下面对S-M模型参数进行确定。

对式(10)“=”左右两边同时取对数得

lnε=lnB+bσ+βlnν0-

lnΓ(2-β)+(1-β)lnt

(11)

式(11)中:σ=ησs,η为0.3或0.5。

图5为0.3σs和0.5σs加载等级下-5、-10、-15 ℃三个温度下的lnε-lnt拟合关系。可以发现在同一加载等级下,不同冻结温度的蠕变与时间在取对数的情况下具有明显的线性关系;同一温度、不同加载系数下蠕变也呈线性变化规律。

图5 不同温度、各应力水平下的蠕变与 时间双对数关系曲线Fig.5 Double logarithmic relationship curve between creep and time at different temperatures and stress levels

2.3.1 参数β的确定

温度为-5 ℃时,0.3、0.5倍等级加载时的线性表达式为

(12)

温度为-10 ℃时,0.3、0.5倍等级加载时的线性表达式为

(13)

温度为-15 ℃时,0.3、0.5倍等级加载时的线性表达式为

(14)

2.3.2 参数B和b的确定

lnB1+0.777b1=0.089 2

(15)

lnB1+1.295b1=0.387 5

(16)

解得B1=0.70,b1=0.58。

同理可得T=-10 ℃和T=-15 ℃时对应的参数B和b的值,整理后各温度下的参数情况如表3所示。由表3可知,随着温度降低,参数B的值不断减小,而b值不断增大,通过拟合发现温度随B值与b值均呈线性变化,如图6所示。

温度T与B和b拟合关系式为

B=0.945 8+0.051 6T

(17)

b=0.413 1-0.030 8T

(18)

由此可以得出与温度有关的分数阶S-M蠕变模型为

表3 各温度下参数情况Table 3 Parameters at various temperatures

图6 冻结黏土参数B和b与温度关系曲线Fig.6 Relations of parameter B and b of frozen clay against temperature

(19)

根据推导出的冻结黏土分数阶S-M蠕变模型计算出该黏土冻结状态下的蠕变值,将计算得出的结果与试验值进行对比,结果如图7所示。由图7可知,相同温度下,0.3σs和0.5σs荷载等级下均与曲线吻合情况较好,其中-10 ℃下0.3σs荷载等级下的试验值与计算值出现偏差,但曲线总体走势基本一致,表明该模型可以很好地计算温度效应下黏土的蠕变,验证了该模型对冻结黏土的适用性。

3 考虑损伤的S-M分数阶蠕变模型

3.1 损伤蠕变效应

在应力水平达到0.7σs时,蠕变开始进入加速阶段,此时的蠕变损伤值较大(图4)。岩土的蠕变损伤是在外加荷载的作用下,岩土内部微裂纹不断扩大造成土体破坏的过程。建立的S-M分数阶蠕变模型对加速阶段将不再适用。

假设岩土微元强度服从Weibull分布,则其概率密度函数可表示为

(20)

式(20)中:ε为土的轴向应变值;N和n、α为Weibull分布参数。

损伤变量D可以定义为岩土在某一级荷载作用下,岩土内部被破坏的微元体数目Vt与岩土内部总微元体数目V之比,即

D=Vt/V

(21)

式(21)中:D为损伤变量;Vt为岩土内部被破坏的微元体数目;V为岩土内部总微元体数目。

其中岩土内部被破坏的微元体数目Vt为

(22)

得

(23)

式中:P(x)可作为微元体概率密度函数。

损伤变量D也称损伤因子,其取值范围为(0,1)。D=0时表示岩土内部结构完全没有裂纹的一种理想状态,D=1时表示岩土完全破坏从而失去承载能力的状态。

3.2 S-M分数阶损伤蠕变模型的建立

将损伤变量D引入上文建立的S-M分数阶蠕变模型中,得到S-M分数阶损伤蠕变模型,表达式为

(24)

针对人工冻结黏土加速变形特点,将损伤变量D进行修正,表达式为

(25)

当n=1时,模型方程为

(26)

式(26)中:ε为蠕变应变;σ为根据各温度下单轴抗压强度平均值确定的施加在试样上不同加载等级的恒定荷载,MPa;t为蠕变时间,h;B为试验有关常数;b为确定应力有关常数;ν0为常数,取ν0=1;α为确定损伤变量参数。

对图4中-5、-10、-15 ℃三个温度下的0.7σs荷载等级蠕变曲线数据计算,得到分数阶损伤蠕变模型中参数α,结果如表4所示。

由表4可知,参数α的数值随着温度的降低不断增大,通过拟合发现温度随参数α的数值呈线性变化且拟合优度为0.99,如图8所示。

参数α与温度T间的关系式为

α=-0.47T+2.67

(27)

人工冻结黏土S-M分数阶蠕变损伤模型为

表4 各温度下参数αTable 4 parameter α at various temperatures

(28)

根据拟合参数计算出冻结黏土的蠕变数据,将计算得出的结果与试验值进行对比,结果如图9所示。

由图9可知,参数α越大,冻结黏土抵抗变形破坏的能力越强,因此,参数α反映了冻结黏土的强度特性。在-5 ℃时计算值与试验值基本趋于一致,在-10 ℃和-15 ℃中计算值与试验值出现偏差并在某一时间两线相交,但曲线的总体走势基本一致,吻合程度较高,故能够较好地反映冻结黏土整个蠕变损伤过程,验证了该模型对冻结黏土的适用性。

图8 冻结黏土参数α与温度关系曲线Fig.8 Relation curve between frozen clay parameter αand temperature

图9 试验值与计算值的对比Fig.9 Comparison between experimental values and calculated values

4 结论

针对传统Singh-Mitchell模型在确定参数、描述人工冻土蠕变特性等方面的不足,基于分数阶导数理论,建立反映人工冻土蠕变稳定阶段的S-M分数阶蠕变模型。将损伤因子引入建立的S-M分数阶蠕变模型,进而建立能反映蠕变加速阶段的人工冻土S-M分数阶损伤蠕变模型,得到如下结论。

(1)利用分数阶导数理论对S-M模型中的时间参数求分数阶导数,克服了传统S-M模型未能反映冻结黏土蠕变介于理想固体和理想流体之间特性的不足,建立的人工冻结黏土分数阶蠕变模型概念清晰、推导严谨。

(2)有别于一般的经验模型参数完全由最小二乘法拟合得到,建立的人工冻结黏土分数阶蠕变模型参数由不同冻结温度下,对不同加载应力产生的蠕变与时间取对数联立方程组求得,因此所建模型参数具有一定的物理意义且易于确定,便于工程应用。

(3)为在所建模型中能反映由于损伤导致蠕变加速发展的现象,基于Weibull分布假设,将损伤因子引入S-M分数阶蠕变模型中,得到用以描述冻结黏土的S-M分数阶损伤蠕变模型。模型拟合的结果与试验数值较为吻合,验证了该模型的适用性,并为描述冻土蠕变损伤提供了一种方法。

文中所建立的模型方法在人工冻结黏土蠕变稳定阶段、加速阶段得到了很好的应用,但该建模方法在其他类型的人工冻结土质如砂土蠕变的应用有待进一步验证。同时为计算不同蠕变阶段分别建立了分数阶蠕变模型,而在实际工程中,很难准确确定发生加速阶段的荷载,在不同的应力边界下进行数值分析时需要对模型进行切换。因此,建立能同时计算稳定阶段、加速阶段蠕变的统一模型是下一步的研究目标。