CN103942387A - 一种基于变分数阶导数建立岩石蠕变本构模型的新方法 - Google Patents

Abstract

本发明公开了一种基于变分数阶导数建立岩石蠕变本构模型的新方法。本方法它包括以下步骤：建立变分数阶导数蠕变模型；获取实验数据；参数拟合。基于本方法建立的模型可以描述全过程蠕变曲线，尤其是非线性加速蠕变特征。本发明中的变分数阶导数可以充分体现岩石在蠕变过程中性质会随时间改变这一特征，并且应用变分数阶导数建立的岩石蠕变本构模型更加简单有效，可以准确的判断和预测岩体的加速变形，做出合理的防护措施，减少工程损失，具有巨大的工程应用价值和经济效益。

Description

一种基于变分数阶导数建立岩石蠕变本构模型的新方法

技术领域：

本发明属于岩石工程技术领域，具体涉及一种基于变分数阶导数建立岩石蠕变本构模型的新方法。

背景技术：

岩石蠕变是指在常值应力持续作用下，岩体变形随时间而持续增长发展的过程。

岩石的蠕变性能是决定工程围岩长期稳定性的重要因素。如软岩流变常常会引起如岩体滑坡、地基失稳、坝基开裂以及矿柱岩爆等自然灾害。而通过岩石蠕变本构模型可准确判断和预测岩体的变形，做出合理的防护措施。因此，对岩石蠕变本构模型的研究是大型岩石工程必不可缺的重要内容，其重要性不言而喻。

由于分数阶导数算子具有记忆效应，近年来不少研究者提出了基于分数阶导数的岩石蠕变本构模型。但是常见的分数阶模型往往都很难反映岩石的全过程蠕变曲线，尤其是非线性加速蠕变阶段。

发明内容：

本发明的目的在于克服现有技术的不足，提供一种基于变分数阶导数建立岩石蠕变本构模型的新方法。

为了解决背景技术所存在的问题，本发明采用以下技术方案：

一种基于变分数阶导数建立岩石蠕变本构模型的新方法，它包括以下步骤：

步骤一：建立变分数阶导数蠕变模型:假设分数阶蠕变模型中的分数阶阶数是一个关于时间的函数，建立起变分数阶导数蠕变模型；

步骤二：获取实验数据:通过室内试验获得岩石出现加速蠕变特征的蠕变数据，进而根据蠕变实验数据绘出蠕变率曲线，从曲线上找出岩石进入加速蠕变阶段的时刻;

步骤三：参数拟合:根据岩石进入加速蠕变阶段的时刻，对变分数阶蠕变模型进行分段拟合，确定出蠕变本构模型的参数值。

进一步的，在步骤一中，分数阶阶数是一个关于时间的函数，可按以下步骤确定函数的形式：

（1）通过室内试验等技术手段得到岩石的蠕变实验曲线；

（2）再根据岩石蠕变实验曲线来确定函数的具体形式：假设分数阶阶数在出现加速蠕变的时刻之前是一个特定的值，出现加速蠕变的时刻之后变成另一个特定的值，这样就确定了函数的形式为分段阶跃函数；

基于Caputo导数定义，分数阶Maxwell由胡克体和Abel粘壶串联而成，

Abel粘壶（A)的应力应变关系为

$\mathrm{\σ}={{\mathrm{\η}}_0}_0^C{D}_t^{\mathrm{\β}}\mathrm{\ϵ}\left(t\right)$

其中表示对ε(t)的β阶分数阶导数，定义如下

式中：n为大于β的最小正整数；f⁽ⁿ⁾(τ)为函数f(τ)的n阶导数；

则式中σ为Abel粘壶的应力，即总应力；η₀为其粘性系数；

认为Abel粘壶的分数阶阶数是一个关于时间的函数，即β=α(t)，那么式中α(t)中为关于时间t的分数阶函数，η_α(t)为对应的粘性系数；

胡克体(H)的应力应变关系为式中σ为胡克体的应力，即总应力；E₀为胡克体中弹簧的弹性模量；

综合考虑两部分应变，则变分数阶导数蠕变模型的本构方程可表示为：

$\mathrm{\ϵ}\left(t\right)=\frac{\mathrm{\σ}}{E_0}+\frac{\mathrm{\σ}}{{\mathrm{\η}}_{\mathrm{\α}\left(t\right)}}\frac{t^{\mathrm{\α}\left(t\right)}}{\mathrm{\Γ}(\mathrm{\α}\left(t\right)+1)}.$

进一步的，在步骤二中，假设α(t)在t_c点以前是个特定的值，在t_c点以后变成另一个特定的值，即：

$\mathrm{\α}\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\α}}_1,& t\≤t_c\\ {\mathrm{\α}}_2,& t>t_c\end{array}\right.,$

式中t_c为岩石进入加速蠕变的时刻；

其蠕变本构方程为： $\mathrm{\ϵ}\left(t\right)=\frac{\mathrm{\σ}}{E_0}+\frac{\mathrm{\σ}}{{\mathrm{\η}}_{\mathrm{\α}\left(t\right)}}\frac{t^{\mathrm{\α}\left(t\right)}}{\mathrm{\Γ}(\mathrm{\α}\left(t\right)+1)},$ 其中 $\mathrm{\α}\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\α}}_1,& t\≤t_c\\ {\mathrm{\α}}_2,& t>t_c\end{array}\right..$

进一步的，在步骤三中，利用非线性最小二乘法拟合出蠕变本构方程：中的E₀、η_α1、η_α2、α₁、α₂5个参数。

本发明对比现有技术，有如下的有益效果：基于本方法建立的模型可以描述全过程蠕变曲线，尤其是非线性加速蠕变特征。本发明中的变分数阶导数可以充分体现岩石在蠕变过程中性质会随时间改变这一特征，并且应用变分数阶导数建立的岩石蠕变本构模型更加简单有效，可以准确的判断和预测岩体的加速变形，做出合理的防护措施，减少工程损失，具有巨大的工程应用价值和经济效益。

附图说明：

图1是为分数阶Maxwell模型图。

图2是岩石的典型蠕变曲线。

图3是通过室内试验得到的盐岩蠕变曲线。

图4是根据蠕变数据得到的蠕变率曲线。

图5是利用变分数阶导数蠕变模型对试验数据的拟合结果。

图6是利用分数阶导数蠕变模型对试验数据的拟合结果。

图7是本发明方法流程图。

具体实施方式：

下面结合附图和具体实施方式对本发明作进一步描述：

图7是本发明方法流程图。一种基于变分数阶导数建立岩石蠕变本构模型的新方法，它包括以下步骤：

步骤一：建立变分数阶导数蠕变模型:假设分数阶蠕变模型中的分数阶阶数是一个关于时间的函数，建立起变分数阶导数蠕变模型；

步骤二：获取实验数据:通过室内试验获得岩石出现加速蠕变特征的蠕变数据，进而根据蠕变实验数据绘出蠕变率曲线，从曲线上找出岩石进入加速蠕变阶段的时刻;

步骤三：参数拟合:根据岩石进入加速蠕变阶段的时刻，对变分数阶蠕变模型进行分段拟合，确定出蠕变本构模型的参数值。

方法原理：

由于常见的分数阶导数蠕变模型能很好的描述减速蠕变和稳定蠕变阶段特征，却不能反映加速蠕变特征。那么我们就可以假设分数阶阶数在出现加速蠕变的时刻之前是一个特定的值，出现加速蠕变的时刻之后变成另一个特定的值，这样就确定了函数的形式为分段阶跃函数。

确定求导阶数的函数形式关键在于判断岩石进入加速蠕变的时间点。因为在稳定蠕变阶段其蠕变率是几乎不变的，到了加速阶段其蠕变率会突然增大，所以可根据蠕变实验数据绘出蠕变率曲线，从曲线上找出岩石进入加速蠕变阶段的时刻。

进一步的，在步骤一中，分数阶阶数是一个关于时间的函数，可按以下步骤确定函数的形式：

（1）通过室内试验等技术手段得到岩石的蠕变实验曲线；

（2）再根据岩石蠕变实验曲线来确定函数的具体形式：假设分数阶阶数在出现加速蠕变的时刻之前是一个特定的值，出现加速蠕变的时刻之后变成另一个特定的值，这样就确定了函数的形式为分段阶跃函数；

基于Caputo导数定义，分数阶Maxwell由胡克体和Abel粘壶串联而成，

Abel粘壶（A)的应力应变关系为

$\mathrm{\σ}={{\mathrm{\η}}_0}_0^C{D}_t^{\mathrm{\β}}\mathrm{\ϵ}\left(t\right)$

其中表示对ε(t)的β阶分数阶导数，定义如下

式中：n为大于β的最小正整数；f⁽ⁿ⁾(τ)为函数f(τ)的n阶导数；

则式中σ为Abel粘壶的应力，即总应力；η₀为其粘性系数；

认为Abel粘壶的分数阶阶数是一个关于时间的函数，即β=α(t)，那么式中α(t)中为关于时间t的分数阶函数，η_α(t)为对应的粘性系数；

胡克体(H)的应力应变关系为式中σ为胡克体的应力，即总应力；E₀为胡克体中弹簧的弹性模量；

综合考虑两部分应变，则变分数阶导数蠕变模型的本构方程可表示为：

$\mathrm{\ϵ}\left(t\right)=\frac{\mathrm{\σ}}{E_0}+\frac{\mathrm{\σ}}{{\mathrm{\η}}_{\mathrm{\α}\left(t\right)}}\frac{t^{\mathrm{\α}\left(t\right)}}{\mathrm{\Γ}(\mathrm{\α}\left(t\right)+1)}.$

进一步的，在步骤二中，假设α(t)在t_c点以前是个特定的值，在t_c点以后变成另一个特定的值，即：

$\mathrm{\α}\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\α}}_1,& t\≤t_c\\ {\mathrm{\α}}_2,& t>t_c\end{array}\right.,$

式中t_c为岩石进入加速蠕变的时刻；

其蠕变本构方程为： $\mathrm{\ϵ}\left(t\right)=\frac{\mathrm{\σ}}{E_0}+\frac{\mathrm{\σ}}{{\mathrm{\η}}_{\mathrm{\α}\left(t\right)}}\frac{t^{\mathrm{\α}\left(t\right)}}{\mathrm{\Γ}(\mathrm{\α}\left(t\right)+1)},$ 其中 $\mathrm{\α}\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\α}}_1,& t\≤t_c\\ {\mathrm{\α}}_2,& t>t_c\end{array}\right..$

进一步的，在步骤三中，利用非线性最小二乘法拟合出蠕变本构方程：中的E₀、η_α1、η_α2、α₁、α₂5个参数。

图1是为分数阶Maxwell模型图。将整数阶Maxwell模型中的牛顿粘壶变为Abel粘壶，就是分数阶Maxwell模型。

图2是岩石的典型蠕变曲线。岩石的典型蠕变曲线一般分为三阶段：1、减速蠕变阶段；2、稳定蠕变阶段；3、加速蠕变阶段

图3是通过室内试验得到的盐岩蠕变曲线。图3是通过室内试验获得盐岩蠕变数据，并绘制成而成的曲线。

图4是根据蠕变数据得到的蠕变率曲线。蠕变数据进行差分求导得到蠕变率曲线。

图5是利用变分数阶导数蠕变模型对试验数据的拟合结果,是根据表1中得到的参数值代入蠕变本构方程中绘出的蠕变理论曲线。从图5可以看出，由于采用了本文提出的变分数阶导数建立的蠕变模型能很好的反映盐岩加速蠕变特征。

表1是基于盐岩蠕变实验的参数拟合结果

(表1中σ表示实验中的应力大小、E₀表示试件的弹性模量、η_α1表示蠕变第一阶段Abel粘壶的粘性系数、η_α2表示第二阶段Abel粘壶的粘性系数、α₁表示第一阶段的分数阶导数阶数、α₂表示第二阶段的分数阶导数阶数。)

图6是利用分数阶导数蠕变模型对试验数据的拟合结果。从图6可以看出，常用的分数阶蠕变模型是不能描述加速蠕变特征的。

本发明首次提出基于变分数阶导数建立岩石蠕变本构模型的新方法，并利用室内试验等技术手段找到岩石进入加速蠕变的时刻，确定了分数阶阶数的分段函数形式，从而建立了变分数阶导数蠕变模型，和常见的分数阶导数蠕变模型相比，基于变分数阶导数建立的蠕变模型能够很好的反映蠕变过程中的加速阶段；本发明中的变分数阶导数可以充分体现岩石在蠕变过程中性质会随时间改变这一特征；本发明提出基于变分数阶导数建立岩石蠕变本构模型的新方法，使得所建模型简单，参数确定容易，可以更准确地判断和预测岩体的变形，做出合理的防护措施，减少工程损失。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于变分数阶导数建立岩石蠕变本构模型的新方法，其特征在于，它包括以下步骤：

步骤一：建立变分数阶导数蠕变模型:假设分数阶蠕变模型中的分数阶阶数是一个关于时间的函数，建立起变分数阶导数蠕变模型；

步骤二：获取实验数据:通过室内试验获得岩石出现加速蠕变特征的蠕变数据，进而根据蠕变实验数据绘出蠕变率曲线，从曲线上找出岩石进入加速蠕变阶段的时刻;

步骤三：参数拟合:根据岩石进入加速蠕变阶段的时刻，对变分数阶蠕变模型进行分段拟合，确定出蠕变本构模型的参数值。

2.根据权利要求1所述的一种基于变分数阶导数建立岩石蠕变本构模型的新方法，其特征在于，在步骤一中，分数阶阶数是一个关于时间的函数，可按以下步骤确定函数的形式：

（1）通过室内试验等技术手段得到岩石的蠕变实验曲线；

（2）再根据岩石蠕变实验曲线来确定函数的具体形式：假设分数阶阶数在出现加速蠕变的时刻之前是一个特定的值，出现加速蠕变的时刻之后变成另一个特定的值，这样就确定了函数的形式为分段阶跃函数；

基于Caputo导数定义，分数阶Maxwell由胡克体和Abel粘壶串联而成，

Abel粘壶（A)的应力应变关系为

$\mathrm{\σ}={{\mathrm{\η}}_0}_0^C{D}_t^{\mathrm{\β}}\mathrm{\ϵ}\left(t\right)$

其中表示对ε(t)的β阶分数阶导数，定义如下

式中：n为大于β的最小正整数；f⁽ⁿ⁾(τ)为函数f(τ)的n阶导数；

则式中σ为Abel粘壶的应力，即总应力；η₀为其粘性系数；

认为Abel粘壶的分数阶阶数是一个关于时间的函数，即β=α(t)，那么式中α(t)中为关于时间t的分数阶函数，η_α(t)为对应的粘性系数；

胡克体(H)的应力应变关系为式中σ为胡克体的应力，即总应力；E₀为胡克体中弹簧的弹性模量；

综合考虑两部分应变，则变分数阶导数蠕变模型的本构方程可表示为：

$\mathrm{\ϵ}\left(t\right)=\frac{\mathrm{\σ}}{E_0}+\frac{\mathrm{\σ}}{{\mathrm{\η}}_{\mathrm{\α}\left(t\right)}}\frac{t^{\mathrm{\α}\left(t\right)}}{\mathrm{\Γ}(\mathrm{\α}\left(t\right)+1)}.$

3.根据权利要求1所述的一种基于变分数阶导数建立岩石蠕变本构模型的新方法，其特征在于，在步骤二中，

假设α(t)在t_c点以前是个特定的值，在t_c点以后变成另一个特定的值，即：

$\mathrm{\α}\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\α}}_1,& t\≤t_c\\ {\mathrm{\α}}_2,& t>t_c\end{array}\right.,$

式中t_c为岩石进入加速蠕变的时刻；

其蠕变本构方程为： $\mathrm{\ϵ}\left(t\right)=\frac{\mathrm{\σ}}{E_0}+\frac{\mathrm{\σ}}{{\mathrm{\η}}_{\mathrm{\α}\left(t\right)}}\frac{t^{\mathrm{\α}\left(t\right)}}{\mathrm{\Γ}(\mathrm{\α}\left(t\right)+1)},$ 其中 $\mathrm{\α}\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\α}}_1,& t\≤t_c\\ {\mathrm{\α}}_2,& t>t_c\end{array}\right..$

4.根据权利要求1所述的一种基于变分数阶导数建立岩石蠕变本构模型的新方法，其特征在于，在步骤三中，

利用非线性最小二乘法拟合出蠕变本构方程：中的E₀、η_α1、η_α2、α₁、α₂5个参数。