周期荷载下岩石分数阶黏弹塑性本构模型研究

，，,， ，，

(1.贵州大学 土木工程学院，贵阳 550025；2.北京交通大学 基建与规划处，北京 100044)

1 研究背景

在工程地质和岩土工程领域中，经常会涉及周期荷载下岩体变形问题，如采动应力下煤矿井中煤壁、抽水蓄能电站上下游库水位的周期性变动、交通隧道及地震荷载下的地基岩体等。长期周期荷载作用下岩体强度会降低而破坏，造成工程事故。因此，开展周期荷载下岩石的变形特性研究具有一定的理论价值和实用意义。

目前，国内对周期荷载下岩石的变形特性的研究主要为理论推导、模型计算、试验验证3个方面。

(1)理论推导方面：郭建强等[1]借鉴蠕变理论，提出弹性、黏性、塑性3个疲劳元件，建立了单轴周期荷载下岩石Maxwell体、Kelvin体及非线性Bingham体串联而成的非线性黏弹塑性疲劳模型；张平阳等[2]基于Weibull分布对岩石损伤软化模型进行拓展，提出岩石循环加卸载本构模型；王春等[3]基于损伤力学原理，采用组合模型方法，建立半正弦波荷载下的一维高应力及重复冲击共同作用下岩石的本构模型。

(2)模型计算方面：刘文韬等[4]通过研究含孔隙和微裂纹的岩石剪胀效应和软化效应的应力-应变关系，提出岩石弹塑性连续损伤本构模型；赵怡晴等[5]基于岩土结构力学观点，运用流变力学理论，通过组合计算，建立了三轴压缩荷载条件下节理岩体损伤本构模型；孟红霞等[6]基于动量和能量守衡原理及岩石冲击开裂模拟试验装置的工作原理，通过模拟计算，建立了岩石动态损伤峰值压力和加压速率计算模型，并通过岩心冲击开裂试验验证模型的适用性。

(3)试验验证方面：杨永杰等[7]对循环荷载作用下煤岩强度及变形特征开展了试验研究，认为煤岩损伤是煤岩中微裂纹扩展和产生后汇聚为宏观裂隙的过程，循环荷载下煤岩的损伤可以分为3个阶段，并基于应变等效性假设提出煤岩损伤本构方程；唐礼忠等[8]对在高应力状态下受小幅循环动力的大理岩进行了扰动力学试验，从损伤力学角度研究岩石变形与循环荷载的关系，提出大理岩损伤疲劳本构模型；王军保等[9]为研究低频周期荷载下岩盐的蠕变特性，对岩盐试件进行恒轴压、三角波循环围压下三轴压缩试验，建立低频循环荷载下Burgers模型，该模型拟合曲线与试验结果相似，但模型仅能描述岩石变形的减速、等速2个变形阶段。

从目前的研究成果来看，周期荷载下岩石疲劳本构方面取得很多成果。但元件疲劳模型方面的研究相对较少，相关方面的研究仍需加强。为此，基于流变力学相关理论，试图建立一种周期荷载下岩石疲劳蠕变本构模型，以描述周期荷载下岩石变形规律。研究成果不仅可丰富岩石力学理论，亦可给相关方面的研究提供一些参考。

2 岩石疲劳本构模型

为便于研究，在围压应力状态下，可将等幅周期荷载作用下岩石所受的偏应力等效于应力函数,即

式中：σav为周期荷载的不考围压平均值；σ1max为周期荷载应力上限值；σs为复杂应力状态下岩石的临界强度值[用(σ1-σ3)表示]，由周期荷载下岩石刚好发生疲劳破坏时(σ1-σ3)-ε曲线所确定的最大应力值；σc为复杂应力状态下静载的屈服强度值，即为静力荷载下(σ1-σ3)-ε曲线上的峰值应力；σ1,σ3分别为岩石所受的大、小主应力；f为周期荷载的频率。

由式(1)可知，当周期荷载的峰值偏应力(σ1max-σ3)>σs时等效偏应力(σ1-σ3)>σav，岩石可发生疲劳破坏；当周期荷载的峰值偏应力(σ1max-σ3)≤σs时，等效偏应力(σ1-σ3)≤σav，岩石不发生疲劳破坏；且σ1-σ3的值随周期荷载的频率增大而增大，说明式(1)可反映周期荷载作用。

当围压应力σ3=0时，周期荷载可等效成式(2)所示的应力函数[1]。

式中：σmax为周期荷载应力上限值；σs为单向应力状态下岩石破坏门槛值(即临界值)；σc为岩石单轴抗压强度。

前人的研究结果表明[7,9-10]周期荷载下岩石疲劳变形可划分为减速、等速、加速3个阶段(见图1)。当岩石所受的周期荷载σ1max-σ3小于岩石发生疲劳破坏时的临界强度值σs时，岩石的变形只有减速、等速2个变形阶段；反之，岩石发生疲劳破坏，变形具有减速、等速、加速3个阶段，体现岩石黏弹塑性的特征。

图1 岩石疲劳变形曲线Fig.1 Fatigue deformation curve of rock

西原模型[11]由广义Kevin体及黏塑性体串联而成，可以描述岩石黏弹特性，却无法描述岩石的非线性变形特征。为能够描述岩石在周期荷载下疲劳变形规律，基于流变力学相关理论，将西原模型中黏塑性体中的常值黏壶替换分数阶黏壶，并串联一个黏性元件，建立一种可描述不同周期荷载下岩石变形规律的分数阶黏弹塑性本构模型，并将此模型推广到复杂应力状态下，推导出三轴条件下岩石疲劳本构方程。

图2 分数阶黏壶Fig.2 Fractional-order dashpot

(3)

式中：η为分数阶黏壶黏滞系数，其物理量纲为[应力·时间]；n为分数阶微分的阶数，反映岩石加速阶段应变速率的材料参数；N为荷载循环系数，N=ft；η,n为材料参数，可通过试验测定。

当σ(t)=σ时，即所受力保持不变，根据分数阶微积分基本理论，可得分数阶黏壶的本构方程式为

(4)

式中：ε(t)为t时刻该黏壶的应变量；ε(N)为周期荷载作用了N时该黏壶的应变量。

对于式(4)，该本构方程式可控制岩石加速变形阶段应变速率。当n≥1时，该元件应变本构ε(N)为n的增函数，n值越大，黏壶的应变速率越大。

2.1 一维岩石疲劳本构方程

一维岩石分数阶黏弹塑性本构模型如图3所示。当等效应力σ>σav时，分数阶黏壶参与变形，模型为描述岩石3个变形阶段的分数阶黏弹塑性本构模型；当σ≤σav时，分阶黏壶不参与变形，模型为描述岩石等速、减速变形的Burgers模型。

图3 一维岩石分数阶黏弹塑性本构模型Fig.3 One-dimensional fractional-order visco-elasto- plastic constitutive model(FVEPCM) for rock

2.1.1 本构方程分解

岩石分数阶黏弹塑性模型由Maxwell、Kelvin及分数阶Bingham模型组成，分别如图3的Ⅰ,Ⅱ,Ⅲ所示。Maxwell模型由黏性和弹性元件串联而成；Kelvin由弹性和黏性元件并联组成；分数阶Bingham模型由分数阶黏壶和塑性元件并联组成。对于组合元件,并联各组件的应变相等，总应力等于各个元件应力之和；串联组件的各元件上的应力相等，总应变等于各组件应变之和。

由于弹性、黏性元件的状态方程为：

σ=Eε(N) ；

(5)

(6)

则各模型的本构方程如下所示。

(7)

(8)

Ⅲ：

(9)

式中：η1，η2，η3分别为Ⅰ，Ⅱ，Ⅲ中黏壶的黏性系数，量纲为[应力·时间]；EM，EK分别为Ⅰ、Ⅱ中弹性系数，量纲为[应力]。

2.1.2 岩石疲劳本构方程

(1)当σ≤σav时，此时岩石变形规律特性由Ⅰ，Ⅱ 2部分组成的Burgers模型描述，由式(7)、式(8)可得基于图3所示模型的本构方程,即

(2)当σ>σav时，此时岩石变形规律由Ⅰ，Ⅱ，Ⅲ 3部分组成的分数阶黏弹塑性本构模型描述，由式(7)、式(8)、式(9)可得基于图3所示模型的岩石疲劳本构方程，即

(11)

2.2 三维岩石疲劳本构方程

在三维受力σij状态下，岩石力学参数不同于一维受力状态。三维模型(见图4)中的参数的意义与一维不同。模型中G1，G2分别为弹性与黏弹性剪切弹性系数，反映岩石的剪切变形特性的指标，量纲为[应力]；η1,η2,η3分别为黏性、黏弹性、分数阶黏塑性剪切黏滞系数。

图4 三维岩石分数阶黏弹塑性本构模型Fig.4 Three-dimensional FVEPCM for rock

2.2.1 基本元件

2.2.1.1 弹性元件

根据弹塑性力学基本理论，可得弹性元件(图4中①)的状态方程为

(12)

(13)

式中σkk为应力张量第一不变量。

2.2.1.2 黏性元件

对于黏性元件(图4中②)，不考虑应力球张量下黏性元件的体积变形，则其状态方程为

(15)

对式(15)积分则可得

(16)

2.2.1.3 黏弹性模型

对于黏弹性模型(图4中③)，不考虑应力球张量下黏弹性模型的体积变形，其状态方程式为

(17)

由式(17)可得黏弹性体本构方程为

2.2.1.4 分数阶黏塑模型

(19)

2.2.2 分数阶黏弹塑性疲劳本构方程

(1)当σ1-σ3≤σav时，岩石不发生破坏，由式(14)—式(18)可得基于图4所示模型所建立的本构方程式为

(20)

式中εij(N)为模型总应变张量。

(2)当σ1-σ3>σav时，岩石发生破坏，由式(14)—式(19)可得基于图4所示模型所建立的岩石疲劳本构方程式为

(21)

2.2.3 三轴条件下岩石疲劳本构方程

三轴循环压缩条件下岩石等效应力状态为

(22)

式(22)的等效偏应力为

(23)

式中：σm为平均应力，σm=(σ1+2σ3)/3；δij为Kronecker等号。

只考虑岩石弹性体积应变，不考虑黏性、黏弹性及分数阶黏塑性体积应变。则三轴压缩状态下，岩石总应变为

(24)

将式(21)代入式(24)得三轴压缩条件下岩石发生变形时，试件的岩石轴向疲劳本构方程为

(25)

3 模型参数的求取及验证

3.1 一维状态下模型拟合及参数的确定

肖建清等[13]利用岩石试验机对频率为0.2 Hz正弦波循环荷载下花岗岩变形特性进行了研究，试件F1-7的循环荷载上、下限应力分别为136.26 MPa和43.03 MPa，通过试验岩石单轴抗压强度为143.43 MPa；赵凯等[14]利用动力试验机对不同频率荷载下石灰岩疲劳特性进行了研究，其中DD-53试件的单轴抗压强度为104.32 MPa,所受正弦波的循环荷载频率为1 Hz，荷载上、下限应力分别为31.30 MPa和93.89 MPa；杨永杰等[7]采用电液伺服岩石试验系统对煤岩试件进行试验研究，得到单轴状态下煤岩临界值不大于单轴抗压强度的81%，本文临界值取为单轴抗压强度的81%。其中BDC2试件的单轴抗压强度为28.4 MPa，所受正弦波的循环荷载频率为0.5 Hz，荷载上、下限应力分别为0.568 MPa和23.04 MPa。提出岩石疲劳本构方程式(11)对以上3个试验数据进行拟合和模型参数的求取，拟合曲线见图5，参数见表1。

如图5及表1所示，由分数阶黏弹塑性疲劳本构方程在对岩石疲劳变形曲线的拟合结果可知：①模型拟合参数均为正，无一负值，且对大部分岩石的拟合相似度达96%以上，但对煤岩疲劳曲线的拟合系数较低，可能因为数据的获取误差造成的；②各岩石的模型参数n从大到小依次是煤岩>石灰岩>花岗岩，表明分数阶黏壶参数n随岩石的单轴抗压强度的增大而减小；③图5中，煤岩的加速阶段应变速率较其他岩石大，其原因为煤岩的拟合参数n较其他岩石大，导致岩石的加速变形阶段应变速率较其他岩石更陡。

图5 一维FVEPCM模型对花岗岩、石灰岩、煤岩疲劳变形的拟合曲线Fig.5 Fitting of fatigue deformation of granite limestone,and coal rock by one-dimensional FVEPCM

表1 一维FVEPCM模型对文献[13]、文献[14]、文献[7]的拟合参数Table 1 Fitted parameters of one-dimensional FVEPCM from literature[13], [14], and [7]

表2 三维FVEPCM模型对文献[15]、文献[16]的拟合参数Table 2 Fitted parameters of three-dimensional FVEPCM from literature[15] and [16]

3.2 三维模型拟合及参数的确定

章清叙等[15]利用岩石多功能试验机对受围压应力下红砂岩的变形特性进行了循环荷载试验研究，其中试件RS-4-3#的加载频率为0.2 Hz，所受的循环荷载上、下限应力分别为51.8 MPa和99.1 MPa，利用Origin软件可得15 MPa围压下红砂岩的静载屈服强度σc=98.22 MPa及屈服临界值σs=104.653 MPa。提出岩石疲劳本构方程式(25)对试件RS-4-3#疲劳曲线进行拟合和模型参数的求取，拟合曲线见图6(a)，参数见表2。

丁祖德等[16]通过自制三轴围压系统及MTS单轴疲劳试验机对富水砂质泥岩在正弦波周期荷载下变形的研究，试验得到围压σ3=200 kPa下软岩的临界动应力幅值在180～300 kPa之间，本文取临界值σs=300 kPa。葛修润等[17]根据其试验得出在三围状态下岩石的临界应力为岩石静载屈服值的0.85～0.9，本文取比值为0.85，则软岩静载屈服值为σc=353 kPa。提出岩石疲劳本构方程式(25)对围压σ3=200 kPa、加载频率f=3 Hz、静偏应力σsp=180 kPa及所受循环荷载动应力σd=300 kPa下岩石曲线进行拟合和模型参数的求取，拟合曲线见图6(b)，参数见表2。此种情况下岩石轴向荷载为

(26)

图6 三维FVEPCM模型对红砂岩、富水砂质泥岩疲劳变形的拟合曲线Fig.6 Fitting of fatigue deformation of red sandstone and water rich sandy mudstone by three-dimensional FVEPCM

由表2、图6可知，分数阶黏弹塑性本构模型在三维受力状态下对岩石的变形有很好的拟合结果。由拟合结果(表2)可知：① 2种岩石的拟合相关系数都在0.995以上；②拟合结果n值富水砂质泥岩>红砂岩表明分数阶黏壶参数n随岩石三轴静载强度的增大而减小。

4 结 论

(1)基于流变力学理论，提出可描述周期荷载下不同岩石完整疲劳变形规律的一维岩石分数阶黏弹塑性本构模型，一维岩石疲劳本构模型推广到复杂应力状态下岩石疲劳本构模型，并推导出三轴条件下分数阶黏弹塑性本构模型。

(2)单轴条件下：当σ>σav时，模型为反映岩石减速、等速、加速3个变形阶段变形规律的分数阶黏弹塑性模型；反之，则为反映岩石减速、加速变形规律的Burgers模型。

(3)三轴条件下：当σ1-σ3>σav时，模型为反映岩石减速、等速、加速3个变形阶段变形规律的分数阶黏弹塑性模型；反之，则为反映岩石减速、加速变形规律的Burgers模型。

(4)用相关文献的岩石疲劳变形试验结果对分数阶黏弹塑性本构模型的合理性进行了验证。结果表明：该模型能够很好地描述岩石疲劳变形全过程的3个阶段，且拟合曲线和试验曲线吻合良好，误差较小，其拟合相关系数在0.96以上。

本文模型仅描述了周期荷载循环次数与单轴轴向应变、三轴轴向应变之间的关系，并没有对其他方面因素进行考虑；模型所描述的曲线呈倒“S”形，描述效果较好，因减速阶段反应时间短，反应迅速，该模型对一些岩石初始减速阶段的拟合不是太好，以及对其他岩石疲劳变形的异形曲线的拟合效果需要进一步研究。