考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法与流程

本发明涉及数值模拟技术领域，尤其涉及一种考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法。

背景技术：

在深部采矿工程或隧道工程中，高地应力条件下地下工程围岩在开挖前处于三维应力状态，且三个方向的主应力大小不一，此时，我们称岩石处于真三轴应力状态。众多国内外学者采用多种手段研究真三轴应力条件下岩石强度准则以及破坏模式，取得了丰硕的成果。室内试验，被认为是一种直观地、可靠地方式，已经广泛应用于岩石力学与岩石工程之中。因为常规三轴试验无法考虑中间主应力对岩石破坏的影响，mogi于1970年最早设计出了适用于各种岩石类材料的并可以独立施加三个方向主应力的真三轴试验装置。通过一系列的真三轴试验，他发现中间主应力对于岩石的强度以及破坏模式确实存在着非常重要的影响。当中间主应力逐渐增大时，岩石的峰值强度会呈现先增大后降低的趋势；同时，岩石的破裂角(宏观裂纹与最小主应力方向夹角)随着中间主应力的增大而增大，最终由剪切破坏转变为劈裂(板裂化)破坏。

尽管国内外学者在室内真三轴试验装置及其试验方面做出了许多努力和尝试，但是应用数值模拟手段进行真三轴应力条件下岩石类材料破坏特性的研究却较为少见。事实上，采用数值模拟开展相关研究不仅可以节省大量的人力物力，还可以避免试验中所产生的高额费用。当前真三轴数值模拟研究相对匮乏的一个重要原因就是现有商业软件(这里主要指基于有限差分、有限元法原理的数值模拟软件，不包括离散元软件，如pfc等)中内置的本构模型不能够很好反映真三轴条件下岩石类材料的破坏特性。岩石强度准则经过数十年发展，目前已出现多种不同形式的强度准则，比如最大拉应力准则、griffith准则、mohr-coulomb准则、hoek-brown准则等，尽管这些准则在岩石强度理论体系中已成为常用准则，但这些准则都没考虑中间主应力对岩石强度的影响。为了考虑中间主应力影响，不少学者通过理论推导发展了真三轴理论强度准则，比如drucker-prager准则、lade准则以及由此发展起来的一系列理论公式，但是这些准则并不能准确预测岩石的强度，因此并不被广泛的应用于实际之中。

1972年，mogi提出了一个基于广义von-mises准则的真三轴破坏准则。由于该公式是八面体剪应力关于平均有效应力的函数，因此mogi准则充分考虑了中间主应力的作用。尽管这一准则可以很好的预测岩石类材料在真三轴应力条件下的强度特性，但是公式中包含的一些经验常数并没有能够与岩石类材料参数或其固有属性相关联，因此该准则必须通过室内真三轴试验数据进行反演拟合得到，实际意义不大。为了将真三轴破坏准则与常规三轴实验数据相关联，al-ajmi和zimmerman于2005年提出了mogi-coulomb(线性版的mogi准则)破坏准则。该准则不仅考虑了中间主应力的作用，其所包含的参数也可以通过岩石粘聚力和内摩擦角获得。因此，即使不进行室内真三轴试验，这一准则仍然适用。基于此，将mogi-coulomb破坏准则引入flac3d数值模拟软件之中开展真三轴应力条件下岩石类材料破坏特性的研究，是非常有科研和实际意义的。

技术实现要素：

本发明目的在于公开一种考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法，以真实反映真三轴应力条件下岩石类材料强度特性与破坏模式。

为达上述目的，本发明公开的考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法，包括：

在flac3d数值模拟软件之中，用c++语言二次开发基于mogi-coulomb破坏准则的有限差分增量迭代计算格式，得到mogi-coulomb本构模型；

建立岩石类材料真三轴应力条件下数值模型，对该数值模型赋以通过单轴压缩、巴西劈裂试验以及常规三轴压缩试验分别得到基本物理力学参数，所述参数包括弹性模量e、泊松比μ、抗压强度σc、抗拉强度σt、粘聚力c以及内摩擦角

根据得到的所述基本物理力学参数，加载和调用所述mogi-coulomb本构模型进行模拟运算已得到相应的数值模拟计算结果，在该模拟运算过程中(具体的模拟运算过程可参照附图1)，对试样施加三个方向的应力，水平方向的中间主应力以及最小主应力采用应力加载的方式加以控制，轴压方向的最大主应力采用位移加载的方式施加于试样顶部，并对相关应力、应变值进行实时监测与记录。

可选的，在加载和调用所述mogi-coulomb本构模型之前，还包括：

将c++语言二次开发的基于mogi-coulomb破坏准则的有限差分增量迭代计算格式导入至上述flac3d数值模拟软件的头文件和源文件中，并生成了基于mogi-coulomb本构模型的动态链接库文件，并将该动态链接库文件复制在上述flac3d数值模拟软件的安装目录下以供加载和调用。优选地，本发明采用平面截断(tension-cutoff)方式进行调用，可以有效避免对于抗拉强度的过高估计。

优选地，上述mogi-coulomb本构模型与上述flac3d数值模拟软件自带的本构模型继承的是同一个基类，以确保两者的执行效率处于同一水平。必要时，还可以对上述数值模型设置边界条件以对相应方向的端部进行约束。

本发明中，上述mogi-coulomb本构模型主要包括弹性增量法则，破坏(屈服)准则和流动法则，塑性修正以及应变软化参数模型四大部分。分述如下：

mogi-coulomb强度准则是基于广义von-mises准则而建立的，表示为：

τoct＝a+bσm,2

式中，为八面体剪应力；σm,2＝(σ1+σ3)/2，为平均有效应力；σ1，σ2，σ3分别为最大主应力，中间主应力以及最小主应力。

优选地，上述mogi-coulomb本构模型还包括：

不考虑试样的应变硬化行为，将应力应变曲线中峰前阶段全部视为线弹性行为，上述的mogi-coulomb破坏准则即为屈服准则。

优选地，上述mogi-coulomb本构模型还包括：

弹性增量法则用通常的应力应变表示的胡克定律增量表达式为：

δτ＝gδγ^e

δσ＝kδe^e＝kδε^e

其中，g和k分别为剪切模量和体积模量；δτ，δσ分别表示胡克定律的剪应力，主应力增量；δγ^e、δe^e/δε^e分别表示胡克定律的剪应变、主应变增量，上标e表示弹性部分。

优选地，上述mogi-coulomb本构模型还包括：

定义剪切破坏包络线f^s为：

定义拉伸破坏包络线f^t为：

f^t＝σm,2-σ^t

式中，和为与粘聚力和内摩擦角有关的常数，为粘聚剪切抗剪强度，为压力相关系数；σ^t为岩石抗拉强度，σm,2为平均有效应力；

其中，当材料参数不为0时，抗拉强度最大值为：

剪切势函数g^s在采用非关联流动法则时，将内摩擦角转换为剪胀角，则有：

τoct为八面体剪应力，拉伸势函数g^t在采用关联流动法则时，内摩擦角等于剪胀角，则有：

g^t＝σm,2。

优选地，上述mogi-coulomb本构模型还包括：

塑性修正为：

首先考虑剪切破坏，根据流动法则可知：

式中，λ^s为待定义的参数，g^s为基于非关联流动法则的塑性剪切势函数，δγ^p、δe^p分别表示塑性阶段剪应变、主应变增量，上标p表示塑性部分；将以及替代δγ^e,δε^e，由公式：

导出：

式中，δγ^e、δe^e/δε^e分别表示胡克定律的剪应变、主应变增量，δσi为弹性应力增量；si为弹性应变增量的线性函数；为弹性应变增量；s1、s2为由弹性增量法则中获得的应力增量；上标s无实际意义，仅用于表示剪切塑性势函数；

当f＝f^s时，则有：

τoctⁿ＝τoctⁱ-λ^sg

式中，τoct为八面体剪应力，上标n表示新的应力状态，i表示由弹性假设得到的应变和原来的应变之和，bψ表示与剪胀角有关的材料常数，g和k分别为剪切模量和体积模量，σm,2为平均有效应力；参照flac手册，推导出：

新的应力张量分量可以有如下公式表示：

τoctⁿ＝μτoctⁱ

式中μ＝1-λ^sg/τⁱ,由偏应力张量sij表达式，可以得到：

sijⁿ＝μsijⁱ

将μ由τoctⁿ＝μτoctⁱ替代，可得：

最终，新的应力分量可以有以下公式进行表示：

σijⁿ＝sijⁿ+σm,2ⁿδij

式中，δij为克罗内克符号；

考虑拉伸破坏时，由流动法则可知：

式中，上标t无实际意义，仅用于表示拉伸塑性势函数，以区别于剪切塑性势函数；

则有即可导出：

当f＝f^t时，则有：

τoctⁿ＝τoctⁱ

σm,2ⁿ＝σm,2ⁱ-λ^tk

并有代入上式可得：

τoctⁿ＝τoctⁱ

σm,2ⁿ＝σm,2^t

针对拉伸破坏模式，sijⁿ＝sijⁱ，根据对于偏应力得定义，可以得出：

σijⁿ＝sijⁿ+σm,2^tδij。

优选地，上述mogi-coulomb本构模型还包括：

应变软化模型为：

采用塑性剪切应变(剪切硬化参数)以及张拉塑性应变(抗拉硬化参数)两个参数作为塑性内变量，塑性剪切应变增量表达式为：

式中，为塑性主应变增量，分别表示最大，中间以及最小塑性主应变增量；为塑性剪切应变主增量，且有：

抗拉硬化参数e^pt用于计算累计得张拉塑性应变，其增量定义为：

本发明具有以下有益效果：

1、本发明将mogi-coulomb破坏准则嵌入有限差分软件flac3d之中，在进行真三轴应力条件下岩石类材料破坏特性分析时，充分考虑了中间主应力对于试样峰值强度以及破坏模式的影响，使得模拟结果更加符合真实情况，克服了数值模拟软件中既有本构模型(基于有限差分理论的模拟软件，尤其针对于flac3d数值模拟软件)不能够准确描述岩石类材料在真三轴条件下破坏特性的弊端。

2、本发明操作简便，计算精确，仅借助单轴压缩、巴西劈裂以及假三轴(常规三轴)等常规性试验便可以获得本模拟所需全部参数，极大地节省了人力、物力以及实验成本，具有较强的实用性。

3、采用本发明获得的数值模拟结果可以为深部采矿工程、隧道工程等真三轴应力状态下岩体失稳机理、破坏特性分析以及稳定性控制对策研究等提供有效可靠的指导依据。

下面将参照附图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为基于mogi-coulomb本构模型二次开发程序流程图；

图2a为solenhofen石灰岩，dunham白云岩以及mizuho粗面岩试样三维计算模型示意图；

图2b为solenhofen石灰岩，dunham白云岩以及mizuho粗面岩试样二维计算模型示意图；

图3为最小主应力为25mpa时，不同中间主应力条件下dunham白云岩试样全应力-应变曲线图；

图4a为不同最小主应力情况下dunham白云岩试样峰值应力-中间主应力变化图(基于mogi试验结果)；

图4b为不同最小主应力情况下dunham白云岩试样峰值应力-中间主应力变化图(基于本数值模拟结果)；

图5a为不同最小主应力情况下solenhofen石灰岩试样峰值应力-中间主应力变化图(基于mogi试验结果)；

图5b为不同最小主应力情况下solenhofen石灰岩试样峰值应力-中间主应力变化图(基于本数值模拟结果)；

图6a为不同最小主应力情况下mizuho粗面岩试样峰值应力-中间主应力变化图(基于mogi试验结果)；

图6b为不同最小主应力情况下mizuho粗面岩试样峰值应力-中间主应力变化图(基于本数值模拟结果)；

图7a为中间主应力为3mpa，最小主应力为0mpa时汨罗花岗岩破坏模式；

图7b为中间主应力为30mpa，最小主应力为0mpa时汨罗花岗岩破坏模式；

图7c为中间主应力为60mpa，最小主应力为0mpa时汨罗花岗岩破坏模式。

图中：σ1，σ2，σ3分别为最大主应力、中间主应力以及最小主应力，ε1为最大主应力方向所对应的应变值。

具体实施方式

基于上述《发明内容》中的技术方案(相关mogi-coulomb本构模型的推导及开发参照上述发明内容部分，不再赘述)，以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。

实施例1

参见附图1至6可知，根据mogi-coulomb强度准则，将推导出的有限差分增量迭代计算格式嵌入有限差分软件flac3d之中，以实现mogi-coulomb本构模型在flac3d软件中的二次开发及应用；涉及的主要成员如：根据应变增量计算应力张量函数run()函数，其是整个模型开发过程中最重要的函数，在flac3d单元计算中的每一个循环、每一子单元都要调用；其中，程序编写主要包括五部分的内容：①基类的描述；②成员函数的描述；③模型的注册；④模型与flac3d之间的信息交换；⑤模型状态指示器的描述；具体的程序开发的流程图请参照附图1。本发明的具体实施例以mogik在2007年出版的专著《experimentalrockmechanics》以及al-ajmiam和zimmermanrw于2005年发表在internationaljournalofrockmechanicsandminingsciences上的文章“relationbetweenthemogiandthecoulombfailurecriteria”为基础，以其中的三类岩石(dunham白云岩，solenhofen石灰岩以及mizuho粗面岩)在真三轴试验中所获得的数据为依据，与本发明的数值模拟结果进行对比验证。具体步骤如下：

(1)建立长方体岩样三维数值模型

根据文献中试样尺寸，对上述三中岩样建立六面体网格数值模型。试样尺寸均为3cm×1.5cm×1.5cm，模型中网格数与节点数分别为31250以及34476，详见附图2a，2b所示。

(2)赋予mogi-coulomb本构模型相应力学参数

设置三种岩样三维数值模型中本构模型为mogi-coulomb模型，相关物理力学参数可由相关参考文献获得，详见表1。

表1.三种岩样尺寸以及相关物理力学参数

注释:e为弹性模量,μ为泊松比,c、和ψ分别为粘聚力、内摩擦角以及剪胀角，σt为岩石抗拉强度。

3、设置岩样三维数值模型边界条件，应力路径

对试样底部z轴方向(竖直方向，也是最大主应力方向)进行约束。对试样顶部施加均匀荷载，采用位移控制的方式进行加载，对试样x和y方向端面同样施加均匀荷载，采用应力控制的方式进行加载。试样六个面需同时加载，调试好加载速率，保证试样在达到峰值强度以前围压即已达到预定值，当最小主应力(x方向)以及中间主应力(y方向)达到预定值时，保持不变，继续施加z轴方向力，至试样破坏以后停止计算。

4、全程监测、记录试样轴向应力、应变值

采用fish编程语言对轴向应力、轴向应变值以及水平应变值等数据进行实时监测与记录。其中，轴向应力可以通过计算底部平面上所有节点的方向的不平衡力之和，再除以底面积而得到。轴向应变即试样的总的轴向应变量，水平应变可以监测每一个单元的水平应变随计算时步的变化。

5、提取监测数据并分析数值模拟结果

输出最小主应力为25mpa时，不同中间主应力条件下dunham白云岩试样全应力-应变数据，并绘制成图，如图3所示，其中横坐标为应变，纵坐标为应力。由图3可知，试样的应力-应变曲线明显受到了中间主应力值的影响。在加载初期，各曲线基本重合。值得注意的是，由于没有考虑岩石应变硬化行为，因此在峰值以前曲线全部为线弹性阶段。图4a，5a，6a为不同最小主应力情况下dunham白云岩试样，solenhofen石灰岩试样以及mizuho粗面岩试样峰值应力-中间主应力变化图(基于mogik以及al-ajmiam和zimmermanrw文献中的数据绘制)，图4b，5b，6b为不同最小主应力情况下dunham白云岩试样，solenhofen石灰岩试样以及mizuho粗面岩试样峰值应力-中间主应力变化图(基于本发明的数值模拟结果)。通过分析对比可知，本发明的提出一种考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法，所得结果与之前实验数据吻合度较高，真三轴应力条件下岩石类材料峰值强度明显受到了中间主应力的影响，其变化规律呈现出先增大后降低的趋势。然而，基于mohr-coulomb或strain-softeningmohr-coulomb本构模型所得结果没有体现出中间主应力值的影响，即当最小主应力值一定时，不同中间主应力值所对应的峰值强度相同。通过以上分析表明：基于本发明的数值模拟方法能够保证真三轴应力条件下岩石类材料峰值强度数据的可靠性与科学性，同时具有较强的适用性。

实施例2

参见附图1，附图7可知，根据mogi-coulomb强度准则，将推导出的有限差分增量迭代计算格式嵌入有限差分软件flac3d之中，以实现mogi-coulomb本构模型在flac3d软件中的二次开发及应用。本发明的具体实施以双轴压缩下(事实上也是真三轴应力条件的一种，即另最小主应力σ3为0，目的是为了更好的体现出中间主应力对于试样破坏模式的影响)汨罗花岗岩试样破坏模式为例，具体步骤如下：

(1)建立长方体岩样三维数值模型

对汨罗花岗岩试样建立六面体网格数值模型。试样尺寸为5cm×5cm×5cm，模型中网格数与节点数分别为125000以及132651，详见附图7a，7b以及7c所示。

(2)赋予mogi-coulomb本构模型相应力学参数

设置汨罗花岗岩试样三维数值模型中本构模型为mogi-coulomg模型，相关物理力学参数可由单轴压缩试验，巴西劈裂试验以及常规三轴试验获得，详见表2。其中，弹性模量e，泊松比μ由单轴压缩试验获得，粘聚力c以及内摩擦角由常规三轴实验计算得到，抗拉强度由巴西劈裂法获取，剪胀角可通过查阅相关文献或手册得知。残余粘聚力、内摩擦角以及剪胀角均设置为初始粘聚力、内摩擦角以及剪胀角的20％。

表2.汨罗花岗岩试样尺寸以及相关物理力学参数

注释：和ψres分别为最终残余粘聚力、内摩擦角以及剪胀角。

3、设置岩样三维数值模型边界条件，应力路径

对试样底部z轴方向(竖直方向，也是最大主应力方向)进行约束。对试样顶部施加均匀荷载，采用位移控制的方式进行加载，对试样x和y方向端面同样施加均匀荷载，采用应力控制的方式进行加载。试样六个面需同时加载，调试好加载速率，保证试样在达到峰值强度以前围压即已达到预定值，当最小主应力(x方向)以及中间主应力(y方向)达到预定值时(其中，最小主应力始终为0mpa，中间主应力预定值为3，30以及60mpa三种情况)，保持不变，继续施加z轴方向力，至试样破坏并完全达到其残余强度以后方可停止计算。

4、全程监测、记录试样轴向应力、应变值

采用fish编程语言输出xx方向(平行于x方向，即垂直于y方向)应变云图，如图7a，7b，7c所示。

5、提取监测数据并分析数值模拟结果

图7a，7b，7c为最小主应力为0mpa时不同中间主应力汨罗花岗岩破坏模式。为了更加直观清晰地描述出试样的破坏模式，将试样内xx方向应变值小于0.005时的区域设置为白色。由图4a可以看出，当中间主应力均为3mpa时(即接近于单轴压缩条件下)，汨罗花岗岩试验的破坏模式为剪切破坏，或者称之为对顶锥形破坏；随着中间主应力值的提高(30mpa时)，试样宏观裂纹开始增多，且破坏面倾角(宏观裂纹与最小主应力方向夹角)开始有所增大；当中间主应力值增大至60mpa时，破坏面倾角进一步增大，此种情况下的破坏模式类似于实验室所观察到的板裂化破坏。通过分析对比可知，本发明的提出一种考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法，所得到的试样破坏模明显受到了中间主应力的影响，且变化规律与先前一些学者所观察到的试验现象基本一致。通过以上分析表明：基于本发明的数值模拟方法能够保证真三轴应力条件下岩石类材料破坏模式的可靠性与科学性。同时，采用数值模拟的方式可节约大量的人力物力以及试验费用，具有较强的适用性，值得推广。

综上，本实施例公开的考虑中间主应力效应的岩石类材料真三轴试验数值模拟方法，具有以下有益效果：

1、将mogi-coulomb破坏准则嵌入有限差分软件flac3d之中，在进行真三轴应力条件下岩石类材料破坏特性分析时，充分考虑了中间主应力对于试样峰值强度以及破坏模式的影响，使得模拟结果更加符合真实情况，克服了数值模拟软件中既有本构模型(基于有限差分原理)不能够准确描述岩石类材料在真三轴条件下破坏特性的弊端。

2、操作简便，计算精确，仅借助单轴压缩、巴西劈裂以及假三轴(常规三轴)等常规性试验便可以获得本模拟所需全部参数，极大地节省了人力、物力以及实验成本，具有较强的实用性。

3、获得的数值模拟结果可以为深部采矿工程、隧道工程等真三轴应力状态下岩体失稳机理、破坏特性分析以及稳定性控制对策研究等提供有效可靠的指导依据。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。