横观各向同性岩体三维蠕变特性的表征及其表征模型的构建的制作方法

本发明涉及岩石力学与工程技术领域，具体涉及一种横观各向同性岩体三维蠕变特性的表征及其表征模型的构建。

背景技术：

层状岩体在地壳中广泛存在，其在力学上具有横观各向同性的特点，岩体的破坏和失稳与其内部的损伤演化有着密切的关系。

随着地下工程的开发建设，层状岩体的力学性质（包括蠕变特性）对隧道、地下厂房、水电等工程围岩变形和应力分布，尤其是长期稳定性有重要的影响。由于层状岩体特殊的构造性质，其力学性质与各向同性岩石相比具有较大的差异，如随着层理倾角的不同，其蠕变特性具有显著的各向异性。如溪洛渡工程坝肩两岸分别存在n200w和近ew向两组基体结构面，这两组结构面在岩体中大量发育，构成岩体的基本层次，两岸岩体沿结构面发生滑移错动挤压坝体，对坝体安全性构成严重威胁；锦屏一级水电站坝址区地层岩性、层厚变化大，岩体组合复杂，其中大理岩夹杂工程地质特性较差的绿泥石组成的层状岩体，主要分布在两岸坝肩、水垫塘及地下厂房等部位，影响到这些工程部位的稳定性，尤其对坝区右岸而言，大理岩中呈透镜状、团块状分布的绿片岩，其片理发育且倾向坡外，形成了层间弱面，成为影响右岸抗滑稳定的主要因素；内蒙古平庄西露天矿某处边坡，其下部的缓倾角砂岩—泥岩软弱夹层在降水的长期浸泡下，岩体软化，沿软弱面发生滑移，造成了滑坡事故。由此可以看出，岩体蠕变特性对工程建筑物的安全有着重要影响。从理论上分析评价这些建筑物的安全性是非常必要的，而构建科学、合理的蠕变本构方程是进行理论分析的基础。这其中就需要选用合理的蠕变参数进行分析预测，并利用分析结果指导工程设计施工，达到确保工程结构的安全性的目的。

现有的针对岩石材料的各向异性蠕变本构模型多是基于塑性理论，其更适用于金属等发生塑性破坏的材料。众所周知，岩体的破坏展示出准脆性性质，而非完全塑性。而现有模型研究成果多根据试验现象，通过构建一维蠕变本构模型，将不同方向蠕变参数简单辨识为不同的相互独立的值，而并未建立真正意义上的横观各向同性三维蠕变本构模型。

此外，由于层理弱面的存在，现有的各向同性蠕变模型难以全面反映横观各向同性岩体的蠕变力学特性。

因此，亟需发展适用于可描述岩体准脆性断裂破坏特性的各向异性黏弹塑性元件模型，以期能够更为真实地表征或反映岩体断裂破坏本质，进而为各项与岩石相关的工程建设提供更精准可靠的决策或设计依据。

技术实现要素：

本发明要解决的技术问题是提供一种横观各向同性岩体三维蠕变特性的表征方法及其表征模型的构建，以期解决现有技术中的各向异性蠕变本构模型难以全面表征反映横观各向同性岩体的蠕变力学特性的技术问题。

为解决上述技术问题，本发明的技术思路如下（参见图1）：

横观各向同性岩体由于层理弱面的存在，表现出各向异性的特点。因此，不同于各向同性岩石，在确定横观各向同性岩体五个弹性常数时，根据弹性本构方程的特点，至少需要根据平行组、垂直组两个试件的单轴或三轴压缩实验结果来计算弹性常数。蠕变参数也是如此，依据横观各向同性岩体三维蠕变本构方程的特点，根据平行组、垂直组两个试件的蠕变试验结果，首次分析得出横观各向同性岩体的蠕变参数。

具体采用的技术方案如下：

设计一种表征横观各向同性岩体三维蠕变特性的方法，即横观各向同性岩体三维蠕变模型，三维空间（三维坐标系）中，设xoy为岩体层理面所在平面，对于层理水平岩样，，采用如下蠕变损伤本构方程进行表征：

——式（vii）；

对于层理垂直岩样，采用如下蠕变损伤本构方程进行表征：

——式（viii）；

以上各式中，ε90轴为x方向的正应变；ε0轴为z方向的正应变；σ1为轴向大主应力；σ3为环向小主应力；μhh、μhv分别为平行和垂直层理方向的泊松比；εhm、ηhm、εhk、ηhk，及对应的n1、n2、n3、n4分别为垂直层理面方向的独立拉压蠕变参数；t为时间。

所述横观各向同性岩体三维蠕变模型的构建方法，包括以下步骤：

（1）将岩石试样制成层理水平岩样和层理垂直岩样，进行横观各向同性岩体的压缩蠕变试验；

（2）采用陈氏加载法进行处理步骤（1）的试验结果，得层理水平岩样和层理垂直岩样的分级蠕变试验曲线；

（3）三维空间中，设xoy为层理面所在平面，采用以下层状岩体的三维蠕变本构方程，进行相应应变计算：

——式（i）；

——式（ii）；

——式（iii）；

——式（iv）；

——式（v）；

——式（vi）；

式（i）至式（vi）中：εx为x方向的正应变；εy为y方向的正应变；εz为z方向的正应变；γyz为yoz平面的剪应变；γxz为xoz平面的剪应变；γxy为xoy平面的剪应变；μhh、μhv分别为平行和垂直层理方向的泊松比；εhm、ηhm、εhk、ηhk，及对应的n1、n2、n3、n4分别为垂直层理面方向的独立拉压蠕变参数；t为时间σx为方向的正应力，σy为方向的正应力，σz为方向的正应力；

（4）对于层理水平岩样，xoy平面为各向同性平面，则σx=σy=σ3，σz=σ1，εz=ε0轴，其应变采用式（iii）计算，式（iii）变换写为式（vii）：

——式（vii）；

对于层理垂直岩样，σx=σy=σ3，σx=σ1，εx=ε90轴，其应变采用式（i）计算，式（i）变换写为式（viii）：

——式（viii）；

式中，σ1为轴向大主应力；σ3为环向小主应力；

即得。

优选的，在步骤（1）中，所述压缩蠕变试验采用单轴或三轴压缩。

与现有技术相比，本发明的有益技术效果在于：

1.本发明首次建立了真正意义上的可反映横观各向同性岩体特性的三维蠕变本构模型。

2.本发明建立模型综合考虑了平行和垂直层理面方向蠕变力学行为的差异性，所提的横观各向同性岩体的三维蠕变本构模型，由八个独立的拉压蠕变参数所组成。

3.本发明建立的模型可实现采用一套蠕变参数同时表征出横观各向同性岩体不同层理方向的蠕变特征的目的。

4.采用本发明建立的模型，可采用传统的蠕变试验设计方案（即根据不同层理倾角千枚岩荷载峰值，分级加载）进行不同应力水平的蠕变试验；也可采用采用相同的应力水平进行蠕变试验，应用范围更为广泛。

5.本发明建立的模型通过选用合理的蠕变参数进行分析，并利用分析结果指导工程设计施工，可确保工程结构的安全性（如对四川杜家山千枚岩隧道进行稳定性分析），还可对工程的长期安全性进行评价，对岩石力学理论以及地下工程建设具有重大的指导意义。

附图说明

图1为横观各向同性岩体三维模型蠕变参数识别流程图；

图2为蠕变模型示意图；

图3为坐标系示意图；

图4为千枚岩层理水平岩样的轴向蠕变分级曲线图；

图5为千枚岩层理垂直岩样的轴向蠕变分级曲线图；

图6为千枚岩蠕变参数回归分析图；

图7为千枚岩层理水平岩样的蠕变试验及发展模型的拟合曲线图；

图8为千枚岩层理垂直岩样的蠕变试验及发展模型的拟合曲线图；

图9为云母片岩层理水平岩样的轴向蠕变分级曲线图；

图10为云母片岩层理垂直岩样的轴向蠕变分级曲线图；

图11为云母片岩蠕变参数回归分析图；

图12为云母片岩层理水平岩样的蠕变试验及发展模型的拟合曲线图；

图13为云母片岩层理垂直岩样的蠕变试验及发展模型的拟合曲线图。

具体实施方式

下面结合附图和实施例来说明本发明的具体实施方式，但以下实施例只是用来详细说明本发明，并不以任何方式限制本发明的范围。

在以下实施例中所涉及的材料如无特别说明，均为市售常规材料；所涉及的步骤或试验方法，如无特别说明，均为常规方法；所涉及的名称参数，如无特别说明，均为本领域常规术语。

实施例一：层状岩体的三维蠕变本构方程的构建

1.层状岩体的蠕变特性具有以下特征：

（1）加载初期具有瞬时弹性变形，模型中应含有弹性元件；

（2）随着时间的增加，蠕变速率逐渐降低，并逐渐近似趋近于某个常数，故模型中应具有粘性元件；

（3）在低应力水平下，蠕变曲线达到某个极限值后，其应变值保持不变，因此模型中应包含弹性元件与粘性元件的并联；

可得描述层岩岩体的蠕变元件模型如图2所示。

2.对图2元件模型进行力学分析。

（1）maxwell体中弹性元件与粘性元件串联，具有以下关系：

——式（a）；

——式（b）；

式中，d为算子，d=d/dt。

（2）kelvin体中弹性元件与粘性元件并联，具有以下关系：

——式（c）；

——式（d）；

式（a）～式（d）中，em，ek分别为maxwell体和kelvin体的弹性参数；ηm，ηk分别为maxwell体和kelvin体的粘性参数。

图2中的一维蠕变模型由maxwell体和kelvin体串联而成，可得：

——式（e）；

将式（b）、（d）代入式（e），得：

——式（f）；

计蠕变柔量：

，则式（f）可写为：

——式（g）；

通过对式（f）中的算子d进行积分，可得到最终的一维蠕变本构方程：

——式（h）；

3.各向同性材料蠕变本构方程

根据一维条件下建立的蠕变本构方程，通过常泊松比方法将其推广到三维应力状态。

思路为：物体蠕变过程中假设泊松比不随时间和应力改变，并等于弹性阶段的泊松比，μ(σ,t)=μ为不变的实常数。基于此假定将蠕变方程从一维状态推广到三维状态。

各向同性材料在三维应力空间的应力应变关系为：

——式（i）；

式中：{ε^e}为弹性应变矩阵，{σ^e}为三维应力矩阵，{h}为柔度矩阵，分别可表示为：

——式（j）；

矩阵（j）中e、μ分别为各向同性材料的弹性模量和泊松比。

式（i）也可写为：

——式（k）；

式中：

。

保持泊松比矩阵[u]不发生变化，采用微分算子替换的方法，将蠕变模型中的柔度模量j(t)替换线弹性本构，即式（k）中的柔度系数1/e，得到蠕变本构方程的基本式：

——式（l）；

根据前文的推导，代入式（l），得

——式（m）；

将式（m）采用张量形式表达，可得：

——式（n）；

式中：i=x、y、z，j=x、y、z，σm=1/3(σx+σy+σz)，δij为kronecker符号。

4.横观各向同性岩体蠕变本构方程

在横观各向同性条件下，假设横观各向同性岩体的自然坐标系为oxyz，横观各向同性面所在平面为xoy平面，z轴为xoy平面法向方向，x、y、z服从右手法则，如图3所示。

基于图3，可得三维应力状态下横观各向同性岩体的弹性应力应变本构关系为：

——式（o）；

式中：

其中，eh、ev、μhh、μhv、g’为5个弹性常数参数，eh、ev分别为平行和垂直层理方向的弹性模量，μhh、μhv分别为平行和垂直层理方向的泊松比，g’为垂直层理方向的剪切模量，{ε}={εx、εy、εz、γyz、γxz、γxy}^t，{σ}={σx、σy、σz、τyz、τxz、τxy}^t。

横观各向同性岩体具有五个独立的弹性参数。saint—venant提出了近似弹性参数表达式（p），其为四个弹性参数。

——式（p）；

假设ev=neh，并结合式（p），则式（o）中的柔度矩阵[s]可写为：

——式（q）；

式中：

针对横观各向同性岩体的特点，基于常泊松比思想，假定μhh和μhv不随应力和时间改变，均等于弹性阶段的常数。采用与各向同性材料推广到三维相同的思想，保持泊松比矩阵不发生变化，将式（l）中的各向同性材料的泊松比矩阵[u]采用式（q）所示的横观各向同性材料泊松比矩阵[utra]代替，则可得层状岩体的三维蠕变本构方程如下：

——式（r）。

实施例二：横观各向同性岩体三维蠕变模型的应用

基于上述模型，可实现采用一套蠕变参数同时表现出横观各向同性岩体不同层理方向的蠕变特征的目的，具体包括以下步骤：

（1）基于横观各向同性岩体的层理面与水平面的不同夹角，岩石试样分别加工制作成层理水平岩样和层理垂直岩样两类，进行横观各向同性岩体的三轴压缩蠕变试验；

（2）将层理水平岩样与层理垂直岩样的试验结果，采用陈氏加载法（参考文献：陈宗基，康文法．岩石的封闭应力、蠕变和扩容及本构方程）进行处理，进一步得到相应的层理水平岩样与层理垂直岩样的分级蠕变试验曲线；

（3）三维空间中，设xoy为层理面所在平面，将实施例一中的式（r）中蠕变柔量进行积分计算，并考虑平行和垂直层理方面蠕变特性的差异，横观各向同性材料三维蠕变模型由垂直层理方向的εhm、ηhm、εhk、ηhk，及对应的n1、n2、n3、n4共8个独立拉压蠕变参数，将式（r）展开，并采用类比的思想，可得横观各向同性材料的蠕变方程表达式，采用如下层状岩体的三维蠕变本构方程，进行相应应变计算：

——式（i）；

——式（ii）；

——式（iii）；

——式（iv）；

——式（v）；

——式（vi）；

式（i）至式（vi）中：εx为x方向的正应变；εy为y方向的正应变；εz为z方向的正应变；γyz为yoz平面的剪应变；γxz为xoz平面的剪应变；γxy为xoy平面的剪应变；μhh、μhv分别为平行和垂直层理方向的泊松比；εhm、ηhm、εhk、ηhk，及对应的n1、n2、n3、n4分别为垂直层理面方向的八个相互独立的拉压蠕变参数；t为时间。

（4）对于层理水平岩样，xoy平面为各向同性平面，则σx=σy=σ3，σz=σ1，εz=ε0轴，其应变采用式（iii）计算，式（iii）通过变换可写为式（vii）：

——式（vii）；

式中，σ1为轴向大主应力；σ3为环向小主应力。

对于层理垂直岩样，σx=σy=σ3，σx=σ1，εx=ε90轴，其应变采用式（i）计算，式（i）通过变换可写为式（viii）：

——式（viii）；

（5）基于式（vii），结合层理水平岩样的蠕变试验曲线，通过levenberg-marquardt非线性最小二乘法进行拟合可分别得到参数组合εhmn1、ηhmn2、εhkn3和ηhkn4；

（6）若层理水平岩样和垂直岩样采用相同的应力水平进行蠕变试验，即层理水平岩样应力水平σh等于垂直岩样的应力水平σv，则将第（5）步拟合得到的εhmn1、ηhmn2、εhkn3和ηhkn4结果代入式（viii）中；再根据层理垂直岩样的蠕变试验曲线结果，拟合得到εhm、ηhm、εhk和ηhk；

（7）最后将εhmn1、ηhmn2、εhkn3、ηhkn4和对应εhm、ηhm、εhk、ηhk的计算结果，将二者相除，从而得出n1、n2、n3和n4。最终得到横观各向同性岩体的八个独立的拉压蠕变参数εhm、ηhm、εhk、ηhk、n1、n2、n3、n4。

试验例1：

取汶川县桃坪隧道进洞口的千枚岩作为试验材料，该岩石材料层状构造显著，为一种典型的横观各向同性岩体。采用实施例的方法验证所构建的横观各向同性岩体三维蠕变模型。

具体步骤如下：

（1）基于层理面与水平面的不同夹角，岩石试样被别加工制作成层理水平岩样和层理垂直岩样两类。

（2）进行三轴压缩蠕变试验：

采用5mpa的恒定围压，分级加载轴向压力的加载制度。其中，层理水平岩样轴向压力依次为9.4mpa、15.1mpa、20.8mpa、26.5mpa、32.2mpa；层理垂直岩样轴向压力依次为26.0mpa、35.6mpa、45.6mpa、55.5mpa、65.4mpa、75.1mpa、79.3mpa。

（3）将层理水平岩样与层理垂直岩样的试验结果采用陈氏加载法进行处理，室内试验为每间隔一段时间测量一个应变值，将不同的应变值连起来，进一步可得层理水平岩样与层理垂直岩样的分级蠕变试验曲线，如图4和图5所示。

由图4和图5可以看出，平行组和垂直组千枚岩在加载瞬间都会产生一定量的瞬时应变，此后在轴向应力作用下，产生随时间增长的蠕变变形。应力水平较低时，蠕变速率短时间内趋于零值；随着应力水平的提高，蠕变速率随时间的增长而衰减，减速蠕变时间增加，最后趋于一个恒定的非零值，如果蠕变时间足够长，蠕变变形达到一定值后最终使岩样发生破坏。

由此可见，千枚岩具有三轴应力作用下产生瞬时应变、衰减蠕变和稳定蠕变的典型蠕变行为特性。

（4）采用实施例中的蠕变参数识别方法与流程，基于步骤（3）所得层理水平岩样的蠕变试验结果，采用式（vii）进行拟合分析，可得出层理水平岩样在9.4mpa、15.1mpa、20.8mpa、26.5mpa、32.2mpa等不同应力水平下参数组合εhmn1、ηhmn2、εhkn3和ηhkn4的具体数据，确定结果如表1所示：

表1层理水平岩样蠕变参数

。

（5）根据层理水平岩样的试验结果，通过回归分析似计算得到层理垂直岩样应力水平对应的εhmn1、ηhmn2、εhkn3和ηhkn4。具体回归分析时，为避免出现过大值、过小值甚至负值而使结果失真，回归曲线一般情况下不选用线性和多项式方程，而采用指数方程、幂函数方程、对数方程等非线性方程。层理水平岩样的蠕变参数回归如图6所示。

由图6可以看出，εhmn1随着应力水平的提高有线性增加的趋势；ηhmn2随着应力水平的提高呈指数函数下降；εhkn3和ηhkn4随着应力水平呈非线性先升高后下降的趋势。

根据回归方程计算得到层理垂直岩样对应应力水平下的εhmn1、ηhmn2、εhkn3和ηhkn4的具体结果如表2所示。

表2参数一览表

。

（6）将表2中的分析结果代入式（viii），并利用式（viii）对层理垂直岩样的蠕变曲线进行拟合得到的蠕变参数如表3所示。

表3参数一览表

。

（7）根据表2和表3中的参数预测结果，将二者相除，可计算得到n1、n2、n3和n4。最终可得到横观各向同性岩体八个独立的蠕变参数，如表4所示。

表4层理垂直岩样蠕变参数

。

（8）将千枚岩蠕变试验所得值（注：图4、图5为试验值）与基于实施例的模型计算所得值拟合（注：拟合过程为采用与室内试验相同的时间t，即图中的横坐标值，代入蠕变方程式（vii）和式（viii）中，方程中的蠕变参数如表4所示，从而对应每个时间t，计算出一个应变值，然后与试验值进行对比），结果如图7和图8所示。

由图7和图8可以看出，实施例二的三维蠕变损伤本构方程得出的计算值与试验值吻合良好，该模型可较好地反映了千枚岩的蠕变特性。证明了本发明横观各向同性岩体三维蠕变模型的表征有效、合理，能够反映出其真实的蠕变特性。

试验例2：

取国家重点工程丹巴水电站的云母片岩作为试验材料，该岩石材料层状构造显著，为一种典型的横观各向同性软岩体，对其蠕变损伤的分析直接影响丹巴水电站安全运营的理论分析结果。采用实施例的方法验证所构建的横观各向同性岩体三维蠕变模型。

具体步骤如下：

（1）基于层理面与水平面的不同夹角，岩石试样被别加工制作成层理水平岩样和层理垂直岩样两类。

（2）进行单轴压缩蠕变试验：

采用0mpa的围压（即无围压），分级加载轴向压力的加载制度。其中，层理水平岩样轴向压力依次为6.0mpa、10.0mpa、16.0mpa；层理垂直岩样轴向压力依次为3.0mpa、5.0mpa、7.0mpa、9.0mpa。

（3）将层理水平岩样与层理垂直岩样的试验结果采用陈氏加载法进行处理，室内试验为每间隔一段时间测量一个应变值，将不同的应变值连起来，进一步可得层理水平岩样与层理垂直岩样的分级蠕变试验曲线，如图9和图10所示。

由图9和图10可以看出，平行组和垂直组云母片岩在加载瞬间都会产生一定量的瞬时应变，此后在轴向应力作用下，产生随时间增长的蠕变变形。应力水平较低时，蠕变速率短时间内趋于零值；随着应力水平的提高，蠕变速率随时间的增长而衰减，减速蠕变时间增加，最后趋于一个恒定的非零值，如果蠕变时间足够长，蠕变变形达到一定值后最终使岩样发生破坏，应力越高，稳态蠕变速率越高。

由此可见，云母片岩在单轴应力作用下会产生瞬时应变、衰减蠕变和稳定蠕变的典型蠕变行为特性。

（4）采用实施例中的蠕变参数识别方法与流程，基于步骤（3）所得层理水平岩样的蠕变试验结果，采用式（vii）进行拟合分析，可得出层理水平岩样在6.0mpa、10.0mpa、16.0mpa等不同应力水平下参数组合εhmn1、ηhmn2、εhkn3和ηhkn4的具体数据，确定结果如表5所示：

表5层理水平岩样蠕变参数

。

（5）根据层理水平岩样的试验结果，通过回归分析似计算得到层理垂直岩样应力水平对应的εhmn1、ηhmn2、εhkn3和ηhkn4。具体回归分析时，为避免出现过大值、过小值甚至负值使结果失真，回归曲线一般情况下不选用线性和多项式方程，而采用指数方程、幂函数方程、对数方程等非线性方程。层理水平岩样的蠕变参数回归如图11所示。

由图11可以看出，εhmn1随着应力水平的提高有线性增加的趋势；ηhmn2随着应力水平的提高呈指数函数先下降后提高的趋势；εhkn3和ηhkn4随着应力水平呈指数函数先升高后下降的趋势。

根据回归方程计算得到层理垂直岩样对应应力水平下的εhmn1、ηhmn2、εhkn3和ηhkn4的具体结果如表6所示。

表6参数一览表

。

（6）将表6中的分析结果代入式（viii），并利用式（viii）对层理垂直岩样的蠕变曲线进行拟合得到的蠕变参数如表7所示。

表7参数一览表

。

（7）根据表6和表7中的参数预测结果，将二者相除，可计算得到n1、n2、n3和n4。

最终可得到横观各向同性岩体八个独立的蠕变参数，如表8所示。

表8层理垂直岩样蠕变参数

。

（8）将云母片岩蠕变试验所得值（注：图9、图10为试验值）与基于实施例的模型计算所得值拟合（注：拟合过程为采用与室内试验相同的时间t，即图中的横坐标值，代入蠕变方程式（vii）和式（viii）中，方程中的蠕变参数如表8所示，从而对应每个时间t，计算出一个应变值，然后与试验值进行对比），结果如图12和图13所示。

由图12和图13可以看出，实施例二的三维蠕变损伤本构方程得出的计算值与试验值吻合良好，该模型可较好地反映了云母片岩的蠕变特性。证明了本发明所构建的横观各向同性岩体三维蠕变模型的的表征有效、合理，能够反映出其真实的蠕变特性。

上面结合附图和实施例对本发明作了详细的说明，但是，所属技术领域的技术人员能够理解，在不脱离本发明宗旨的前提下，还可以对上述实施例中的各个具体参数进行变更，形成多个具体的实施例，均为本发明的常见变化范围，在此不再一一详述。