干燥与饱和岩石裂隙受压变形与声发射特性研究

李博; 叶鹏进; 黄林; 王丁; 赵程; 邹良超

doi:10.11779/CJGE202112011

干燥与饱和岩石裂隙受压变形与声发射特性研究 English Version

李博^{1, 2,},
叶鹏进¹,
黄林²,
王丁^1, ,,
赵程²,
邹良超³

1.
绍兴文理学院岩石力学与地质灾害浙江省重点实验室,浙江　绍兴　312000
2.
同济大学地下建筑与工程系,上海　200092
3.
瑞典皇家理工学院可持续发展与环境工程系,斯德哥尔摩　SE-100　44,瑞典

基金项目:

浙江省自然科学基金项目 LR19E090001

国家自然科学基金项目 42077252

国家自然科学基金项目 41804127

详细信息

作者简介:
李博(1981— ),男,博士,教授,主要从事岩体水力耦合特性方面的研究工作。E-mail：libo@usx.edu.cn

通讯作者:
王丁, E-mail：wdgeophysics@hotmail.com

中图分类号: TU443
计量
- 文章访问数: 262
- HTML全文浏览量: 32
- PDF下载量: 223
出版历程
- 收稿日期: 2021-03-25
- 网络出版日期: 2022-11-30
- 刊出日期: 2021-11-30

Deformation and acoustic emission characteristics of dry and saturated rock fractures

1.
Key Laboratory of Rock Mechanics and Geohazards of Zhejiang Province, Shaoxing University, Shaoxing 312000, China
2.
Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China
3.
Royal Institute of Technology, Department of Sustainable Development and Environmental Engineering, Stockholm SE-100 44, Sweden

摘要

摘要: 天然岩体处于各类复杂的地质环境中,定量描述岩石裂隙在各种条件下的变形破坏行为是开展各类与岩石裂隙力学特性相关研究的基础。通过单轴压缩试验和弹塑性接触数值计算研究了两种具有不同岩性的岩石裂隙在干燥与饱和状态下的受压变形破坏行为,并对试验过程中的声发射信号进行了采集分析。结果表明：试验与数值计算得到的法向应力-位移曲线和塑性破坏区域皆基本吻合,验证了弹塑性接触计算方法的可靠性；塑性变形均值增量u_p随法向应力σ_n的增长呈现变化速率降低的非线性减小关系,利用力学及几何参数建立的拟合公式可较好地预测该关系；声发射源位置与试验及数值模拟得到的损伤破坏区域基本匹配,干燥状态下两种岩石声发射振铃计数、累计计数均高于饱和状态；声发射振铃计数与塑性变形均值增量的变化趋势基本一致。研究结果定量地揭示了岩石基本力学参数和裂隙表面特征对其受压变形破坏行为的控制作用。
- 岩石裂隙 /
- 变形 /
- 弹塑性接触 /
- 声发射 /
- 法向应力
Abstract: The natural rock masses are situated in various complex geological conditions. Quantitative description of the deformation and failure behaviors of rock fractures under these conditions is of fundamental importance for the studies related to their mechanical behaviors. In this study, the unconfined compression tests and the elastic-plastic contact numerical simulations are implemented to study the compressive deformation and failure behavior of two kinds of rock fractures under dry and saturated conditions, together with acoustic emission detection and analysis. The results show that the normal stress-displacement curves and the plastic failure areas obtained from the experiment and the numerical simulation agree well with each other, which verifies the reliability of the contact method. The mean increment of plastic deformation decreases nonlinearly with the increasing normal stress with a decreasing rate, and a fitting formula is established using mechanical and geometric parameters. The position of acoustic emission (AE) sources matches with the damage area obtained from the experiment and the numerical simulation. Both the AE ringing count and the cumulative count are higher in dry rocks than those in saturated rocks. The AE ringing count and the mean increment of plastic deformation follow an identical changing trend. These results reveal the controlling role played by the fundamental mechanical parameters and geometric properties in the deformation and failure behaviors of rock fractures.
- rock fracture /
- deformation /
- elastic-plastic contact /
- acoustic emission /
- normal stress

HTML全文

0. 引言

传统弹塑性模型由于隐含了共轴性假设,无法模拟土体变形的非共轴特性^[1-3]。针对这一问题,国内外学者从不同角度入手,提出了多种形式的非共轴模型。一些研究将土体的非共轴行为归因于其细观结构的各向异性,并基于土体细观结构的演化规律,采用多尺度的方法建立非共轴模型^[4-9]。另一方面,通过直接采用数学公式描述宏观试验现象来构建土体本构模型也是一种非常有效的方法,因此,一些研究者也采用引入非共轴塑性理论的方式建立土体本构模型^[10-13]。Rudnicki等^[14]在研究裂隙岩体分叉行为时,为了提高理论预测的精确性,引入了非共轴应力率引发的塑性变形,由此提出了原始的非共轴理论。但这一早期理论是在二维应力空间内所建立,并不适用于三维化的模型。后来,Hashiguchi等^[15]通过重新定义非共轴应力率表达式,建立了三维形式的非共轴理论。钱建固等^[16]采用Gram-Schimidt张量正交化理论,将第三应力不变量的影响加入到了非共轴理论中。基于钱建固等^[16]的研究,李学丰等^[17],陈洲泉等^[18]分别采用非共轴模型对空心圆柱扭剪试验中砂土的力学响应进行了模拟。刘元雪等^[13]研究了三维条件下主应力轴旋转加载过程中应力率的数学描述问题,并在此基础上建立了考虑非共轴塑性变形的本构模型。这些研究表明,通过将非共轴理论引入已有的共轴模型中,可以一定程度上再现单剪试验与空心圆柱试验中的非共轴现象。

分析已有研究结果可以发现,传统非共轴理论均采用从总应力率减去所谓共轴应力率的方式定义非共轴应力率,不同理论中定义的非共轴应力率之间的区别还没有被明确地探讨过。最近,Li等^[19]从纯数学角度对两个张量之间的非共轴性进行了研究,但与传统非共轴理论同样地采用从总量减去所谓共轴量的方式定义非共轴量,这种定义方式一定程度上给非共轴应力率的理解带来了困难。另一方面,虽然传统非共轴理论在模拟土体变形方面取得了一些成功,但非共轴塑性应变率与非共轴应力率之间通常采用线性假设并且二者始终同向,这与土体的实际变形特性不符。针对上述问题,本文在总结已有理论基础上,通过建立一种针对对称二阶张量的分解法将应力率进行了正交分解,由此可从另一角度对非共轴应力率进行表述。进一步地,在Zienkiewicz^[20]建立的广义塑性力学基础上,针对诱发非共轴塑性变形的每一分量应力率建立相应的非线性加载机制,提出了广义形式的非共轴理论。然后,基于该理论建立了一个非共轴模型并采用两种数值试验对模型的合理性进行了检验,以期广义非共轴理论能够为土体非共轴模型的构建提供更为广阔的理论基础。

1. 传统的非共轴理论

首先需要说明的是,本文对于非共轴理论的介绍均是基于小变形假设。因此,对于采用了Jaumann应力率的非共轴理论,本文在介绍时均采用普通应力率代替。

在Rudnicki等^[14]提出的早期非共轴理论中,总应变率 $\dot{ε}$ 被分解为弹性应变率 ${\dot{ε}}^{e}$ ,共轴塑性应变率 ${\dot{ε}}^{cp}$ 与非共轴塑性应变率 ${\dot{ε}}^{np}$ 三者之和：

$\dot{ε} = {\dot{ε}}^{e} + {\dot{ε}}^{cp} + {\dot{ε}}^{np}$ ,

(1)

Rudnicki等^[14]认为非共轴塑性应变率 ${\dot{ε}}^{np}$ 与静水压力部分无关,其方向与 $s$ 正交,由此将 ${\dot{ε}}^{np}$ 表示为

${\dot{s}}_{nc} = \dot{s} - \frac{s : \dot{s}}{s : s} s$ ,

(2a)

${\dot{ε}}^{np} = \frac{1}{H_{t}} {\dot{s}}_{nc}$ ,

(2b)

式中, $H_{t}$ 为非共轴塑性模量, $s$ 为偏应力, $\dot{s}$ 为偏应力率, ${\dot{s}}_{nc}$ 为非共轴应力率。然而式（2a）,（2b）由于没有考虑第三应力不变量的影响,并不适用于一般化的三维模型。

Hashiguchi等^[15]通过假设非共轴塑性应变率与屈服面偏切线方向的应力率同向且二者之间为线性关系,将非共轴塑性应变率 ${\dot{ε}}^{np}$ 的表达式修正为

${\dot{ε}}^{np} = T E^{- 1} : {\dot{s}}_{nc}$ 。

(3)

式中 $T$ 为一个无量纲的标量系数； $E$ 为弹性刚度张量； ${\dot{s}}_{nc}$ 为

${\dot{s}}_{nc} = Ι_{t} : \dot{σ}$ ,

(4a)

式中

$Ι_{t} = I - \frac{1}{3} 1 \otimes 1 - \frac{{f^{'}}_{σ}}{‖ {f^{'}}_{σ} ‖} \otimes \frac{{f^{'}}_{σ}}{‖ {f^{'}}_{σ} ‖}$ ,

(4b)

${f^{'}}_{σ} = f_{σ} - \frac{1}{3} tr(f_{σ})1$ ,

(4c)

其中, $σ$ 为应力, $\dot{σ}$ 为应力率, $Ι_{t}$ 为偏切线投影张量, $I$ 为对称转换张量,其分量形式为 $I_{i j k l} =$ $(δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) / 2$ , $δ_{i j}$ 为Kronecker符号, ${f^{'}}_{σ}$ 为屈服函数梯度的偏量部分, $f_{σ} = \partial f / \partial σ$ 表示屈服函数的梯度（本文中 $f_{a} = \partial f / \partial a$ ,下文不再说明）, $1$ 为二阶单位张量,其分量为 $δ_{i j}$ 。

钱建固等^[16]为考虑第三应力不变量的影响提出了另一种三维化的非共轴理论。在该理论中,非共轴塑性应变率的表达式与式（2b）相同,但非共轴应力率 ${\dot{s}}_{nc}$ 的表达式修改为

${\dot{s}}_{nc} = \dot{s} - \frac{s : \dot{s}}{s : s} s - \frac{S : \dot{s}}{S : S} S$ ,

(5a)

$S = s^{2} - \frac{2}{3} J_{2} 1 - \frac{3}{2} \frac{J_{3}}{J_{2}} s$ 。

(5b)

式中 $S$ 为二阶张量,分别正交于 $1$ 和 $s$ ； $J_{2} = tr(s^{2})/2$ 为偏应力第二不变量； $J_{3} = tr(s^{3})/3$ 为偏应力第三不变量。

分析上述非共轴理论可知,非共轴应力率的定义均采用从总应力率中减去所谓共轴应力率的方式,这与Li等^[19]对于两个张量之间非共轴性的定义方式相类似。但这种定义方式不能说明非共轴应力率是由多项正交应力率分量所构成这一事实。为了对此问题进行更加清楚的解释,有必要对应力率的分量表达形式进行研究。

2. 非共轴应力率的正交分解

令3个单位基矢量 ${\hat{e}}_{1}$ , ${\hat{e}}_{2}$ , ${\hat{e}}_{3}$ 分别代表3个主应力 $σ_{1}$ , $σ_{2}$ , $σ_{3}$ 所在的方向。则在这一坐标系中,应力张量 $σ$ 可为

$σ = σ_{1} {\hat{e}}_{1} \otimes {\hat{e}}_{1} + σ_{2} {\hat{e}}_{2} \otimes {\hat{e}}_{2} + σ_{3} {\hat{e}}_{3} \otimes {\hat{e}}_{3}$ 。

(6)

3个主应力 $σ_{1}$ , $σ_{2}$ , $σ_{3}$ 分别为

$σ_{1} = p + \frac{2}{3} q \sin (θ + \frac{2}{3} π)$ ,

(7a)

$σ_{2} = p + \frac{2}{3} q \sin (θ)$ ,

(7b)

$σ_{3} = p + \frac{2}{3} q \sin (θ - \frac{2}{3} π)$ 。

(7c)

式中 $p$ 为平均应力, $p = tr(σ)/3$ ； $q$ 为广义剪应力, $q = \sqrt{3 s : s / 2}$ ； $θ$ 为罗德角, $θ = - [\sin^{- 1} (27 J_{3} / 2 q^{3})] / 3$ 。

现假设 $σ$ 经过微小的增量 $Δ σ$ 变化为 $σ^{*}$ ,则可以采用如下取极限的方式定义一个与 $σ$ 相关的方向量：

$n = \lim_{σ^{*} \to σ} \frac{σ^{*} - σ}{‖ σ^{*} - σ ‖}$ 。

(8)

进一步地,假设 $Δ σ$ 由平均应力变化 $Δ p$ 引发,令 $Δ p \to 0^{+}$ ,则根据式（8）可得到方向 ${\hat{n}}_{p}$ 为

${\hat{n}}_{p} = \frac{\sqrt{3}}{3} {\hat{e}}_{1} \otimes {\hat{e}}_{1} + \frac{\sqrt{3}}{3} {\hat{e}}_{2} \otimes {\hat{e}}_{2} + \frac{\sqrt{3}}{3} {\hat{e}}_{3} \otimes {\hat{e}}_{3}$ 。

(9)

同理,当 $Δ σ$ 分别是由广义剪应力 $q$ 和罗德角 $θ$ 的独立变化所引发时,令 $Δ q \to 0^{+}$ 和 $Δ θ \to 0^{+}$ 可分别得到方向 ${\hat{n}}_{q}$ 和 ${\hat{n}}_{θ}$ 为

$\begin{matrix} {\hat{n}}_{q} = \frac{\sqrt{6}}{3} \sin (θ + \frac{2}{3} π) {\hat{e}}_{1} \otimes {\hat{e}}_{1} + \\ \frac{\sqrt{6}}{3} \sin (θ) {\hat{e}}_{2} \otimes {\hat{e}}_{2} + \frac{\sqrt{6}}{3} \sin (θ - \frac{2}{3} π) {\hat{e}}_{3} \otimes {\hat{e}}_{3} \end{matrix}$ ,

(10)

$\begin{matrix} {\hat{n}}_{θ} = \frac{\sqrt{6}}{3} \cos (θ + \frac{2}{3} π) {\hat{e}}_{1} \otimes {\hat{e}}_{1} + \\ \frac{\sqrt{6}}{3} \cos (θ) {\hat{e}}_{2} \otimes {\hat{e}}_{2} + \frac{\sqrt{6}}{3} \cos (θ - \frac{2}{3} π) {\hat{e}}_{3} \otimes {\hat{e}}_{3} \end{matrix}$ 。

(11)

进一步地,假设 $σ^{*}$ 是由 $σ$ 经过绕 $σ_{1}$ 轴旋转 $α$ 角所得,则 $σ^{*}$ 可表示为 $σ^{*} = Q \cdot σ \cdot Q^{T}$ ,相应的正交变换张量 $Q$ 为

$Q_{i j} = [\begin{matrix} 1 & 0 & 0 \\ 0 & \cos (α) & \sin (α) \\ 0 & - \sin (α) & \cos (α) \end{matrix}]$ 。

(12)

将 $σ^{*} = Q \cdot σ \cdot Q^{T}$ 和式（12）代入式（8）,令 $Δ α \to 0^{+}$ 可得方向 ${\hat{n}}_{r1}$ 为

${\hat{n}}_{r1} = \frac{\sqrt{2}}{2} {\hat{e}}_{2} \otimes {\hat{e}}_{3} + \frac{\sqrt{2}}{2} {\hat{e}}_{3} \otimes {\hat{e}}_{2}$ 。

(13)

假设 $σ^{*}$ 分别由 $σ$ 绕 $σ_{2}$ 轴和 $σ_{3}$ 轴旋转所得,经过与前述相同的运算,可以得到方向 ${\hat{n}}_{r2}$ , ${\hat{n}}_{r3}$ 为

${\hat{n}}_{r2} = \frac{\sqrt{2}}{2} {\hat{e}}_{1} \otimes {\hat{e}}_{3} + \frac{\sqrt{2}}{2} {\hat{e}}_{3} \otimes {\hat{e}}_{1}$ ,

(14)

${\hat{n}}_{r3} = \frac{\sqrt{2}}{2} {\hat{e}}_{1} \otimes {\hat{e}}_{2} + \frac{\sqrt{2}}{2} {\hat{e}}_{2} \otimes {\hat{e}}_{1}$ ,

(15)

可以证明,这样得到的6个方向 ${\hat{n}}_{p}$ , ${\hat{n}}_{q}$ , ${\hat{n}}_{θ}$ , ${\hat{n}}_{r1}$ , ${\hat{n}}_{r2}$ , ${\hat{n}}_{r3}$ 是相互正交的（本文中任意两个方向量 $a$ 与 $b$ 正交意味着 $a : b = 0$ ）并且6个方向量与对称转换张量 $I$ 存在如下关系：

$\begin{matrix} I = {\hat{n}}_{p} \otimes {\hat{n}}_{p} + {\hat{n}}_{q} \otimes {\hat{n}}_{q} + {\hat{n}}_{θ} \otimes {\hat{n}}_{θ} + \\ {\hat{n}}_{r1} \otimes {\hat{n}}_{r1} + {\hat{n}}_{r2} \otimes {\hat{n}}_{r2} + {\hat{n}}_{r3} \otimes {\hat{n}}_{r3} \end{matrix}$ 。

(16)

当应力 $σ$ 处于屈服面上时,并假设屈服函数具有 $f (p, q, θ)$ =0的形式,则其偏法线方向 ${\hat{n}}_{f}$ 可以表示为

${\hat{n}}_{f} = (f_{q} q n_{q} + f_{θ} n_{θ})/ \sqrt{{(f_{q} q)}^{2} + f_{θ}^{2}}$ 。

(17)

屈服面的偏切线方向 ${\hat{n}}_{t}$ 为

${\hat{n}}_{t} = (f_{q} q {\hat{n}}_{θ} - f_{θ} {\hat{n}}_{q}) / \sqrt{{(f_{q} q)}^{2} + f_{θ}^{2}}$ 。

(18)

容易证明,6个方向量 ${\hat{n}}_{p}$ , ${\hat{n}}_{f}$ , ${\hat{n}}_{t}$ , ${\hat{n}}_{r1}$ , ${\hat{n}}_{r2}$ , ${\hat{n}}_{r3}$ 之间也是相互正交的并满足关系：

$\begin{matrix} I = {\hat{n}}_{p} \otimes {\hat{n}}_{p} + {\hat{n}}_{f} \otimes {\hat{n}}_{f} + {\hat{n}}_{t} \otimes {\hat{n}}_{t} + \\ {\hat{n}}_{r1} \otimes {\hat{n}}_{r1} + {\hat{n}}_{r2} \otimes {\hat{n}}_{r2} + {\hat{n}}_{r3} \otimes {\hat{n}}_{r3} \end{matrix}$ 。

(19)

需注意的是,上述讨论都是限制在主应力空间内进行的。为将式（19）推广至一般应力空间,现假设 $e_{1}$ , $e_{2}$ , $e_{3}$ 为一般应力空间的3个单位基矢量,令单位正交张量 $Q$ 表示主应力空间与一般应力空间之间的转换,则 $Q$ 可表示为

$Q = Q_{i j} {\hat{e}}_{i} \otimes {\hat{e}}_{j}$ ,

(20a)

$Q_{i j} = {\hat{e}}_{i} • e_{j}$ 。

(20b)

在一般应力空间中,定义方向 $n_{p} = Q^{T} \cdot {\hat{n}}_{p} \cdot Q$ , $n_{q} = Q^{T} \cdot {\hat{n}}_{q} \cdot Q$ , $n_{θ} = Q^{T} \cdot {\hat{n}}_{θ} \cdot Q$ , $n_{f} = Q^{T} \cdot {\hat{n}}_{f} \cdot Q$ , $n_{t} = Q^{T} \cdot {\hat{n}}_{t} \cdot Q$ , $n_{r1} = Q^{T} \cdot {\hat{n}}_{r1} \cdot Q$ , $n_{r2} = Q^{T} \cdot {\hat{n}}_{r2} \cdot Q$ , $n_{r3} = Q^{T} \cdot {\hat{n}}_{r3} \cdot Q$ 。则可以证明, $n_{p}$ , $n_{q}$ , $n_{θ}$ , $n_{f}$ , $n_{t}$ , $n_{r1}$ , $n_{r2}$ , $n_{r3}$ 与对称转换张量 $I$ 之间存在关系：

$\begin{array}{l} I = n_{p} \otimes n_{p} + n_{q} \otimes n_{q} + n_{θ} \otimes n_{θ} + \\ n_{r1} \otimes n_{r1} + n_{r2} \otimes n_{r2} + n_{r3} \otimes n_{r3} \\ = n_{p} \otimes n_{p} + n_{f} \otimes n_{f} + n_{t} \otimes n_{t} + \\ n_{r1} \otimes n_{r1} + n_{r2} \otimes n_{r2} + n_{r3} \otimes n_{r3} \end{array}$ 。

(21)

利用式（21）,可将应力率 $\dot{σ}$ 表示为

$\begin{matrix} \dot{σ} = I : \dot{σ} = {\dot{σ}}_{p} + {\dot{σ}}_{q} + {\dot{σ}}_{θ} + {\dot{σ}}_{r1} + {\dot{σ}}_{r2} + {\dot{σ}}_{r3} \\ = {\dot{σ}}_{p} + {\dot{σ}}_{f} + {\dot{σ}}_{t} + {\dot{σ}}_{r1} + {\dot{σ}}_{r2} + {\dot{σ}}_{r3} \end{matrix}$ 。

(22)

式中 ${\dot{σ}}_{p} = (\dot{σ} : n_{p}) n_{p}$ , ${\dot{σ}}_{q} = (\dot{σ} : n_{q}) n_{q}$ , ${\dot{σ}}_{θ} = (\dot{σ} :$ $n_{θ}) n_{θ}$ , ${\dot{σ}}_{f} = (\dot{σ} : n_{f}) n_{f}$ , ${\dot{σ}}_{t} = (\dot{σ} : n_{t}) n_{t}$ , ${\dot{σ}}_{r1} = (\dot{σ} :$ $n_{r1}) n_{r1}$ , ${\dot{σ}}_{r2} = (\dot{σ} : n_{r2}) n_{r2}$ , ${\dot{σ}}_{r3} = (\dot{σ} : n_{r3}) n_{r3}$ 。式（22）表明,总应力率 $\dot{σ}$ 既可以分解为 ${\dot{σ}}_{p}$ , ${\dot{σ}}_{q}$ , ${\dot{σ}}_{θ}$ , ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 之和,也可以分解为 ${\dot{σ}}_{p}$ , ${\dot{σ}}_{f}$ , ${\dot{σ}}_{t}$ , ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 之和,并且上述分解方式适用于任意的对称二阶张量。此外,从上述推导过程可知,应力率分量 ${\dot{σ}}_{p}$ , ${\dot{σ}}_{q}$ , ${\dot{σ}}_{θ}$ , ${\dot{σ}}_{f}$ , ${\dot{σ}}_{t}$ 都是由主应力大小的改变所引发,而 ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 3个应力率分量是由主应力方向的改变所引发。进一步地,

$n_{q} = \lim_{Δ q \to 0^{+}} \frac{σ^{*} - σ}{‖ σ^{*} - σ ‖} = \frac{\partial q / \partial σ}{‖ \partial q / \partial σ ‖} = \frac{s}{\sqrt{s : s}}$ ,

(23)

$n_{θ} = \lim_{Δ θ \to 0^{+}} \frac{σ^{*} - σ}{‖ σ^{*} - σ ‖} = \frac{\partial θ / \partial σ}{‖ \partial θ / \partial σ ‖} = \frac{S}{\sqrt{S : S}}$ 。

(24)

利用式（22）～（24）,可将Rudnicki等^[14]定义的非共轴应力率表示为

${\dot{s}}_{nc} = \dot{σ} - {\dot{σ}}_{p} - {\dot{σ}}_{q} = {\dot{σ}}_{θ} + {\dot{σ}}_{r1} + {\dot{σ}}_{r2} + {\dot{σ}}_{r3}$ 。

(25)

式（25）的几何解释如图1所示。

图 1 传统的非共轴理论^[14]

Figure 1. Traditional non-coaxial theory by Ref. [14]

下载: 全尺寸图片幻灯片

同理,可将Hashiguchi等^[15]定义的非共轴应力率表示为

${\dot{s}}_{nc} = \dot{σ} - {\dot{σ}}_{p} - {\dot{σ}}_{f} = {\dot{σ}}_{t} + {\dot{σ}}_{r1} + {\dot{σ}}_{r2} + {\dot{σ}}_{r3}$ 。

(26)

式（26）的几何解释如图2所示。

图 2 传统的非共轴理论^[15]

Figure 2. Traditional non-coaxial theory by Ref. [15]

下载: 全尺寸图片幻灯片

而钱建固等^[16]定义的非共轴应力率可重新表示为

${\dot{s}}_{nc} = \dot{σ} - {\dot{σ}}_{p} - {\dot{σ}}_{q} - {\dot{σ}}_{θ} = {\dot{σ}}_{r1} + {\dot{σ}}_{r2} + {\dot{σ}}_{r3}$ 。

(27)

式（27）的几何解释如图3所示。

图 3 传统的非共轴理论^[16]

Figure 3. Traditional non-coaxial theory by Ref. [16]

下载: 全尺寸图片幻灯片

通过上述分析可清楚得知,Rudnicki等^[14]定义的非共轴应力率包含 ${\dot{σ}}_{θ}$ , ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 共4项正交分量,Hashiguchi等^[15]定义的非共轴应力率包含 ${\dot{σ}}_{t}$ , ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 共4项正交分量,而钱建固等^[16]给出非共轴应力率包含 ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 共3项正交分量。其中,Rudnicki等^[14]理论与Hashiguchi等^[15]理论的主要区别体现在 ${\dot{σ}}_{θ}$ 和 ${\dot{σ}}_{t}$ 的不同。而钱建固等^[16]定义的非共轴应力率仅包含主应力方向改变引发的应力率 ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ ,从而完全排除了主应力大小改变的影响,这是其区另外两种非共轴理论的标志。

需要指出的是,上述结论是基于式（8）定义的方向量和式（22）表示的应力率的正交分解基础上推出的,而基于传统非共轴理论以及Li等^[19]研究给出的非共轴应力率的表述方式难以得出类似结论。

3. 广义非共轴理论

传统弹塑性理论^[21]中通常将加载强度指标 $L$ 定义为

$L = \frac{1}{K_{p}} \frac{\partial f}{\partial σ} : \dot{σ}$ ,

(28)

式中, $K_{p}$ 为塑性模量。利用应力率的分解式（22）并根据分量之间的正交关系,可将式（28）表示为

$L = \frac{1}{K_{p}} \frac{\partial f}{\partial σ} : \dot{σ} = \frac{1}{K_{p}} \frac{\partial f}{\partial σ} ({\dot{σ}}_{p} + {\dot{σ}}_{f})$ 。

(29)

由式（29）可知,传统弹塑性模型中仅考虑了应力率分量 ${\dot{σ}}_{p}$ 与 ${\dot{σ}}_{f}$ 诱发的塑性变形,而忽略了非共轴应力率分量 ${\dot{σ}}_{t}$ , ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 的影响。尽管Hashiguchi等^[15]非共轴理论可以反映这4项的影响,但该理论将非共轴塑性应变率与非共轴应力率之间表示为了线性关系并且二者始终同向,这一假设与岩土材料的试验结果不相符。例如,Gutierrez等^[22]在对Toyoura砂进行空心圆柱扭剪试验时发现,在主应力幅值保持不变而方向连续改变的旋转加载条件下,非共轴塑性应变率显著偏离了非共轴应力率的方向。此外,砂土试验结果^[23-26]表明,非共轴应力率所诱发的塑性变形不仅包括偏应变也包括体应变,而且砂土的非共轴塑性应变率还受当前孔隙比、应力状态和各向异性等因素的影响。此外,需要注意的是, ${\dot{σ}}_{t}$ 表示主应力大小的改变,而 ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 三项代表了主应力方向的改变,因而 ${\dot{σ}}_{t}$ 与 ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 代表了两类不同的加载机制。

为了精确描述非共轴应力率的影响,本文建议针对每一分量应力率建立独立的非线性加载机制。首先以分量 ${\dot{σ}}_{t}$ 为例,本文将与其相关的加载称为第一类非共轴加载。参考Zienkiewicz等^[20]提出的广义塑性力学,通过直接定义加载强度指标、塑性模量和塑性流动方向,可将相应的加载机制表示为

$L_{t} = \frac{1}{K_{p}^{t}} (n_{t} : \dot{σ}) = \frac{1}{K_{p}^{t}} (n_{t} : {\dot{σ}}_{t})$ ,

(30a)

${\dot{ε}}_{t}^{p} = L_{1} (n_{t}^{*} + \frac{1}{3} D_{t} \sqrt{2 / 3})$ ,

(30b)

$K_{p}^{t} = K_{p}^{t} (σ, ψ, C)$ ,

(30c)

$D_{t} = D_{t} (σ, ψ, C)$ 。

(30d)

式中 $L_{t}$ 为 $n_{t}$ 方向的加载强度指标； ${\dot{ε}}_{t}^{p}$ 为分量 ${\dot{σ}}_{t}$ 诱发的塑性应变率； $n_{t}^{*}$ 为 ${\dot{σ}}_{t}$ 诱发的塑性偏应变率的方向； $K_{p}^{t}$ , $D_{t}$ 为与 ${\dot{σ}}_{t}$ 对应的塑性模量和剪胀系数,二者可表示为前应力 $σ$ 、状态参数 $ψ$ 和代表材料特性的一组参数 $C$ 的函数。考虑到当前的研究普遍认为非共轴性源于材料的各向异性,因此参数 $C$ 中应当包括代表各向异性程度的幅值参量,并且当材料为各向同性时,需满足塑性模量 $K_{p}^{t} = + \infty$ ,以保证 $n_{t}$ 方向的应力率不会诱发塑性变形。注意到,在Gutierrez等^[22]进行的主应力轴纯旋转加载试验中,随着应力比的增加,砂土变形的非共轴性会相应减弱,且当应力状态接近临界状态时,非共轴性几乎消失。基于这一规律,Gutierrez等^[22]建议用应力率与边界面交点处的外法线方向表示塑性流动方向,如图4所示。

图 4 塑性流动方向^[22]

Figure 4. Plastic flow direction by Ref [22]

下载: 全尺寸图片幻灯片

关于 $n_{t}^{*}$ 方向的确定,这里采用与Gutierrez等^[22]相似的方法,将 ${\dot{σ}}_{t}$ 诱发的塑性偏应变率方向定义为 $n_{t}$ 与边界面交点处的外法线方向,该方法的几何解释如图5所示。

图 5 几何映射方法示意图

Figure 5. Schematic diagram of geometric mapping method

下载: 全尺寸图片幻灯片

此外,也可以采用如下的方式定义 $n_{t}^{*}$ 为

$n_{t}^{*} = \sqrt{1 - χ} n_{t} + \sqrt{χ} n_{f}$ ,

(31)

式中, $χ$ 称为方向系数且 $0 \leq χ \leq 1$ 。理论上,通过将 $χ$ 表示为随应力状态变化的函数也可达到与上述几何映射方法相似的描述效果,如图6所示。

图 6 插值系数法示意图

Figure 6. Schematic diagram of interpolation coefficient method

下载: 全尺寸图片幻灯片

另一方面,从几何角度可知, $n_{t}$ 和 $- n_{t}$ 均表示偏平面上屈服函数的切线方向。关于加载过程中 $n_{t}$ 和 $- n_{t}$ 的选择,可类比于传统弹塑性模型中加卸载的概念,引入加载与反向加载的概念,如图7所示。

图 7 加载与反向加载

Figure 7. Loading and reverse loading

下载: 全尺寸图片幻灯片

在弹塑性加载过程中,首先按照弹性变形计算试探应力率 ${\dot{σ}}^{trial}$ ,若 ${\dot{σ}}^{trial} : n_{t} \geq 0$ ,选用 $n_{t}$ 计算；相反,若 ${\dot{σ}}^{trial} : n_{t} < 0$ ,选用 $- n_{t}$ 代替式（30）中的 $n_{t}$ 进行计算。

至此,完成了对应于分量 ${\dot{σ}}_{t}$ 的加载机制的描述。采用同样的方式,可以将对应于分量 ${\dot{σ}}_{r1}$ , ${\dot{σ}}_{r2}$ , ${\dot{σ}}_{r3}$ 的加载称为第二类非共轴加载,并将相应的加载机制表示为

$L_{r i} = \frac{1}{K_{p}^{r i}} (n_{r i} : \dot{σ}) = \frac{1}{K_{p}^{r i}} (n_{r i} : {\dot{σ}}_{t})$ ,

(32a)

${\dot{ε}}_{r i}^{p} = L_{r i} (n_{r i}^{*} + \frac{1}{3} D_{r i} \sqrt{2 / 3})$ ,

(32b)

$K_{p}^{r i} = K_{p}^{r i} (σ, ψ, C)$ ,

(32c)

$D_{r i} = D_{r i} (σ, ψ, C)$ 。

(32d)

式中指标 $r i$ 分别取 $r1$ , $r2$ , $r3$ ； $L_{r i}$ 为 $n_{r i}$ 方向的加载强度指标； ${\dot{ε}}_{r i}^{p}$ 为分量 ${\dot{σ}}_{r i}$ 诱发的塑性应变率； $n_{r i}^{*}$ 为 ${\dot{σ}}_{r i}$ 诱发的塑性偏应变率的方向； $K_{p}^{r i}$ , $D_{r i}$ 为与 ${\dot{σ}}_{r i}$ 对应的塑性模量和剪胀系数,二者同样地都可以表示为当前应力 $σ$ 、状态参数 $ψ$ 和代表材料特性参数 $C$ 的函数。3个加载机制的示意图如图8所示。

图 8 第二类非共轴加载

Figure 8. Second kind of non-coaxial loading

下载: 全尺寸图片幻灯片

其中,非共轴塑性流动方向 $n_{r i}^{*}$ 可为

$n_{r i}^{*} = \sqrt{1 - χ} n_{r i} + \sqrt{χ} n_{f}$ 。

(33)

也可采用图5所示的映射方法确定,只需将图中的 $n_{t}$ 理解为 $n_{r i}$ 即可。

将非共轴应力率进行正交分解并参考Zienkiewicz等^[20]的广义塑性力学,通过直接定义塑性指标变量建立了多重非线性加载机制,用于描述每一非共轴应力率分量诱发的塑性变形。通过采用边界面理论中的几何映射法则或采用方向插值系数法可描述非共轴塑性流动方向在屈服面偏切线方向与偏法线方向之间连续变化。因而相较于传统理论,本理论在描述土体的非共轴塑性变形方面具有更强的灵活性。

4. 数值试验

为了评估本文建立的非共轴理论在分析岩土材料非共轴变形方面的合理性,结合该理论与Li等^[27]发展的状态相关模型,建立砂土的非共轴本构模型。模型中的状态变量定义为

$ψ = e - e_{c} = e - [e_{Γ} - λ_{c} {(p / p_{a})}^{ξ}]$ 。

(34)

式中 $e$ 为当前孔隙比； $e_{c}$ 为临界孔隙比； $e_{Γ}$ , $λ_{c}$ , $ξ$ 均为材料常数； $p_{a}$ 为大气压力。

剪切模量 $G$ 和体积模量 $K$ 为

$G = G_{0} \frac{{(2.97 - e)}^{2}}{1 + e} \sqrt{p p_{a}}$ ,

(35)

$K = G \frac{2 (1 + ν)}{3 (1 - 2 ν)}$ ,

(36)

式中, $G_{0}$ 为材料常数, $ν$ 为泊松比。

屈服函数 $f$ 为

$f = q - M g (θ) p = 0$ ,

(37)

式中, $M$ 为硬化内变量, $g (θ)$ 为角隅函数^[10],

$g (θ) = \frac{\sqrt{{(1 + c^{2})}^{2} + 4 (1 - c^{2}) \sin 3 θ} - (1 + c^{2})}{2 (1 - c) \sin 3 θ}$ ,

(38)

其中, $c$ 为三轴拉伸与压缩条件下极限应力比的比值。

塑性模量 $K_{p}$ 为

$K_{p} = h G e^{n ψ} [\frac{M_{c} g (θ) e^{- n ψ}}{η} - 1]$ 。

(39)

式中 $M_{c}$ 为临界状态应力比；应力比 $η = q / p$ ； $n$ 为材料常数； $h = h_{1} - h_{2} e$ , $h_{1}$ , $h_{2}$ 为2个材料常数。剪胀系数 $D$ 和塑性势函数 $g$ 为

$D = \frac{d_{0}}{M_{c} g (θ)} [M_{c} g (θ) e^{m ψ} - η]$ ,

(40a)

$g = q - \frac{d_{0} M_{c} g (θ) e^{m ψ} p}{d_{0} - M_{c} g (θ)} [1 - {(\frac{p}{p_{0}})}^{\frac{d_{0}}{M_{c} g (θ)} - 1}]$ ,

(40b)

其中, $d_{0}$ , $m$ 为2个材料常数, $p_{0}$ 为塑性势面与 $p$ 轴的交点。

基于小变形假设,将总应变率 $\dot{ε}$ 分解为弹性应变率 ${\dot{ε}}^{e}$ 、共轴塑性应变率 ${\dot{ε}}^{cp}$ 与非共轴塑性应变率 ${\dot{ε}}_{t}^{p}$ , ${\dot{ε}}_{r1}^{p}$ , ${\dot{ε}}_{r2}^{p}$ , ${\dot{ε}}_{r3}^{p}$ 之和：

$\dot{ε} = {\dot{ε}}^{e} + {\dot{ε}}^{cp} + {\dot{ε}}_{t}^{p} + {\dot{ε}}_{r1}^{p} + {\dot{ε}}_{r2}^{p} + {\dot{ε}}_{r3}^{p}$ 。

(41)

式（41）中的共轴塑性应变率 ${\dot{ε}}^{cp}$ 由传统位势理论确定,而非共轴塑性应变率 ${\dot{ε}}_{t}^{p}$ , ${\dot{ε}}_{r1}^{p}$ , ${\dot{ε}}_{r2}^{p}$ , ${\dot{ε}}_{r3}^{p}$ 由多重非共轴加载机制确定。令 $L_{1} = L$ , $L_{2} = L_{t}$ , $L_{3} = L_{r 1}$ , $L_{4} = L_{r2}$ , $L_{5} = L_{r3}$ , $N_{1} = \partial f / \partial σ$ , $N_{2} =$ $n_{t}$ , $N_{3} = n_{r1}$ , $N_{4} = n_{r2}$ , $N_{5} = n_{r3}$ , $N_{1}^{*} = \partial g / \partial σ$ , $N_{2}^{*} = n_{t}^{*} + D_{t} \sqrt{2 / 3} 1 / 3$ , $N_{3}^{*} = n_{r1}^{*} + D_{r1} \sqrt{2 / 3} 1 / 3$ , $N_{4}^{*} =$ $n_{r2}^{*} + D_{r2} \sqrt{2 / 3} 1 / 3$ , $N_{5}^{*} = n_{r3}^{*} + D_{r3} \sqrt{2 / 3} 1 / 3$ 。利用弹性本构关系与式（30）,（32）,（41）可得

$\dot{σ} = E : \dot{ε} - \sum_{I = 1}^{5} L_{I} E : N_{I}^{*}$ 。

(42)

令 $K_{p}^{1} = K_{p}$ , $K_{p}^{2} = K_{p}^{t}$ , $K_{p}^{3} = K_{p}^{r1}$ , $K_{p}^{4} = K_{p}^{r2}$ , $K_{p}^{5} = K_{p}^{r3}$ ,分别用 $N_{1}$ , $N_{2}$ , $N_{3}$ , $N_{4}$ , $N_{5}$ 双点积式（42）,则可得

$\begin{matrix} (K_{p}^{I} + N_{I} : E : N_{I}^{*}) L_{I} + \sum_{J = 1, J \neq I}^{5} N_{I} : E : N_{J}^{*} L_{J} \\ = N_{I} : E : {\dot{ε}}^{e} \end{matrix}$ 。

(43)

式（43）是关于加载强度指标 $L_{I}$ 的一组线性方程,令 $Λ$ 表示 $L_{I}$ 的系数矩阵,求解该方程组可得

$L_{J} = \sum_{I = 1}^{5} b_{I} Λ_{I J}^{*} / \det (Λ)$ 。

(44)

式中 $b_{I} = N_{I} : E : {\dot{ε}}^{e}$ ； $Λ^{*}$ 为 $Λ$ 的伴随矩阵； $\det (Λ)$ 为 $Λ$ 的行列式。将式（44）代入式（42）,可得非共轴本构关系式：

$\dot{σ} = [E^{e} - E^{e} : (\sum_{I = 1}^{5} N_{I}^{*} \otimes M_{I}) : E^{e}] : \dot{ε}$ ,

(45a)

$M_{J} = \sum_{I = 1}^{5} N_{I} Λ_{I J}^{*} / \det (Λ)$ 。

(45b)

为了评估模型描述砂土本构反应的合理性,本文设计了应力与应变两种探测试验。如图8（b）所示,应力探测试验模拟主应力绕中主应力 $σ_{2}$ 的纯旋转加载。

材料假设为Toyoura砂,共轴模型参数取自文献[2],并假设 $K_{p}^{t} = K_{p}^{r1} = K_{p}^{r2} = K_{p}^{r3} = K_{p}$ , $D_{t} = D_{r1} =$ $D_{r2} =$ $D_{r3} = D$ 。初始时刻平均应力 $p =$ 98 kPa,剪应力 $q =$ 50 kPa,罗德角 $θ =$ 0.0,初始孔隙比 $e_{0}$ =0.78。加载过程中,大主应力方向角 $α$ 从0°连续增加至180°。非共轴塑性流动方向采用式（31）确定,方向系数 $χ$ 分别为0.0,0.25,0.5,0.75,模拟结果如图9所示。

图 9 应力探测示意图

Figure 9. Schematic diagram of stress probe tests

下载: 全尺寸图片幻灯片

图9表明,当 $χ = 0.0$ 时,计算的非共轴塑性应变率方向与应力路径正切；当 $χ \neq 0.0$ 时,计算的非共轴塑性应变率方向处于应力路径的切线方向与法线方向之间,更符合试验结果。随着 $χ$ 的增大,非共轴塑性增量方向由切线方向逐渐向法线方向靠近。因此,基于广义非共轴理论建立的模型可以更加真实的反映砂土的变形特点。

如图10所示,在应变探测试验中,初始时刻的应力状态用A点表示,并假设与前述应力探测试验中的初始状态相同。然后在二维空间内,沿任意方向施加应变率幅值 $‖ \dot{ε} ‖$ 为0.0001的应变荷载。Gudehus^[28]首先提出了应变探测方法,并认为合理的应力率响应包络线应满足的基本条件是连续性。

图 10 应变探测示意图

Figure 10. Schematic diagram of stress probe tests

下载: 全尺寸图片幻灯片

图11（a）展示了弹性模型、共轴弹塑性模型和非共轴弹塑性模型（ $χ = 0.0$ ）计算的应力率响应包络线。显然,3种模型预测的包络线满足连续性要求。此外,由图11（a）可知,沿着0°方向和180°方向加载,非共轴模型（ $χ = 0.0$ ）与共轴模型计算的应力率相同,而沿其他方向加载,非共轴模型（ $χ = 0.0$ ）的计算结果显示出更加软化的力学响应。图11（b）展示了不同 $χ$ 条件下,非共轴模型计算的应力率包络线。由图11可知,随着 $χ$ 的增加,在90°～180°～270°范围内,模型展示出了相对硬化的力学响应,而在270°～0°～90°范围内, $χ$ 对包络线的影响并不明显。此外,不同 $χ$ 条件下,模型计算的包络线均满足连续性要求。应变探测试验结果初步表明,基于广义非共轴理论建立的模型满足数值计算的稳定性要求。

图 11 应力率响应包络线

Figure 11. Response envelopes of stress rate

下载: 全尺寸图片幻灯片

5. 结论

（1）传统非共轴理论将非共轴应力率定义为总应力率减去所谓共轴应力率的剩余部分,无法说明非共轴应力率由多项正交分量构成的事实。本文通过严格的数学推导证明了总应力率可以分解为6项正交应力率分量之和。在此基础上,对于已有的非共轴理论进行了总结和比较,指出,Rudnicki等与Hashiguchi等定义的非共轴应力率包含4项正交的应力率分量,而钱建固等定义的非共轴应力率包含3项正交的应力率分量。

（2）传统非共轴理论假设非共轴塑性应变率与非共轴应力率之间为线性关系并且二者始终同向,这与试验结果不相符。将非共轴应力率正交分解为4项,然后采用多重非线性加载机制描述每一非共轴应力率分量诱发的塑性变形,而总的非共轴塑性变形视为分量非共轴塑性变形之和,从而建立了一种广义形式的非共轴理论。新建理论能够描述非共轴应力率与非共轴塑性应变率之间的非线性关系,同时也能描述二者方向之间的非一致性。

（3）推导了基于广义非共轴理论建立弹塑性模型时的应力-应变关系式。然后用两种数值试验评估了模型的合理性,结果表明,基于广义非共轴理论建立的模型预测的应力率响应包络线满足连续性要求,在描述砂土非共轴变形方面具有更强的灵活性。数值试验结果初步显示了广义非共轴理论在砂土本构建模方面的优越性和合理性,表明该理论能够为土体非共轴模型的构建提供更为广阔的理论基础。

图 1 岩石试样和裂隙错位情况

Figure 1. Rock samples and dislocated fractures

下载: 全尺寸图片幻灯片

图 2 试验设备及传感器设置

Figure 2. Diagram of experimental equipment and settings of sensors

下载: 全尺寸图片幻灯片

图 3 裂隙面几何形貌

Figure 3. Morphological characteristics of fracture samples

下载: 全尺寸图片幻灯片

图 4 试验与数值计算得到的法向应力-法向位移曲线对比

Figure 4. Comparison of experimentally and numerically obtained normal stress-normal displacement curves

下载: 全尺寸图片幻灯片

图 5 试验测量与数值计算得到的表面破坏区域对比

Figure 5. Comparison of experimentally measured and numerically calculated asperity damage areas

下载: 全尺寸图片幻灯片

图 6 塑性变形均值增量-法向应力关系图

Figure 6. Relationship between mean plastic deformation increment and normal stress

下载: 全尺寸图片幻灯片

图 7 u_p/e₀与力学及几何参数之间的关系

Figure 7. Relationship among u_p/e₀, mechanical and geometric parameters

下载: 全尺寸图片幻灯片

图 8 声发射定位点与试验测量得到的表面破坏区域对比图

注：红色为三维定位点，青色为三维定位点的X-Y投影，蓝色为三维定位点的X-Z投影，绿色为三维定位点的Y-Z投影。

Figure 8. Comparison between AE positioning points and experimentally measured asperity damage areas

下载: 全尺寸图片幻灯片

图 9 岩石裂隙声发射响应特征

Figure 9. AE response characteristics of rock fractures

下载: 全尺寸图片幻灯片

表 1 岩样的基本性质

Table 1 Fundamental properties of rock samples

试样编号	UCS/MPa	弹性模量/GPa	泊松比	P波波速/(m/s)	旋转角度/(°)
G-D	148.94	47.33	0.24	3179	6.11
G-S	133.20	44.99	0.23	4292	6.80
R-D	71.45	15.10	0.40	2689	3.49
R-S	39.98	11.20	0.35	3022	14.77
注：G为花岗岩,R为红砂岩,D为干燥,S为饱和。

下载: 导出CSV

表 2 试样加卸载方案

Table 2 Loading and unloading plans

试样编号	加载速率/(MPa·s^-1)	目标应力/MPa	卸载速率/(MPa·s^-1)	保载时间/s
G-D	0.1	10	0.1	45
G-W	0.1	10	0.1	45
R-D	0.1	10	0.1	45
R-W	0.05	5	0.05	45

下载: 导出CSV

参考文献(29)

[1]	ERGULER Z A, ULUSAY R. Water-induced variations in mechanical properties of clay-bearing rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2009, 46(2): 355-370. doi: 10.1016/j.ijrmms.2008.07.002
[2]	RAJABZADEH M A, MOOSAVINASAB Z, RAKHSHANDEHROO G. Effects of rock classes and porosity on the relation between uniaxial compressive strength and some rock properties for carbonate rocks[J]. Rock Mechanics and Rock Engineering, 2012, 45(1): 113-122. doi: 10.1007/s00603-011-0169-y
[3]	WONG L N Y, MARUVANCHERY V, LIU G. Water effects on rock strength and stiffness degradation[J]. Acta Geotechnica, 2016, 11(4): 713-737. doi: 10.1007/s11440-015-0407-7
[4]	BANDIS S C, LUMSDEN A C, BARTON N R. Fundamentals of rock joint deformation[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1983, 20(6): 249-268.
[5]	PEI L, HYUN S, MOLINARI J F, et al. Finite element modeling of elasto-plastic contact between rough surfaces[J]. Journal of the Mechanics & Physics of Solids, 2005, 53(11): 2385-2409.
[6]	LAVROV A. Fracture permeability under normal stress: a fully computational approach[J]. Journal of Petroleum Exploration and Production Technology, 2017, 7(1): 181-194. doi: 10.1007/s13202-016-0254-6
[7]	TIAN X F, BHUSHAN B. A numerical three-dimensional model for the contact of rough surfaces by variational principle[J]. Journal of Tribology, 1996, 118(1): 33-42. doi: 10.1115/1.2837089
[8]	HOPKINS D L. The Effect of Surface Roughness on Joint Stiffness, Aperture, and Acoustic Wave Propagation[D]. Berkeley: University of California, 1991.
[9]	LI B, ZHAO Z H, JIANG Y J, et al. Contact mechanism of a rock fracture subjected to normal loading and its impact on fast closure behavior during initial stage of fluid flow experiment[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(13): 1431-1449. doi: 10.1002/nag.2365
[10]	GREENWOOD J A, WILLIAMSON J B P P. Contact of nominally flat surfaces[J]. Proceedings of the Royal Society of London, 1966, 295(1442): 300-319.
[11]	KLING T, VOGLER D, PASTEWKA L, et al. Numerical simulations and validation of contact mechanics in a granodiorite fracture[J]. Rock Mechanics and Rock Engineering, 2018, 51(9): 2805-2824. doi: 10.1007/s00603-018-1498-x
[12]	ZOU L C, LI B, MO Y Y, et al. A high-resolution contact analysis of rough-walled crystalline rock fractures subject to normal stress[J]. Rock Mechanics and Rock Engineering, 2020, 53(5): 2141-2155. doi: 10.1007/s00603-019-02034-w
[13]	RUDAJEV V, VILHELM J, LOKAJÍČEK T. Laboratory studies of acoustic emission prior to uniaxial compressive rock failure[J]. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(4): 699-704. doi: 10.1016/S1365-1609(99)00126-4
[14]	陈国庆, 陈毅, 孙祥, 等. 开放型岩桥裂纹贯通机理及脆性破坏特征研究[J]. 岩土工程学报, 2020, 42(5): 908-915. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202005018.htm CHEN Guo-qing, CHEN Yi, SUN Xiang, et al. Crack coalescence and brittle failure characteristics of open rock bridges[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(5): 908-915. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202005018.htm
[15]	LIU S M, LI X L, WANG D K, et al. Mechanical and acoustic emission characteristics of coal at temperature impact[J]. Natural Resources Research, 2019, 29(4): 1755-17772.
[16]	龚囱, 李长洪, 赵奎. 红砂岩短时蠕变声发射b值特征[J]. 煤炭学报, 2015, 40(增刊1): 85-92. https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB2015S1013.htm GONG Cong, LI Chang-hong, ZHAO Kui. Study on b-value characteristics of acoustic emission of red sandstone during short-time creep process[J]. Journal of China Coal Society, 2015, 40(S1): 85-92. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB2015S1013.htm
[17]	王春来, 廖泽锋, 李长峰, 等. 花岗岩岩爆声发射时空熵值动态特征实验研究[J]. 采矿与安全工程学报, 2019, 36(3): 626-633. https://www.cnki.com.cn/Article/CJFDTOTAL-KSYL201903027.htm WANG Chun-lai, LIAO Ze-feng, LI Chang-feng, et al. Experimental investigation of dynamic characteristics of AE spatio-temporal entropy for granitic rockburst[J]. Journal of Mining & Safety Engineering, 2019, 36(3): 626-633. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-KSYL201903027.htm
[18]	LI L R, DENG J H, ZHENG L, et al. Dominant frequency characteristics of acoustic emissions in white marble during direct tensile tests[J]. Rock Mechanics and Rock Engineering, 2017, 50(5): 1337-1346. doi: 10.1007/s00603-016-1162-2
[19]	赵兴东, 李元辉, 袁瑞甫, 等. 基于声发射定位的岩石裂纹动态演化过程研究[J]. 岩石力学与工程学报, 2007, 26(5): 944-950. doi: 10.3321/j.issn:1000-6915.2007.05.011 ZHAO Xing-dong, LI Yuan-hui, YUAN Rui-pu, et al. Study on crack dynamic propagation process of rock samples based on acoustic emission location[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(5): 944-950. (in Chinese) doi: 10.3321/j.issn:1000-6915.2007.05.011
[20]	刘建坡, 徐世达, 李元辉, 等. 预制孔岩石破坏过程中的声发射时空演化特征研究[J]. 岩石力学与工程学报, 2012, 31(12): 2538-2547. doi: 10.3969/j.issn.1000-6915.2012.12.018 LIU Jian-po, XU Shi-da, LI Yuan-hui, et al. Studies of AE time-space evolution characteristics during failure process of rock specimens with prefabricated holes[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(12): 2538-2547. (in Chinese) doi: 10.3969/j.issn.1000-6915.2012.12.018
[21]	DONG L J, HU Q C, TONG X J, et al. Velocity-free MS/AE source location method for three-dimensional hole-containing structures[J]. Engineering, 2020, 6(7): 827-834. doi: 10.1016/j.eng.2019.12.016
[22]	GEIGER L. Probability method for determination of earthquake epicenters form arrival time only[J]. Bulletin of Saint Louis University, 1912, 8: 60-71.
[23]	赵兴东, 刘建坡, 李元辉, 等. 岩石声发射定位技术及其实验验证[J]. 岩土工程学报, 2008, 30(10): 1472-1476. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200810012.htm ZHAO Xing-dong, LIU Jian-po, LI Yuan-hui, et al. Experimental verification of rock locating technique with acoustic emission[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(10): 1472-1476. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200810012.htm
[24]	BOUSSINESQ J. Application des potentiels: à l'étude de l'équilibre et du mouvement des solides élastiques[J]. Gauthier-Villard, Paris, 1885.
[25]	MOMBER A W. Fracture toughness effects in geomaterial solid particle erosion[J]. Rock Mechanics and Rock Engineering, 2015, 48(4): 1573-1588.
[26]	PRADHAN P S, KING R L, YOUNAN N H, et al. Estimation of the number of decomposition levels for a wavelet-based multiresolution multisensor image fusion[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(12): 3674-3686.
[27]	BROWN S R, SCHOLZ C H. Closure of rock joints[J]. Journal of Geophysical Research: Solid Earth, 1986, 91(B5): 4939-4948.
[28]	王笑然, 李楠, 王恩元, 等. 岩石裂纹扩展微观机制声发射定量反演[J]. 地球物理学报, 2020, 63(7): 2627-2643. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX202007013.htm WANG Xiao-ran, LI Nan, WANG En-yuan, et al. Microcracking mechanisms of sandstone from acoustic emission source inversion[J]. Chinese Journal of Geophysics, 2020, 63(7): 2627-2643. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX202007013.htm
[29]	YAO Q L, CHEN T, JU M H, et al. Effects of water intrusion on mechanical properties of and crack propagation in coal[J]. Rock Mechanics and Rock Engineering, 2016, 49(12): 4699-4709.