基于分数阶导数模型的巷道围岩变形失稳规律研究

肖书文

（山西泽州天泰锦辰煤业有限公司，山西 晋城 048000）

深部巷道围岩长期处于高地应力的环境中会发生蠕变变形[1-3]，对于蠕变本构模型的建模研究主要包括经验模型与组合元件模型，其中组合元件模型由于物理意义清晰而得到了广泛应用，传统的元件模型包括Kelvin 模型、Burgers 模型、Nishihara 模型等，它们仅能描述前两个蠕变变形阶段，具有较多局限性[3-5]。本文通过分数阶导数的建模方法建立组合元件模型，给出了本构模型的解析解。经数据验证，基于分数阶导数的本构模型可以更好地反映蠕变变形特征，尤其是蠕变的第三阶段。

1 基于分数阶导数的元件模型

1.1 分数阶导数定义

分数阶微积分相比于整数阶微积分的区别在于积分阶次由整数变为了分数，分数阶微积分具有多种数学定义方式，其中Riemann-Liouville 分数阶微积分定义方法的应用最为广泛。设f在(0,+∞)上连续分布，并且在[0,+∞] 的任何有限子区间可积，当t＞0，Re(γ)＞0时，称式（1）为函数f(t)的γ阶Riemann-Liouville分数阶导数：

基于以上分数阶微积分的定义，将传统的蠕变元件模型进行改进，得到分数阶导数的蠕变元件模型。

1.2 分数阶导数蠕变元件模型

在元件模型中，常运用Newton 黏壶表示黏性变形结构，其本构方程如式（2）所示：

式中，σ为应力，MPa；η为黏性系数；ε为应变，%；t为时间，h；γ为导数的阶次。当γ=1 时，此方程为传统的Newton 黏壶，用于表示黏性变形状态。当γ=0 时，此方程为弹簧元件，用于表示理想弹性变形状态。根据物理意义可知，当0＜γ＜1时，此方程表示分数阶黏壶（Abel黏壶），用于表示黏弹性变形状态。

在蠕变过程中，应力σ为常数，根据Riemann- Liouville 分数阶微积分算子理论，将式（2）进行变换得到以时间t为自变量的应变方程如式（3）所示：

在岩石的蠕变过程中通常伴随着损伤过程，随着损伤的加剧岩石的变形程度也逐渐增加，最终衍生出宏观裂隙使岩石发生破坏，在工程中常表现为巷道的变形与失稳的破坏过程。因此，考虑将损伤变量引入到黏性系数中，用于表示黏性系数的劣化过程，将损伤变量定义为负指数形式如式（4）所示：

式中，D为损伤变量；α 为损伤劣化系数，用于表示损伤变化程度。

联立式（2）、式（3）与式（4），保持应力σ为常数，根据Riemann-Liouville 分数阶微积分算子理论，得到以时间t为自变量、包含黏性损伤劣化过程的应变方程如式（5）所示：

2 分数阶导数蠕变本构模型

2.1 蠕变本构模型

蠕变过程通常分为初始蠕变阶段、稳定蠕变阶段与快速蠕变阶段，组合原件模型的实际意义是通过不同元件模型的组合搭配可以建立本构方程，以期对上述阶段进行描述。其中，Nishihara 模型被广泛用于描述初始蠕变阶段与稳定蠕变阶段，如图1 所示。模型中包括描述弹性应变εe的胡克体、描述黏弹性应变εve的黏弹性体、描述黏塑性应变εvp的黏塑性体。

图1 Nishihara蠕变模型

2.2 蠕变本构方程

在Nishihara 模型中，将Newton 黏壶替换为分数阶Abel黏壶得到图2所示的分数阶Nishihara蠕变模型。

图2 分数阶Nishihara蠕变模型

对图2 中的模型本构方程进行分析，其中胡克体的本构方程如式（6）所示：

式中，E0为胡克体弹簧元件的弹性模量，GPa。

1.4 统计学处理 采用 SPSS 22.0 软件进行数据分析。呈正态分布的计量资料以 ±s 表示，两组间比较采用两独立样本 t 检验；呈偏态分布的计量资料以中位数（下四分位数，上四分位数）表示，两组间比较采用 Mann-Whitney U 检验；计数资料以例数和百分数表示，组间比较采用 χ2 检验。以是否合并糖尿病微血管病变为因变量，年龄、体质量、腰围、糖尿病病程、血脂、HbA1c 水平、GA水平、AIR 指数、ACR 指数、INSAUC 和 HOMA-IR指数等为自变量进行 logistic 回归分析。检验水准（α）为 0.05。

黏弹性体的本构模型如式（7）所示：

式中，E1为黏弹性体弹簧元件的弹性模量，GPa；为黏弹性体Abel黏壶的粘性系数。

根据分数阶微积分的数学理论，为简化运算，将Riemann-Liouville 分数阶导数替换为Caputo 分数阶导数。考虑到初始状态下，即t=0 时，εVe=0，经过Laplace双重变换得到式（8）所示的黏弹性体分数阶本构方程：

黏塑性体中包含摩擦滑块，摩擦滑块主要表示为应力阈值点，如式（9）所示：

式中，σp为摩擦滑块应力，MPa；σs为阈值应力，MPa。当摩擦滑块的应力大于阈值应力时，摩擦滑块发生变形。结合考虑黏性系数劣化过程，并经过 Laplace双重变化最终得到式（10）所示的黏塑性体本构方程：

综合考虑图2 所示的组合元件模型结构，最终得到三阶段分数阶蠕变模型的本构方程如式（11）所示：

式中，α、β为函数参数；z为函数自变量；k为函数阶次。

3 围岩岩石变形失稳规律分析

为验证式（13）所示的本构方程的合理性，分析岩石的变形失稳规律，采用山西省晋城市锦辰煤业3 采区巷道围岩岩石进行单轴应力加载实验。3 采区整体为单斜构造，水文地质条件中等，地压正常约为5~9 MPa，3 采区开采条件较好，采区储量为18.04 Mt，预计巷道的投入使用时间为4~5年，因此，将直接面临围岩蠕变变形的影响。本次实验采用三轴流变仪如图3所示。

图3 三轴流变仪

实验的岩石岩性为粉砂岩，仪器加载的轴向应力为5 MPa、6 MPa、7 MPa、8 MPa、9 MPa，最终得到图4所示的岩石加载蠕变曲线。

图4 不同轴向应力加载下蠕变曲线

在实验中考虑到巷道围岩所受的水平应力即围压通常较小，所以仅进行轴向应力的持续加载。根据图4 可知，在长时间的轴向应力作用下，岩石发生了蠕变变形。在5~9 MPa 的不同应力加载试验中，均发生了三阶段的蠕变变形，随着应力的逐渐增加，岩石的变形趋势更加明显。其中第一阶段的初始蠕变阶段主要发生在0 h~25 h 之间，此时岩石发生蠕变变形的速率逐渐减小；第二阶段的稳定蠕变阶段主要发生在25 ~225 h 之间，岩石发生蠕变的速度保持不变，变形速率为0，岩石进入稳定的变形阶段；第三阶段的快速蠕变阶段主要发生在225 h以后，岩石的蠕变变形速度快速增加，并伴随着出现宏观裂隙的产生。在实际工程中，主要原因是在应力的长时间作用下，岩石内部的损伤逐渐积累，在达到破坏的阈值点后发生宏观破坏，此时巷道围岩即进入了快速失稳破化阶段。相应地在所提出的本构方程中，负指数形式的损伤变量可以很好地表征此过程。最终试验数据与本构模型具有很好的拟合效果，如图5所示。

图5 本构模型与实验数据的拟合曲线

本构模型的拟合力学参数拟合结果如表1所示，仅改变应力水平即可对不同应力水平下的岩石蠕变变形曲线进行拟合。根据数据可知，随着应力水平的增加，岩石进入第三阶段（快速蠕变阶段）的时间逐渐缩短。

表1 本构方程参数表

4 本构方程参数分析

根据式（13）所示的本构模型可知，影响蠕变变形的三个主要参数为应力σ、分数阶导数的阶次γ、损伤系数α。其中应力σ主要影响岩石不同阶段的持续时间，当应力越大时，岩石进入不同变形阶段的速度就会加快；分数阶导数的阶次γ主要影响岩石的变形程度，而不影响岩石进入不同变形阶段的速度，当阶次逐渐降低时，岩石变形幅度越小，当分数阶导数阶次为1时，式（13）可以退化为Nishihara 模型。损伤系数α 主要影响岩石第三阶段的变形幅度，当损伤系数越大时，岩石进入第三阶段的损伤累计速度越快，黏性系数的劣化越逐渐明显。

5 结语

1）基于分数阶微积分理论，将传统Newton 黏壶改进为分数阶Abel 黏壶，并引入负指数形式的黏性系数损伤变量，建立了分数阶Nishihara模型的本构方程。

2）对砂岩岩石样品进行单轴应力加载实验，得到5~9 MPa 应力下的蠕变变形曲线，随着应力水平的增加，均发生了三阶段蠕变变形，岩石进入第三阶段快速蠕变阶段的时间逐渐缩短。

3）提出的考虑损伤演化的分数阶Nishihara 模型的本构方程可以很好地表征岩石的蠕变变形。应力、导数阶次、损伤系数是影响蠕变变形的主要参数，当导数阶次为1时，分数阶模型可退化为Nishihara模型。