CN104361211B - 岩石统计损伤本构模型的构建和应用方法 - Google Patents

Abstract

本发明提出了一种岩石统计损伤本构模型，采用Mohr‑Coulomb准则作为岩石微元强度，和采用工程可靠度领域中常用的分布形式，即对数正态分布，作为微元强度的分布状态，利用比较成熟且应用比较广泛的Mohr‑Coulomb准则作为微元强度破坏准则，同时给出该本构模型的损伤参数A、B、C、H与围岩压力σ₃之间的函数关系，从而使得该本构模型可以直接采用数值计算软件应用于岩土工程的计算中，并且在实际的工程应用中验证了效果。

Description

岩石统计损伤本构模型的构建和应用方法

技术领域

本发明涉及一种岩石本构模型的构建方法，属于岩土工程本构模型研究领域，适用于采矿、边坡、隧道、公路、坝基等各类岩土工程本构模型的研究。

具体可以指导(1)地下采矿的安全生产，即在确保地面已有建构筑物设施安全的前提下，最大限度地科学开采地下矿产资源。(2)露天矿山最终境界的合理圈定，即确定合适的坡面角和最终边坡角，以确保在安全的边坡环境下最大限度的回收矿产资源，减少剥岩量，降低采矿成本。(3)岩土工程领域中其它涉及该本构模型的相关问题。

技术背景

岩石是一种非均质的地质材料，本构关系复杂，其研究一直是岩石力学和岩土工程的一个重要方面，受到工程界的普遍重视。近年来许多学者从损伤力学的角度出发，考虑岩石的损伤过程，试图找到符合岩石全过程应力-应变关系曲线的本构模型，不失为一种较好的方法。

本构模型是岩土工程力学计算的基础，由于本构模型构建不合理而导致岩土工程事故、矿山资源的浪费和岩土工程经济损失等时有发生。因此，深入开展本构模型的研究是保障岩土工程合理设计、安全施工、确保效益的基础，具有重要的科学研究意义和工程应用价值。

目前，以损伤力学作为深入研究本构模型的一种有效的技术手段。损伤力学的研究方法大致分为三种：金属物理学方法、维象学方法和统计学方法。本发明是在损伤力学的统计学方法下进行研究的，统计损伤力学是研究岩石损伤本构模型的一种有效方法。

目前在应用统计损伤力学研究岩石损伤本构模型中，D.Krajcinovic、唐春安等提出从岩石内部缺陷分布随机性出发，采用岩石的轴向应变作为微元强度，利用微元强度服从Weibull分布建立岩石统计损伤本构模型，由于轴向应变无法准确表示微元强度，因而该本构模型具有一定的局限性；曹文贵、赵明华等探讨了基于Weibull分布的岩石损伤软化参数与围压的关系，提出新的岩石微元强度表示法，同时采用Drucker-Prager准则建立了岩石统计损伤本构模型，由于Drucker-Prager准则较为保守，使得采用该准则建立起来的岩石微元强度缺少合理性，采用该准则构建起来的岩石统计损伤本构模型应用于工程计算将会导致计算结果偏于安全。

张明、王菲、杨强等人在论文“基于三轴压缩试验的岩石统计损伤本构模型”(参见《岩土工程学报》，2013年11月，第35卷第11期)中以Mohr-Coulomb准则作为岩石微元强度，以服从对数正态分布为基础构建了常规三轴压缩条件下岩石统计损伤本构模型，

式中，c₂＝2vσ₃，a＝1/ξ，b＝-(ln S₀)/ξ，S₀和ξ为分布参数。

φ(内摩擦角)、v(泊松比，v相当于本申请的μ)、E(弹性模量)为已知，σ₃(围岩压力)为试验条件(可事先给定)，c₁、c₂可由上述公式计算得出，σ₁(最大主应力)、ε₁(最大主应力对应的应变)之间的函数关系可以通过试验获得，ξ、S₀为分布参数，可以由一个试验条件(围岩压力σ₃)根据试验结果获得的σ₁、ε₁之间的函数关系曲线反推确定，从而确定式(1)(该试验曲线)的一个σ₃给定条件下的具体数学表达式。

论文给出了不同围岩压力σ₃下的岩石统计损伤本构模型表达式：

当σ₃＝3.45MPa时，

σ₁＝95000ε₁+1.73，S≤0

当σ₃＝6.9MPa时，

σ₁＝95000ε₁+3.45，S≤0

当σ₃＝13.8MPa时，

σ₁＝95000ε₁+6.90，S≤0

当σ₃＝27.6MPa时，

σ₁＝95000ε₁+13.80，S≤0

当σ₃＝23.5MPa时，

σ₁＝29000ε₁+11.75，S≤0

当σ₃＝50MPa时，

σ₁＝29000ε₁+25.00，S≤0

论文通过上述六组试验数据，得到了在六个不同围岩压力σ₃条件下的六组不同的分布参数ξ、S₀，从而确定了a和b₁，并最终确定了六个围岩压力σ₃下的岩石统计损伤本构模型表达式，而该论文没有阐述ξ、S₀是否与围岩压力σ₃存在必然联系，也没有用一种确定的函数方程形式表达他们之间的关系。而为了使该本构模型能够应用到工程计算中，就要求该本构模型在连续变化的不同围岩压力σ₃的条件下都要有相对应的函数关系表达式。虽然该本构模型为表达应力与应变之间的关系方面提供了理论基础，但是将该本构模型直接应用于数值计算软件，以解决生产中的实际问题时，由于该本构模型没有建立起在连续变化的不同围岩压力σ₃条件下相对应的函数关系表达式，导致将其应用于数值计算软件计算围岩压力σ₃参数变化过程中需要连续不断地输入与连续变化的围岩压力σ₃参数相对应的该本构模型的参数值，这就决定了该本构模型不能直接在数值计算软件中得以应用推广，更不能在岩土工程领域的工程应用计算方面发挥有效的指导作用。

总之，现有技术中，各岩石统计损伤本构模型的应用范围具有一定的局限性的困难，如：D.Krajcinovic、唐春安等提出从岩石内部缺陷分布随机性出发，采用岩石的轴向应变作为微元强度，利用微元强度服从Weibull分布建立岩石统计损伤本构模型，由于轴向应变无法准确表示微元强度，因此，应用轴向应变表示微元强度构建起来的岩石统计损伤本构模型也无法准确地表示岩石应力与应变之间的关系。曹文贵、赵明华等探讨了基于Weibull分布的岩石损伤软化参数与围压的关系，提出新的岩石微元强度表示法，同时采用Drucker-Prager准则建立了岩石统计损伤本构模型，Drucker-Prager准则较为保守，因此，采用该准则构建起来的岩石统计损伤本构模型应用于工程计算将会导致计算结果偏于安全，具体到矿山开采领域，将会导致计算开采范围比可开采的实际范围要小，造成资源的浪费。张明、王菲、杨强以Mohr-Coulomb准则作为岩石微元强度，服从对数正态分布为基础构建了常规三轴压缩条件下岩石统计损伤本构模型，该模型只能通过试验获得随围岩压力σ₃变化的有限个离散的应力与应变关系表达式，不能获得随围岩压力σ₃连续变化的应力与应变函数关系表达式，距离实际应用于数值计算软件还存在一定差距。

发明内容

本发明的目的在于针对现有技术中所存在的技术问题，提出一种新型的岩石统计损伤本构模型，采用Mohr-Coulomb准则作为岩石微元强度，克服了采用Drucker-Prager准则较为保守的缺点；采用工程可靠度领域中常用的分布形式(对数正态分布)作为微元强度的分布状态，比较成熟且应用比较广泛的Mohr-Coulomb准则作为微元强度破坏准则，同时给出该本构模型参数与围岩压力σ₃之间的函数关系，具体由损伤参数A、B、C、H表示，由此弥补了上述张明等人提出的岩石统计损伤本构模型的不足，使得该本构模型可以直接采用数值计算软件应用于岩土工程的计算中，并且在实际的工程应用中验证了其效果。

根据本发明的第一个实施方案，提供一种构建岩石统计损伤本构模型的方法，该方法包括以下步骤：

1)对于岩石样品进行岩石室内力学试验，以获取岩石的弹性模量E、泊松比μ，和内摩擦角φ试验参数；

2)构建本构模型，该步骤包括以下的子步骤：

2.1)利用摩尔库伦准则作为正六面体岩石微元强度破坏准则，和微元强度服从对数正态分布准则，建立岩石统计损伤本构模型表达式(1)：

式(1)中：σ₁、σ₂、σ₃为正六面体微元各面上的主应力，其中σ₁为最大主应力；σ₂为中间主应力；σ₃为最小主应力；ε₁为最大主应力对应的应变；D为岩石损伤变量；Φ为标准正态分布函数；

式(1)中：其中σ₁ ^*为有效最大主应力，σ₃ ^*为有效最小主应力；F为有关有效最大主应力σ₁ ^*、有效最小主应力σ₃ ^*、内摩擦角φ的函数，F表示微元强度；

根据J.Lemaitre应变等价性假说，有效应力与主应力及岩石损伤变量之间具有如下关系：σ_n ^*＝σ_n/[1-D]，n＝1或3，σ_n ^*表示有效应力，σ_n表示主应力；

2.2)为了简化试验和方便计算，忽略中间主应力σ₂的影响，将正六面体微元简化为圆柱体微元，此时，σ₃转化为围岩压力，将式(1)变成下式(2)：

在以上式(2)中，中间参数F₀是微元强度F满足对数正态分布条件下的均值，它与围岩压力σ₃之间存在线性关系(3)：

F₀＝Aσ₃+B (3)

式(3)中，损伤参数A、B是F₀与σ₃之间的统计概率的拟合线性关系系数；和

在以上式(2)中，中间参数S₀是微元强度F满足对数正态分布条件下的方差，它与围岩压力σ₃之间存在线性关系(4)：

S₀＝Cσ₃+H (4)

式(4)中，损伤参数C、H是S₀与σ₃之间的统计概率的拟合线性关系系数；

2.3)根据上述关系得出岩石统计损伤本构模型的最终表达式(5)：

优选的是，以上所述的步骤1)是：

针对岩石样品，通过室内单轴压缩试验，获得岩石的弹性模量E、泊松比μ试验参数；通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t试验参数；通过室内三轴压缩试验，获得岩石的粘结力c和内摩擦角φ；以及(例如通过室内试验，如采用天平法或密度仪法)测定岩石的密度ρ。

一般来说，式(3)中的A、B是通过不同围岩压力σ₃条件下的不同F₀的多组数据拟合获得的，一个围岩压力σ₃对应一个F₀，不同的围岩压力σ₃对应的F₀不同。

一般来说，式(4)中的C、H是通过不同围岩压力σ₃条件下的不同S₀的多组数据拟合获得的，一个围岩压力σ₃对应一个S₀，不同的围岩压力σ₃对应的S₀不同。

当岩石损伤变量D满足0≤D≤1的条件时，从而有下式(6)：

当岩石损伤变量D不在0≤D≤1的范围内时，岩石没有发生损伤，公式(5)不适用，此时采用σ₁＝Eε₁+2μσ₃。

通常，在实验室做出各地质岩层的岩石样品在不同的围岩压力σ₃(例如0Mpa、3.45Mpa、6.9Mpa、13.8Mpa、27.6Mpa)情况下的全应力应变曲线，利用通过线性拟合的方法得出中间参数F₀、S₀，然后根据式(3)F₀与围岩压力σ₃的关系采用线性拟合方法确定损伤参数A、B，和根据式(4)S₀与围岩压力σ₃的关系采用线性拟合方法确定损伤参数C、H。

根据本发明的第二个实施方案，提供一种岩石统计损伤本构模型在岩土工程中的应用方法，该应用方法包括以下步骤：

(A)针对岩石样品，通过室内单轴压缩试验，获得岩石的弹性模量E、泊松比μ试验参数；通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t试验参数；通过室内三轴压缩试验，获得岩石的粘结力c和内摩擦角φ；以及测定岩石的密度ρ；

(B)将以上所述方法获得的岩石统计损伤本构模型(5)应用于数值计算软件，获得本构模型的二次开发文件；

(C)根据岩石室内力学试验所测得的试验参数，按下述公式计算岩石体积模量K：

同时，计算岩石切变模量G，按下述公式计算：

在数值计算软件中输入岩石室内力学试验所测得的单轴抗拉强度σ_t、粘结力c、密度ρ试验参数；

然后应用步骤(B)中获得的二次开发文件进行工程计算，获得工程参数，和

(D)按照步骤(C)中所获得的工程参数进行工程应用。

这里所述的工程包括：矿石开采(尤其地下矿产资源开采)；或水利工程；或市政地下工程等。所述的工程参数包括：受力岩石在x方向位移，y方向位移，z方向位移即沉降量。此外，所述的工程参数还包括受力岩石在x方向的应力、应变量，y方向的应力、应变量，z方向的应力、应变量。

优选的是，数值计算软件是FLAC^3D、ANSYS或MIDAS，更优选是FLAC^3D。

当数值计算软件是FLAC^3D时，将本发明的岩石统计损伤本构模型应用于该软件，获得该本构模型二次开发文件usermohr.h和usermohr.cpp，然后编译生成usermohr.dll文件并在步骤(B)的工程计算时加载。

根据本发明的第三个实施方案，提供一种通过构建岩石统计损伤本构模型在岩土工程中进行实际工程应用的方法，该方法包括：

I)针对岩石样品，通过室内单轴压缩试验，获得岩石的弹性模量E、泊松比μ试验参数；通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t试验参数；通过室内三轴压缩试验，获得岩石的粘结力c和内摩擦角φ；以及测定岩石的密度ρ；

II)根据以上所述的本发明的第一个实施方案所述的方法，构建以上所述的岩石统计损伤本构模型(5)；

III)将以上II)所述的方法构建的岩石统计损伤本构模型(5)应用于数值计算软件，获得该本构模型的二次开发文件；

IV)然后应用该文件进行工程计算，获得工程参数；和

V)按照步骤IV)中所获得的工程参数进行实际工程应用。

一般，在以上步骤(IV)的计算过程中输入岩石室内力学试验所测得的单轴抗拉强度σ_t、粘结力c、密度ρ试验参数。

本发明的岩石统计损伤本构模型及其在岩土工程中的应用方法，能够在确保地面已有建构筑物设施安全的前提下，最大限度地科学开采地下矿产资源，以及露天矿山最终境界的合理圈定，即确定合适的坡面角和最终边坡角，以确保在安全的边坡环境下最大限度的回收矿产资源。

总之，本发明的最终的岩石统计损伤本构模型表达式如下：

σ₁ ^*＝σ₁/[1-D]

σ₃ ^*＝σ₃/[1-D]

F₀＝Aσ₃+B

S₀＝Cσ₃+H

如以上所述，本发明的岩石统计损伤本构模型中的参数有弹性模量E、泊松比μ、内摩擦角φ，上述各参数对于本领域的技术人员来说均可以在实验室内按照本领域的常规方法通过试验获得；该本构模型中的损伤参数A、B和C、H可以根据式(3)F₀与围岩压力σ₃的关系采用线性拟合方法确定损伤参数A、B，和根据式(4)S₀与围岩压力σ₃的关系采用线性拟合方法确定损伤参数C、H，其中中间参数F₀、S₀利用在实验室做出各地质岩层的岩石样品试样在不同的围岩压力σ₃情况下的全应力应变曲线，通过线性拟合的方法获得。因此，这使得本发明的岩石统计损伤本构模型具有实用价值。

在本申请中，采用式(5)计算σ₁的条件是：0≤D≤1。当D不在0≤D≤1的范围内时，岩石没有发生损伤，公式(5)不适用，此时采用σ₁＝Eε₁+2μσ₃。

本发明的一个目的就是找出岩石随围岩压力σ₃连续变化的应力与应变之间的函数关系表达式。该本构模型具体可以解决如下问题：(1)地下采矿的安全生产与资源回收，即在确保地面已有建构筑物设施安全的前提下，最大限度地科学开采地下矿产资源。(2)露天矿山最终境界的合理圈定，即确定合适的坡面角和最终边坡角，以确保在安全的边坡环境下最大限度的回收矿产资源。(3)岩土工程领域中其它有关应力与应变关系的相关问题，如市政地下工程。

另外，值得指出的是，本发明对岩石取样没有任何限制，符合各行业领域相关采样设计规范要求即可。不同的采样设计方法只会影响岩石样品的代表性，并不会影响本发明所述的岩石统计损伤本构模型中各参数间的相关关系。

本发明对岩石室内力学试验的测试过程和条件及测试设备没有任何限制，比如岩石的密度ρ测试方法可以采用天平法或密度仪法，天平仪器可以采用全自动、半自动或手动天平仪。密度仪可以采用静压式、振动式、浮子式和放射性同位素式等类型的密度仪。使用不同的仪器设备和方法只会影响岩石的弹性模量E、泊松比μ，岩石的单轴抗拉强度σ_t，岩石的粘结力c、内摩擦角φ，岩石的密度ρ基础参数的精度，并不会影响本发明所述的岩石统计损伤本构模型中各参数间的相关关系。

上述基础参数属于岩土工程领域通用型的基础力学参数，本领域技术人员可以通过相关规范规定获取上述基础参数。

本发明的优点或有益技术效果

本发明的岩石统计损伤本构模型成功地解决了岩石随围岩压力σ₃连续变化的应力与应变之间的函数关系问题，即公式(5)是本发明的核心。本发明是将公式(5)应用于现有的数值计算软件，如：FLAC^3D数值计算软件，非常理想地解决了以下问题：(1)地下采矿的安全生产与资源回收，即在确保地面已有建构筑物设施安全的前提下，最大限度地科学开采地下矿产资源；(2)露天矿山最终境界的合理圈定，即确定合适的坡面角和最终边坡角，以确保在安全的边坡环境下最大限度的回收矿产资源；和(3)岩土工程领域中其它有关受力与应变关系的相关问题。因此，本发明将公式(5)应用于现有的数值计算软件，如：FLAC^3D数值计算软件，就能够合理圈定地下采矿范围。

本发明提供一种原理简单、操作方便、成本低廉、结果合理可靠的岩石统计损伤本构模型，应用该本构模型在地下采矿工程设计时，可以合理圈定开采范围，在保障地面已有建构筑物设施安全的前提下，最大限度地回收地下矿产资源；在露天矿山采矿设计中，能够合理确定坡面角和最终边坡角，在边坡安全的前提下，降低剥岩量，从而降低露天采矿成本，且最大限度地回收矿产资源。

附图说明

图1是不同围岩压力σ₃下X与Y的散点图。

图2是不同围岩压力σ₃下三轴试验的拟合曲线。

图3是F₀与岩石统计损伤本构模型的关系曲线。

图4是S₀与岩石统计损伤本构模型的关系曲线。

图5是参数S₀与围岩压力σ₃关系曲线。

图6是参数F₀与围岩压力σ₃关系曲线。

图7是岩石统计损伤本构模型计算流程图。

图8是岩石统计损伤本构模型校验结果对比图。

图9是采矿范围以及地面构筑物关系图。

图10是坝体X方向位移比较图。

图11是坝体Y方向位移比较图。

图12是坝体Z方向位移比较图。

图13是FLAC^3D计算中应用本申请中所述的岩石统计损伤本构模型的计算流程图。

图14是地质体、保安矿柱及地表构筑物关系图。

图15是不同围岩压力σ₃下全应力应变曲线。

图16是S₀与围岩压力σ₃关系曲线。

图17是F₀与围岩压力σ₃关系曲线。

具体的实施方式

下面相对详细地描述本发明的具体实施方式：

1、岩石统计损伤本构模型的建立

1.1)损伤变量与本构关系的建立

根据J.Lemaitre应变等价性假说(参见Lemaitre J.A continuous damagemechanics model for ductile fracture[J].J.Engng.Master.Tech.，1985，107(1)：83-89)，建立岩石本构关系如下：

[σ^*]＝[σ]/[1-D]＝[C][ε]/(1-D) (1’)

式中：[C]为岩石材料的弹性矩阵，[σ^*]为有效应力矩阵，[σ]为名义应力矩阵，[ε]为应变矩阵，D为岩石损伤变量。建立岩石损伤本构模型的关键在于岩石损伤变量D的确定，采用轴线应变表示岩石微元强度存在严重的不足，依据曹文贵等提出的新的岩石微元强度表示法(参见曹文贵，赵明华，刘成学.基于Weibull分布的岩石损伤软化模型及其修正方法研究[J].岩石力学与工程学报，2004，23(19)：3226-3231)，岩石破坏准则如下：

f(σ^*)-k₀＝0 (2’)

式中：k₀为岩石内摩擦角和粘聚力有关的常数，F＝f(σ^*)反应了岩石破坏的危险程度，可作为岩石的微元强度，破坏概率为P[f(σ^*)]，定义损伤变量与破坏概率的关系如下：

1.2岩石微元强度的确定

由以上1.1可知，F＝f(σ^*)表示岩石微元强度，直接取决于岩石的破坏准则，鉴于Mohr-Coulomb准则来自于试验，参数简单、应用普遍。本申请采用该准则作为岩石的微元强度，Mohr-Coulomb准则如下：

τ＝c+σtanφ (4’)

式中：τ为剪应力，c为粘结力，σ为名义正应力，φ为内摩擦角。

因此，准则可采用σ₁、σ₃表示如下：

式中：σ₁为名义最大主应力，σ₃为名义最小主应力，其余同前。

根据公式(5’)得出

式中：σ₁ ^*为有效最大主应力，σ₃ ^*为有效最小主应力。

1.3基于对数正态分布的岩石损伤模型

在以往的岩石损伤本构模型研究中，微元强度通常假设服从Weibull分布，但是在结构可靠度分析中应用较普遍的是以对数正态分布假设作为结构抗力的概率分布，本申请中尝试采用岩石微元强度服从对数正态分布的假设建立岩石统计损伤本构模型。微元强度分布函数如下：

式中：中间参数F₀、S₀为对数正态分布参数，其余同前。

根据式(3’)和(7’)可得岩石损伤变量计算式如下：

式中：为标准正态分布函数。

岩石三轴试验中可以测得名义应力σ₁、σ₃(名义应力σ₁即为最大主应力，名义应力σ₃即为正六面体微元条件下的最小主应力，亦即圆柱微元体条件下的围岩压力)，名义应变ε(即各名义应力所对应的应变)，弹性模量E和泊松比μ，由(1’)式和胡克定律可得岩石三维损伤统计本构模型如下：

上述三维损伤本构模型的确定关键在于中间参数F₀和S₀的确定，由(9’)式可得：

由式(12’)可得，

式中σ₁、σ₃、ε₁可由试验确定，因此式(13’)对于一组σ₁、σ₃、ε₁试验值可以得到一个对应的Φ值，由标准正态分布表可查出与之对应的的值，设Y＝lnF，如果以Y为纵坐标，X为横坐标，则Y＝F₀-S₀X，其中，Y可由式(6’)求出，根据多组试验数据便可求出中间参数F₀和S₀。

为了求出中间参数F₀和S₀，本申请中引用现有技术的试验数据(例如参见：曹文贵，赵明华，刘成学.岩石破裂过程的统计损伤模拟研究[J].岩石力学与工程学报，2003，25(2)：184-187)，例如，其中岩石的弹性模量E＝90MPa、泊松比μ＝0.25，内摩擦角φ＝31.3039°，得出X与Y的散点图，见图1(Scatter diagram of X and Y with differentconfining pressure)。采用曲线拟合的方法得到不同围岩压力σ₃下的岩石统计损伤本构模型，其曲线及参数见图2和下表1。

表1三轴试验曲线拟合参数值表

由此可以看出该曲线具有以下特点：(1)该曲线能够反映岩石应力与应变的全过程，尤其是岩石软化性质。(2)随着围岩压力σ₃的增大，岩石的强度逐渐升高，与实际吻合。

2岩石统计损伤本构模型的修正

2.1)中间参数F₀与S₀对模型的影响

为了探讨F₀与S₀对岩石统计损伤本构模型的影响，本申请选用上述实验数据中围岩压力σ₃为6.9MPa时进行分析，分析时保持其它参数不变。F₀对模型的影响见图3，S₀对模型的影响见图4。

从图3和图4可以得出如下结论：中间参数F₀与岩石统计损伤本构模型曲线的峰值强度有关，F₀越大峰值强度越大，但对峰值强度前的弹性变形没有任何影响。中间参数S₀与岩石统计损伤本构模型峰后曲线斜率和峰值强度均有关，S₀越大峰后曲线越缓，表明岩石峰后破裂速度越慢，延性增强，S₀越小峰后曲线越陡，表明岩石峰后破裂速度越快，脆性增强，同时，S₀越小峰值强度越大，S₀越大峰值强度越小，同样S₀对峰值强度前的弹性变形没有任何影响。

2.2)中间参数F₀与S₀与围岩压力σ₃的关系

岩石全应力应变曲线普遍存在峰值强度和峰后曲线随着围压变化而变化(参见：韦立德，徐卫亚，杨春和，杨圣奇.具有统计损伤的岩石弹塑性本构模型研究[J].岩石力学与工程学报，2004，23(12):1971-1975)，而峰前弹性变形部分并不随围压变化的规律，本申请前述中间参数F₀与S₀与峰值强度和峰后曲线有一定的关系，同时也不影响曲线峰前弹性变形部分，显然围岩压力σ₃与中间参数F₀、S₀对岩石全应力应变曲线的影响是相似的，因此，不难看出中间参数F₀、S₀与围岩压力σ₃之间存在某种联系。

采用表1中的数据，以围岩压力σ₃为横坐标，分别以中间参数S₀、F₀为纵坐标做出图5和图6的散点图，进行线性拟合得出如下关系式：

S₀＝0.0014σ₃+0.3153 (14’)

F₀＝0.0136σ₃+5.7449 (15’)

线性拟合相关系数分别为0.9424、0.9768，将(14’)(15’)代入式(9’)可得岩石统计损伤本构模型修正后的表达式：

式中损伤参数A、B、C、H为材料特性参数，每种材料的损伤参数一般不相同，A、B是通过一定围岩压力σ₃条件下的不同F₀的多组数据拟合获得的，C、H是通过一定围岩压力σ₃条件下的不同S₀的多组数据拟合获得的。式(17’)为岩石统计损伤本构模型的通式。

2.3)本构模型的构建

2.3.1)以上利用摩尔库伦准则作为正六面体岩石微元强度破坏准则，和微元强度服从对数正态分布准则，获得了岩石统计损伤本构模型通式(17’)，该表达式(17’)在下面的修正或演变过程中作为起始的表达式(1)：

式(1)中：σ₁、σ₂、σ₃为正六面体微元各面上的主应力，其中σ₁为最大主应力；σ₂为中间主应力；σ₃为最小主应力；ε₁为最大主应力对应的应变；D为岩石损伤变量；Φ为标准正态分布函数；

式(1)中：其中σ₁ ^*为有效最大主应力，σ₃ ^*为有效最小主应力；F为有关有效最大主应力σ₁ ^*、有效最小主应力σ₃ ^*、内摩擦角φ的函数，F表示微元强度；

根据J.Lemaitre应变等价性假说，有效应力与主应力及岩石损伤变量之间具有如下关系：σ_n ^*＝σ_n/[1-D]，n＝1或3，σ_n ^*表示有效应力，σ_n表示主应力；

2.3.2)为了简化试验和方便计算，忽略中间主应力σ₂的影响，将正六面体微元简化为圆柱体微元，此时，σ₃转化为围岩压力，将式(1)变成下式(2)：

在以上式(2)中，F₀是微元强度F满足对数正态分布条件下的均值，它与围岩压力σ₃之间存在线性关系(3)：

F₀＝Aσ₃+B (3)

式(3)中，损伤参数A、B是F₀与σ₃之间的统计概率的拟合线性关系系数；和

在以上式(2)中，S₀是微元强度F满足对数正态分布条件下的方差，它与围岩压力σ₃之间存在线性关系(4)：

S₀＝Cσ₃+H (4)

式(4)中，损伤参数C、H是S₀与σ₃之间的统计概率的拟合线性关系系数；

2.3.3)根据上述关系得出岩石统计损伤本构模型的最终表达式(5)：

3岩石统计损伤本构模型的FLAC^3D实现

3.1)模型的应用过程及参数

根据前述岩石统计损伤本构模型的原理及实验参数，采用FLAC^3D数值计算软件提供的二次开发udm.zip软件包，在Visual Studio 2005编译环境下编写动态链接库文件(.dll)实现本发明的岩石统计损伤本构模型。(参见：蓝航，姚建国，张华兴，徐乃忠.基于FLAC^3D的节理岩体采动损伤本构模型的开发及应用[J].岩石力学与工程学报，2008，27(3):572-579；以及，褚卫江，徐卫亚，杨圣奇，周维垣.基于FLAC^3D岩石粘弹性流变模型的二次开发研究[J].岩土力学，2001，11(27):2005-2010)。

生成本发明的岩石统计损伤本构模型的过程包括：基类的继承、成员变量及成员函数的定义、所述本构模型的注册、所述本构模型与FLAC^3D之间的信息交换、所述本构模型状态指示器生成。将udm.zip软件包导入Visual Studio 2005编译环境下，修改usermohr.h和usermohr.cpp文件，通过编译连接调试生成usermohr.dll二次开发文件，将该文件拷贝至FLAC^3D安装文件下，通过命令config usermohr和model loadusermohr.dll将本构模型加载注册到程序中。(参见：杨文东，张强勇，张建国，贺如平，曾纪全.基于FLAC^3D的改进Burgers蠕变损伤模型的二次开发研究[J].岩土力学，2010，6(31):1956-1964.)。

本发明的岩石统计损伤本构模型的参数包括外部输入参数和内部计算中间参数。外部输入参数：内摩擦角φ、粘结力c、密度ρ、单轴抗拉强度σ_t、损伤参数A、B、C、H以及体积模量K、切变模量G，其中岩石弹性模量E、泊松比μ由体积模量和切变模量计算求得。内部计算中间参数：岩石损伤变量D，取值范围为(0～1)。FLAC^3D数值计算软件的外部参数输入语句为：例prop bulk 42.25e9 shear 37.15e9 coh 0.3e6 fri 38ten 4.5 e6 A 0.4769 B-6.4238 C0.1976 H-6.3824；ini de 3830。

或者，本发明的岩石统计损伤本构模型的参数包括外部输入参数和内部计算中间参数。外部输入参数：岩石弹性模量E、泊松比μ、内摩擦角φ、粘结力c、密度ρ、单轴抗拉强度σ_t、损伤参数A、B、C、H。内部计算中间参数：岩石损伤变量D，取值范围为(0～1)。FLAC^3D数值计算软件的外部参数输入语句为：例prop el 86.18e9 nu 0.16 coh 0.3e6 fri 38 ten4.5 e6 A0.4769 B-6.4238 C 0.1976 H-6.3824；ini de 3830。

3.2)模型的实现

本发明的岩石统计损伤本构模型的应用计算关键在于求解岩石损伤变量D，微元强度F满足对数正态分布关系。岩石受外力作用时，对于每一对有效应力σ₁ ^*、σ₃ ^*，总有一个F值与之对应，由F值可计算出新的D和新的应变ε₁，如此循环，直至有效应力σ₁ ^*、σ₃ ^*满足终止条件(最大不平衡力≤50N)为止，其计算流程见图7。

3.3)模型的验证

为了验证FLAC^3D二次开发出的本发明的岩石统计损伤本构模型的正确性，仍然采用现有技术(例如参见：曹文贵，赵明华，刘成学.岩石破裂过程的统计损伤模拟研究[J].岩石力学与工程学报，2003，25(2)：184-187)的试验数据作为对比数据，在FLAC^3D数值计算软件中创建圆柱形试样的几何模型，直径5cm，高10cm，两端施加压力从0开始逐渐增加直至应变达到设定值6×10^-6，围压设定为6.90Mpa，岩石力学室内试验的参数如前文所述，对比结果见图8。由图8可以看出数值计算模拟曲线与本发明的岩石统计损伤本构模型理论曲线吻合程度很高，且与试验曲线基本吻合，能够表达试验曲线，由此验证了模型的正确性。

3.4)工程应用

本申请的发明人将本申请开发的岩石统计损伤本构模型应用于本申请的发明人之前所研究的工程实例(参见：房智恒，王李管，熊张友.基于Micromine-FLAC^3D耦合技术的金属矿采矿扰动影响分析[J].采矿与安全工程学报，2012，6(29):870-875)中进行对比。该论文已经详细描述了各地质体的建模方法，本申请不再赘述。

(一)工程概况

花岭沟铁矿属于本溪市歪头山矿床的一个矿区，其中四号矿体地质赋存条件复杂，空间上四号矿体北部延伸至露天坑附近，南端深入尾矿库底部，尾矿库初期坝顶轴线端点坐标为(41551800，4593375和41552250,4593375)，四号矿体埋深标高范围为-490m～-60m。采用无底柱分段崩落法采矿。矿区地层从上至下依次为第四系表土层，厚度为0.5～20m、斜长角闪岩地层厚度为0～30m、铁矿层和条带状阳起石英岩层延伸数十米。按图9中圈定的采矿范围计算采矿后地表尾矿坝体的位移等参数。

(二)力学参数

四号矿体围岩主要为斜长角闪岩、石英岩，在进行力学计算时，考虑到岩体的各向异性以及裂隙、节理等弱面和水的影响，对实验室力学试验中获得的力学参数进行折减，根据花岭沟矿区现场地质情况以及相关经验选取折减系数，折减后的力学参数见表2。本申请中重力加速度选为9.81m/s²。

表2矿体及围岩(ore and rock)物理力学参数

(三)计算结果

采用本申请所开发的岩石统计损伤本构模型，对四号矿体进行开采模拟，计算过程中监测点仍然选用尾矿库初期坝顶上的同样位置，图10～12为尾矿库初期坝顶轴线处不同位置监测点随东坐标变化时x、y、z方向发生的位移分量曲线图。坝体x方向位移从西向东依次减小至零后反向并依次增大，最大值为-5.765×10^-4m，坝体y方向位移均为正值，最大位移发生在东坐标52042.5处，值为8.674×10^-4m，坝体z方向位移即沉降量均为负值，最大沉降量发生在坐标52012.5处，值为-9.366×10^-4m。与采用Mohr-Coulomb准则计算的结果较为吻合，因此，此次工程应用也从另一方面验证了本申请所开发的岩石统计损伤本构模型的可信性。

岩石统计损伤本构模型的FLAC^3D实现结果总结：

(1)本申请提出的岩石统计损伤本构模型能够反应岩石全应力应变曲线，与理论和试验数据基本吻合。

(2)本申请结合FLAC^3D数值计算软件进行了岩石统计损伤本构模型的工程应用开发，工程应用也从另一方面验证了本申请所开发的岩石统计损伤本构模型的可信性。同时此次工程应用也为本申请的本构模型在实际工程中的应用奠定了基础。

实施例1

以位于本溪市歪头山矿床的一个矿区“花岭沟铁矿”的矿山地下开采合理圈定保安矿柱范围为例(工程概况如前所述)，介绍一种岩石统计损伤本构模型构建方法的具体应用，包含如下步骤：

(1)进行岩石室内力学试验，获取相关的基础参数：首先，在所应用工程的现场进行岩石取样，分别采取铁矿石、第四系、斜长角闪岩、石英岩各岩土层的样品。通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t；通过室内单轴压缩试验，获得弹性模量E、泊松比μ；通过室内三轴压缩试验，获得岩石的粘结力c，内摩擦角φ；采用天平法或密度仪法测定岩石的密度ρ。

表1为在现场取样后，进行岩石室内试验获得的各岩土层的力学试验参数。

表1各地层实验室所测各岩土层的力学试验参数

(2)根据各岩土层的已有力学试验参数计算岩石体积模量K，按下述公式计算：计算岩石切变模量G，按下述公式计算：

表2为各岩土层计算所得的岩石体积模量K和切变模量G。

表2各岩土层计算所得岩石体积模量K和切变模量G

(3)在实验室做出各岩土层试样在0Mpa、3.45Mpa、6.9Mpa、13.8Mpa、27.6Mpa等不同围岩压力σ₃情况下全应力应变曲线(见图15)，利用线性拟合的方法得出中间参数F₀、S₀与围岩压力σ₃的关系曲线(见图16、17)。根据F₀、S₀与围岩压力σ₃的关系曲线，采用线性拟合方法确定各岩土层的损伤参数A、B、C、H。

表3各岩土层的本构模型损伤参数A、B、C、H表

Claims (13)

1.构建岩石统计损伤本构模型的方法，该方法包括以下步骤：

1)对于岩石样品进行岩石室内力学试验，以获取岩石的弹性模量E、泊松比μ，和内摩擦角φ试验参数；

2)构建岩石统计损伤本构模型，该步骤包括以下的子步骤：

2.1)利用摩尔库伦准则作为正六面体岩石微元强度破坏准则，和微元强度服从对数正态分布准则，建立岩石统计损伤本构模型表达式(1)：

式(1)中：σ₁、σ₂、σ₃为正六面体微元各面上的主应力，其中σ₁为最大主应力；σ₂为中间主应力；σ₃为最小主应力；ε₁为最大主应力对应的应变；D为岩石损伤变量；Φ为标准正态分布函数；岩石损伤变量F₀是微元强度F满足对数正态分布条件下的均值；S₀是微元强度F满足对数正态分布条件下的方差；

式(1)中：其中σ₁ ^*为有效最大主应力，σ₃ ^*为有效最小主应力；F为有关有效最大主应力σ₁ ^*、有效最小主应力σ₃ ^*、内摩擦角φ的函数，F表示微元强度；

有效应力与主应力及岩石损伤变量之间具有如下关系：σ_n ^*＝σ_n/[1-D]，n＝1或3，σ_n ^*表示有效应力，σ_n表示主应力；

2.2)为了简化试验和方便计算，忽略中间主应力σ₂的影响，将正六面体微元简化为圆柱体微元，此时，σ₃转化为围岩压力，将式(1)变成下式(2)：

在以上式(2)中，中间参数F₀是微元强度F满足对数正态分布条件下的均值，它与围岩压力σ₃之间存在线性关系(3)：

F₀＝Aσ₃+B (3)

式(3)中，损伤参数A、B是F₀与σ₃之间的统计概率的拟合线性关系系数；式(3)中的A、B是通过不同围岩压力σ₃条件下的不同F₀的多组数据拟合获得的，一个围岩压力σ₃对应一个F₀，不同的围岩压力σ₃对应的F₀不同；和

在以上式(2)中，中间参数S₀是微元强度F满足对数正态分布条件下的方差，它与围岩压力σ₃之间存在线性关系(4)：

S₀＝Cσ₃+H (4)

式(4)中，损伤参数C、H是S₀与σ₃之间的统计概率的拟合线性关系系数；式(4)中的C、H是通过不同围岩压力σ₃条件下的不同S₀的多组数据拟合获得的，一个围岩压力σ₃对应一个S₀，不同的围岩压力σ₃对应的S₀不同；

2.3)根据上述关系得出岩石统计损伤本构模型的最终表达式(5)：

2.根据权利要求1所述的方法，其中步骤1)是：

针对岩石样品，通过室内单轴压缩试验，获得岩石的弹性模量E、泊松比μ试验参数；通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t试验参数；通过室内三轴压缩试验，获得岩石的粘结力c和内摩擦角φ；以及测定岩石的密度ρ。

3.根据权利要求1或2所述的方法，其中当岩石损伤变量D满足0≤D≤1的条件时，有下式(6)：

4.根据权利要求1或2所述的方法，其中在实验室做出各地质岩层的岩石样品在不同的围岩压力σ₃情况下的全应力应变曲线，利用通过线性拟合的方法得出中间参数F₀、S₀，然后根据式(3)F₀与围岩压力σ₃的关系采用线性拟合方法确定损伤参数A、B，和根据式(4)S₀与围岩压力σ₃的关系采用线性拟合方法确定损伤参数C、H。

5.根据权利要求4所述的方法，其特征在于：围岩压力σ₃为0Mpa、3.45Mpa、6.9Mpa、13.8Mpa或27.6Mpa。

6.岩石统计损伤本构模型在岩土工程中的应用方法，该应用方法包括以下步骤：

(A)针对岩石样品，通过室内单轴压缩试验，获得岩石的弹性模量E、泊松比μ试验参数；通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t试验参数；通过室内三轴压缩试验，获得岩石的粘结力c和内摩擦角φ；以及测定岩石的密度ρ；

(B)将权利要求1-5中任何一项的方法所获得的岩石统计损伤本构模型(5)应用于数值计算软件，获得所述本构模型的二次开发文件；

(C)根据岩石室内力学试验所测得的试验参数，按下述公式计算岩石体积模量K：

同时，计算岩石切变模量G，按下述公式计算：

在数值计算软件中输入岩石室内力学试验所测得的单轴抗拉强度σ_t、粘结力c、密度ρ试验参数；

然后应用步骤(B)中获得的二次开发文件进行工程计算，获得工程参数，和

(D)按照步骤(C)中所获得的工程参数进行工程应用。

7.根据权利要求6所述的应用方法，其中数值计算软件是FLAC^3D、ANSYS或MIDAS。

8.根据权利要求6所述的应用方法，其中数值计算软件是FLAC^3D并且在步骤(B)中获得所述本构模型的二次开发文件usermohr.h和usermohr.cpp，编译生成usermohr.dll文件并在步骤(C)的工程计算时加载。

9.根据权利要求8所述的应用方法，其中所述的工程包括：矿石开采；或水利工程；或市政地下工程。

10.根据权利要求9所述的应用方法，其中所述的矿石开采是地下矿产资源开采。

11.根据权利要求8所述的应用方法，其中所述的工程参数包括：受力岩石在x方向的位移，y方向的位移，z方向的位移。

12.根据权利要求9所述的应用方法，其中所述的工程参数包括：受力岩石在x方向的位移，y方向的位移，z方向的位移。

13.岩石统计损伤本构模型在岩土工程中进行实际工程应用的方法，该方法包括：

I)针对岩石样品，通过室内单轴压缩试验，获得岩石的弹性模量E、泊松比μ试验参数；通过室内拉伸试验，获得岩石的单轴抗拉强度σ_t试验参数；通过室内三轴压缩试验，获得岩石的粘结力c和内摩擦角φ；以及测定岩石的密度ρ；

II)根据权利要求1-5中任何一项所述的方法，构建以上所述的岩石统计损伤本构模型(5)；

III)将以上II)所述的方法构建的所述本构模型(5)应用于数值计算软件，获得所述本构模型的二次开发文件；

IV)然后应用该二次开发文件进行工程计算，获得工程参数；和

V)按照步骤IV)中所获得的工程参数进行实际工程应用。