基于分数阶导数的岩石蠕变本构模型研究

胡其志,王芝超，丁志刚,2

(1.湖北工业大学 土木建筑与环境学院,湖北 武汉 430068;2.中交路桥南方工程有限公司，北京 101149)

0 引 言

岩石自身的蠕变特性对工程稳定性有重要影响。近年来，国内外学者对岩石蠕变特性作了细致而深入研究，提出很多岩石蠕变本构模型，主要分为组合元件模型、积分型本构模型和经验模型[1-4]。组合元件模型结构简单是，物理意义清晰明确，在实际工程中应用广泛。元件模型通常是线性模型，而实际工程中岩石蠕变往往呈现明显的非线性特点。针对元件模型这一劣势，研究发现运用分数阶微积分理论推导的蠕变模型在对蠕变加速阶段非线性曲线刻画上具有线性元件无法比拟的优势，因此得到了广泛关注和研究。

徐国文等[5]用分数阶黏壶和NVPB模型分别替换了西原模型中的牛顿黏壶和黏塑性体，提出改进西原模型并推导其三维蠕变本构方程；吴斐等[6]提出一种反映加速阶段分数阶非线性黏壶元件，并和传统的Abel黏壶和弹簧元件串联形成新的非线性分数阶蠕变模型；宋勇军等[7]将软体元件与胡克元件串联，引入一个具有幂函数特征的黏塑性体，推导一个新的四元件黏弹塑性蠕变模型，并给出该模型的本构和蠕变方程；何利军等[8]按照遗传流变理论，利用修正的Burgers模型导数形式作为蠕变核函数，得到一个反映软黏土蠕变的经验模型；何志磊等[9]基于分数阶微积分理论，将西原模型中的牛顿黏壶替换为Abel黏壶，并考虑蠕变参数的非定常性，得到基于分数阶导数的非定常蠕变本构模型；苏腾等[10]结合Scott-Blair 分数阶元件和变系数分数阶元件，提出一种三维变阶分数阶非线性黏弹塑性蠕变模型；许多等[11]对锦屏深部大理岩进行三轴蠕变试验，得到其轴向和环向变形规律，运用分数阶理论改进了大理岩黏弹塑性损伤蠕变模型；丁靖洋等[12]引入分数阶和损伤理论对经典三元件模型进行改进，得到基于损伤演变的盐岩分数阶三元件模型；王立业等[13]在线性元件模型上串联Abel黏壶，推导出渗透吸力条件下的饱和盐渍土分数阶蠕变模型，并进行验证。

借鉴经典元件模型的组合思路，将基于分数阶理论的软体元件并与一个胡克体和一个黏塑性体串联形成一个四元件分数阶模型，并推导其本构方程。运用绿片岩、砂岩和盐岩的蠕变试验数据进行拟合分析，发现该模型可很好模拟岩石蠕变的3阶段，且结构简单，力学意义明确。

1 分数阶软体元件

岩石作为一种实际的工程材料具有比较明显的蠕变特性，既像固体一样具有一定刚度，也可类似液体具有一定的流动性。岩石蠕变在低应力状态下往往呈现出黏弹性特征，因此可将弹簧体和黏性体进行串并联等方式组合，经典的Maxwell模型、Bergers模型和Kevin模型就是此类形式。然而这些模型都是线性模型，无法反映岩石加速蠕变阶段，因此经典元件模型的应用存在一定的局限性。研究发现运用分数阶理论在描述岩石蠕变理论曲线上具有明显优势，尤其是在刻画岩石的等速和加速蠕变阶段上。

分数阶微积分是微积分理论一个分支，它研究函数任意阶次的微分、积分算子特性，在数学和力学建模方面应用广泛。相对于整数阶微积分，分数阶微积分全局关联度好，试验结果和理论模拟值吻合度高[14]。同时在对复杂物理力学问题模拟中，分数阶模型比经典模型表述也更简洁，物理意义更明确。因此，分数阶理论在众多领域得到广泛应用。

对于分数阶微积分定义，理论上有各种方法。岩石蠕变本构模型领域中Riemann-Liouville 型分数阶微积分算子理论应用最普遍。由该理论可知，函数f(t)在可积空间[0,t]上的β阶积分定义为

(1)

相应函数f(t)β阶微分定义为

(2)

即

(3)

式中：β为分数阶数，0≤β≤1；Γ(β)为伽马函数，当β>0时，Γ(β)定义为

(4)

以上分数阶微积分的Laplace变换公式为

(5)

式中，f(p)为f(t)Laplace变换。

在经典元件模型中，胡克定律σ(t)=Eε(t)，可以准确反映理想弹性体的应力-应变关系。其中σ(t)为应力，ε(t)为应变；牛顿定律即σ(t)=ηdε(t)/d(t)可以准确反映牛顿流体的应力-应变关系。引入R-L理论，则存在一种软体元件，它的力学性质介于弹性体和牛顿流体之间，并且本构关系为

(6)

式中：σ为应力；ε为应变；t为加载时间；ηβ为黏滞系数。

当应力σ保持不变时，根据分数阶R-L积分定义对式(6)两边积分得

(7)

式中，Γ(β+1)为Gamma函数。同理，当σ(t)为常数时，即应力恒定情况下，分数阶元件将可描述应力松弛，经过推导可得松弛方程为

(8)

其中，ε0为对应于初始应力的初始应变。

2 非线性黏塑性体

在岩石的蠕变全过程中，低应力下岩石主要表现为弹性和黏弹性变形。随着应力提高,岩石所受荷载逐渐大于长期强度,蠕变进入加速阶段。此时岩石内部损伤加剧,导致黏塑性变形大大增加,从而产生破坏。因此长期强度即为岩石在长时间应力作用下维持稳定所能承受的最大荷载。长期强度主要根据岩石蠕变试验间接获得,确定长期强度最普遍的方法是等时应力-应变曲线法[15]，相等时间内岩石在不同应力作用下得到的蠕变变形与应力关系曲线即为等时应力-应变曲线。经过时间推移和应力增大，等时曲线会逐渐弯曲,非线性程度会愈加明显,曲线簇上会出现拐点，拐点对应的应力值即为岩石长期强度。

当岩石经过等速蠕变阶段进入加速蠕变阶段时，其内部损伤和变形加速发展，呈现出明显的黏塑性变形特征,少量可恢复的弹性变形和大量不可恢复的塑性变形同时出现。因此可将一个应力触发形式的非线性黏壶元件和摩擦键并联[16]，组成一个可以描述岩石加速蠕变阶段特性的非线性黏塑性元件，如图1所示。

图1 非线性黏塑元件

在岩石所受应力σ作用下，非线性黏塑性元件的蠕变方程为

(9)

式中：σs为应力阈值；a，b，η1为蠕变参数,由试验确定。

H(σ)为Heaviside单位阶跃函数，其表达式为

(10)

3 分数阶黏弹塑性模型

一般情况下,岩石的蠕变变形可简化为应力初加载时的瞬时变形、加载持续期间内的黏弹性变形、蠕变破坏期间的黏塑性变形3部分,因此用胡克体表征岩石的瞬时蠕变εe,分数阶软体元件表征岩石的黏弹性变形特征εve,带应力触发的非线性黏塑性元件表征岩石黏塑性变形εvp。将四元件串联形成一个新的分数阶黏弹塑模型(图2)。

图2 岩石分数阶黏弹塑性模型

当岩石所受应力未达到其屈服强度,即σ<σs时，蠕变变形处于衰减蠕变和等速蠕变阶段。此时岩石蠕变变形由瞬时弹性变形εe和黏弹性变形εve组成,蠕变总应变为

ε=εe+εve。

(11)

(1)胡克体对应的应力-应变关系为

(12)

式中：σ为胡克体的应力，即总应力；E1为胡克体中弹簧的弹性模量。

(2)分数阶软体元件的应力-应变关系为

(13)

式中：η1为其黏性系数:β为其求导阶数。

当σ≥σs时，岩石蠕变变形进入加速阶段，岩石除了出现少部分弹性和黏弹性变形外，还有比较明显的黏塑性变形εvp，此时蠕变总变形为

ε=εe+εve+εvp。

(14)

(3)黏塑性元件对应的应力-应变关系为

(15)

因此，综合3部分应变建立的岩石分数阶蠕变模型的表达式为

(16)

4 模型拟合分析

模型参数拟合分析主要有最小二乘法、退火算法、粒子群算法和全局优化算法等，其中以最小二乘法应用最为普遍。最小二乘法通过最小化误差的平方和寻找数据,从而得到匹配程度最优的函数。本次拟合采用基于最小二乘算法的MATLAB软件识别参数并绘制拟合曲线。

采用徐卫亚等[17]对锦屏一级水电站绿片岩进行三轴蠕变试验的试验数据和齐亚静等[18]对万州红层砂岩进行三轴蠕变试验的试验数据以及邱贤德等[19]对长山盐岩进行单轴压缩的试验数据进行参数拟合。为说明拟合过程，现以绿片岩试验数据为研究对象，砂岩和盐岩的参数拟合方法和绿片岩一致。绿片岩蠕变试验中其试验围压为15 MPa,主应力为100 MPa，岩石进入加速蠕变阶段的时间点为1.5 h，砂岩蠕变试验的围压为6 MPa,主应力为61 MPa，盐岩所受应力为14.41 MPa。

t<1.5 h时，岩石未到达加速蠕变阶段，先求解瞬时弹性模量E1。E1=σ/ε1,ε1为t=0时岩石瞬时加载产生的形变，可从蠕变曲线中求得。通过式(16)σ<σs时的公式进行数据拟合，求解到η1和β。

t>1.5 h时，岩石进入到加速蠕变阶段，此时利用式(16)中σ≥σs时的公式进行数据拟合，求解得到余下的参数。

σ<σs时，岩石处于衰减和等速蠕变阶段，岩石变形趋于稳定，黏塑性体此时不起作用。根据式(7)，将应力和黏滞系数设为定值，β为变量，得到岩石蠕变曲线，如图3所示。由于β不断增大，岩石蠕变特性逐渐增强。当β趋近于0时，分数阶软体元件类似理想弹性体；当β趋近于1时，分数阶软体元件类似理想流体。具体在实际工程中岩石蠕变特性表现各异。一般状态下，硬岩蠕变不明显，其性质和理想弹性体相近，适合用较小的β进行模拟；软岩在荷载应力作用下具有明显蠕变特点，适合用较大的β进行模拟。

图3 等速阶段下不同β值的蠕变-时间曲线

表1 绿片岩蠕变拟合参数

表2 红砂岩蠕变拟合参数

表3 长山盐岩蠕变拟合参数

由图4～6可知，该分数阶模型可以较好反映岩石衰减蠕变阶段随时间增长应变逐渐减小的特点,并且分数阶软体元件对岩石等速蠕变阶段内蠕变平滑缓慢增加的线性形变发展趋势描绘准确。当进入加速蠕变阶段时,岩石蠕变曲线曲率增大，同时岩石变形急剧扩展，轴向蠕变明显增加。而非线性黏塑性元件可准确描绘岩石应变发展随时间变化呈现指数型函数增长的特点。模型理论值和实际岩石蠕变值基本一致，拟合程度较高。说明推导的分数阶黏弹塑性模型能有效描述岩石的3阶段蠕变全过程。该模型参数数量少，求解简单，在实际工程中有一定应用性。

图4 绿片岩试验值与模型拟合值对比曲线

图5 砂岩试验值与模型拟合值对比曲线

图6 盐岩试验值与拟合值对比曲线

5 结 语

(1)基于分数阶理论，引入分数阶软体元件。并将该软体元件和弹簧元件串联，引入一个表征非线性加速阶段的黏塑性体，形成一个非线性黏弹塑性蠕变模型，从而更好反映岩石蠕变全过程。

(2)通过调整分数阶数β值大小，发现β值对于岩石等速阶段蠕变速率有较大影响。一般情况下，硬岩用较大β值进行模拟，软岩用较小β值进行模拟。

(3)结合绿片岩、万州红岩和长山盐岩的蠕变实验数据，运用MATLAB软件进行拟合，并对提出的基于分数阶导数的岩石蠕变模型参数进行识别，拟合结果表明理论值与实际值比较贴合，相关系数较高，具有一定的适用性。