基于冻土蠕变的分数阶广义Kelvin模型

长江大学石油工程学院 湖北 武汉 430100

0 引言

我国的多年冻土面积占国土陆地面积的21.5%,季节冻土分布面积占53.5%,是世界第三冻土大国。[1]我国许多重大管道工程不可避免地要穿越冻土区,例如已建成的国家重点工程—西气东输管线。其中轮南-武威段季节性冻土冻深为0.8～1.4m,武威-靖边段季节性冻土冻深为1.0～1.3m。冻土蠕变是一种较为复杂的现象,土壤的冻结与融化作用会引起管道变形。近年来,随着冻土区一系列重要基础设施和重大工程的逐步开展以及人工冻结法施工技术在隧道开挖中的广泛应用,冻土在复杂环境下(尤其是高应力水平下)的蠕变变形计算成为一个非常重要的研究课题。所以准确建立冻土蠕变模型,对保证冻土区管道安全具有重要意义。

1 研究背景

冻土蠕变研究最早是在国外开展的,苏联始于1930年代初期,北美也在1950年代初开始对冻土强度进行研究,逐步建立了冻土强度与蠕变理论。[2]自20世纪60年代以来,冻土的流变本构理论研究取得了较大的发展。其中经验模型和整数阶元件模型是目前最常用的两类冻土蠕变本构模型。经验模型通过对人工冻土蠕变试验数据进行最小二乘法拟合,得到应变与时间的函数关系。该模型公式简单,针对性强,可以很好地对蠕变曲线进行拟合,而模型中的参数缺乏明确的物理意义[3]。整数阶元件模型是由弹簧元件、粘壶等元件经过各种组合而得,整数阶元件蠕变模型的构成具有明确的物理力学意义,直观易懂[4]。整数阶伯格斯元件模型、宾厄姆元件模型及西原元件模型等是常用的人工冻土经典蠕变模型[5]。然而,模拟一个完整的人工冻土蠕变模型,往往需要很多元件进行组合,进而导致所建模型参数繁多。而分数阶导数理论是数学分析的一个分支,是专门研究任意阶积分和微分的数学性质及其应用的领域。分数阶导数模型可以精确的拟合冻土的蠕变过程,而且形式简单、统一,在计算过程中需要调整的参数很少,仅采用很少的几个元件的组合就可以获得很好的结果。相比较于整数阶导数以及积分,分数阶微积分可以更简洁准确的描述具有历史记忆和空间全域相关性等复杂力学和物理过程。[6]目前对于冻土蠕变模型的构建,建立考虑冻土应力-应变-强度特性随时间变化的蠕变本构模型,有助于解决工程应用中遇到的蠕变问题。研究人员从不同方向出发,提出了一系列用于预估冻土蠕变变形的模型。例如:Hou等[7]同时考虑了应力引起的弱化及强化效应对冻土蠕变过程的影响,提出了一种分数阶蠕变本构模型。李鑫等[8]通过引入硬化因子与损伤因子来考虑温度、应力造成的冻土材料强化与弱化,进而提出了适用于冻土的改进Nishihara模型。Yang等[9]针对冻结粉土进行了三轴压缩及蠕变试验,发现围压的增长能抑制冻土的变形,并提出了扰动状态的黏弹塑性本构模型及一种新的元件蠕变本构模型。姚亚峰等[10]通过改进遗传算法,获得各温度下蠕变模型,其能准确反映冻土蠕变的整体规律。周志伟等[11]进行不同温度和围压条件下的三轴压缩和蠕变试验,基于超塑性理论和过应力原理建立了冻结黄土的速率相关本构模型。李东伟等[12]在Nishihara模型的基础上,通过分析深层人工冻土在不同围压和温度下的三轴蠕变和剪切试验结果,建立了深层人工冻土的分数阶本构模型。

经典Kelvin模型是由弹簧与牛顿粘壶串联组合而成,能描述冻土的衰减蠕变过程,而不能完整描述非衰减蠕变,局限性较多。因此我们提出将Kelvin模型中的整数阶牛顿粘壶替换为分数阶Abel粘壶,并串联一个应力-应变双控的非线性黏壶元件,从而得到改进后的分数阶导数Kelvin模型。改进后的模型实用性及合理性通过已有的蠕变试验结果进行验证。

2 模型改进

传统Kelvin模型是普通的元件模型,通过虎克弹性体和牛顿黏性体的串并联组合来模拟冻土蠕变的力学行为。由于这些元件均是线性元件,只能描述线性变化,所以我们引入非线性元件,即分数阶元件,以便更好拟合冻土蠕变过程。如图1、图2,将Kelvin模型中原本的牛顿粘壶替换成Abel粘壶,使得改进后的模型参数少且物理意义明确,与实验结果也能很好的进行拟合,为人工冻土领域的计算提供了一种更精确的新方法。进一步通过分析发现,改进后的Abel粘壶虽能够描述非线性变化,但该模型却始终无法较好的拟合冻土加速蠕变阶段的蠕变情况。实际上,冻土蠕变大部分时间处于非稳定蠕变阶段和稳定蠕变阶段,只有当应力应变达到某一临界值时,才会出现加速蠕变阶段,且经历时间很短,即可导致冻土结构物破坏。分数阶Kelvin模型无法更好的模拟加速蠕变的原因也在此,而采取继续串联一个应力-应变触发的Abel粘壶即可将两种情况分开,如图3。

2.1 Abel粘壶 Abel粘壶的应变满足以下关系:

式中:η2是分数阶元件的粘滞系数,r是分数阶元件求导阶数,t为时间。

图4 Abel粘壶Figure4 Abel sticky pot

2.2 应力-应变触发的Abel粘壶 复杂应力状态下,冻土蠕变与时间关系曲线分为两类:衰减型蠕变和非衰减型蠕变。图-5中横轴表示时间,纵轴表示应变偏量对时间的变化率。当应力偏量较小时,蠕变呈衰减型,如图5a所示,在衰减过程中,变形速率逐渐趋近于零,蠕变变形收敛于某一个变形水平。而当应变超过某一界限值,其蠕变呈现非衰减型,如图5b所示,非衰减型蠕变一般可分为3个阶段:非稳定蠕变阶段AB,稳定蠕变阶段BC和加速蠕变或破坏阶段CD。实际工程中,加速蠕变阶段CD往往经历很短时间,即可导致冻土结构物破坏。因此,对实际工程有研究意义的通常是衰减型蠕变以及非衰减型的蠕变阶段AB段和BC段

图5 变形随时间变化的蠕变曲线Fig.5 creep curves of deformation with time

在高应力情况下,冻土内部受到破坏而发生加速蠕变。吴斐等研究盐岩时引入了一种带有应变触发的粘壶,通过应变的条件触发该元件的开关,既可以判断蠕变过程是否进入加速阶段,又能清楚反应衰减型蠕变与非衰减性蠕变的区别。当应力满足σ≤σs(σs为屈服应力)且应变满足ε≤εa(εa为冻土进入加速蠕变的临界值),该元件开始触发,冻土显示加速蠕变特性,如图6.

图6 应力-应变触发的Abel粘壶Figure6 Abel sticky pot triggered by stress-strain

应力-应变触发的Abel粘壶的本构关系为:

3 本构方程推导

如图3,改进后的模型是采用线性元件与非线性元件串并联组合而成的元件模型,既可以针对低应力水平下的线性变化作出相应的表示,又可以描述超过屈服应力和屈服应变情况下的加速蠕变阶段。该模型的应力应变本构方程满足元件串并联关系,根据串并联公式,可得:

冻土的蠕变特性与应力大小相关,当冻土处于低应力水平下,蠕变只表现衰减型蠕变的特征。元件没有触发最后一个装置,简化为图2所示的模型,本构方程为(1)式。当冻土处于高应力水平下,冻土开始进入加速蠕变阶段,应力-应变双触发粘壶开关打开,进入工作状态,此时可用来拟合衰减型模型的冻土蠕变特性。

4 数据拟合

为了验证修改后模型的适用性,本文数据采用单轴应力状态下冻结的兰州沙土蠕变蠕变实验数据进行验证与分析。[13]模型参数如下表:

表1 模型参数Table 1 model parameter

式(1)、(2)为改进后的分数阶Kelvin模型,根据实际数据与计算值可以看出,模型对试验数据有较好的描述。在低应力下的试样开始表现出衰减蠕变,随后进入稳定蠕变阶段。

图8 蠕变阶段应变随时间变化图Fig.8 variation of strain in creep stage

5 结论

引入应力-应变双触发粘壶对分数阶Kelvin模型进行改进,得到应变与时间的函数,可以描述冻土蠕变过程。相对之前的模型,改进后的模型能很好地描述衰减型蠕变与非衰减性蠕变两种情况,且互不影响,既保证了之前模型的优点,又能很好拟合加速蠕变阶段的蠕变特性,为冻土的蠕变特征计算提供了一种新的方法。

引入应力-应变双触发粘壶元件,仅是拟合冻土蠕变众多方法中的一种,模型中对于某些数据的拟合还不够精确,还具有改进空间。或许还有更好的方法,需要我们继续去探索研究,从而更好解决冻土区管道工程问题。