基于改进分数阶黏滞体的岩石非线性蠕变模型

王军保，刘新荣，王铁行



基于改进分数阶黏滞体的岩石非线性蠕变模型

王军保1,2，刘新荣1,2，王铁行1

(1. 西安建筑科技大学土木工程学院，陕西西安，710055；2. 重庆大学土木工程学院，重庆，400045)

通过将表征应力水平对岩石非线性蠕变特性影响的函数引入到常规分数阶黏滞体的本构关系中，提出一种改进的分数阶黏滞体；基于岩石非线性流变力学理论，提出一个能够描述岩石加速蠕变的非线性黏滞体，其黏滞系数随应力水平提高和蠕变时间延长而逐渐减小；将改进的分数阶黏滞体、非线性黏滞体与基本弹性体和塑性体进行组合，建立一个新的4元件非线性黏弹塑性蠕变模型，给出模型的蠕变方程，并利用岩石蠕变试验结果对模型合理性进行验证。研究结果表明：该模型能够较好地描述岩石蠕变全过程的3个阶段，且拟合曲线和试验曲线吻合良好，误差较小。

岩石；蠕变模型；改进的分数阶黏滞体；非线性黏滞体

岩石元件组合蠕变模型有助于从概念上认识岩石变形的弹性、黏性和塑性分量，可以将岩石复杂的力学性质直观的表现出来，且其本构方程形式简单、参数物理意义明确，因而得到了广泛应用。元件组合模型的本构关系本质上是一种整数阶微分型本构关系，分数阶微积分理论是研究任意阶微分和积分的理论，是整数阶微积分向任意阶的推广。因此，采用分数阶微积分理论来解决岩石蠕变本构建模问题，其在曲线形式的描述方面将具有明显优势[1]。鉴于分数阶微积分在岩石蠕变本构建模方面的优越性，近年来一些学者借鉴元件组合模型的建模思路，以分数阶黏滞体代替经典黏滞体来描述岩土材料的蠕变特性[1−6]，为岩土蠕变模型研究提供了一种新思路。但应当看到，用分数阶黏滞体代替经典黏滞体来建立蠕变模型，实质上相当于考虑了黏滞系数的时间相关性[6]。而岩石蠕变过程中，其黏滞系数的变化不仅与时间有关，还与所受荷载水平有关[7−11]；此外，利用分数阶黏滞体和经典元件进行组合得到的蠕变模型虽然可以很好地拟合岩石衰减蠕变和稳态蠕变阶段的试验曲线，但仍然无法描述岩石的加速蠕变。基于此，本文作者对常规分数阶黏滞体的本构关系进行了改进，使之能够反映应力水平对岩石非线性蠕变特性的影响；在此基础上，提出了一个能够描述岩石加速蠕变的非线性黏滞体，其黏滞系数将随应力水平提高和蠕变时间延长而逐渐减小；之后，将改进的分数阶黏滞体、非线性黏滞体与经典弹性元件和塑形元件进行组合，建立了一个能够反映岩石蠕变全过程的4元件非线性黏弹塑性蠕变模型，并用相关文献岩石蠕变试验结果对模型的合理性进行了验证。

1 分数阶黏滞体的改进

1.1 分数阶微积分基本定义

分数阶微积分有多种不同的定义方法，其中在岩石流变本构模型研究中应用较为广泛的是Riemann-Liouville型分数阶微积分算子理论[2−6]。按照该理论，函数()在可积空间[0，]上的阶积分定义为[12−13]

相应的，函数()的阶微分定义为

即

式中：为微分阶数，＞0，且−1＜≤(为正整数)；为伽马函数，当＞0时，定义为

当()在=0附近可积，且0≤≤1时，以上分数阶微积分的Laplace变换公式为

式中：()为()Laplace变换。

1.2 分数阶黏滞体

对于分数阶黏滞体，尽管不同文献有不同的命名方式，但其本构关系是一致的，即

式中：为应力；为应变；为时间；为分数阶微分阶数，0≤≤1；η为与黏滞系数类似的参数，本文称之为类黏滞系数，单位为“应力单位∙时间单位”。

由式(6)可见：当=0时，令η­=，则有=，为弹性模量，此时分数阶黏滞体代表理想固体；当=1时，令η­=，则有，为黏滞系数，此时分数阶黏滞体代表理想流体；当0＜＜1时，分数阶黏滞体可以模拟处于理想固体和理想流体之间如岩土材料等的变形情况。

当应力保持不变时，根据Riemann−Liouville型分数阶微积分算子理论，将式(6)两侧进行分数阶积分，可得

式(7)即为分数阶黏滞体的蠕变方程。

取=10 MPa，η­=4 000 MPa∙h，根据式(7)可得到分数阶黏滞体在取不同值时的一组蠕变曲线，如图1(a)所示。由图1(a)可见：单个分数阶黏滞体元件就可以表现出黏弹性力学性质。当较小时，模型蠕变变形量也较小；随着增大，相同时刻蠕变量和蠕变速率均逐渐增大，同时衰减蠕变阶段曲线的曲率半径逐渐增大，经历的时间逐渐延长。对于蠕变性不同的岩石材料，可通过调整来改变蠕变曲线的形状，从而达到精确拟合试验结果的目的，这正是用分数阶黏滞体描述岩石蠕变特性的优势所在。

(a) γ的影响；(b) σ的影响

取η­=4 000 MPa∙h，=0.15，根据式(7)可得到分数阶黏滞体在不同应力水平下的蠕变曲线，如图1(b)所示。由图1(b)可以看出：随着应力水平提高，蠕变变形量和蠕变速率等均逐渐增大。但是，相同时刻蠕变量随应力水平提高而呈线性规律增大，即对于某一时刻来说，在相同应力增量情况下，图1(b)中。这一情况通常在岩石所受应力水平较低时才近似成立，当应力水平相对较高时，根据同一岩石试样在不同应力水平下的蠕变试验结果[14−16]，在相同应力增量情况下，对于同一时刻，通常有。因此，常规分数阶黏滞体无法反映应力水平对岩石非线性蠕变特性的影响。

1.3 改进的分数阶黏滞体

利用分数阶黏滞体代替常规黏滞体来建立岩石蠕变模型，实际上相当于考虑了黏滞系数的时间相关性；而在岩石蠕变过程中，其黏滞系数不仅与蠕变时间有关，还与荷载水平有关。熊良宵等[10−11]研究表明：岩石黏滞系数与应力的幂次方成倒数关系。因此，为了使分数阶黏滞体能够反映应力水平对岩石非线性蠕变特性的影响，本文假定分数阶黏滞体的类黏滞系数η与应力之间存在如下关系

式中：η0为类黏滞系数的初始值；为反映应力水平对岩石非线性蠕变特性影响的非线性系数，且≥1；0为单位应力，取为1，量纲与相同。

将式(8)代入式(6)并整理可得改进后的分数阶黏滞体的本构关系为

当应力保持不变时，将式(9)两侧进行分数阶积分，根据Riemann−Liouville型分数阶微积分算子理论，可得

式(10)即为改进分数阶黏滞体的蠕变方程。对比式(10)和式(7)可以看到：当=1时，本文改进的分数阶黏滞体可退化为常规分数阶黏滞体。

取=10 MPa，=0.15，η0=4 000 MPa∙，利用式(10)可得到改进分数阶黏滞体在系数取不同值时的一组蠕变曲线，如图2(a)所示；取=0.15，=1.8，η0= 4 000 MPa∙，利用式(10)可得到改进分数阶黏滞体在不同应力水平下的蠕变曲线，如图2(b)所示。

(a) n的影响；(b) σ的影响

由图2可以看出：随着增大和应力水平提高，蠕变曲线变化规律是一致的，即越大，应力水平越高，蠕变变形量和蠕变速率越大。并且，相同时刻的蠕变量随应力水平提高而呈非线性加速增大趋势，即对于某一时刻来说，在相同应力增量情况下，图2(b)中，这与实际情况较为符合。因此，改进后的分数阶黏滞体能够更好地反映应力水平对岩石非线性蠕变特性的影响。

2 可描述加速蠕变的非线性黏滞体

研究表明[17]，在较高应力水平下，荷载持续作用会导致岩石内部裂纹演化扩展和损伤逐步积累，从而导致其黏滞系数不断减小，宏观上表现为岩石蠕变变形随时间增加而非线性加速增大；并且，黏滞系数的减小过程不仅与蠕变时间有关，还与应力水平有关，应力水平越高，减小越快。从目前对于岩石非线性蠕变模型的研究来看，元件非线性化过程主要是采用经验公式的方法。借鉴这一建模思路，经过仔细分析，本文提出了一个新的非线性黏滞体，其黏滞系数是应力及其作用时间的衰减函数，即

式中：0为黏滞系数初始值；和分别为反映应力和时间影响的非线性系数，且满足≥1，/0≤1；0为单位应力，取为1，量纲与相同；0为单位时间，取为1，量纲与相同。

在应力恒定的情况下，将式(11)对时间求一阶导数并整理可得

由于/0≤1，则′(，)≤0，说明随时间增加而单调递减。

结合常规黏滞体的本构关系及式(11)，可得非线性黏滞体的本构方程为

在应力保持不变的情况下，将式(14)对时间积分，并考虑初始条件：=0，=0，则可解得本文非线性黏滞体的蠕变方程为

由式(14)可以看出：当=1和=0时，本文非线性黏滞体实际上退化为常规黏滞体；当=1和＞0时，非线性黏滞体可反映蠕变时间对岩石非线性蠕变特性的影响；当＞1，=0时，非线性黏滞体可反映应力水平对岩石非线性蠕变特性的影响；当＞1和＞0时，非线性黏滞体可综合反映应力水平和蠕变时间的影响。

取=10 MPa，=1.5，0=1 000GPa·h，利用式(14)可得到系数取不同值时非线性黏滞体的一组蠕变曲线，如图3(a)所示；取=10 MPa，=0.005，0=1 000GPa·h，利用式(14)可得到系数取不同值时非线性黏滞体的一组蠕变曲线，如图3(b)所示。

(a) m的影响；(b) p的影响

由图3可以看到：参数和对岩石加速蠕变曲线的影响非常大，且两者对蠕变曲线的影响规律是一致的。在同一值(或值)情况下，随着时间延长，蠕变量和蠕变速率均呈非线性加速增大；随着(或)增大，相同时刻岩石蠕变变形量和蠕变速率均增大，岩石蠕变曲线的加速特征越来越明显。这表明系数和充分反映了岩石非线性加速蠕变特性，(或)越大，岩石黏滞系数衰减越快，破坏也就越快。但刚开始一段时间内(约40 h之内)，参数对蠕变曲线的影响不如参数对蠕变曲线的影响大。

3 岩石非线性蠕变模型及验证

3.1 非线性蠕变模型的建立

为了较好地描述岩石蠕变全过程，本文将前面提出的改进分数阶黏滞体、非线性黏滞体与基本弹性体和塑性体进行组合，建立了一个4元件的非线性黏弹塑性蠕变模型。考虑到模型参数过多会给参数求解增加难度，本文取改进分数阶黏滞体中反映应力水平影响的非线性系数和非线性黏滞体中反映应力水平影响的非线性系数为同一值，且以来表示。非线性黏弹塑性蠕变模型示意图如图4所示。由图4可见：非线性模型由弹性体(第I部分)、改进的分数阶黏滞体(第II部分)和由非线性黏滞体与塑性体并联组成的非线性黏塑性体(第III部分)等3部分串联组成。其中，弹性体描述岩石的瞬时变形，改进的分数阶黏滞体反映岩石的衰减及稳态蠕变变形，而非线性黏塑性体则模拟岩石加速蠕变变形。

图4 非线性黏弹塑性蠕变模型

假定在=0时刻对模型施加恒定应力，当≤s时，模型中第III部分不参与蠕变，此时有

式中：和分别为模型总应力和总应变；和分别为第I部分和第II部分的应力；和分别为第I部分和第II部分的应变。

对于第I部分，其本构关系为

式中：为弹性模量。

对于第II部分，结合式(10)可得其蠕变方程为

当0＞s时，模型中所有部分均参与蠕变，则有

此时模型中前2部分的蠕变方程与≤s时相同。对于第III部分，结合常规黏塑性体的蠕变方程及式(14)可得其蠕变方程为

当≤s时，由式(15)~(17)可得非线性模型的蠕变方程为

当＞s时，由式(16)~(19)可得非线性模型的蠕变方程为

3.2 模型验证

为了说明本文模型的合理性和有效性，利用文献[15]不同应力水平下冻结软岩的单轴蠕变试验结果对模型进行了验证。软岩蠕变试验曲线见图5散点曲线，图5中数字为轴向应力。由试验曲线可以看出：当应力为2.8 MPa和3.2 MPa时，软岩仅出现了前2个蠕变阶段，即衰减蠕变和稳态蠕变；当应力为4.5 MPa时，软岩蠕变过程出现了完整的3个阶段。根据软岩蠕变试验结果，基于最小二乘法原理，利用1stOpt数学优化软件，采用自定义函数的方法对本文模型参数进行了反演拟合，软岩强度参数s=3.24 MPa。参数反演结果见表1，拟合曲线和试验曲线对比情况见图5。

图5 拟合曲线与试验曲线对比

表1 蠕变参数反演结果

由表1可见：本文非线性蠕变模型参数拟合相关性较高，相关系数均在0.99以上，接近于1。另外，尽管在不同应力水平下参数，η0，和取值有一定的变化，但总体来看变化不大，这说明对于某一种岩石而言，这些参数可近似认为是常量。同时，由图5可以看出：本文非线性模型能够很好地描述岩石蠕变全过程的3个阶段，且拟合曲线和试验曲线吻合的相当理想，误差较小，从而也证明了其合理性和有 效性。

4 结论

1) 在常规分数阶黏滞体本构关系中引入了表征应力水平对岩石非线性蠕变特性影响的函数，提出了一种改进的分数阶黏滞体。

2) 基于非线性流变力学理论，提出了一个能够模拟岩石加速蠕变新的非线性黏滞体，其黏滞系数是应力及其作用时间的衰减函数。

3) 将改进的分数阶黏滞体、非线性黏滞体与基本弹性体和塑性体进行组合，建立了一个新的4元件非线性黏弹塑性蠕变模型。

4) 用相关文献岩石蠕变试验结果对非线性模型的合理性进行了验证。结果表明：该模型能够很好的描述岩石蠕变全过程的3个阶段，且拟合曲线和试验曲线吻合良好，误差较小。

[1] 何利军, 孔令伟, 吴文军, 等. 采用分数阶导数描述软黏土蠕变的模型[J]. 岩土力学, 2011, 32(S2): 239−245.HE Lijun, KONG Lingwei, WU Wenjun, et al. A description of creep model for soft soil with fractional derivative[J]. Rock and Soil Mechanics, 2011, 32(S2): 239−245.

[2] 殷德顺, 任俊娟, 和成亮, 等. 一种新的岩土流变模型元件[J]. 岩石力学与工程学报, 2007, 26(9): 1899−1903.YIN Deshun, REN Junjuan, HE Chengliang, et al. A new rheological model element for geomaterials[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(9): 1899−1903.

[3] 孙海忠, 张卫. 一种分析软土黏弹性的分数导数开尔文模型[J]. 岩土力学, 2007, 28(9): 1983−1986.SUN Haizhong, ZHANG Wei. Analysis of soft soil with viscoelastic fractional derivative Kelvin model[J]. Rock and Soil Mechanics, 2007, 28(9): 1983−1986.

[4] 刘林超, 闫启方, 孙海忠. 软土流变特性的模型研究[J]. 岩土力学, 2006, 27(S1): 214−217.LIU Linchao, YAN Qifang, SUN Haizhong. Study on model of rheological property of soft clay[J]. Rock and Soil Mechanics, 2006, 27(S1): 214−217.

[5] 陈亮, 陈寿根, 张恒, 等. 基于分数阶微积分的非线性黏弹塑性蠕变模型[J]. 四川大学学报(工程科学版), 2013, 45(3): 7−11. CHEN Liang, ZHANG Shougen, ZHANG Heng, et al. A nonlinear viscoelasto−plastic creep model based on fractional calculus[J]. Journal of Sichuan University (Engineering Science Edition), 2013, 45(3): 7−11.

[6] 康永刚, 张秀娥. 岩石蠕变的非定常分数伯格斯模型[J]. 岩土力学, 2011, 32(11): 3237−3241.KANG Yonggang, ZHANG Xiue. Nonstationary parameter fractional Burgers model of rock creep[J]. Rock and Soil Mechanics, 2011, 32(11): 3237−3241.

[7] 熊良宵, 杨林德, 张尧. 岩石的非定常Burgers 模型[J]. 中南大学学报(自然科学版), 2010, 41(2): 679−684.XIONG Liangxiao, YANG Linde, ZHANG Yao. Non-stationary Burgers model for rock[J]. Journal of Central South University (Science and Technology), 2010, 41(2): 679−684.

[8] 孙钧. 岩土材料流变及其工程应用[M]. 北京: 中国建筑工业出版社, 1999: 545−550.SUN Jun. Rheological behavior of geomaterials and its engineering applications[M]. Beijing: China Architecture and Building Press, 1999: 545−550.

[9] 曹平, 刘业科, 蒲成志, 等. 一种改进的岩石黏弹塑性加速蠕变力学模型[J]. 中南大学学报(自然科学版), 2011, 42(1): 142−146.CAO Ping, LIU Yeke, PU Chengzhi, et al. An improved accelerated creep mechanical model of viscoelasto-plastic rock[J]. Journal of Central South University (Science and Technology), 2011, 42(1): 142−146.

[10] 熊良宵, 杨林德. 硬脆岩的非线性黏弹塑性流变模型[J]. 同济大学学报(自然科学版), 2010, 38(2): 188−193.XIONG Liangxiao, YANG Linde. Viscoelasto-plastic rheological model for brittle hard rock[J]. Journal of Tongji University (Natural Sciences), 2010, 38(2): 188−193.

[11] 赵延林, 曹平, 文有道, 等. 岩石弹黏塑性流变试验和非线性流变模型研究[J]. 岩石力学与工程学报, 2008, 27(3): 477−486.ZHAO Yanlin, CAO Ping, WEN Youdao, et al. Viscoelasto-plastic rheological experimental and nonlinear rheological model of rocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(3): 477−486.

[12] Yang P, Lam Y C, Zhu K Q. Constitutive equation with fractional derivatives for the generalized UCM model[J]. Journal of Non-Newtonian Fluid Mechanics, 2010, 165(3/4): 88−97.

[13] Metzlera R, Theo F N. Fractional relaxation processes and fractional rheological models for the description of a class of viscoelastic materials[J]. International Journal of Plasticity, 2003, 19(7): 941−959.

[14] 宋飞, 赵法锁, 李亚兰. 一种岩石损伤流变模型及数值分析[J]. 水利水电科技进展, 2008, 28(1): 12−15.SONG Fei, ZHAO Fasuo, LI Yalan. Damage rheological constitutive model for numerical simulation to gypsum breccias[J]. Advances in Science and Technology of Water Resources, 2008, 28(1): 12−15.

[15] 李栋伟, 汪仁和, 范菊红. 白垩系冻结软岩非线性流变模型试验研究[J]. 岩土工程学报, 2011, 33(3): 398−403.LI Dongwei, WANG Renhe, FAN Juhong. Nonlinear rheological model for frozen soft rock during Cretaceous period[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(3): 398−403.

[16] 于怀昌, 李亚丽, 刘汉东. 粉砂质泥岩常规力学、蠕变以及应力松弛特性的对比研究[J]. 岩石力学与工程学报, 2012, 31(1): 60−70.YU Huaichang, LI Yali, LIU Handong. Comparative study of conventional mechanical, creep and stress relaxation properties of silty mudstone under triaxial compression[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(1): 60−70.

[17] 蒋昱州, 张明鸣, 李良权. 岩石非线性黏弹塑性蠕变模型研究及其参数识别[J]. 岩石力学与工程学报, 2008, 27(4): 832−839.JIANG Yuzhou, ZHANG Mingming, LI Liangquan. Study on nonlinear viscoelasto-plastic creep model of rock and its parameter identification[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(4): 832−839.

(编辑 杨幼平)

A nonlinear creep model for rocks based on modified fractional viscous body

WANG Junbao1, 2, LIU Xinrong1, 2, WANG Tiehang1

(1. College of Civil Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China; 2. College of Civil Engineering, Chongqing University, Chongqing 400045, China)

By introducing a stress influence function to the constitutive relation of the conventional fractional viscous body, a modified fractional viscous body which could reflect the influence of stress level on the creep property of rocks was proposed. Moreover, a new nonlinear viscous body which could describe the tertiary creep of rocks was put forward according to the nonlinear rheology theory, the viscous coefficient of which decreased with the increase of stress and time.Then, through connecting the modified fractional viscous body and the nonlinear viscous body with the traditional elastic and plastic bodies in forms of series and parallel, a new nonlinear viscoelasto-plastic creep model with 4 components for rocks was established, the rationality of which was verified by using creep test data of rocks. The results indicate that the proposed model can well reflect the three stages of rock creep process, and the fitting curves and test curves conform very well to each other.

rock; creep model; modified fractional viscous body; nonlinear viscous body

10.11817/j.issn.1672-7207.2015.04.037

TU45

A

1672−7207(2015)04−1461−07

2014−04−13；

2014−06−20

国家自然科学基金资助项目(51404184)；陕西省教育厅专项科研计划项目(14JK1401)(Project (51404184) supported by the National Natural Science Foundation of China; Project (14JK1401) supported by the Specialized Scientific Research Program of Education Department of Shaanxi Provincial Government)

王军保，博士，讲师，从事岩石力学与地下工程方面研究；E-mail：xajdwangjunbao@163.com