考虑围压效应和强度脆性跌落的岩石全应力–应变曲线统计损伤模拟方法

张超; 曹文贵; 赵衡; 贺敏

doi:10.11779/CJGE202205017

考虑围压效应和强度脆性跌落的岩石全应力–应变曲线统计损伤模拟方法 English Version

张超^1,,
曹文贵²,
赵衡²,
贺敏³

1.
湖南科技大学岩土工程稳定控制与健康监测湖南省重点实验室, 湖南湘潭 411201
2.
湖南大学岩土工程研究所, 湖南长沙 410082
3.
湖南工业大学土木工程学院, 湖南株洲 412007

基金项目:

国家自然科学基金项目 51909086

湖南省自然科学基金项目 2018JJ3163

湖南省自然科学基金项目 2020JJ5131

湖南省教育厅科学研究项目 21B0489

湖南省教育厅科学研究项目 18C0505

详细信息

作者简介:
张超(1985—)，男，博士，主要从事岩土工程教学与研究工作。E-mail: flyheartzc@hnust.edu.cn

中图分类号: TU452
计量
- 文章访问数: 221
- HTML全文浏览量: 24
- PDF下载量: 107
出版历程
- 收稿日期: 2020-10-19
- 网络出版日期: 2022-09-22
- 刊出日期: 2022-04-30

Statistical damage simulation method for complete stress-strain path of rocks considering confining pressure effect and strength brittle drop

1.
Hunan Provincial Key Laboratory of Geotechnical Engineering for Stability Control and Health Monitoring, Hunan University of Science and Technology, Xiangtan 411201, China
2.
Geotechnical Engineering Institute of Hunan University, Changsha 410082, China
3.
College of Civil Engineering, Hunan University of Technology, Zhuzhou 412007, China

摘要

摘要: 为了建立能够较准确地描述岩石全应力–应变曲线的统计损伤模拟方法，针对目前统计损伤本构模型在模拟围压效应和强度脆性跌落方面存在的严重缺陷，首先，基于三轴压缩条件下脆性岩石基本变形特征，将岩石视为由软物质和硬物质两种细观化材料串联组成，建立宏细观物质变形分析模型；然后，基于瞬时应变分析方法和统计损伤理论，分别提出软物质和硬物质的变形分析方法，从而建立能够考虑围压效应和强度脆性跌落的岩石全应力–应变曲线统计损伤本构模型，给出参数的确定方法；最后，通过模型验证与分析，模型不但能够描述岩石初始宏观变形特征与围压之间的关系，也能够消除岩石轴向峰后理论变形曲线存在的随机性，从而能够较准确地模拟出岩石全应力–应变曲线，表明方法具有一定的合理性和可行性。
- 脆性岩石 /
- 统计损伤 /
- 围压效应 /
- 强度脆性跌落 /
- 全应力–应变曲线
Abstract: To establish a statistical damage simulation method which can accurately describe the full stress-strain curve of rocks, aiming at the serious defects of the existing statistical damage constitutive models in simulating confining pressure effect and strength brittle drop, firstly, based on the basic deformation characteristics of brittle rocks under triaxial compression, the rocks are regarded as a series of two meso-materials: soft and hard materials, and the deformation analysis model for macro-and meso-materials is established. Then, based on the instantaneous strain analysis method and statistical damage theory, the deformation analysis methods for the soft and hard materials are proposed respectively, so as to establish the statistical damage constitutive model for the full stress-strain curve of rocks considering confining pressure effect and strength brittle drop, and the determination methods for the model parameters are given. Finally, through model verification and analysis, the proposed model can describe the relationship between the initial macro-deformation characteristics of rocks and the confining pressure and eliminate the randomness of the axial post-peak theoretical deformation curve, so it can accurately simulate the full stress-strain curve of rocks, which shows that the above methods are reasonable and feasible.
- brittle rock /
- statistical damage /
- confining pressure effect /
- strength brittle drop /
- full stress-strain curve

HTML全文

0. 引言

岩石作为天然地质材料，含有大量多尺度且随机分布的初始空隙，在荷载作用下产生的力学响应极为复杂^[1]。三轴压缩试验已表明，荷载作用初期，初始空隙的闭合将引起初始宏观非线性变形和弹性模量的变化，随着荷载持续作用，脆性岩石将表现出屈服硬化、应变软化和残余强度变形特征。由于脆性岩石强度破坏过程无明显征兆，产生的破坏后果严重，因此，建立能够较准确地预测其全应力–应变过程的模拟方法已成为岩石力学特性研究的长期热点。

自伺服试验机得出岩石全应力–应变曲线以后，基于经典弹塑性理论建立的本构模型已不再适用于模拟岩石变形破坏全过程。近年来，基于应变等效性假设和有效应力原理，将强度统计方法和连续介质损伤理论相结合，提出了统计损伤概念^[2]，在此基础上，众多学者建立了不同表达形式的统计损伤本构模型，并对其合理性和可行性进行了广泛研究，其内容主要包括是否考虑损伤阈值^[3-4]、应力状态^[5-6]和微元体强度随机分布^[7-8]的影响，强度破坏准则^[9-10]和模型参数的确定方法^[11-12]是否合理，以及损伤模型的局限性^[13-14]等，这些问题的研究较好地完善了统计损伤本构理论，然而，该类模型在模拟围压效应和强度脆性跌落方面存在严重缺陷，使其不能够较准确地模拟脆性岩石全应力–应变曲线，这正是本文研究的核心内容。

为此，本文在探讨脆性岩石及其组成物质基本变形特征基础上，将岩石细观化为由软物质和硬物质两部分材料串联组成，提出宏细观物质变形分析模型以及软物质和硬物质的变形分析方法，建立能够考虑围压效应和强度脆性跌落的统计损伤本构模型，给出参数的确定方法，使其能够较准确地预测脆性岩石三轴压缩全应力–应变曲线，以期提高该类模型在岩石工程结构变形分析与计算方面的应用价值。

1. 脆性岩石三轴压缩基本变形特征

MTS815电液伺服岩石力学实验系统广泛使用以来，获得了大量岩石三轴压缩偏应力与轴向应变全过程试验曲线，其典型特征曲线^[15]，见图 1，具有明显阶段性变形特征和围压效应，尽管并非所有岩石都能够反映这些特征，但这不意味着岩石全应力–应变过程客观上就缺少这些特征^[16]。

图 1 大理岩三轴压缩试验曲线^[15]

Figure 1. Triaxial compression test curves of marble^[15]

下载: 全尺寸图片幻灯片

岩石基本变形特征与其围压水平密不可分，随着围压的增大，应力特征值不断增大，峰后强度脆性跌落程度不断减弱，与此同时，发生闭合的初始空隙数量不断增多，初始宏观非线性上凹程度不断降低，弹性模量逐渐增大，当所有初始空隙能够仅在围压作用下完成闭合时，初始宏观变形由非线性转为线性，弹性模量不再随围压增大而变化，而是基本保持不变，称此特征为“围压效应”。孔隙率低的岩石，围压效应微弱，可忽略之，但孔隙率高的岩石，围压效应显著，不能忽略。

岩石先进行仅围压作用，待变形稳定以后再进行恒围压增轴压作用，在此过程岩石与其各组成物质变形协调，其中，初始孔隙率较大或变形力学参数较小的部分物质能够在低应力作用下发生较大变形甚至完成变形，使岩石发生不可恢复的非线性变形，将此类基本变形特征的组成物质细观化为“软物质”；初始孔隙率较小或变形力学参数较大的部分物质在低应力作用下产生弹性微变形，当其所受应力未达到损伤阈值时仅发生弹性变形，超过损伤阈值时发生的变形部分不可恢复，将此类基本变形特征的组成物质细观化为“硬物质”。岩石基本变形特征受其内部细观结构的控制，将其视为由软物质和硬物质两部分细观化材料串联组成，因此需对岩石宏细观物质之间的变形关系进行分析。

2. 岩石宏细观物质变形分析模型

沿某主应力 ${\sigma _i}$ （i=1，2，3）作用方向取出岩石代表性柱状单元（RCU）并进行细观化处理，提出宏细观物质变形分析模型，见。假设RCU、硬物质和软物质在 ${\sigma _i}$ 加载前的长度分别为 $l_i^{}$ ， $l_i^{\text{r}}$ 和 $l_i^{\text{v}}$ ， ${\sigma _i}$ 加载后的变形量分别为 $\Delta l_i^{}$ ， $\Delta l_i^{\text{r}}$ 和 $\Delta l_i^{\text{v}}$ ，软物质初始长度占比（ $l_i^{\text{v}}$ / $l_i^{}$ ）为参数 $\gamma$ 。根据应变分析基本方法，可得

，

${\varepsilon _i} = (1 - \gamma )\varepsilon _i^{\text{r}} + \gamma \varepsilon _i^{\text{v}} \text{，}$

(1)

图 2 宏细观物质变形分析模型

Figure 2. Deformation analysis model for macro- and meso-materials

下载: 全尺寸图片幻灯片

式中， $\varepsilon _i^{\text{v}}$ ， $\varepsilon _i^{\text{r}}$ 为软物质和硬物质的应变。

2.1 软物质变形分析方法

由于软物质持续产生的变形源于持续增大的应力，基于瞬时应变分析方法^[17]可得软物质变形分析表达式，即

，

$\varepsilon _i^{\text{v}} = 1 - \exp ( - {\sigma _i}/{E_{\text{v}}}) \text{，}$

(2)

式中， ${E_{\text{v}}}$ 为软物质弹性模量。 $\varepsilon _i^{\text{v}}$ 随 ${\sigma _i}$ 的增大而增大并表现为负指数变化规律，由于 ${E_{\text{v}}}$ 数量级较小，软物质能够在低应力作用下产生较大变形甚至完成变形。

2.2 硬物质变形分析方法

硬物质变形分析与计算取决于所受应力是否达到损伤阈值，需引入损伤力学理论进行研究。

（1）已有损伤模型局限性分析

损伤模型是建立硬物质变形分析方法的关键。目前，已有损伤模型使统计损伤本构模型在模拟强度脆性跌落方面存在严重缺陷，要弄清楚该问题，有必要对模型构建基本方法进行简要梳理。基于应变等效性假设可得损伤模型的基本表达式^[2]，即

。

${\sigma _i} = {\sigma '_i}(1 - D) 。$

(3)

由于未损伤材料有效应力 ${\sigma '_i}$ 与变形服从广义虎克定律，于是，损伤本构模型基本表达式可为

，

$\{ \sigma \} = [C]\{ \varepsilon \} (1 - D) \text{，}$

(4)

式中， $\{ \sigma \}$ ， $\{ \varepsilon \}$ 为名义应力张量和名义应变张量，[C]为弹性矩阵，D为损伤变量。由于岩石属于非均质材料，微元体强度F各不相同，因此，假设F服从Weibull分布，D为F破坏的累积效应^[8]，

，

$D = \int_0^F {P(F){\text{d}}F} = 1 - \exp [ - {(F/{F_0})^m}] \text{，}$

(5)

式中，P(F)为概率密度函数，m，F₀为分布参数。于是，统计损伤本构模型基本表达式可为

。

$\{ \sigma \} = [C]\{ \varepsilon \} \exp [ - {(F/{F_0})^m}] 。$

(6)

该模型参数仅包括[C]，m和F₀。[C]为基本力学参数，它们不因模型的变化而改变，而是取决于岩石物理性质，可由岩石力学试验获得；m，F₀的确定方法主要有峰值点法^[11]、反演分析法^[12]和线性拟合法^[18]等。峰值点法物理意义明确且计算简单，得到了广泛应用，即

，

${\left. {{\sigma _1}({\varepsilon _1})} \right|_{{\varepsilon _1} = {\varepsilon _{1{\text{c}}}}}} = {\sigma _{1{\text{c}}}} \text{，}$

(7)

，

${\left. {\frac{{\partial {\sigma _1}}}{{\partial {\varepsilon _1}}}} \right|_{{\sigma _1} = {\sigma _{1{\text{c}}}},\;{\varepsilon _1} = {\varepsilon _{1{\text{c}}}}}} = 0 \text{，}$

(8)

式中， ${\sigma _{1{\text{c}}}}$ ， ${\varepsilon _{1{\text{c}}}}$ 为试验曲线峰值处最大主应力及其对应的应变。

联立式（7），（8）可得m和F₀的确定方法，然后将所求参数代入式（6）可得理论模型曲线，这就是构建统计损伤本构模型的基本方法。

由此可知，m和F₀的变化仅取决于 ${\sigma _{1{\text{c}}}}$ 和 ${\varepsilon _{1{\text{c}}}}$ ，特定围压下m和F₀具有唯一性，这就导致特定围压下理论曲线有且仅有一条。为了分析参数变化对理论曲线产生的影响，根据式（4），（5）可知强度脆性跌落仅取决于D，而D的变化仅与m和F₀有关，见图 3，4。随着F的增大，D由0递增至1，D–F曲线呈S型，其中，随着m的增大，D–F曲线形状不变并绕点A作逆时针转动，损伤变化率（∂D/∂F）和强度脆性跌落程度不断增大；随着F₀的增大，D–F曲线形状不变并向F增大方向移动，∂D/∂F和强度脆性跌落程度存在不断减小的趋势。因此，当前理论曲线峰后变化规律本质上取决于m和F₀，然而，特定围压下有且仅有一条的理论曲线无法起到模拟具有随机性的试验曲线峰后变化的作用，究其原因，峰值点法无法对理论曲线峰后变化规律起到约束作用，尽管有文献出现理论曲线与试验曲线峰后吻合良好的实例，但这仅仅表明具有唯一性和随机性的理论曲线与试验曲线之间峰后恰巧接近而已，实际上不具有普适性。

图 3 m和F₀对D的影响

Figure 3. Effects of m and F₀ on D

下载: 全尺寸图片幻灯片

图 4 m和F₀对理论模型曲线的影响

Figure 4. Effects of m and F₀ on theoretical model curves

下载: 全尺寸图片幻灯片

（2）考虑初始损伤的双参数型硬物质损伤模型

由于两类细观化材料属于串联结构，硬物质弹性模量远大于软物质，因此可知软物质变形基本发生在初始空隙压密变形阶段，岩石宏观线弹性及其以后变形过程仅为硬物质产生，而且硬物质弹性模量与岩石相等。根据损伤理论可将硬物质视为由未损伤、初始损伤和后继损伤3种材料并联组成，见，它们所受应力分别为有效应力 ${\sigma '_i}$ 、残余应力 ${\sigma _{\text{r}}}$ 和0，作用面积分别为A_u，A_s和A_d，总作用面积为A_t。

图 5 硬物质损伤模型

Figure 5. Damage model for hard materials

下载: 全尺寸图片幻灯片

根据竖向力学平衡分析方法，可建立考虑初始损伤的硬物质损伤模型，即

，

${\sigma _i} = {\sigma '_i}(1 - {D_{\text{d}}} - {D_{\text{s}}}) + {\sigma _{\text{r}}}{D_{\text{s}}} \text{，}$

(9)

式中， ${D_{\text{d}}}$ ， ${D_{\text{s}}}$ 为硬物质的初始损伤变量和后继损伤变量。由于目前统计损伤本构理论框架已渐趋完善，因此，延续已有损伤模型表达范式与峰值点法优越性，使拟建理论模型能够较准确地模拟试验曲线峰后变化过程，是解决模型局限性的最佳途径。为此，在式（9）基础上引入双参数λ和η，建立考虑初始损伤和强度脆性跌落的双参数型硬物质损伤模型，即

。

${\sigma _i} = {\sigma '_i}(1 - {D_{\text{d}}} - \lambda {D_{\text{s}}}) + \eta {\sigma _{\text{r}}}{D_{\text{s}}} 。$

(10)

由于弹性模量的变化是损伤典型的宏观表现^[3-4]，因此，根据岩石类材料损伤变量与弹性模量之间的关系，可得 ${D_{\text{d}}}$ 的确定方法，即

，

${D_{\text{d}}} = 1 - {E_{\text{d}}}/{E_{\text{u}}} \text{，}$

(11)

式中， ${E_{\text{d}}}$ ， ${E_{\text{u}}}$ 为损伤和无损状态下硬物质的弹性模量。由于硬物质含有大量多尺度且不同形态的微裂纹和孔隙，而微裂纹，尤其是孔隙的完全闭合要极高压应力，很难测得 ${E_{\text{u}}}$ ，因此，基于大量三轴压缩试验资料对 ${E_{\text{d}}}$ 和 ${\sigma _3}$ 之间的关系进行统计表达，即

，

${E_{\text{d}}}/{E_{\text{c}}} = a + b\exp [ - ({\sigma _3}/{\sigma _{\text{c}}})/k] \text{，}$

(12)

式中，a，b，k为试验参数， ${\sigma _{\text{c}}}$ ， ${E_{\text{c}}}$ 为单轴压缩下岩石峰值应力和弹性模量。当 ${\sigma _3}$ 趋于极大值时，基于极限方法可认为空隙近似完全闭合， ${E_{\text{d}}}$ 在数值上趋于 ${E_{\text{u}}}$ ，则 ${E_{\text{u}}}$ 的计算表达式可为

。

${E_{\text{u}}} = a{E_{\text{c}}} 。$

(13)

（3）硬物质统计损伤本构模型

假设未损伤材料应力–应变关系服从广义虎克定律，即

，

${\sigma '_i} = {E_{\text{u}}}\varepsilon _i^{\text{r}} + \mu ({\sigma '_j} + {\sigma '_k}) \text{，}$

(14)

式中，μ为硬物质的泊松比。由于三轴压缩下岩石轴向发生破裂并丧失承载力，侧向约束变形并提高承载力，因此可忽略侧向损伤即认为有效应力 ${\sigma '_j}$ 和 ${\sigma '_k}$ 分别与名义应力 ${\sigma _j}$ 和 ${\sigma _k}$ 相等，于是式（14）可表达为

，

${\sigma '_i} = {E_{\text{u}}}\varepsilon _i^{\text{r}} + \mu ({\sigma _j} + {\sigma _k}) \text{，}$

(15)

于是，将式（15）代入式（10）可得硬物质变形分析表达式，

，

$\varepsilon _i^{\text{r}} = \frac{{{\sigma _i} - \xi {D_{\text{s}}} - \mu ({\sigma _j} + {\sigma _k})(1 - {D_{\text{d}}})}}{{(1 - {D_{\text{t}}}){E_{\text{u}}}}} \text{，}$

(16)

其中，

，

$\xi = \eta {\sigma _{\text{r}}} - \lambda \mu ({\sigma _j} + {\sigma _k}) \text{，}$

(17)

。

${D_{\text{t}}} = {D_{\text{d}}} + \lambda {D_{\text{s}}} 。$

(18)

由于微元体力学特性受应力状态和损伤阈值的双重影响，因此可认为三轴压缩下微元体服从M-C强度破坏准则^[14]，

，

$F = {E_{\text{d}}}\varepsilon _1^{\text{r}} + (2\mu - \alpha ){\sigma _3} - \kappa \text{，}$

(19)

其中，

，

$\alpha = (1 + \sin {\varphi _{\text{y}}})/(1 - \sin {\varphi _{\text{y}}}) \text{，}$

(20)

。

$\kappa = 2{c_{\text{y}}}\cos {\varphi _{\text{y}}}/(1 - \sin {\varphi _{\text{y}}}) 。$

(21)

式中， ${\varphi _{\text{y}}}$ ， ${c_{\text{y}}}$ 为临界损伤状态下硬物质的内摩擦角与黏聚力。当硬物质处于临界损伤状态以下时，微元体未发生破坏，D_s=0；当硬物质处于临界损伤状态以上时，微元体不断发生破坏，D_s由0递增直至硬物质发生完全损伤破坏，D_t=1。于是，D_s可表述为

（ ） （ ） （ ） ，

${D}_{\text{s}}=\left\{\begin{array}{l}0&（F < 0）\\ 1-\mathrm{exp}[-{(}{F/F_0)^m}]&（0\le F < {F}_{\text{s}}）\\ (1-{D}_{\text{d}})/\lambda \text{}&（{F}_{\text{s}}\le F）\end{array} \right.\text{，}$

(22)

式中， ${F_{\text{s}}}$ 为 ${D_{\text{s}}}$ 恰增至(1- ${D_{\text{d}}}$ )/ $\lambda$ 时的微元体强度值。由于硬物质发生损伤可视为连续过程，因此，根据分段函数端点的连续性可得 ${F_{\text{s}}}$ 的计算表达式，

。

${F_{\text{s}}} = {F_0}{\left( { - \ln \frac{{{D_{\text{d}}} + \lambda - 1}}{\lambda }} \right)^{1/m}} 。$

(23)

3. 岩石全应力–应变曲线模拟方法

当硬物质处于弹性状态（F < 0）时，岩石变形由两类细观化物质组成，将式（2）代入式（1）可得

。

$\varepsilon _i^{\text{r}} = \frac{{{\varepsilon _i} - \gamma [1 - \exp ( - {\sigma _i}/{E_{\text{v}}})]}}{{1 - \gamma }} 。$

(24)

当硬物质处于损伤变形阶段（0≤F < ${F_{\text{s}}}$ ）时，软物质变形基本结束，exp( $- {\sigma _i}/{E_{\text{v}}}$ )→0，岩石变形与硬物质变形等同，由式（24）可得

。

$\varepsilon _i^{\text{r}} = ({\varepsilon _i} - \gamma )/(1 - \gamma ) 。$

(25)

当硬物质处于完全损伤变形阶段（ ${F_{\text{s}}}$ ≤F）时，岩石进入残余强度变形阶段。由式（16），（22）可得

。

${\sigma }_{1}=\left\{\begin{array}{l}\frac{{E}_{\text{d}}\{{\varepsilon }_{1}-\gamma [1-\mathrm{exp}(-{\sigma }_{1}/{E}_{\text{v}})]\}}{1-\gamma }+2\mu {\sigma }_{3}(1-{D}_{\text{d}})\text{(}F < 0\text{)}\\ \frac{{E}_{\text{u}}(1-{D}_{\text{t}})({\varepsilon }_{1}-\gamma )}{1-\gamma }+\xi {D}_{\text{s}}+2\mu {\sigma }_{3}(1-{D}_{\text{d}})\text{(}0\le F < {F}_{\text{s}}\text{)}\\ \frac{(1-{D}_{\text{d}})\eta {\sigma }_{\text{r}}}{\lambda }\text{(}{F}_{\text{s}}\le F\text{)}\end{array}\right.。$

(26)

由于三轴压缩下岩石实测偏应力 ${\sigma _{{\text{1t}}}}$ –实测轴向应变 ${\varepsilon _{{\text{1t}}}}$ 试验曲线与 ${\sigma _{\text{1}}}$ - ${\varepsilon _{\text{1}}}$ 理论曲线存在差异，因此需对 ${\sigma _{{\text{1t}}}}$ - ${\varepsilon _{{\text{1t}}}}$ 与 ${\sigma _{\text{1}}}$ - ${\varepsilon _{\text{1}}}$ 之间的关系进行转化：

，

${\sigma _{\text{1}}} = {\sigma _{{\text{1t}}}} + {\sigma _{\text{3}}} \text{，}$

(27)

，

${\varepsilon _{\text{1}}} = {\varepsilon _{{\text{1t}}}} + {\varepsilon _{\text{0}}} \text{，}$

(28)

其中，

。

${\varepsilon _0} = \frac{{(1 - \gamma )[1 - 2\mu (1 - {D_{\text{d}}})]{\sigma _3}}}{{{E_{\text{u}}}(1 - {D_{\text{d}}})}} + \gamma [1 - \exp ( - {\sigma _3}/{E_{\text{v}}})] 。$

(29)

式中， ${\varepsilon _0}$ 为仅围压作用下岩石发生的应变。于是，将式（27），（28）代入式（26）可得三轴压缩下脆性岩石全应力–应变曲线模拟方法，即

${\sigma _{{\text{1t}}}} = \left\{ \begin{gathered} \frac{{{E_{\text{d}}}}}{{1 - \gamma }}\{ {\varepsilon _{{\text{1t}}}} - \gamma '[1 - \exp ( - \frac{{{\sigma _{{\text{1t}}}}}}{{{E_{\text{v}}}}})]\} \;\,\,\;\quad {\text{ (}}F < 0{\text{)}} \hfill \\ \frac{{{E_{\text{u}}}(1 - {D_{\text{t}}})}}{{1 - \gamma }}({\varepsilon _{{\text{1t}}}} - \gamma ') + (\xi - \xi '){D_{\text{s}}}\quad {\text{(}}0 \leqslant F < {F_{\text{s}}}{\text{)}} \hfill \\ \frac{{(1 - {D_{\text{d}}})\eta {\sigma _{\text{r}}}}}{\lambda } - {\sigma _3}\quad \;\quad \;\quad \quad \quad \quad \;{\text{(}}{F_{\text{s}}} \leqslant F{\text{)}} \hfill \\ \end{gathered} \right. \text{，}$

(30)

其中，

$\gamma ' = \gamma \exp ( - {\sigma _3}/{E_{\text{v}}}) \text{，}$

(31)

$\xi ' = \lambda [1 - 2\mu (1 - {D_{\text{d}}}){\sigma _3}/(1 - {D_{\text{d}}})] 。$

(32)

4. 参数的确定方法

由前述分析可知，m和F₀仍然由峰值点法进行确定。根据式（26），联立式（7），（8）可得m和F₀的确定方法，

$m = \frac{{{F_{{\text{sc}}}}(1 - {D_{\text{d}}} - \lambda {D_{{\text{sc}}}})}}{{(1 - {D_{\text{d}}})(1 - {D_{{\text{sc}}}})[\xi (1 - \gamma ) - \lambda {E_{\text{u}}}({\varepsilon _{1{\text{c}}}} - \gamma )]\ln (1 - {D_{{\text{sc}}}})}} \text{，}$

(33)

${F_0} = {F_{{\text{sc}}}}{[ - \ln (1 - {D_{{\text{sc}}}})]^{ - 1/m}} \text{，}$

(34)

其中，

${F_{{\text{sc}}}} = {E_{\text{d}}}({\varepsilon _{1{\text{c}}}} - \gamma )/(1 - \gamma ) + (2\mu - \alpha ){\sigma _3} - \kappa \text{，}$

(35)

${D_{{\text{sc}}}} = \frac{{(1 - \gamma )[{\sigma _{1{\text{c}}}} - 2\mu {\sigma _3}(1 - {D_{\text{d}}})] - {E_{\text{u}}}(1 - {D_{\text{d}}})({\varepsilon _{{\text{1c}}}} - \gamma )}}{{\xi (1 - \gamma ) - \lambda {E_{\text{u}}}({\varepsilon _{1{\text{c}}}} - \gamma )}} 。$

(36)

此外，参数还包括E_v，E_d， $\gamma$ ，D_d， $\lambda$ 和 $\eta$ ，下面简要介绍其确定方法。首先，当轴向应力低于损伤阈值时，初始空隙压密变形阶段试验数据需满足式（30）第一式，可根据拟合分析方法获得E_v；然后，线弹性及其以后变形过程为硬物质所产生，根据式（30）第一式可知线弹性变形反向延长线在 ${\varepsilon _{{\text{1t}}}}$ 轴上的截距为 $\gamma '$ ，将E_v代入式（30）可得 $\gamma$ ；最后，线弹性直线斜率理论表达式为E_d/(1- $\gamma$ )，基于已知 $\gamma$ 和斜率试验值可得E_d，再将E_d代入式（11）可得D_d。 $\lambda$ 和 $\eta$ 的作用未知，待参数分析以后给出其确定方法。

5. 模型验证

为了表明本文模型和方法能够较准确地模拟出围压效应和强度脆性跌落过程，下面基于大理岩试验资料^[15]分别进行验证。

5.1 围压效应验证

由于围压效应发生在岩石初始宏观变形阶段，因此，根据式（30）第一式进行验证。该理论表达式不包含 $\lambda$ 和 $\eta$ ，根据本文模型和参数确定方法可得相关参数值，见表 1，将其代入理论表达式可获得三轴压缩下大理岩初始宏观变形理论曲线，见图 6，与此同时，也可获得此变形阶段内大理岩宏观与其细观化物质理论曲线之间的关系，见图 7。

表 1 围压效应模型参数

Table 1. Model parameters of confining pressure effect

${\sigma _{\text{3}}}$ /MPa	$\gamma$ /10^-2	${E_{\text{d}}}$ /GPa	${E_{\text{v}}}$ /MPa	D_d
1.38	0.54	10.37	4.73	0.42
2.76	0.47	12.42	6.00	0.31
5.52	0.53	13.69	5.26	0.24
6.89	0.64	16.54	4.27	0.08

下载: 导出CSV

| 显示表格

图 6 初始变形阶段理论曲线与试验数据的比较

Figure 6. Comparison between theoretical curves and test data at initial deformation stage

下载: 全尺寸图片幻灯片

图 7 初始变形阶段宏细观物质理论曲线的关系(

${\sigma _{\text{3}}}$ =2.76 MPa)

Figure 7. Relationship between macro- and meso-theoretical curves at initial deformation stage (

${\sigma _{\text{3}}}$ =2.76 MPa)

下载: 全尺寸图片幻灯片

（1）不同围压作用下轴向荷载作用初期变形阶段理论曲线与大理岩初始空隙压密和线弹性等变形阶段试验数据吻合良好，表明本文模型和方法在反映围压效应方面具有较强的适用性。

（2）大理岩轴向宏观变形可由软物质和硬物质两部分细观变形加权构成，权重分别为 $\gamma$ 与1- $\gamma$ ，克服了将岩石等同颗粒骨架且忽略初始空隙闭合给岩石宏观变形造成的不利影响，具有明显的优越性。

（3）软物质变形是引起大理岩初始宏观变形为非线性的主因，其变形的结束即为初始空隙压密变形阶段的终点，由此可确定出空隙闭合应力，表明本文模型和方法能够阐释出岩石宏观与其细观化物质之间的理论变形关系。

5.2 强度脆性跌落验证

由于强度脆性跌落发生在岩石损伤变形破坏阶段，因此根据式（30）第二、三式进行验证。为了表明仅考虑初始损伤的硬物质损伤模型即式（9）的缺陷，以及在此基础上考虑强度脆性跌落的双参数型硬物质损伤模型即式（10）的优越性，下面分别进行探讨。

（1）考虑初始损伤的本构模型

由于分别基于式（9），（10）建立的统计损伤本构模型的区别仅在于理论表达式是否包含 $\lambda$ 和 $\eta$ ，因此，将 $\lambda$ 和 $\eta$ 均赋值为1，两者模型相同，根据本文模型和方法可获得仅考虑初始损伤的大理岩轴向变形理论曲线及其损伤演化曲线，见图 8。

图 8 理论曲线与试验数据的比较(λ=η=1)

Figure 8. Comparison between theoretical curves and test data (λ=η=1)

下载: 全尺寸图片幻灯片

a）该模型能够较准确地模拟出大理岩峰前轴向变形试验曲线，包括初始空隙压密、线弹性和屈服硬化等阶段，但在模拟应变软化和残余强度变形等阶段存在缺陷，尤其在强度脆性跌落过程中理论曲线与试验曲线之间的应变差值较大，无法将两者近似等同。

b）D_s– ${\varepsilon _{{\text{1t}}}}$ 和D_t– ${\varepsilon _{{\text{1t}}}}$ 曲线形状基本呈S型，两者纵向间距为D_d。当轴向应力未达到损伤阈值时，D_s，D_t分别等于为0，D_d并保持不变；当轴向应力达到损伤阈值时，其内部损伤开始累积，D_s，D_t不断增大；当轴向应力达到残余应力时，D_s，D_t分别等于1-D_d，1并保持不变，表明该损伤演化曲线能够较准确地阐释出大理岩阶段性变形特征。

c）随着 ${\sigma _3}$ 的增大，D_d不断减小，损伤阈值和残余应力不断增大，增大 ${\sigma _3}$ 对大理岩发生轴向损伤变形破坏具有有效抑制作用。

（2）考虑强度脆性跌落的本构模型

为了解决峰后理论曲线严重偏离于试验曲线的不足，建立了考虑强度脆性跌落的本构模型，但该模型含有未知 $\lambda$ ， $\eta$ ，无法给出大理岩轴向变形峰后理论曲线，因此将在参数分析基础上探讨其作用并给出其确定方法，然后对该模型进行验证，具体过程如下。

为了探讨 $\lambda$ 的作用，将 $\eta$ 赋值为1，根据本文模型和方法可获得不同 $\lambda$ 下损伤阈值后轴向变形理论曲线以及损伤演化理论曲线，见，可以看出：λ的变化对不同 ${\sigma _3}$ 下峰前变形理论曲线变化规律不起作用，仅对峰后理论曲线变化规律起作用，随着 $\gamma$ 的增大，强度脆性跌落程度和损伤变化率不断增大，残余应力不断减小，损伤演化曲线形状基本保持为S型并绕点A作逆时针转动。

图 9 不同λ下理论曲线变化规律(η=1)

Figure 9. Variation of theoretical curves under different λ (η=1)

下载: 全尺寸图片幻灯片

同理，为了探讨 $\eta$ 的作用，将λ赋值为1，根据本文模型和方法可获得不同 $\eta$ 下损伤阈值后轴向变形理论曲线以及损伤演化理论曲线，见，可以看出：η的变化对不同 ${\sigma _3}$ 下峰前变形理论曲线变化规律不起作用，仅对峰后理论曲线变化规律起作用，随着η的增大，强度脆性跌落程度不断减小，残余应力和损伤变化率不断增大，损伤演化曲线形状基本保持为S型并向轴向应变减小方向移动。

图 10 不同η下理论曲线变化规律(λ=1)

Figure 10. Variation of theoretical curves under different η (λ=1)

下载: 全尺寸图片幻灯片

由此可知， $\lambda$ ， $\eta$ 具有m，F₀对轴向变形理论曲线以及损伤演化理论曲线的等效作用，在m，F₀保持不变的前提下可通过改变 $\lambda$ ， $\eta$ 的方式以达到模拟强度脆性跌落过程的目的，表明参数引入量及其引入方式合理且可行。由于 $\lambda$ ， $\eta$ 对峰后理论曲线变化规律起控制作用，而峰后变形特征表征指标主要指脆性指标和残余应力，在此基础上可给出λ和η的确定方法。

根据岩石脆性特征定量评价方法^[19]，峰后应力降绝对速率B_pa和峰后应力降相对速率B_pr可表示为

${B_{{\text{pa}}}} = ({\sigma _{1{\text{c}}}} - {\sigma _{\text{r}}})/({\varepsilon _{1{\text{r}}}} - {\varepsilon _{{\text{1c}}}}) \text{，}$

(37)

${B_{{\text{pr}}}} = ({\sigma _{1{\text{c}}}} - {\sigma _{\text{r}}})/{\sigma _{{\text{1c}}}} \text{，}$

(38)

式中， ${\varepsilon _{1{\text{r}}}}$ 为轴向应力恰降至 ${\sigma _{\text{r}}}$ 时的应变。于是考虑峰前与峰后变化特征的脆性指标B_n可建立为

${B_{\text{n}}} = {B_{{\text{pr}}}} \cdot \frac{2}{{\text{π }}}\arctan {B_{{\text{pa}}}} 。$

(39)

显然，B_n在[0，1]变化，脆性程度随B_n增大而增大，于是根据式（37）~（39）可得

${\varepsilon _{{\text{1r}}}} = {\varepsilon _{{\text{1c}}}} + ({\sigma _{{\text{1c}}}} - {\sigma _{\text{r}}}){\tan ^{ - 1}}\frac{{{\text{π }}{B_{\text{n}}}{\sigma _{{\text{1c}}}}}}{{2({\sigma _{{\text{1c}}}} - {\sigma _{\text{r}}})}} 。$

(40)

由式（1），（19）可得偏应力恰跌至σ_r时F_s，

${F_{\text{s}}} = {E_{\text{d}}}({\varepsilon _{{\text{1r}}}} - \gamma )/(1 - \gamma ) + (2\mu - \alpha ){\sigma _3} - \kappa 。$

(41)

于是结合式（23），（41）可得λ的确定方法，

$\lambda =\frac{1-{D}_{\text{d}}}{1-\mathrm{exp}\{{}-{[({E}_{\text{d}}({\varepsilon }_{\text{1r}}-\gamma )/(1-\gamma )+(2\mu -\alpha ){\sigma }_{3}-\kappa )/{F}_{0}]}^{m}\text{}\}}。$

(42)

由三轴压缩条件下岩石基本变形特征分析可得

$\mathop {\lim }\limits_{{\varepsilon _1} \to {\varepsilon _{1{\text{r}}}}} {\sigma _1} = {\sigma _{\text{r}}} 。$

(43)

于是将式（26）第三式代入式（43）可得η的确定方法，即

$\eta = \lambda /(1 - {D_{\text{d}}}) 。$

(44)

根据λ和η的确定方法可获取相应的参数，见，然后，将其代入式（30）可得不同 ${\sigma _3}$ 作用下大理岩全应力–应变理论曲线，见图 11，可以看出：本文模型和方法不仅能够较准确地反映围压效应，也能够较准确地模拟出强度脆性跌落过程，消除了已有模型峰后理论曲线变化的随机性，较好地缓解了峰后变形理论曲线与试验数据的严重偏离程度，从而表明本文模型和方法具有一定的合理性与可行性。

表 2 不同

${\sigma _3}$ 下模型参数λ和η

Table 2. Model parameters λ and η under different σ₃

${\sigma _3}$ /MPa	λ	η
1.38	0.58	1.002
2.76	0.69	0.996
5.52	0.76	0.995
6.89	0.92	0.997

下载: 导出CSV

| 显示表格

图 11 全应力–应变理论曲线与大理岩试验数据的比较

Figure 11. Comparison between theoretical curves of complete stress-strain and test data of marble

下载: 全尺寸图片幻灯片

6. 结论

（1）岩石可视为由软物质和硬物质两种细观化材料的串联结构，在此基础上建立了宏细观物质变形分析模型，其中软物质服从负指数变形规律，硬物质可根据统计损伤理论构建本构模型。

（2）m，F₀能够使模型较准确地模拟轴向峰前变形过程，但在模拟轴向峰后变形过程时带有强烈的随机性；λ，η对轴向峰前变形理论曲线变化规律不起作用，但对峰后变形理论曲线变化规律起控制作用，具有m，F₀对轴向变形理论曲线以及损伤演化理论曲线的等效作用。

（3）建立了能够考虑围压效应和强度脆性跌落的脆性岩石全应力–应变曲线统计损伤模拟方法，给出了参数的确定方法，通过模型验证与分析表明本文模型和方法具有一定的合理性与可行性，完善了岩石统计损伤本构理论。

图 1 大理岩三轴压缩试验曲线^[15]

Figure 1. Triaxial compression test curves of marble^[15]

下载: 全尺寸图片幻灯片

图 2 宏细观物质变形分析模型

Figure 2. Deformation analysis model for macro- and meso-materials

下载: 全尺寸图片幻灯片

图 3 m和F₀对D的影响

Figure 3. Effects of m and F₀ on D

下载: 全尺寸图片幻灯片

图 4 m和F₀对理论模型曲线的影响

Figure 4. Effects of m and F₀ on theoretical model curves

下载: 全尺寸图片幻灯片

图 5 硬物质损伤模型

Figure 5. Damage model for hard materials

下载: 全尺寸图片幻灯片

图 6 初始变形阶段理论曲线与试验数据的比较

Figure 6. Comparison between theoretical curves and test data at initial deformation stage

下载: 全尺寸图片幻灯片

初始变形阶段宏细观物质理论曲线的关系( ${\sigma _{\text{3}}}$ =2.76 MPa)

Relationship between macro- and meso-theoretical curves at initial deformation stage ( ${\sigma _{\text{3}}}$ =2.76 MPa)

下载: 全尺寸图片幻灯片

图 8 理论曲线与试验数据的比较(λ=η=1)

Figure 8. Comparison between theoretical curves and test data (λ=η=1)

下载: 全尺寸图片幻灯片

图 9 不同λ下理论曲线变化规律(η=1)

Figure 9. Variation of theoretical curves under different λ (η=1)

下载: 全尺寸图片幻灯片

图 10 不同η下理论曲线变化规律(λ=1)

Figure 10. Variation of theoretical curves under different η (λ=1)

下载: 全尺寸图片幻灯片

图 11 全应力–应变理论曲线与大理岩试验数据的比较

Figure 11. Comparison between theoretical curves of complete stress-strain and test data of marble

下载: 全尺寸图片幻灯片

表 1 围压效应模型参数

Table 1 Model parameters of confining pressure effect

${\sigma _{\text{3}}}$ /MPa	$\gamma$ /10^-2	${E_{\text{d}}}$ /GPa	${E_{\text{v}}}$ /MPa	D_d
1.38	0.54	10.37	4.73	0.42
2.76	0.47	12.42	6.00	0.31
5.52	0.53	13.69	5.26	0.24
6.89	0.64	16.54	4.27	0.08

下载: 导出CSV

不同 ${\sigma _3}$ 下模型参数λ和η

Table 2 Model parameters λ and η under different σ₃

${\sigma _3}$ /MPa	λ	η
1.38	0.58	1.002
2.76	0.69	0.996
5.52	0.76	0.995
6.89	0.92	0.997

下载: 导出CSV

参考文献(19)

[1]	陈颙, 黄庭芳, 刘恩儒. 岩石物理学[M]. 合肥: 中国科学技术大学出版社, 2009. CHEN Yong, HUANG Ting-fang, LIU En-ru. Rock Physics[M]. Hefei: Press of University of Science and Technology of China, 2009. (in Chinese)
[2]	KRAJCINOVIC D, SILVA M A G. Statistical aspects of the continuous damage theory[J]. International Journal of Solids and Structures, 1982, 18(7): 551–562. doi: 10.1016/0020-7683(82)90039-7
[3]	ZHAO H, ZHANG C, CAO W G, et al. Statistical mesodamage model for quasi-brittle rocks to account for damage tolerance principle[J]. Environmental Earth Sciences, 2016, 75(10): 1–12.
[4]	JIANG H B, LI K N, HOU X B. Statistical damage model of rocks reflecting strain softening considering the influences of both damage threshold and residual strength[J]. Arabian Journal of Geosciences, 2020, 13(7): 1–8.
[5]	曹文贵, 戴笠, 张超. 深部岩石统计损伤本构模型研究[J]. 水文地质工程地质, 2016, 43(4): 60–65. https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG201604011.htm CAO Wen-gui, DAI Li, ZHANG Chao. A study of statistical damage constitutive models for deep earth rocks[J]. Hydrogeology & Engineering Geology, 2016, 43(4): 60–65. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG201604011.htm
[6]	张慧梅, 谢祥妙, 彭川, 等. 三向应力状态下冻融岩石损伤本构模型[J]. 岩土工程学报, 2017, 39(8): 1444–1452. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16984.shtml ZHANG Hui-mei, XIE Xiang-miao, PENG Chuan, et al. Constitutive model for damage of freeze-thaw rock under three-dimensional stress[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(8): 1444–1452. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16984.shtml
[7]	ZHU Z N, TIAN H, WANG R, et al. Statistical thermal damage constitutive model of rocks based on Weibull distribution[J]. Arabian Journal of Geosciences, 2021, 14(6): 1-14.
[8]	张明, 王菲, 杨强. 基于三轴压缩试验的岩石统计损伤本构模型[J]. 岩土工程学报, 2013, 35(11): 1965–1971. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract15325.shtml ZHANG Ming, WANG Fei, YANG Qiang. Statistical damage constitutive model for rocks based on triaxial compression tests[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(11): 1965–1971. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract15325.shtml
[9]	DENG J, GU D S. On a statistical damage constitutive model for rock materials[J]. Computers & Geosciences, 2011, 37(2): 122–128. https://www.sciencedirect.com/science/article/pii/S0098300410002700
[10]	金俊超, 佘成学, 尚朋阳. 基于Hoek-Brown准则的岩石应变软化模型研究[J]. 岩土力学, 2020, 41(3): 939–951. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202003025.htm JIN Jun-chao, SHE Cheng-xue, SHANG Peng-yang. A strain-softening model of rock based on Hoek-Brown criterion[J]. Rock and Soil Mechanics, 2020, 41(3): 939–951. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202003025.htm
[11]	ZHAO H, SHI C J, ZHAO M H, et al. Statistical damage constitutive model for rocks considering residual strength[J]. International Journal of Geomechanics, 2017, 17(1): 04016033. doi: 10.1061/(ASCE)GM.1943-5622.0000680
[12]	刘齐建, 杨林德, 曹文贵. 岩石统计损伤本构模型及其参数反演[J]. 岩石力学与工程学报, 2005, 24(4): 616–621. doi: 10.3321/j.issn:1000-6915.2005.04.012 LIU Qi-jian, YANG Lin-de, CAO Wen-gui. Statistical damage constitutive model for rock and back analysis of its parameters[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(4): 616–621. (in Chinese) doi: 10.3321/j.issn:1000-6915.2005.04.012
[13]	温韬, 唐辉明, 马俊伟, 等. 考虑初始损伤和残余强度的岩石变形过程模拟[J]. 地球科学, 2019, 44(2): 652–663. https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX201902026.htm WEN Tao, TANG Hui-ming, MA Jun-wei, et al. Deformation simulation for rock in consideration of initial damage and residual strength[J]. Earth Science, 2019, 44(2): 652–663. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX201902026.htm
[14]	CAO W G, TAN X, ZHANG C, et al. Constitutive model to simulate full deformation and failure process for rocks considering initial compression and residual strength behaviors[J]. Canadian Geotechnical Journal, 2019, 56(5): 649–661. doi: 10.1139/cgj-2018-0178
[15]	ROSENGREN K J, JAEGER J C. The mechanical properties of an interlocked low-porosity aggregate[J]. Géotechnique, 1968, 18(3): 317–326. doi: 10.1680/geot.1968.18.3.317
[16]	徐志英. 岩石力学[M]. 北京: 中国水利水电出版社, 1993. XU Zhi-ying. Rock Mechanics[M]. Beijing: China Water Power Press, 1993. (in Chinese)
[17]	ZHAO Y, LIU H H. An elastic stress–strain relationship for porous rock under anisotropic stress conditions[J]. Rock Mechanics and Rock Engineering, 2012, 45(3): 389–399. doi: 10.1007/s00603-011-0193-y
[18]	LI X, CAO W G, SU Y H. A statistical damage constitutive model for softening behavior of rocks[J]. Engineering Geology, 2012, 143/144: 1–17. doi: 10.1016/j.enggeo.2012.05.005
[19]	TARASOV B, POTVIN Y. Universal criteria for rock brittleness estimation under triaxial compression[J]. International Journal of Rock Mechanics and Mining Sciences, 2013, 59: 57–69. doi: 10.1016/j.ijrmms.2012.12.011

施引文献(16)

期刊类型引用(7)

1.	刘婷婷，阳润超，丁鹿阳，李新平，曾乐乐. 含非贯通节理花岗岩的力学特性与细观起裂机制研究. 岩石力学与工程学报. 2023(05): 1070-1082 . 百度学术
2.	罗亚琼，张超，李洁，卢霖，任中俊. 土石混合体大型三轴剪切变形特性及模拟方法. 水利水电技术(中英文). 2023(05): 168-176 . 百度学术
3.	蒋培林，丁书学，南华，乔静，刘江，田永超，李佳赛，李树杰. 中间主应力影响下加锚试样应力-应变曲线特征参数分析. 科学技术与工程. 2023(18): 7725-7733 . 百度学术
4.	戴俊，杨清清，张敏，蒋昊轩，王林涛，任一平. 微波照射下花岗岩单轴压缩损伤本构模型. 科学技术与工程. 2023(19): 8350-8357 . 百度学术
5.	刘兆阳，于开治，熊文彬. 围压作用下的泡沫铝特性研究. 科技与创新. 2023(15): 105-107+111 . 百度学术
6.	辛子朋，柴肇云，孙浩程，李天宇，刘新雨，段碧英. 砂质泥岩峰后破裂承载特征与块体分布规律研究. 岩土力学. 2023(08): 2369-2380 . 百度学术
7.	宋清蔚，朱万成，徐晓冬，张鹏海，陈承桢，刘溪鸽. 岩石破坏过程实时监测预警软件系统——搭建与初步应用. 金属矿山. 2022(10): 155-164 . 百度学术