干湿循环条件下岩石蠕变本构模型研究及应用

侯 迪,许 珍,莫 至 坤,姚 池,叶 志 伟

(1.贵州水利水电勘测设计研究院,贵州 贵阳 550002; 2.江西省尾矿库工程安全重点实验室,江西 南昌 330031; 3.南昌大学 工程建设学院,江西 南昌 330031)

0 引 言

岩石蠕变是水电工程边坡变形失稳的主要原因之一[1-3]。水电工程库水位周期性涨落使库区边坡岩石反复经历干湿循环过程,岩石蠕变特性随即发生显著变化,进而影响库区边坡的蠕变过程及稳定性。研究干湿循环条件下岩石蠕变特性对于保持边坡稳定及保障库坝运行安全具有十分重要的意义。

开展岩石蠕变试验是研究干湿循环条件下岩石蠕变特性的重要手段。Liu等[4]开展了红砂岩单轴蠕变试验研究,得到岩石在加速蠕变阶段的蠕变破坏时间和蠕变速率与干湿循环次数的关系。马芹永等[5]通过单轴蠕变试验,研究了不同干湿循环次数下粉砂岩的蠕变特性,发现随着干湿循环次数增多,粉砂岩轴向蠕变应变和轴向稳态蠕变速率呈非线性增大,瞬时变形模量呈对数降低。李安润等[6]开展了干湿循环条件下泥岩的剪切蠕变试验研究,结果表明相同剪切荷载条件下,泥岩弹性模量、黏弹性模量、黏滞系数等蠕变参数都随着干湿循环次数的增加而减小。试验研究结果虽然能揭示岩石力学性质随干湿循环次数变化规律,但都是针对某种特定岩石,因而具有一定局限性。

建立考虑干湿循环作用的蠕变本构模型有助于描述不同岩体工况下岩石经历干湿循环后的蠕变特性。Li等[7]通过单轴压缩蠕变试验,得到了砂岩长期强度随干湿循环作用次数的关系,并建立了相关蠕变模型。Wang等[8-10]基于水库消落带岩石的蠕变试验结果,提出一种非定常非线性的砂岩蠕变本构模型。王新刚[11]提出了考虑岩石饱水-失水循环次数的非线性黏弹塑性蠕变模型。岩石的蠕变曲线通常可分为衰减蠕变、稳态蠕变和加速蠕变3个阶段[12],明晰加速蠕变阶段岩石的变形特征对于滑坡临滑预报具有十分重要的意义。目前,通常采用非线性蠕变元件代替常规的线性蠕变元件,建立能够描述岩石加速蠕变阶段的蠕变本构模型[13]。蠕变过程是岩石内部损伤累积直至破坏的过程,当满足一定条件时岩石才会进入加速蠕变阶段[14]。现有本构模型在描述岩石由稳态蠕变进入加速蠕变过程时并没有设置判别条件。对此,本文将考虑判别条件的非线性阈值元件与西原模型串联,建立了可以描述干湿循环条件下包含加速蠕变在内3个阶段的岩石蠕变本构模型,并将该模型嵌入FLAC3D分析程序中,以促进新岩石蠕变模型在工程实践中的应用,为保持水电工程边坡稳定及保障库坝运行安全提供理论依据与数值模拟分析方法。

1 考虑干湿循环作用的岩石蠕变模型

现有岩石蠕变本构模型中,西原模型因能较好地描述岩石线性蠕变特性而被广泛应用,但其无法描述岩石加速蠕变特性[15]。岩石进入加速蠕变阶段后岩石应力与应变通常不满足一一对应的关系,且岩石变形是不可逆的,因此本文将应变控制的非线性阈值元件与西原模型串联,并充分考虑干湿循环次数对岩石力学参数的影响,建立可以描述干湿循环条件下包含加速蠕变在内的3个阶段的岩石蠕变本构模型,如图1所示。图中带有(n)的物理量代表第n次干湿循环后岩石的相关力学参数。本文通过改变干湿循环后岩石物理力学参数来考虑干湿循环作用,即岩石物理力学参数是n的函数。

图1 考虑干湿循环作用的岩石非线性蠕变模型Fig.1 Nonlinear rock creep model considering dry-wet cycle

引入非线性阈值元件后,当岩石蠕变位移小于元件阈值时,图1所示岩石非线性蠕变模型退化为传统的西原模型;反之,岩石进入加速蠕变阶段。非线性阈值元件满足如下关系:

(1)

非线性阈值元件的黏滞系数与岩石蠕变损伤满足如下关系:

ηnl(t)=ηnl,0(n)(1-D)

(2)

式中:ηnl,0为加速蠕变初始时刻非线性阈值元件的黏滞系数;D为加速蠕变过程岩石的损伤变量,0≤D<1。

岩石加速蠕变阶段损伤变量随时间的变化规律如式(3)所示:

D=1-e-at

(3)

式中:a为与岩石材料性质相关的常数,受干湿循环次数的影响。

联立式(1)～(3),可得到加速蠕变过程中非线性元件的蠕变位移:

(4)

基于非线性阈值元件的应力应变关系,结合西原模型的蠕变方程,得到干湿循环条件下岩石非线性蠕变位移方程为

(5)

式中:Sij(t)为t时刻蠕变位移;带下标0,1,2及nl的G(n)、η(n)分别代表Hooke体、Kelvin体、Bingham体及非线性体在干湿循环次数为n时的弹性模量和黏滞系数;K为体积模量;σmδij为球应力张量,其中δij为Kronecker符号;sij、eij分别为应力偏张量和应变偏张量;F为岩石屈服函数;Q为塑性势函数。

2 岩石蠕变模型验证

基于本文建立的非线性岩石蠕变本构模型及文献[5]、文献[8]中的试验结果,采用MATLAB中1stOpt进行参数辨识,得到不同干湿循环次数下岩石蠕变全过程曲线,如图2～3所示。对比分析可知,岩石蠕变过程明显分为3个阶段:即减速蠕变阶段、稳态蠕变阶段及加速蠕变阶段,且稳态蠕变阶段与加速蠕变阶段之间存在明显的分界点,这从侧面说明了本文引入非线性阈值元件是与实际吻合的。新本构模型描述的岩石蠕变曲线与试验数据吻合较好,验证了本文提出模型的有效性和可靠性。

图2 砂岩试验结果与模型拟合结果Fig.2 Experimental results and model fitting results of sandstone

图3 粉砂岩实测结果与模型拟合结果Fig.3 Experimental results and model fitting results of siltstone

3 蠕变本构模型计算程序开发

3.1 模型有限差分形式

为了促进新岩石蠕变本构模型在工程实践中的应用,对蠕变本构模型满足的应力应变关系进行推导,并将其改写成有限差分形式,以便对边坡蠕变进行数值模拟分析。本文模型分为4个部分,第1部分为Hooke体,第2部分为Kelvin模型,第3部分为Bingham模型,第4部分为非线性体。根据各部分模型的串联关系可知,各部分模型的应力相等,总应变为各部分的应变之和,即有:

Sij=SijE=SijK=SijB=Sijnl

(6)

式中:Sij为总应力;带上标E,K,B及nl的Sij分别代表Hooke体,Kelvin体,Bingham体及非线性体的应力。

εij=εijE+εijK+εijB+εijnl

(7)

式中:εij为总应变;带上标E,K,B及NL的εij分别代表Hooke体,Kelvin体,Bingham体与非线性体的应变。

三维应力状态下,Hooke体的应变增量可表示为

(8)

Kelvin体的本构关系可表示为

(9)

Bingham体的应变关系中包含了塑性本构关系,其表示为

(10)

根据FLAC3D的计算方法,在循环计算过程中,各部分的增量关系为

(11)

(12)

式中:上标N、O分别表示一个时间步长的两次迭代过程中的新老变量值。

(13)

Kelvin体的总应变值为

(14)

其中:

(15)

结合上式,得到非线性黏弹塑性蠕变模型的新的偏应力为

(16)

其中:

(17)

在非线性阈值元件触发条件下,非线性阈值的蠕变位移为

(18)

根据式(18),得到非线性元件的位移速率为

(19)

因此,非线性元件的应变增量为

(20)

结合各分部模型的位移和应变计算公式,可实现岩石非线性蠕变本构模型的应力更新,并计算得到不同蠕变阶段应变增量。

3.2 模型验证

将提出的考虑干湿循环条件的蠕变本构模型的有限差分形式嵌入FLAC3D分析软件中,以模拟分析干湿循环作用下岩石的蠕变特性。为了验证新模型二次开发的正确性,在FLAC3D中建立了一个单元体模型,如图4所示,对单元体的上表面施加1 MPa的恒定应力,分别对不同阶段蠕变变形进行验证。为了验证自定义蠕变模型在黏-弹阶段的正确性,可将单元材料的黏聚力与抗拉强度设置为无穷大,以保证单元体材料不会发生塑性屈服。开启蠕变计算模块后得到未进入屈服状态时单元体模型的蠕变过程曲线,试验中将单元体顶部节点设置为位移监测点,对比分析监测得到的轴向位移与对应的蠕变模型理论位移,验证新蠕变本构模型的有限差分形式在黏弹性蠕变阶段的正确性。单元体蠕变参数设置见表1。

表1 单元体蠕变参数Tab.1 Monomer creep parameters of element model

图4 单元体模型Fig.4 Diagram of element model

单元体荷载的矢量图如图4所示,本次数值试验中下底面设置为固定边界条件,在应力的作用下单元体将产生侧向位移与轴向位移。将监测点的轴向位移过程曲线与理论计算得到的曲线进行对比,如图5所示。发现数值模拟试验中得到的蠕变位移时程曲线与理论计算得到的结果十分吻合,验证了新本构模型的有限差分形式在黏弹性蠕变阶段的有效性。

图5 黏弹性蠕变阶段验证Fig.5 Visco-elasticity creep stage verification

为了验证新本构模型的有限差分形式在黏塑性蠕变阶段的有效性,调整单元体黏聚力的大小,保证单元体在相同应力的作用下能进入屈服阶段,得到的屈服状态下监测点的蠕变位移时程曲线如图6所示。由图可知,当黏聚力为1 MPa时,单元进入剪切屈服状态,单元体监测点处蠕变位移时程曲线也表现出黏塑性蠕变行为。将监测点的轴向位移时程曲线与理论计算得到的曲线进行对比分析,发现数值模拟中得到的蠕变位移时程曲线与理论计算得到曲线的十分吻合,验证了新蠕变本构模型的有限差分形式在黏塑性蠕变阶段的有效性。

图6 黏塑性蠕变阶段验证Fig.6 Visco-plasticity creep stage verification

为了验证新本构模型的有限差分形式在描述岩石非线性蠕变性质方面的有效性,在数值试验中设置单元体的蠕变位移阈值。通过主程序判断主应变值是否大于元件应变阈值,若是,则非线性元件触发,单元体的应力与应变呈现非线性关系,单元体进入加速蠕变阶段。试验中逐级设置轴向应力得到蠕变过程如图7所示。由图7可知,随着荷载的递增,蠕变位移逐渐增大,当蠕变位移大于应变阈值时,岩石进入加速蠕变阶段。

图7 单元体分级加载蠕变变形验证Fig.7 Creep deformation verification of element under graded loading

综合上述结果可知,基于FLAC3D软件二次开发得到的非线性蠕变模型可以准确可靠地描述包含加速蠕变在内的蠕变全过程。

3.3 案例计算

锦屏一级水电站位于四川省凉山彝族自治州盐源县与木里县交界的洼里乡灯盏窝,是雅砻江水能资源最富集的中、下游河段5级水电开发中的第一级。本文以锦屏一级水电站左岸边坡工程为例开展蠕变数值模拟研究工作。首先建立左岸边坡有限元模型,如图8所示。边坡开挖施工时间段为2009年8月31日至2013年7月31日,以此时间段作为蠕变时长,利用嵌入FLAC3D的非线性蠕变本构模型,计算边坡横河向位移,如图9所示。由图9可知,边坡Ⅱ-Ⅱ断面的横河向最大位移达到了29 mm。根据监测资料,1 885 m坝顶平台以下边坡监测点TPL64的横河向位移为21.8 mm,与数值模拟对应位置处位移相差较小,如图10(a)所示;对比分析图10(b)中监测点TPL48处的现场监测结果与相应的数值模拟结果,发现二者随时间变化的趋势和幅度基本吻合,这验证了该工程案例数值模拟计算结果的可靠性。

图8 左岸边坡计算模型Fig.8 Calculation model of left bank slope

图9 边坡位移云图(单位:m)Fig.9 Displacement nephogram of slope

图10 锦屏Ⅰ级水电站左岸边坡蠕变曲线Fig.10 Creep curves of left bank slope at Jinping Ⅰ Hydropower Station

运行期水位升降一次意味着边坡岩石经历一次干湿循环。岩石经历干湿循环过程中,后一次干湿循环岩石蠕变参数降低幅度为前一次的40%。基于此,利用前述二次开发得到的新蠕变本构模型计算得到不同干湿循环次数下边坡岩石蠕变结果,如图11所示。由图可知,随着干湿循环次数的增加,边坡蠕变变形的程度逐渐增大,且由图容易得出蠕变变形显著增大区域的位置和范围,这为后续采取合理工程措施以保持边坡稳定提供了重要参考和依据。

图11 干湿循环条件下边坡蠕变位移(单位:m)Fig.11 Creep displacement of slope under dry-wet cycles

4 结 论

本文将非线性阈值元件与西原模型串联,建立了考虑干湿循环作用的岩石蠕变本构模型,获得的主要结论如下:

(1) 考虑加速蠕变判别条件的新蠕变本构模型能准确描述干湿循环条件下岩石蠕变的全过程曲线。

(2) 成功将新蠕变本构模型嵌入FLAC3D分析软件,模拟结果与室内试验结果及边坡工程监测结果吻合较好,验证了程序的有效性。

(3) 数值模拟结果表明随着干湿循环次数的增加,边坡蠕变变形的程度逐渐增大,且可直观给出蠕变变形显著增大区域的位置和范围。

本文研究基于各向同性的假定,没有考虑实际工程中岩体存在的软弱结构面对蠕变的影响,对于实际工程岩体的各向异性特征没有体现。总体而言,本文建立的蠕变本构模型对工程研究具有一定的参考价值和进步意义,但还需进一步修正和完善。