岩石轴压CT图像的分区损伤信息与临界破坏识别

刘树新; 邢杰; 郑旭; 陈世杰

doi:10.11779/CJGE202103005

岩石轴压CT图像的分区损伤信息与临界破坏识别 English Version

1.
内蒙古科技大学矿业研究院,内蒙古　包头　014010
2.
中国科学院寒区旱区环境与工程研究所,甘肃　兰州　730000

基金项目:

国家自然科学基金项目 51264028

内蒙古自治区高等学校科学研究项目 NJZY11147

详细信息

作者简介:
刘树新（1971— ）,男,博士,教授,从事岩石力学与采矿工程方面的教学科研工作。E-mail：sxbt@163.com

中图分类号: TU452
计量
- 文章访问数: 320
- HTML全文浏览量: 24
- PDF下载量: 208
出版历程
- 收稿日期: 2020-06-18
- 网络出版日期: 2022-12-04
- 刊出日期: 2021-02-28

Zonal damage information and critical failure identification of CT images of rock under triaxial compression

1.
Mining Institute of Inner Mongolia University of Science and Technology, Baotou 014010, China
2.
Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou 730000, China

摘要

摘要: 岩石三轴压缩试验过程中虽然损伤不断积累,但宏观破坏前其CT图像却没有明显的变化,因此在峰值前仅通过CT数难以准确分析其损伤破坏信息。为了识别岩石的损伤信息及临界破坏特征,将轴压过程中的砂岩CT图像进行分区,使用位置线上灰度的变化代表分区内灰度的变化。CT图像中某一位置线上灰度值的变化曲线是一种起伏、粗糙的曲线,利用描述粗糙曲线的分形R_d指标定量分析砂岩损伤规律,并基于CT图像上全部位置线的R_d值总和来定义损伤变量,用其表征砂岩的损伤程度。计算结果表明荷载水平在0.84左右时,R_d值波动大,损伤变量加速增大,表明损伤进入了加速发展阶段,岩石处于损坏破坏的临界状态。由此可以合理地预测岩石的损伤破坏。
- 砂岩 /
- CT图像 /
- 分形指标 /
- 临界损伤信息
Abstract: The CT image and value do not change obviously before the macro-failure under triaxial compression experiments on rock during damage accumulation, so it is difficult to analyze accurately the damage information only through the change of CT number. In order to identify the damage information and critical failure characteristics of sandstone, the CT image is partitioned, and the change of grey level on the position line is used to represent the change of the gray level in the subarea. The grey value curve is a rough curve representing grey change of a certain location online in the CT image. It is used to describe quantitatively the rules of sandstone damage by fractal R_d index, the damage variable is defined as the sum of R_d value of all the position line in the CT image, and it reflects the damage degree of rock. The calculated results show that when the load is about 0.84, the value of R_d has obvious fluctuations and the damage variable increases rapidly, which denotes that the damage has entered a phase of accelerated development, and the rock is in a critical damage state. The proposed method can be used to predict reasonably the damage and failure of rock.
- sandstone /
- CT image /
- fractal index /
- critical damage information

HTML全文

0. 引言

水力压裂技术被广泛应用于解决非常规油气资源开发中的增透提渗、矿山治理中的注浆加固及地应力测量等问题^[1-4]。水力压裂破裂压力是指井内液体压力不断增大使地层产生水力裂缝或张开原有裂缝时所达到的峰值压力，与围岩物理力学特性、地应力状态、流体压力及天然裂缝发育情况密切相关，是控制水力压裂设计和施工的重要参数^[5-7]。现有水力压裂破裂压力计算模型多基于应用广泛的Hubbert-Willi(H-W)模型和Haimson-Fairhurst(H-F)模型^[8-9]，其中H-W模型未考虑围岩的渗透效应，可得到破裂压力的上限值；H-F模型考虑了围岩的高渗透性，得到破裂压力的下限值^[10]。王六鹏等^[11]基于最优化理论扩展了H-W模型，利用井周围岩应力分布和张性破裂准则，分别建立了裸眼和射孔完井条件下破裂压力非线性约束优化计算模型；郭天魁等^[12]在H-F模型基础上考虑了渗流引起的应力，提出了新的水平井射孔孔壁的应力分布，建立了考虑页岩射孔水平井水力压裂在围岩内破裂、沿天然裂缝剪切破裂和沿天然裂缝张性破裂3种破裂压力计算模型。李传亮等^[13-14]基于双重有效应力概念，集合了H-W模型和H-F模型，建立了裸眼井完井条件下围岩破裂压力计算模型，由此提出了射孔完井条件下产生垂向裂缝和水平裂缝情况下围岩破裂压力计算模型。杨兆中等^[15]在井筒围岩应力场计算中，考虑套管和水泥环的影响，建立了煤层直井井周应力场解析模型。上述模型虽考虑了诸多影响因素，但水力压裂过程中注入流体渗流和围岩应力的耦合作用尚未涉及。

水力压裂是井内流体压力、围岩内流体渗流和围岩应力相互耦合的过程^[16-17]；水力压裂过程中，随着流体注入，井内流体压力逐渐升高；压力驱动作用下流体持续注入井壁围岩内，导致围岩流体压力改变，围岩有效应力场亦改变，影响围岩内部结构的演化，引起流体应力场再次改变；当围岩应力状态达到破裂准则，射孔围岩破裂，此时井内流体压力即为水力压裂破裂压力。故水力压裂过程中，射孔围岩应力场和流体场处于相互联系、相互影响和相互制约的水力耦合状态。渗透率和孔隙度的应力敏感性是研究水力耦合问题的关键。黄远智等^[18]基于渗流–应力耦合试验，认为渗透率演化在于有效围压作用下岩石孔隙变形，且渗透率会随有效围压的增加呈指数下降。卢家亭等^[19]研究了低渗砂岩微观结构与渗透率应力敏感性的关系，认为乘幂函数可表述渗透率和有效应力的关系。薛永超等^[20]探讨了不同渗透率岩石的应力敏感性规律，认为中高渗储层岩石应力敏感性符合“平缓型”模式，而低渗透储层岩石应力敏感性呈“先陡后缓”的规律，且低渗透岩石具有更强的应力敏感性。上述研究表明，岩石渗透率和孔隙度的应力敏感性研究取得了较多的成果，但将其考虑至水力压裂耦合的研究中鲜有成果。

综上所述，现存水力压裂射孔围岩破裂模拟研究多考虑单纯的应力作用，与流体渗流耦合研究成果较少。射孔围岩流体渗流具有非稳态性，导致与围岩应力的水力耦合作用变得异常复杂，且相关围岩应力状态亦产生较大改变。本文基于渗透率和孔隙度的应力敏感性，构建了考虑流体渗流和围岩应力耦合的射孔围岩水力压裂破裂力学模型，并基于有限容积法提出了水力耦合作用下射孔围岩水力压裂破裂数值模拟方法，进而验证了方法的正确性。

1. 射孔围岩应力分布

1.1 物理模型及假设

储层岩石处于垂向地应力 ${\sigma _{\text{V}}}$ 、最大水平应力 ${\sigma _{\text{H}}}$ 和最小水平应力 ${\sigma _{\text{h}}}$ 状态下。假设射孔井是由井筒和射孔组成的正交系统，在远场地应力与井内流体压力作用下，近井区域围岩应力重分布，引起双重应力集中。

为了准确分析射孔围岩应力状态，提出以下假设：①储层岩石为均匀各向同性、线弹性多孔材料，且井壁围岩处于平面应变状态；②井筒与井壁胶结良好时，水力压裂过程中光滑接触；井筒与井壁胶结不好时，不考虑井筒与井壁的相互作用，认为井壁光滑且与井筒间存在微环隙，流体较容易流入微环隙；③在渗流计算时不考虑流体的压缩性；④井筒内和射孔内流体压力相等，储层内均产生渗流。

1.2 井筒围岩应力分布

考虑围岩在远场地应力和井筒内流体压力共同作用下(图 1)，结合井壁渗流状态，利用叠加原理得到极坐标系下围岩的应力分布^[21-22]：

图 1 井筒围岩力学模型

Figure 1. Mechanical model for surrounding rock of well

下载: 全尺寸图片幻灯片

${\sigma _r} = \frac{{({\sigma _{\text{H}}} + {\sigma _{\text{h}}})}}{2}\left( {1 - \frac{{r_{\text{w}}^{\text{2}}}}{{{r^2}}}} \right) +\\ \frac{{({\sigma _{\text{H}}} - {\sigma _{\text{h}}})}}{2}\left( {1 - \frac{{4r_{\text{w}}^{\text{2}}}}{{{r^2}}} + \frac{{3r_{\text{w}}^{\text{4}}}}{{{r^4}}}} \right)\cos 2\theta + \\ {\delta _1} \cdot {P_{\text{w}}}\frac{{2(1 - \nu ')n\frac{{r_{\text{b}}^2}}{{{r^2}}}}}{{\left[ {1 + (1 - 2\nu ')n} \right]\frac{{r_{\text{b}}^2}}{{r_{\text{a}}^2}} - (1 - n)}} + (1 - {\delta _1}){P_{\text{w}}}\frac{{{r_{\text{w}}}^2}}{{{r^2}}} + \\ {\delta _2}\left[ {\frac{{\alpha (1 - 2\nu )}}{{(1 - \nu )}}\frac{1}{{{r^2}}}\int_{{r_{\text{w}}}}^r {(P(r) - {P_0})r{\text{d}}r - \phi } (P(r) - {P_0})} \right] \text{，}$

(1)

$\begin{gathered} {\sigma _\theta } = \frac{{({\sigma _{\text{H}}} + {\sigma _{\text{h}}})}}{2}\left( {1 + \frac{{r_{\text{w}}^2}}{{{r^2}}}} \right) - \frac{{({\sigma _{\text{H}}} - {\sigma _{\text{h}}})}}{2}\left( {1 + \frac{{3r_{\text{w}}^4}}{{{r^4}}}} \right)\cos 2\theta - \hfill \\ {\text{ }}{\delta _1} \cdot {P_{\text{w}}}\frac{{2(1 - \nu ')n\frac{{r_{\text{b}}^2}}{{{r^2}}}}}{{\left[ {1 + (1 - 2\nu ')n} \right]\frac{{r_{\text{b}}^2}}{{r_{\text{a}}^2}} - (1 - n)}} - (1 - {\delta _1}){P_{\text{w}}}\frac{{r_{\text{w}}^2}}{{{r^2}}} - \hfill \\ \end{gathered} \\ {\delta _2}\left[ {\frac{{\alpha (1 - 2\nu )}}{{(1 - \nu )}}\frac{1}{{{r^2}}}\int_{{r_{\text{w}}}}^r {(P(r) - {P_0})r{\text{d}}r} } \right. -\\ \left. {\left( {\alpha \frac{{1 - 2\nu }}{{1 - \nu }} - \phi } \right)(P(r) - {P_0})} \right] \text{，}$

(2)

${\sigma _z} = {\sigma _{\text{v}}} - 2\nu ({\sigma _{\text{H}}} - {\sigma _{\text{h}}})\cos 2\theta + \\ \left( {\alpha \frac{{1 - 2\nu }}{{1 - \nu }} - \phi } \right)(P(r) - {P_0}) \text{，}$

(3)

${\tau _{r\theta }} = - \frac{{({\sigma _{\text{H}}} - {\sigma _{\text{h}}})}}{2}\left( {1 + \frac{{2r_{\text{w}}^{\text{2}}}}{{{r^2}}} - \frac{{3r_{\text{w}}^{\text{4}}}}{{{r^4}}}} \right)\sin 2\theta \text{，}$

(4)

${\tau _{r\theta }} = {\tau _{z\theta }} = 0 。$

(5)

式中， ${\sigma _r}$ ， ${\sigma _\theta }$ ， ${\sigma _z}$ 为井筒围岩的三向正应力， ${\tau _{r\theta }}$ ， ${\tau _{rz}}$ ， ${\tau _{z\theta }}$ 为剪切应力， ${r_{\text{w}}}$ 为井筒半径， $\theta$ 是与 ${\sigma _{\text{H}}}$ 的极坐标角， $\nu$ 是泊松比， $n = E(1 + \nu ')/(E'(1 + \nu ))$ 为围岩和井筒剪切模量之比， $E$ 和 $E'$ 分别为围岩和井筒的弹性模量， $\nu$ 和 $\nu '$ 分别为围岩和井筒的泊松比， ${r_{\text{a}}}$ 和 ${r_{\text{b}}}$ 为井筒内径和外径， ${P_{\text{0}}}$ 为初始流体压力， $P{\text{(}}r{\text{)}}$ 为实际流体压力， $\phi$ 为有效孔隙度， $\alpha$ 为有效应力系数。 ${\delta _1}$ 为井筒与井壁胶结状态系数，胶结不好时取0，胶结良好时取1； ${\delta _2}$ 为井壁的渗透性系数，不渗透时取0，可渗透时取1。当井筒与井壁胶结良好时，井筒围岩不考虑渗流应力。

1.3 射孔围岩应力分布

假设射孔孔眼和井眼正交，力学模型及坐标转换如图 2所示，基于井筒围岩应力可得到射孔围岩应力：

图 2 射孔围岩力学模型及坐标转换

Figure 2. Mechanical model and coordinate transformation

下载: 全尺寸图片幻灯片

$\begin{gathered} {\sigma _s} = \frac{{{\text{(}}{\sigma _z} + {\sigma _\theta }{\text{)}}}}{2}\left( {1 - \frac{{s_{\text{w}}^{\text{2}}}}{{{s^2}}}} \right) + \hfill \\ {\text{ }}\frac{{({\sigma _\theta } - {\sigma _z})}}{2}\left( {1 - \frac{{4s_{\text{w}}^{\text{2}}}}{{{s^2}}} + \frac{{3s_{\text{w}}^{\text{4}}}}{{{s^4}}}} \right)\cos 2\varphi + {P_{\text{w}}}\frac{{s_{\text{w}}^{\text{2}}}}{{{s^2}}} + \hfill \\ \end{gathered}\\ \frac{{\alpha (1 - 2\nu )}}{{(1 - \nu )}}\frac{1}{{{r^2}}}\int_{{r_{\text{w}}}}^r {(P(r) - {P_0})r{\text{d}}r - \phi } (P(r) - {P_0}) \text{，}$

(6)

${\sigma _\varphi } = \frac{{({\sigma _z} + {\sigma _\theta })}}{2}\left( {1 + \frac{{s_{\text{w}}^{\text{2}}}}{{{s^2}}}} \right) - \frac{{({\sigma _\theta } - {\sigma _z})}}{2}\left( {1 + \frac{{3s_{\text{w}}^{\text{4}}}}{{{s^4}}}} \right)\cos 2\varphi -\\ {\text{ }}{P_{\text{w}}}\frac{{s_{\text{w}}^{\text{2}}}}{{{s^2}}} - \frac{{\alpha (1 - 2\nu )}}{{(1 - \nu )}}\frac{1}{{{r^2}}}\int_{{r_{\text{w}}}}^r {(P(r) - {P_0})r{\text{d}}r + } {\text{ }}\\ \left( {\alpha \frac{{1 - 2\nu }}{{1 - \nu }} - \phi } \right)(P(r) - {P_0}) \text{，}$

(7)

${\sigma _{zz}} = {\sigma _r} - 2\nu ({\sigma _\theta } - {\sigma _z})\frac{{s_{\text{w}}^{\text{2}}}}{{{s^2}}}\cos 2\theta + \\ \left( {\alpha \frac{{1 - 2\nu }}{{1 - \nu }} - \phi } \right)(P(r) - {P_0}) \text{，}$

(8)

${\tau _{zz\varphi }} = 2{\tau _{r\theta }}\sin \varphi \text{，}$

(9)

${\tau _{szz}} = {\tau _{s\varphi }} = 0 。$

(10)

式中， ${\sigma _s}$ ， ${\sigma _\varphi }$ 分别为射孔横截面上的径向应力和周向应力， ${\sigma _{zz}}$ 为轴向应力， ${\tau _{zz\varphi }}$ ， ${\tau _{szz}}$ ， ${\tau _{s\varphi }}$ 为剪切应力； ${s_{\text{w}}}$ 是射孔半径； $\varphi$ 是 ${\sigma _s}$ 与 ${\sigma _\theta }$ 方向的夹角。

2. 射孔围岩流体压力

2.1 非Darcy渗流方程

射孔围岩不是完全刚性的，流体压力改变引起有效应力变化，导致围岩渗透率 $K$ 和孔隙度 $\phi$ 呈现敏感性。流体在多孔介质中的渗流可用Darcy定律描述；对于非Darcy渗流(雷诺数Re > 10^[23])，Forchheimer ^[24]公式被用于代替Darcy定律：

$\left. \begin{array}{l}-\frac{\partial P}{\partial x}=\frac{\mu }{{K}_{x}}{v}_{x}+\beta \rho {v}_{x}{}^{2}\text{ }\text{，}\hfill \\ -\frac{\partial P}{\partial y}=\frac{\mu }{{K}_{y}}{v}_{y}+\beta \rho {v}_{y}{}^{2}\text{ }\text{，}\hfill \\ -\frac{\partial P}{\partial z}=\frac{\mu }{{K}_{z}}{v}_{z}+\beta \rho {v}_{z}{}^{2}\text{ }\text{，}\hfill \end{array} \right\}$

(11)

式中， ${v_x}$ ， ${v_y}$ ， ${v_z}$ 分别为流体的三向渗流速度， ${K_x}$ ， ${K_y}$ ， ${K_z}$ 分别为介质的三向渗透率， $\rho$ 为流体密度， $\mu$ 为流体黏度。该方程比Darcy定律多一个渗流速度的二次方项，代表惯性力影响， $\beta$ 被称为惯性系数。本文采用适用于致密和松散的砂岩惯性系数公式^[25]：

$\beta {\text{ = }}\frac{{1.59 \times {{10}^3}}}{{\sqrt K {\phi ^{5.5}}}}\text{，}$

(12)

式中，渗透率K的单位为mD。

为了与Darcy定律对比分析，引入非Darcy渗透率 ${K_N}$ ：

${K_N}{\text{ = }}\frac{K}{{1 + \frac{{\beta \rho Kv}}{\mu }}}。$

(13)

式(11)可简化为

$\left. \begin{array}{l}{v}_{x}=-\frac{{K}_{Nx}}{\mu }\frac{\partial P}{\partial x}\text{ }\text{，}\hfill \\ {v}_{y}=-\frac{{K}_{Ny}}{\mu }\frac{\partial P}{\partial y}\text{ }\text{，}\hfill \\ {v}_{z}=-\frac{{K}_{Nz}}{\mu }\frac{\partial P}{\partial z}\text{ }。\hfill \end{array} \right\}$

(14)

2.2 渗流偏微分方程

储层围岩中的流体压力可通过渗流偏微分方程得到；将围岩中的流体视为单相微可压缩牛顿流体，在等温无源非稳态渗流下的连续性方程可表示为

$\frac{{\partial (\rho \phi )}}{{\partial t}} + \nabla \cdot (\rho \overrightarrow v ) = 0。$

(15)

考虑Darcy定律和非Darcy定律，认为流体黏度 $\mu$ 与空间位置无关，得到二维非稳态渗流偏微分方程：

${c_{\text{t}}}\mu \phi \frac{{\partial P}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {{K_x}\frac{{\partial P}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{K_y}\frac{{\partial P}}{{\partial y}}} \right)\text{，}$

(16)

${c_{\text{t}}}\mu \phi \frac{{\partial P}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {{K_{Nx}}\frac{{\partial P}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{K_{Ny}}\frac{{\partial P}}{{\partial y}}} \right)。$

(17)

式中， ${c_{\text{t}}} = {c_{\text{f}}} + {c_\phi }$ 为综合压缩系数， ${c_{\text{f}}}$ 为流体压缩系数， ${c_\phi }$ 为围岩孔隙压缩系数。渗流开始时，认为储层围岩内任意位置的流体压力均为初始孔隙压力，即初始条件为 $P(r, 0) = {P_0}$ ；其边界条件如表 1所示。

表 1 渗流分析的边界条件

Table 1. Boundary conditions of fluid flow analysis

外边界	自由边界	$P({r_{\text{e}}}, t) = {\text{0}}$
内边界	射孔壁和井壁压力	$P({r_{\text{w}}}, t) = {P_{\text{w}}}(t){\text{ }}, {\text{ }}P({s_{\text{w}}}, t) = {P_{\text{w}}}(t)$

下载: 导出CSV

| 显示表格

2.3 渗透率和孔隙度的应力敏感性

渗透率和孔隙度的应力敏感性直接影响围岩内流体压力。现有研究结果表明围岩渗透率随有效应力增加呈指数降低^[26-27]：

$K = {K_0}{{\text{e}}^{ - M\sigma '}} = {K_0}{{\text{e}}^{ - M\left( {\sigma - \alpha P} \right)}}\text{，}$

(18)

式中， ${K_0}$ 为初始渗透率， $M$ 为渗透率应力敏感性系数， $\sigma '$ 为Terzaghi有效正应力， $P$ 为孔隙压力， $\alpha$ 为Biot’s系数。基于多孔介质双重有效应力概念，孔隙度与 $\alpha$ 紧密相关^[28]：

$\left. \begin{array}{l}\alpha ={\left(\frac{\varphi }{{\varphi }_{\text{c}}}\right)}^{\frac{2}{3}}\varphi \le {\varphi }_{\text{c}}\text{ }\text{，}\hfill \\ \alpha =1\text{ }\varphi > {\varphi }_{\text{c}}\text{ }\text{，}\hfill\end{array} \right\}$

(19)

式中， ${\phi _{\text{c}}}$ 为多孔介质材料的临界孔隙度。Nur^[29]基于大量室内试验研究，给出了 ${\phi _{\text{c}}}$ 的建议值为0.4。基于Kozeny-Carman方程^[30]，得到渗透率和孔隙度的关系：

$K = \frac{\phi }{{{k_2}{S_{\text{P}}}^2}}。$

(20)

式中 ${S_{\text{P}}}$ 为固体颗粒表面积与孔隙体积的比值；水力压裂过程中，假设注入流体后围岩变形为小变形， ${S_{\text{P}}}$ 认为不变。 ${k_2}$ 为常数，由此得到渗透率 $K$ 与孔隙度 $\phi$ 成正比：

$\phi = {\phi _0}{{\text{e}}^{ - M\sigma '}} = {\phi _0}{{\text{e}}^{ - M\left( {\sigma - \alpha P} \right)}}。$

(21)

2.4 破裂力学准则

Hossain等^[31]认为最大拉应力准则对于预测破裂压力是准确的。由于井筒和射孔相交的双重应力集中，最大拉应力发生在近交点围岩区域内；当水平方向存在主应力差，射孔方位与最大水平应力方向垂直时，破裂发生在井筒围岩最大水平应力处。针对射孔围岩应力分布，需采用复合应力理论先求出最大拉应力，即第三主应力 ${\sigma _3}$ 对于角度 $\varphi$ 的极值 ${\sigma _{3, \min }}$ ，然后利用最大拉应力准则(22)判定围岩破裂和位置^[32]。

${\sigma _{3, \min }} - \alpha P = - {\sigma _{\text{t}}} \text{，}$

(22)

随着射孔方位角 $\theta$ 的增大， ${\sigma _{3, \min }}$ 逐渐增大，破裂压力亦逐渐增大；当井筒与井壁胶结不好时，最大水平应力方向井壁处围岩应力不会随射孔方位角 $\theta$ 变化；当射孔方位角 $\theta$ 增大时，最大水平应力方向的井壁处围岩会首先达到破裂条件，破裂准则也应该包括如下判据：

$({\sigma _{\text{H}}} + {\sigma _{\text{h}}}) - 2({\sigma _{\text{H}}} - {\sigma _{\text{h}}}) - {P_{\text{w}}} - {\delta _2}\left( {\alpha \frac{{1 - 2\nu }}{{1 - \nu }} - \phi } \right)P - \alpha P = - {\sigma _{\text{t}}} 。$

(23)

3. 水力压裂破裂数值模拟方法

3.1 有限容积法空间区域及控制方程离散化

流体在围岩内渗流属非稳态过程，采用有限容积法(FVM)开展数值模拟计算^[33]。首先，将计算区域离散为互不重叠的控制容积，且使每个网格节点周围有一个控制容积(U，W，E，S，N点)。采用单元中心法对整个计算区域开展离散，即单元位于子区域中心，此时子区域就是控制容积(N₁，N₂，N₃，N₄点组成的阴影区域)，且网格线与界面线(线段w，e，s，n)重合，故采用单元代表整个控制容积上的力学属性(图 3)。其次，将守恒型的控制方程(微分方程)对每个控制容积积分，离散化后转化成各节点上的离散方程组，离散方程中的未知数就是节点上因变量的数值。

图 3 网格单元中心法离散

Figure 3. Discrete mesh by center element method

下载: 全尺寸图片幻灯片

控制方程积分和离散后，需要假定被求函数的型线，求出控制容积积分。本模型选用分段线性分布型线，即假定界面线上被求函数本身及其一阶导数(包括对时间和空间)均在节点间的线性分布上。

3.2 数值计算模型及参数

取垂直井，考虑对称布置的两个射孔，射孔方位角为 $\theta$ ，选取井筒和射孔围岩的四分之一( $\varphi$ =90°)构建模型(图 4)，压裂液体为水，注入流量为15 mL/min，泵容量为200 mL；模型计算参数见表 2。

图 4 模型及网格划分示意图

Figure 4. Meshing of numerical model

下载: 全尺寸图片幻灯片

表 2 FVM模拟计算模型参数

Table 2. Model parameters of numerical simulation in FVM

井筒外径/mm	井筒内径/mm	射孔直径/mm	射孔长度 /mm	模型尺寸/(mm×mm)	网格尺寸/(mm×mm)
6	4	1	4	50×50	0.5×0.5
最大水平主应力/MPa	最小水平主应力/MPa	垂向主应力 /MPa	围岩泊松比	井筒弹性模量 /GPa	井筒泊松比
8.53	6.57	9.85	0.25	200	0.15
弹性模量/MPa	抗拉强度/MPa	初始孔隙度	应力敏感性模数M	初始渗透率 /10^-17m²
368.79	0.715	0.124	0.3	27.194

下载: 导出CSV

| 显示表格

求解过程中，该数值方法亦存在一定的假设；①模型为二维渗流问题，控制容积界面线上的渗流是均匀的，即渗透速度在界面线上处处相同；②模型左边界和下边界为对称边界，对称边界不发生垂直于边界的渗流；射孔井壁存在渗流；③应力边界在模型无穷远处。

3.3 二维非稳态渗流方程离散

渗流偏微分方程是由控制容积的质量守恒推导的连续性方程，对渗流控制方程离散，即把渗流偏微分方程转换为控制容积节点上的离散方程组(图 5)。

图 5 控制容积内质量守恒示意图

Figure 5. Mass conservation in controlling volume

下载: 全尺寸图片幻灯片

时间间隔[ $t, t + \Delta t$ ]内，在控制容积U上(图 6)，对渗流偏微分方程(16)积分：

$\int_s^n {\int_w^e {\int_t^{t + \Delta t} {{c_{\text{t}}}\mu \phi } } } \frac{{\partial P}}{{\partial t}}{\text{d}}x{\text{d}}y{\text{d}}t = \int_t^{t + \Delta t} {\int_s^n {\int_w^e {\frac{\partial }{{\partial x}}\left( {{K_x}\frac{{\partial P}}{{\partial x}}} \right)} } } {\text{d}}x{\text{d}}y{\text{d}}t +\\ \int_t^{t + \Delta t} {\int_s^n {\int_w^e {\frac{\partial }{{\partial y}}\left( {{K_y}\frac{{\partial P}}{{\partial y}}} \right)} } } {\text{d}}x{\text{d}}y{\text{d}}t \text{，}$

(24)

图 6 控制容积示意图

Figure 6. Schematic diagram of controlling volume

下载: 全尺寸图片幻灯片

左侧非稳态项可直接离散为

$\int_s^n {\int_w^e {\int_t^{t + \Delta t} {{c_{\text{t}}}\mu \phi } } } \frac{{\partial P}}{{\partial t}}{\text{d}}x{\text{d}}y{\text{d}}t = {(\mu {\phi ^t}{c_{\text{t}}})_U} \times (P_U^{t + \Delta t} - P_U^t)\Delta x\Delta y \text{，}$

(25)

右侧扩散项可离散为

$\int_t^{t + \Delta t} {\int_s^n {\int_w^e {\frac{\partial }{{\partial x}}\left( {{K_x}\frac{{\partial P}}{{\partial x}}} \right)} } } {\text{d}}x{\text{d}}y{\text{d}}t + \int_t^{t + \Delta t} {\int_s^n {\int_w^e {\frac{\partial }{{\partial y}}\left( {{K_y}\frac{{\partial P}}{{\partial y}}} \right)} } } {\text{d}}x{\text{d}}y{\text{d}}t \\= \left[ {{K_{xe}}\frac{{{P_E} - {P_U}}}{{{{(\delta x)}_e}}} - {K_{xw}}\frac{{{P_U} - {P_W}}}{{{{(\delta x)}_w}}}} \right]\Delta y\Delta t +\\ \left[ {{K_{yn}}\frac{{{P_N} - {P_U}}}{{{{(\delta y)}_n}}} - {K_{ys}}\frac{{{P_U} - {P_S}}}{{{{(\delta y)}_s}}}} \right]\Delta x\Delta t \text{，}$

(26)

由此得到离散方程：

$\left. \begin{array}{l}{a}_{U}{P}_{U}={a}_{E}{P}_{E}+{a}_{W}{P}_{W}+{a}_{N}{P}_{N}+{a}_{S}{P}_{S}+{a}_{U}^{t}{P}_{U}^{t}\text{ }\text{，}\\ {a}_{U}={a}_{E}+{a}_{W}+{a}_{N}+{a}_{S}+{a}_{U}^{t}\text{ }\text{，} \end{array} \right\}$

(27)

式(27)中离散方程系数计算方法如下:

${a}_{E}=\frac{\Delta y}{{(\delta x)}_{e}/{K}_{xe}}，{a}_{W}=\frac{\Delta y}{{(\delta x)}_{w}/{K}_{xw}}，{a}_{N}=\frac{\Delta x}{{(\delta y)}_{n}/{K}_{yn}}\\ {a}_{S}=\frac{\Delta x}{{(\delta y)}_{s}/{K}_{ys}}\text{，}{a}_{U}^{t}=\frac{{(\mu {\varphi }^{t}{c}_{\text{t}})}_{U}\Delta x\Delta y}{\Delta t} 。$

(28)

3.4 应力方程的迭代

对于每个控制容积U()，其应力状态 $\sigma _U^{t + \Delta t}$ 是关于井内流体压力 ${P_{\text{w}}}{\text{(}}t{\text{)}}$ 、孔隙度 ${\phi ^t}$ 和流体压力 $P_U^t$ 的函数，可直接参与每一时间步的迭代计算：

${\sigma }_{U}^{t+\Delta t}=f({P}_{\text{w}}(t)\text{，}{\varphi }^{t}\text{，}{P}_{U}^{t}) 。$

(29)

为耦合围岩应力分布与流体渗流状态，需修正渗透率和孔隙度与有效应力的关系：

$K_U^{t + \Delta t} = {K_0}{{\text{e}}^{ - M\left( {\sigma _U^t - \alpha _U^tP_U^t} \right)}}\text{，}$

(30)

$\phi _U^{t + \Delta t} = {\phi _0}{{\text{e}}^{ - M\left( {\sigma _U^t - \alpha _U^tP_U^t} \right)}}\text{，}$

(31)

$\left. \begin{array}{l}{\alpha }^{t+\Delta t}={\left(\frac{{\phi }^{t+\Delta t}}{{\phi }_{\text{c}}}\right)}^{\frac{2}{3}}~~~~{\phi }^{t+\Delta t}\le {\phi }_{\text{c}}\text{ }\text{，}\hfill \\ \alpha =1~~~~{\phi }^{t+\Delta t} > {\phi }_{\text{c}}\text{ }。\hfill \end{array} \right\}$

(32)

且考虑非Darcy渗透时，应对渗透率予以非Darcy修正：

$v_U^t = \frac{{K_{NU}^t}}{\mu }(P_{U - 1}^t - P_U^t) \text{，}$

(33)

${v}_{U}^{t} > {v}_{L}\text{ }\text{，}Re > 10 。$

(34)

$\left. \begin{array}{l}{\beta }_{U}^{t+\Delta t}=\frac{1.59\times {10}^{3}}{\sqrt{{K}_{U}^{t}}{({K}_{U}^{t})}^{5.5}}\text{ }\text{，}\\ {K}_{NU}^{t+\Delta t}=\frac{{K}_{U}^{t}}{1+\frac{{\beta }_{U}^{t+\Delta t}\rho {K}_{U}^{t}{v}_{U}^{t}}{\mu }}\text{ }\text{，} \end{array} \right\}$

(35)

式中， $v_U^t$ 为控制单元流速， ${v_L}$ 为临界流速，可由 ${v_L} = Re{\phi _{\text{c}}}\mu /(\rho d)$ 得到，其中d为岩石孔隙直径。

3.5 边界条件及初始条件

针对上述离散方程，须考虑相应的边界条件予以求解；包括对称边界条件、内边界条件和外边界条件。其中对称边界处不发生垂直于边界的渗流，故边界节点系数满足： ${a_{{\text{sym}}}} = 0$ ；针对内边界，由于井筒和射孔内节点的整个控制容积内流体压力均为 ${P_{\text{w}}}(t)$ ，且直接作用在界面线上(图 7)。

图 7 内边界网格系统示意图

Figure 7. Mesh system of inner boundary

下载: 全尺寸图片幻灯片

与内边界相邻节点的 ${{\text{(}}\delta x{\text{)}}_w}$ 或 ${{\text{(}}\delta y{\text{)}}_s}$ 为原来的一半，其系数变为

${a}_{\text{In}, W}=\frac{2\Delta y}{{(\delta x)}_{\text{In}, }{}_{w}/{K}_{\text{In}, }{}_{xw}}\text{，}{a}_{\text{In}, S}=\frac{2\Delta x}{{(\delta y)}_{\text{In}, s}/{K}_{\text{In}, ys}} 。$

(36)

外边界节点实际为宏观内部节点，故可补充边界节点代数方程对其处理：

${P}_{M}=\frac{{P}_{M+1}+{P}_{M-1}}{2}\text{，}{a}_{M+1}={a}_{M} 。$

(37)

针对初始条件，当 $t = 0$ ， ${\phi ^t} = {\phi ^0}$ ，那么对于确定的控制容积尺寸， $a_U^t = a_U^0$ , $P_U^t = P_U^0 = {P_0}$ ；对于井内和射孔内流体压力始终满足： $P({r_{\text{w}}}, t) = {P_{\text{w}}}(t)$ ， $P({s_{\text{w}}}, t) =$ ${P_{\text{w}}}{\text{(}}t{\text{)}}$ ，且 ${P_{\text{w}}}(0) = 0$ 。

3.6 数值模拟迭代流程

根据上述离散方程、边界条件和初始条件，结合破裂准则，可迭代计算出井内流体压力、围岩流体压力和应力，由此得到破裂时间和破裂位置。每一时间步的具体迭代流程如下：

(1) 输入围岩力学参数及流体参数，计算临界流速；判断流体是否符合Darcy定律，进而修正非Darcy渗透率 ${K_N}$ 。

(2) 考虑参数应力敏感性，修正渗透率 ${K^t}$ 和孔隙度 ${\phi ^t}$ (若不考虑应力敏感性，可跳过参数修正)；基于井筒内流体压力 ${P_{\text{w}}}$ ，计算围岩流体压力 ${P^t}$ ，得到渗流场分布。

(3) 基于流体压力 ${P^t}$ ，计算围岩应力分布，判断第三有效主应力和最小有效周向应力是否满足破裂条件。

(4) 若满足破裂条件，则输出破裂压力和破裂时间等结果，结束计算；若不满足则重复迭代计算，直至满足破裂条件。

4. 数值模拟结果

4.1 模型验证

考虑达西渗流的情况下，基于广泛使用的Hubbert-Willis(H-W)和Haimson-Fairhurst(H-F)解析模型对本文数值模拟方法予以验证。H-W模型不考虑岩石的渗透性，其结果可作为验证的上限，计算公式如下^[8]：

${p_{\text{b}}} = 3{\sigma _{\text{h}}} - {\sigma _{\text{H}}} + {\sigma _{\text{t}}} - {P_0} 。$

(38)

式中， ${p_{\text{b}}}$ 为破裂压力。H-W模型计算得到的破裂压力为11.90 MPa。

H-F模型考虑岩石的高渗透性，结果可作为验证的下限，计算公式如下^[9]：

${p_{\text{b}}} = \frac{{3{\sigma _{\text{h}}} - {\sigma _{\text{H}}} + {\sigma _{\text{t}}} - 2\eta {P_0}}}{{2{\text{(}}1 - \eta {\text{)}}}} \text{，}$

(39)

$\eta = \frac{{{\phi _0}(1 - 2\nu )}}{{2(1 - \nu )}} 。$

(40)

式中， ${\phi _0}$ 为储层初始孔隙度， $\eta$ 为模型参数。相同条件下，H-F模型计算得到的破裂压力为6.20 MPa。

在射孔方位角为0°，考虑井筒与井壁胶结良好、井筒与井壁胶结不好且不考虑井壁渗透、井筒与井壁胶结不好且考虑井壁渗透条件下，本文数值模拟得到的破裂压力分别为11.56，6.25，6.74 MPa。对比上述结果，当井筒与井壁胶结良好，流体经射孔直接作用于围岩上，破裂压力计算结果接近于H-W模型；当井筒与井壁胶结不好时，井内流体易进入井壁裂隙中，流体直接作用在井壁围岩上，导致破裂压力的降低，计算结果接近于H-F模型。且井筒与井壁胶结不好时，考虑井壁渗透条件下计算得到的破裂压力高于不考虑井壁渗透；主要原因在于井壁渗流条件下考虑了渗透作用引起的附加压应力。3种条件下得到的破裂压力计算结果均处于H-W模型和H-F模型破裂压力计算结果之间，验证了模型的正确性。

4.2 井筒与井壁胶结状态的影响

不考虑渗透率和孔隙度应力敏感性，对井筒与井壁不同胶结状态条件下的射孔围岩水力压裂破裂压力、破裂时间和破裂位置开展模拟研究；主要考虑井与壁胶结良好、井与壁胶结不好但井壁不可渗、井与壁胶结不好且井壁可渗3种工况。3种工况下，射孔围岩破裂压力和破裂时间随射孔方位角增大而增加(图 8)；破裂压力越大，破裂需要的流体压力越高，破裂时间越长。当射孔方位角小于45°，破裂压力和破裂时间增长幅度较大，当射孔方位角大于45°，破裂压力和破裂时间增长速度趋于稳定，由此认为水力压裂宜采用略低的射孔方位角。

图 8 不同射孔方位角下射孔围岩破裂规律

Figure 8. Breakdown of perforated surrounding rock under different perforation azimuths

下载: 全尺寸图片幻灯片

井与壁胶结良好时的破裂压力最大，井筒内流体压力对围岩的影响可以忽略，破裂位置位于射孔顶端。井筒与井壁胶结不好时，不考虑井壁渗透作用，流体通过围岩裂隙直接作用于井壁上，与射孔内的流体在井筒与射孔相交处形成双重应力集中，故在低射孔方位角下，破裂位置会在射孔末端；随着射孔方位角的增大，最大水平应力方向近井处的孔隙流体压力会逐渐减小，破裂位置开始转向最大水平应力方向的井壁处，故高射孔方位角时破裂压力呈现增大趋势。井筒与井壁胶结不好时，考虑井壁渗透作用，最大水平应力方向井壁处的破裂条件与射孔方向角无关，破裂均发生在最大水平主应力方向井壁处；且井壁渗透导致井内流体压力增长缓慢，渗透作用有效延迟了围岩破裂时间。

此外，射孔方位角0°条件下，井筒与井壁胶结良好时，流体通过射孔进入地层，流体压力呈关于射孔段的椭圆形分布(图 9)；随着压裂时间的增长，流体压力分布沿椭圆形扩展，渗流影响范围增加。当井筒与井壁胶结不好且不考虑井壁渗透时，虽然最大流体压力不同，但渗流规律与井壁胶结良好时基本一致。当井筒与井壁胶结不好且考虑井壁渗透时，流体渗透面积增加，不仅在射孔处产生渗流，井壁处亦产生渗流，故流体压力分布是沿井眼向外扩散的，渗流影响范围大幅度增加(图 9)。

图 9 射孔围岩破裂时流体压力分布图

Figure 9. Distribution of fluid pressure during breakdown of perforated surrounding rock

下载: 全尺寸图片幻灯片

监测不同压裂时刻典型截面(图 4，截面1—1，坐标：y=0.5 mm)流体压力的时空演化(图 10)。流体进入射孔围岩后，流体压力大幅度降低；当达到破裂时间t₀，射孔围岩流体压力显著增加，直接影响射孔围岩的破裂。当井筒与井壁胶结不好且考虑井壁渗透时，流体压力时空演化曲线较井筒与井壁胶结良好时平缓，印证了井筒与井壁胶结不好且考虑井壁渗透时渗流影响范围增加显著。

图 10 不同压裂时刻典型截面流体压力时空演化曲线

Figure 10. Time-space evolution curves of fluid pressure of typical section

下载: 全尺寸图片幻灯片

4.3 应力敏感性的影响

渗流影响射孔围岩水力压裂的破裂压力和破裂时间，考虑渗透率和孔隙度应力敏感性的条件下，井筒与井壁胶结良好时，破裂压力和破裂时间均随射孔方位角的增大而增加()。但考虑应力敏感性时，围岩破裂压力和破裂时间均小于不考虑应力敏感性的状态；且考虑应力敏感性时，随着射孔方位角增大，破裂压力增长幅度降低显著，射孔方位角越大，考虑应力敏感性引起的围岩破裂压力和破裂时间降低越明显。其原因在于围岩流体压力增大，有效应力系数 $\alpha$ 增加，井壁围岩有效周向应力减小，产生应力集中，直至达到抗拉强度准则开始破裂，导致破裂压力和破裂时间的降低。井筒与井壁胶结不好且不考虑井壁渗透作用下，应力敏感性对不同射孔方位角下围岩破裂压力和破裂时间基本不存在影响。井筒与井壁胶结不好且考虑井壁渗透作用下，考虑应力敏感性，围岩破裂压力和破裂时间均降低；且在高射孔方位角下，降低明显(图 12)。

图 11 井筒与井壁胶结良好时破裂规律

Figure 11. Breakdown under good cementation of wellbore with well

下载: 全尺寸图片幻灯片

图 12 井筒与井壁胶结不好时破裂规律

Figure 12. Breakdown under bad cementation of wellbore with well

下载: 全尺寸图片幻灯片

射孔方位角0°条件下，考虑渗透率和孔隙度应力敏感性，井筒与井壁胶结良好时，围岩内流体压力亦呈关于射孔段的椭圆形分布(图 13)。但相对于不考虑应力敏感性，渗流影响范围大幅度增加，原因在于围岩流体压力增大，导致有效应力降低，由于渗透率应力敏感性的存在，其渗透率不断增大，流体渗流范围随之增大，且渗透率演化规律和流体压力分布规律基本一致。井筒与井壁胶结不好且考虑井壁渗透时，流体压力分布规律与不考虑应力敏感性时基本一致，但渗流影响范围扩大，渗透率的演化规律亦与流体压力分布规律基本一致，均沿水平方向变化(图 14)。

图 13 井筒与井壁胶结良好时流体压力及渗透率分布

Figure 13. Distribution of fluid pressure and permeability under good cementation of wellbore with well

下载: 全尺寸图片幻灯片

图 14 井筒与井壁胶结不好且井壁可渗流体压力及渗透率分布

Figure 14. Distribution of fluid pressure and permeability under bad cementation of wellbore with well

下载: 全尺寸图片幻灯片

监测不同压裂时刻典型截面(图 4，截面1—1，坐标：y=0.5 mm)流体压力的时空演化。考虑渗透率和孔隙度应力敏感性，井筒与井壁胶结良好和井筒与井壁胶结不好时(图 15)，由于近井围岩渗透率增大，围岩流体压力时空演化曲线均趋于平缓，流体压力分布更加均匀，导致流体压力的降低速率变缓。

图 15 考虑应力敏感性下典型截面流体压力时空演化曲线

Figure 15. Time-space evolution curves of fluid pressure of typical section considering stress sensitivity

下载: 全尺寸图片幻灯片

5. 结论

考虑流体渗流和围岩渗透率、孔隙度的应力敏感性，构建考虑水力耦合的下射孔围岩水力压裂破裂力学模型的基础上，采用有限容积法提出了相关数值模拟方法；主要结论如下：

(1) 考虑初始地应力和流体渗流对射孔围岩的影响，应用力学叠加原理得到了射孔围岩的应力分布；且提出了射孔围岩水力压裂破裂准则不仅应包括满足围岩破裂的最大拉应力准则，亦应考虑最大水平应力方向井壁处围岩的破裂判据。

(2) 考虑围岩渗透率和孔隙度的应力敏感性，确定了射孔围岩的流体压力；由此构建了水力耦合作用下射孔围岩水力压裂破裂力学模型；并基于有限容积法(FVM)提出了水力耦合作用下射孔围岩水力压裂破裂数值模拟方法。

(3) 基于上述方法开展数值模拟计算，射孔围岩破裂压力和破裂时间随射孔方位角的增大而增加；且井筒与井壁胶结良好时破裂压力大于井筒与井壁胶结不好时；考虑渗透率和孔隙度的应力敏感性，围岩破裂压力和破裂时间均降低；考虑渗透率和孔隙度的应力敏感性导致流体压力分布更加均匀，近井区域流体压力增大，围岩渗透率增加，引起渗流影响范围扩大。

图 1 扫描层分布

Figure 1. Distribution of scanning layers

下载: 全尺寸图片幻灯片

图 2 应力应变关系及不同的扫描次序

Figure 2. Relationship between stress and strain

下载: 全尺寸图片幻灯片

图 3 不同轴压岩样的三维图像

Figure 3. Image of 3D rock specimen

下载: 全尺寸图片幻灯片

图 4 不同加载状况下第20扫描层的CT图像

Figure 4. No.1 to No.8 CT images of 20th layer

下载: 全尺寸图片幻灯片

图 5 不同加载状况下第30扫描层的CT图像

Figure 5. No.1 to No.8 CT images of 30th layer

下载: 全尺寸图片幻灯片

图 6 加载过程中H值、SD值和轴向应变间关系

Figure 6. Relationship among H, SD and axial strain

下载: 全尺寸图片幻灯片

图 7 位置线示意图

Figure 7. Schematic diagram of position line

下载: 全尺寸图片幻灯片

图 8 第30层第一次加载损伤后 $s_{i j}$ 分布曲线

Figure 8. Distribution of $s_{i j}$ after first loading

下载: 全尺寸图片幻灯片

图 9 第20层R_d指标随荷载变化关系

Figure 9. Relationship between R_d and load

下载: 全尺寸图片幻灯片

图 10 第30层R_d指标随荷载变化关系

Figure 10. Relationship between R_d and load

下载: 全尺寸图片幻灯片

图 11 第20扫描层D和d_j的变化规律

Figure 11. Relationship among D, d_j and load

下载: 全尺寸图片幻灯片

图 12 第30扫描层D和d_j的变化规律

Figure 12. Relationship among D, d_j and load

下载: 全尺寸图片幻灯片

表 1 加载过程中的H值和SD值

Table 1 Values of H and SD during loading

扫描号	轴向应力/MPa	轴向应变/%	第20扫描层		第30扫描层
扫描号	轴向应力/MPa	轴向应变/%	H	SD	H	SD
1	0	0	2207.4	37.6	2192.0	40.9
2	5.82	0.25	2208.8	40.7	2191.7	39.9
3	12.71	0.43	2208.8	38.9	2193.9	40.4
4	24.53	0.76	2209.2	40.0	2192.4	41.6
5	32.14	1.01	2207.1	44.2	2190.5	45.7
6	38.24	1.26	2204.2	50.2	2182.7	48.3
7	25.61	1.67	2177.6	123.3	2168.4	70.7
8	24.36	1.83	2172.5	142.3	2165.9	79.8

下载: 导出CSV

表 2 第20扫描层各位置线上R_d值

Table 2 Values of R_d of different position lines

扫描次序	荷载水平/%	轴向应力/MPa	R_d
扫描次序	荷载水平/%	轴向应力/MPa	15°	45°	75°	105°	135°	165°	总和
1	0	0	94.05	98.36	100.15	94.77	94.19	95.25	576.77
2	0.15	5.82	90.10	94.52	93.78	93.20	93.03	90.20	554.83
3	0.33	12.72	85.34	95.27	90.81	86.61	92.91	89.70	540.64
4	0.64	24.54	85.97	98.44	91.60	87.90	93.43	89.29	546.63
5	0.84	32.14	81.45	85.23	89.90	85.99	84.26	84.48	511.31
6	1.00	38.24	86.65	88.37	91.69	84.99	90.73	93.89	536.32

下载: 导出CSV

表 3 第30扫描层各位置线上R_d值

Table 3 Value of R_d of different position lines

扫描次序	荷载水平/%	轴向应力/MPa	R_d
扫描次序	荷载水平/%	轴向应力/MPa	15°	45°	75°	105°	135°	165°	总和
1	0	0	98.26	105.05	99.05	100.75	99.77	99.39	602.27
2	0.15	5.82	93.52	98.68	94.29	95.20	91.28	98.03	571.00
3	0.33	12.72	95.07	85.71	90.34	94.70	88.61	96.91	551.34
4	0.64	24.54	95.27	96.50	90.87	94.29	91.21	94.02	562.16
5	0.84	32.14	85.23	81.14	86.35	88.48	89.80	85.26	516.26
6	1.00	38.24	88.18	82.67	73.86	78.20	90.21	87.12	500.24

下载: 导出CSV

表 4 第20,30扫描层的损伤变量值

Table 4 Values of damage variables for 20th and 30th scanning layers

扫描次序	荷载水平/%	损伤变量D		$d_{j}$
扫描次序	荷载水平/%	第20层	第30	第20层	第30层
第1次	0	0	0	576.77	602.27
第2次	0.15	0.23	0.27	554.83	571.00
第3次	0.33	0.38	0.49	540.64	551.34
第4次	0.64	0.40	0.52	546.63	562.16
第5次	0.84	0.84	0.86	511.31	516.26
第6次	1.00	1	1	536.32	500.24

下载: 导出CSV

参考文献(17)

[1]	李树刚, 陈高峰, 双海清, 等. 加载速率和初始损伤对砂岩能量演化影响的试验研究[J]. 采矿与安全工程学报, 2019, 36(2): 373-380. https://www.cnki.com.cn/Article/CJFDTOTAL-KSYL201902022.htm LI Shu-gang, CHEN Gao-feng, SHANG Hai-qing, et al. Experimental study on effect of loading rate and initial damage on energy evolution of sandstone[J]. Journal of Mining & Safety Engineering, 2019, 36(2): 373-380. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-KSYL201902022.htm
[2]	MEES F, SWENNEN R, VAN Geet M, et al. Applications of X-ray computed tomography in the geosciences[J]. Geological Society, London, Special Publications, 2003, 215(1): 1-6. doi: 10.1144/GSL.SP.2003.215.01.01
[3]	MARTÍNEZ-MARTÍNEZ J, FUSI N, GALIANA-MERINO J J, et al. Ultrasonic and X-ray computed tomography characterization of progressive fracture damage in low-porous carbonate rocks[J]. Engineering Geology, 2016, 200: 47-57. doi: 10.1016/j.enggeo.2015.11.009
[4]	杨更社, 刘慧. 基于CT 图像处理技术的岩石损伤特性研究[J]. 煤炭学报, 2007, 32(5): 463-468. doi: 10.3321/j.issn:0253-9993.2007.05.004 YANG Geng-she, LIU Hui. Study on the rock damage characteristics based on the technique of CT image processing[J]. Journal of China Coal Society, 2007, 32(5): 463-468. (in Chinese) doi: 10.3321/j.issn:0253-9993.2007.05.004
[5]	杨更社, 孙钧, 谢定义. 岩石材料损伤变量与CT数间的关系分析[J]. 力学与实践, 1998, 20(4): 47-49. https://www.cnki.com.cn/Article/CJFDTOTAL-LXYS804.015.htm YANG Geng-she, SUN Jun, XIE Ding-yi. Analysis of the relation between damage variable and CT value of rock material[J]. Mechanics in Engineering, 1998, 20(4): 47-49. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXYS804.015.htm
[6]	张国凯, 李海波, 王明洋, 等. 岩石单轴压缩下损伤表征及演化规律对比研究[J]. 岩土工程学报, 2019, 41(6): 1074-1082. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201906013.htm ZHANG Guo-kai, LI Hai-bo, WANG Ming-yang, et al. Comparative study on damage characterization and damage evolution of rock under uniaxial compression[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(6): 1074-1082. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201906013.htm
[7]	任建喜, 惠兴田. 裂隙岩石单轴压缩损伤扩展细观机理CT分析初探[J]. 岩土力学, 2005, 26(增刊1): 18-23. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2005S1010.htm REN Jian-xi, HUI Xing-tian. Primary study on meso-damage propagation mechanism of cracked-sandstone using computerized tomography under uniaxial compression[J]. Rock and Soil Mechanics, 2005, 26(S1): 18-23. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2005S1010.htm
[8]	REN J X, LUO Y, LIU W G. Application of computerized topography testing technology on studying rock failure mechanism under loading and unloading[J]. Journal of Glaciology and Geocryology, 2002, 24(5): 672-675. doi: 10.3969/j.issn.1000-0240.2002.05.036
[9]	YU Q, YANG S, RANJITH P G, et al. Numerical modeling of jointed rock under compressive loading using X-ray computerized tomography[J]. Rock Mechanics and Rock Engineering, 2016, 49(3): 877-891. doi: 10.1007/s00603-015-0800-4
[10]	曹文贵, 林星涛, 张超, 等. 基于非线性动态强度准则的岩石动态变形过程统计损伤模拟方法[J]. 岩石力学与工程学报, 2017, 36(4): 794-802. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201704003.htm CAO Wen-gui, LIN Xing-tao, ZHANG Chao, et al. A statistical damage simulation method of dynamic deformation process for rocks based on nonlinear dynamic strength criterion[J]. Chinese Journal of Rock Mechanics and Engineering, 2017, 36(4): 794-802. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201704003.htm
[11]	丁卫华, 仵彦卿, 蒲毅彬. 受力岩石密度损伤增量及其数字图像[J]. 西安理工大学学报, 2000, 16(1): 61-64. https://www.cnki.com.cn/Article/CJFDTOTAL-XALD200001013.htm DING Wei-hua, WU Yan-qing, PU Yi-bin. The density damage increment and its digital image of rock in compression[J]. Journal of Xi'an University of Technology, 2000, 16(1): 61-64. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-XALD200001013.htm
[12]	丁卫华, 仵彦卿, 蒲毅彬. 低应变率下岩石内部裂纹演化的X射线CT方法[J]. 岩石力学与工程学报, 2003, 22(11): 1793-1797. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200311007.htm DING Wei-hua, WU Yan-qing, PU Yi-bin. X-ray CT approach on rock-interior crack evolution under low strain rate[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(11): 1793-1797. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200311007.htm
[13]	党发宁, 方建银, 丁卫华. 基于CT的混凝土试样静动力单轴拉伸破坏裂纹分形特征比较研究[J]. 岩石力学与工程学报, 2015(增刊1): 2922-2928. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2015S1041.htm DANG Fa-ning, FANG Jian-yin, DING Wei-hua. Fractal comparison research of fracture of concrete samples under static and dynamic uniaxial tensile using CT[J]. Chinese Journal of Rock Mechanics and Engineering, 2015(S1): 2922-2928. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2015S1041.htm
[14]	易成, 张亮, 陈忠辉, 等. 一种新的描述粗糙表面形貌尺度分维参数的研究[J]. 中国矿业大学学报, 2007, 36(1): 75-80. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGKD200701015.htm YI Cheng, ZHANG Liang, CHEN Zhong-hui, et al. A novel description of roughness surface with a modified fractal index R_d[J]. Journal of China University of Mining & Technology, 2007, 36(1): 75-80. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGKD200701015.htm
[15]	李夕兵, 翁磊, 谢晓锋, 等. 动静载荷作用下含孔洞硬岩损伤演化的核磁共振特性试验研究[J]. 岩石力学与工程学报, 2015, 34(10): 1985-1993. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201510005.htm LI Xi-bing, WENG Lei, XIE Xiao-feng, et al. Study on the degradation of hard rock with a pre-existing opening under static-dynamic loadings using nuclear magnetic resonance technique[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(10): 1985-1993. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201510005.htm
[16]	曾鹏, 刘阳军, 纪洪广, 等. 单轴压缩下粗砂岩临界破坏的多频段声发射耦合判据和前兆识别特征[J]. 岩土工程学报, 2017, 39(3): 509-517. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201703021.htm ZENG Peng, LIU Yang-jun, JI Hong-guang, et al. Coupling criteria and precursor identification characteristics of multi-band acoustic emission of gritstone fracture under uniaxial compression[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(3): 509-517. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201703021.htm
[17]	HOUNSFIELD G N. Computerized transverse axial scanning (tomography) part 1: description of system[J]. The British Journal of Radiology, 1973, 46(552): 1016-1022.