考虑细观形貌特征的岩石裂隙受压时变闭合规律

Study on Time⁃Dependent Closure Behavior of Rock Fractures Subject to Normal Stress

摘要

基于不同高径比微凸体的受压时效性试验，根据赫兹接触理论，得到了不同微凸体弹性模量随时间的衰减规律. 对红砂岩和石灰岩新鲜裂隙面开展了不同法向应力条件下的受压时变闭合试验，结合小波分析法、区域生长算法和参考面法，提出了一种岩石裂隙细观尺度微凸体形貌的识别方法，对比了试验前后微凸体数量、高度和高径比的差异. 建立了考虑微凸体间相互作用的接触力学模型，求解Boussinesq方程并考虑弹性模量随时间的衰减，对两种岩石裂隙开展了逐级增加应力条件下的受压时变闭合计算. 通过对比计算和试验得到的损伤面积和蠕变变形验证了模型的有效性，阐明了微凸体的应变、接触面积和接触应力随时间的演化规律，揭示了不同细观形貌特征的微凸体在裂隙受压时变闭合过程中的关键性作用.

Abstract

This study conducted time⁃dependent compression tests on asperities with different height⁃to⁃radius ratios using ultra⁃hard gypsum. According to Hertz contact theory, the attenuation laws of the elastic modulus of different asperities over time were fitted. Time⁃dependent closure tests were performed on fresh fracture surfaces of red sandstone and limestone under varying normal stresses. By integrating wavelet analysis, region growth algorithms, and the reference surface method, a novel approach was developed for identifying the mesoscale asperity morphology of rock fractures, and compared the differences in the number, height, and height⁃to⁃radius ratio of asperities before and after the experiment. Utilizing Boussinesq's solution, an influence matrix was constructed to account for interactions between asperities. Based on the law of the elastic modulus decaying over time, enabling time⁃dependent closure calculations for different rock fractures under variable stress conditions. This approach precisely analyzes the temporal evolution of strain, contact area, and contact stress for individual asperity, with simulation results matching experimental data in terms of damage area and creep deformation. The study reveals the pivotal role of asperities with distinct mesoscale morphological features in the time⁃dependent closure process of rock fractures under compression.

Graphical abstract

关键词

岩石裂隙 / 微凸体 / 形貌特征 / 时变闭合 / 赫兹接触 / 岩石力学.

Key words

rock fracture / asperity / morphological characteristics / time⁃dependent closure / Hertzian contact / rock mechanics

0 引言

在我国“三深”战略目标的引领下，能源结构转型正加速向地球内部和深部拓展演进，形成了涵盖高放废物深地处置、CO₂地质封存、地热能开发及油气资源高效开采的深地工程体系. 此类工程均存在高地应力环境下岩体裂隙变形问题（Viswanathan et al.， 2022）. 裂隙闭合既可以阻滞核素迁移，又可以在碳封存工程中维持盖层结构完整性（Rohmer et al.， 2016；Chen et al.， 2023）；但在增强型地热系统和页岩气储层改造方面，裂隙闭合导致的渗透率时效衰减严重制约着能源开采效率. 因此，裂隙时变闭合行为是控制深地工程长期安全与开发效益的关键科学问题之一（孙钧，2007；张亮亮和王晓健，2020；李博等，2021； Kang et al.， 2022；宋振宇等，2024）.

在细观尺度上，微凸体指岩石裂隙表面在细观尺度上的天然凸起结构（通常为μm⁃mm量级），其几何特征包括高度、曲率半径和空间分布等参数. 岩石裂隙的时变闭合是由不同受压微凸体时效变形行为及相互作用关系控制的问题. 为便于定量分析和简化计算，学者们将微凸体概化成圆柱状、锥状和抛物面等多种几何形式（Brown and Scholz， 1985；Pyrak⁃Nolte et al.， 2016；Ma et al.， 2023）. 其中，椭球面的变形和受力特性可以用Hertz接触理论等解释，其曲率连续变化特性可避免圆锥模型在顶点处应力无穷大的问题，同时比圆柱模型更能反映微凸体的非对称接触特征，在接触载荷下产生的应力场更接近实际，因此被广泛采用（Ciavarella et al.， 2006；Greenwood， 2006； Tang and Zhang， 2021）.

微凸体受压时效性试验与数值研究最初源于材料科学与摩擦学领域，主要是衡量硬质合金等材料的耐用性问题（Brot et al.， 2008）. 常见的建模方法是将裂隙面视为多个半椭球体的组合体，该几何表征能有效反映裂隙面的波状起伏特性，并有助于用赫兹接触理论分析微凸体的力学特性（Greenwood and Williamson， 1966）. 需注意，微凸体并非像传统圆柱体和立方体等规则试样在不同应力水平下均保持相同的面积，其接触面积与应力并非是恒定的，而是随时间演化（Song et al.， 2025）. 受几何形状影响，微凸体主要通过两种途径调节时间依赖性变形和失效模式：（1）在恒定载荷作用下，应力减小而接触面积增大. 在压缩过程中，高应力区域迅速屈服，使微凸体扁平化，导致接触边缘处出现强烈的应力集中（Chung et al.， 2010；Rostami et al.， 2014）；（2）微凸体越扁平，尖端到最大拉应力点的距离会更大；体积越小，就越容易发生局部失效（Serati et al.， 2015）. 早期描述单个微凸体与时间相关闭合的数学模型是从试验数据中得出的. 通过采用基本初等函数，将应力、应变和时间等物理量与拟合参数联系起来，提出了大量经验模型，如幂律模型，指数模型和Garofalo蠕变方程（Malamut et al.， 2009；Ovcharenko et al.， 2009；Goedecke et al.， 2010）. 后来，学者们基于分形理论和元件组合模型，将多个微凸体之间的相互作用推广到裂隙尺度（Matsuki et al.， 2001；Kang et al.， 2020）.

以往对于粗糙裂隙面的接触力学研究多聚焦于法向受压闭合过程中的应力-变形规律和准静态压痕试验（Li et al.， 2015，2022；Kumamoto et al.， 2017；Alamos et al.， 2022），而忽略了裂隙闭合过程中的时间依赖性（Matsuki et al.， 2001；Kang et al.， 2020；Alamos et al.， 2022），主要是因为该研究存在两个关键问题：一个是缺乏岩石等脆性材料表面微凸体的受压时效性试验；另一个是缺乏精确识别和区分裂隙面上不同形貌微凸体的方法. 在裂隙受压闭合过程中，法向变形通常由少数关键微凸体的应力集中和压碎机制控制（Li et al.， 2015），不同岩石的矿物组成和形成机制差异导致了其微凸体形貌与时变闭合机制的分化. 现有研究证实，微凸体形貌在裂隙闭合过程中发挥着关键控制作用，其破坏模式、变形和受力特性等方面与传统规则试样均存在显著差异. 在恒定荷载作用下，微凸体的接触应力和接触面积长期处于动态平衡，受形貌特征影响，接触半径更小的微凸体尖端接触应力会显著增大（Kumamoto et al.， 2017）. 学者们通常将蠕变过程中岩石强度的衰减归因于弹性模量、强度和黏性等参数随时间增长产生的减低（Matsuki et al.， 2001；Sun， 2007；Song et al.， 2025），这种参数的时效性衰减受局部高应力微凸体的渐进破坏控制. 在识别微凸体形貌方面，统计方法忽略了微凸体空间分布和几何多样性（Ciavarella et al.， 2006），分形理论复杂度高、计算量大（Kang et al.， 2020）；分水岭方法易受采样分辨率和数据噪声影响，稳定性较差（Wen et al.， 2022）.

针对上述问题，首先基于不同高径比微凸体的受压时效性试验结果，根据赫兹接触理论，得到了不同微凸体弹性模量随时间的衰减规律；其次，利用小波分析法、区域生长算法和参考面法，提出了一种岩石裂隙细观尺度微凸体的识别方法；然后，将边界元法与共轭梯度法和快速傅里叶变换结合，采用Boussinesq解建立影响矩阵，提出了考虑微凸体之间相互作用的接触算法，开展了不同法向应力作用下红砂岩和石灰岩裂隙受压时变闭合数值计算，获取了损伤面积和蠕变演化特征；最后，结合试验观测和数值模拟，揭示了力学蠕变下，由于细观尺度微凸体接触面积扩展和应力重分布造成的裂隙时变闭合现象的内在机理.

1 微凸体受压时效变形特性

1.1　微凸体受压时效性试验

裂隙面上的微凸体在细观尺度上呈现出不同的形貌特征. 一般将微凸体定义为具有不同初始高径比λ₀的半椭球体，如式（1）：

{\lambda }_{\mathrm{0}}=\frac{L}{R}，（1）

式中：L为微凸体高度（mm）；R为底面半径（mm）.

采用超硬牙科石膏制作了4种不同高径比（A⁃1：0.5，A⁃2：1.0，A⁃3：1.5，A⁃4：2.0）的微凸体试样，开展了受压时效性试验，试验概况如图1所示，具体试验方案见前序研究（Song et al.， 2025）. 试验发现：高径比较小的微凸体（A⁃1和A⁃2）稳定性较高，初始接触面积较大，分担荷载的能力较强，在恒定法向荷载条件下，法向应力不断减小，接触面积不断扩大，在较长的时间发生缓慢变形，从而实现应力重分布而达到平衡，最终突然脆性断裂而破坏；高径比较大的微凸体（A⁃3和A⁃4）初始接触面积较小，尖端发生应力集中，短时间内无法释放高应力，需要通过迅速增大接触面积来平衡，在蠕变过程中容易发生非线性变形而失稳.

根据赫兹接触理论，接触半径和法向变形的关系如式（2）：

a^{\mathrm{2}}(t)=\omega (t)R^{\mathrm{*}}，（2）

式中：ω（t）为t时刻的法向变形（mm）；R^*为刚性平面半径R₁和微凸体曲率半径R₂的等效半径（mm）；1/R^*=1/R₁+1/R₂，当R₁→∞时，R^*=ρ.

已知应变的定义如式（3）：

\epsilon (t)=\omega (t)/L.（3）

将不同形貌微凸体的尖端概化为不同曲率半径的半椭球体，式（2）可改写为：

a^{\mathrm{2}}(t)=\epsilon (t)\rho L，（4）

式中：ρ为微凸体尖端曲率半径（mm）；本研究中，ρ=5/λ₀.

综合式（2）~（4）得：

\omega \left(t\right)=\frac{F_{\mathrm{n}}L}{E\left(t\right)A\left(t\right)}，（5）

式中：F_n为法向荷载（kN）；A（t）为微凸体在不同时刻的接触面积（mm²）.

最终得到恒定法向荷载作用下，微凸体弹性模量随时间的衰减规律，如式（6）：

E\left(t\right)=\frac{F_{\mathrm{n}}L}{\mathrm{\pi }{\omega }^{\mathrm{2}}\left(t\right)\rho }.（6）

然而，式（6）仅适用于描述单个微凸体在初始弹性阶段的模量衰减，在真实岩石裂隙受压时变闭合过程中，微凸体往往经历弹塑性变形并伴随颗粒重排. 通过开展花岗岩裂隙试样的法向受压时变闭合试验（Matsuki et al.， 2001），发现弹性模量、强度和黏性等参数通常都会随时间的增长而降低，岩石力学性质的劣化实质上是力学参数随时间弱化的结果（Sun， 2007）. 因此，需引入式（7）所示的指数衰减模型来表征材料劣化效应，通过E（∞）和τ两个参数，可同时反映弹性模量的渐进衰减趋势与最终稳定状态.

E(t)=E(\mathrm{\infty })+[E(\mathrm{0})-E(\mathrm{\infty })]{\mathrm{e}}^{-\frac{t}{\tau }}，（7）

式中：E（0）为弹性模量初始值（MPa）；E（∞）为弹性模量随时间增长的稳定值（MPa）；τ为拟合参数.

基于上述理论，开展了不同微凸体在统一法向荷载500 N条件下的蠕变试验，保持时间均为12 h，得到了A⁃1~A⁃4的弹性模量-时间曲线和应变-时间曲线，如图2. 传统的弹性模量测试方法保证了法向应力在接触面上均匀分布，而未考虑几何形状对接触特性的影响. 已有文献表明（Kumamoto et al.， 2017），受尺寸效应影响，尖端接触半径越小，硬度越高. 根据赫兹接触理论，微凸体的接触应力呈半椭球型分布，而高径比增大会显著提升最大接触应力. 这种应力集中使得表观弹性模量E_a 远高于材料实际值，高径比越大，尖端应力水平越高，对颗粒重排和应力重分布的驱动力越强，导致表观弹性模量的衰减速率越快. 在相同法向荷载下，高径比越小的微凸体应变越大、弹性模量衰减程度越小、越接近实际弹性模量. 但高径比过高则会引发接触区失稳（A⁃4），从而触发应变快速增大和弹性模量急剧降低. 在前序研究中也观察到了不同微凸体在统一荷载550 N条件下保持7 h，仅有A⁃4发生了加速蠕变现象（Song et al.， 2025）.

通过拟合图2中A⁃1~A⁃3的试验数据，得到了弹性模量随时间衰减规律可用下式描述：

E_a(t)=A{\mathrm{e}}^{-\frac{t}{\rho }}+B.（8）

采用1stopt中的Levenberg⁃Marquardt迭代算法，得到不同高径比微凸体弹性模量随时间演化的拟合参数A、B，随着微凸体高径比增大，A和B呈指数型递增，如表1所示. 本文中不同形貌微凸体的A和B均通过插值方法获得.

1.2　岩石裂隙受压时变闭合试验

选用红砂岩和石灰岩两种不同类型的沉积岩开展岩石裂隙受压时变闭合试验，试验参数如表2，试验照片和变形曲线如图3所示. 裂隙面尺寸均为10×10 cm. 试验前用三维轮廓扫描仪采集裂隙面点云数据（分辨率0.01 mm）. 将巴西劈裂后的试样下表面保持固定，上表面向右错位5 mm，得到非完全咬合的裂隙面，突出微凸体的控制作用（图3a，3b）. 为避免动力冲击效应对微凸体造成的瞬间破坏，同时充分反映裂隙面微凸体的受压时效变形，以0.1 MPa/min的速率施加法向应力，在1、3、5、7和9 MPa时分别保持72 h（王振等， 2017；Zhang et al.， 2019）. 在裂隙闭合过程中，随着法向应力增大，接触凸体之间的法向接触刚度不断增大，法向变形增幅逐渐减小，裂隙闭合速率逐渐减缓. 相比而言，石灰岩整体更密实、强度高，抵抗变形能力强，因此逐级加载的闭合量Δω小于红砂岩（图3c）.

2 数值模拟方法

2.1　微凸体识别方法

由于岩石裂隙面的分形特征，在较大的微凸体上往往分布着尺寸更小的微凸体. 因此，采用小波分析法对裂隙面进行降噪处理（Zou et al.， 2020），获得了裂隙面上所有峰值点，并编号. 利用区域生长算法，从峰值点向下搜索高程小于中心点的网格点，直到高度梯度出现反转时，停止搜索，从而逐个识别裂隙面上的所有微凸体. 传统的分水岭算法仅适用于初始粗糙度快速评估（Wen et al.， 2022），易将微噪声误判为微凸体导致过分割，或合并或忽略多个微凸体导致欠分割. 因此，本文采用参考面法对识别后的每个微凸体逐一划分参考面，选取微凸体边缘高程最高点，划定水平面，将参考面以上部分概化为半椭球体，并计算出该半椭球体的高径比. 识别流程如图4.

运行上述程序，先识别单个微凸体模型，如图5所示. 对比理论轮廓线和程序识别出的表面轮廓线，确认对光滑微凸体A⁃1~A⁃4的识别精度均达到1×10^-4 mm（图5）. 进一步识别真实岩石裂隙面，红砂岩裂隙的上、下表面分别识别出75个和89个微凸体；石灰岩裂隙的上、下表面分别有67个和79个微凸体. 这些微凸体在裂隙受压过程中并非全部接触，随着时间推移和法向应力增大，产生接触的微凸体数量逐渐增长，在法向应力9 MPa条件下维持72 h后，红砂岩和石灰岩最终发生接触的微凸体分别有60个和33个. 真实裂隙面上微凸体的识别效果如图8所示. 图6a展示了表面降噪处理后微凸体的识别效果；图6b展示了不同高径比微凸体的识别和计算方法，所有微凸体在计算应变时均采用统一基准面；图6c展示了传统分水岭方法的局限性，该法对所有微凸体共用一个全局参考面，通过模拟“浸水”过程划分集水盆地边界，但未针对单个微凸体单独定义参考面；图6d绘制了微凸体在裂隙面上的位置和分布特征.

图7和图8分别为红砂岩和石灰岩裂隙面微凸体的高度分布图和高径比分布图. 经统计，两种岩石微凸体高度多数处于0~0.2 mm范围内，高径比多数处于0.1~0.12范围内. 红砂岩裂隙面微凸体的归一化频率和数量随高度增加呈缓慢下降趋势，微凸体呈现出均匀、分布松散的特点，仅在0.6~0.7 mm区间内出现“间断”，最大高度为0.80 mm，最大高径比为0.27；而石灰岩裂隙面微凸体的归一化频率和数量随高度增加均呈快速下降趋势，微凸体呈现出粗糙、棱角分明的特点，在0.48~0.76 mm较长的区间内出现“间断”，最大高度为0.86 mm，最大高径比为0.30. 相比来说，红砂岩中的黏土矿物使其颗粒间结合更紧密，因此呈现出较小的粗糙度，而石灰岩层理发育明显，易沿层理面破裂，因此呈现出较大的粗糙度. 值得注意的是，高度和高径比较大的微凸体一般体积和强度均较大，它们在分担法向应力和维持裂隙长期稳定方面发挥着关键作用.

2.2　受压裂隙时效变形计算方法

为定量描述岩石裂隙长期受压时变闭合特性，基于最小余能原理，采用边界元法、共轭梯度法和快速傅里叶变换，对法向应力增量下产生的法向变形进行迭代求解. 系统内的弹性应变能减去外力对系统做的功等于系统内的最小势能（Tian et al.， 1996；Kling et al.， 2018；Zou et al.， 2020； Li et al.， 2022），如式（9）.

V^{\mathrm{*}}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \underset{\mathrm{\Omega }}{\overset{}{\iint }}}{\sigma }_{\mathrm{n}}{\overline{u}}_{\mathrm{z}}(x,y)\mathrm{d}x\mathrm{d}y-{\displaystyle \underset{\mathrm{\Omega }}{\overset{}{\iint }}}{\sigma }_{\mathrm{n}}{\overline{u}}_{\mathrm{z}}^{\mathrm{*}}(x,y)\mathrm{d}x\mathrm{d}y,（9）

式中：Ω 为接触区域；{\overline{u}}_{\mathrm{z}}（x，y）为接触区域内的复合表面法向位移；{\overline{u}}_{\mathrm{z}}^{\mathrm{*}}（x，y）为基于几何干涉的假定接触区域内两个接触微凸体的总位移.

利用Boussinesq解关联接触应力与半无限体表面上的位移，考虑微凸体自身的变形以及微凸体之间的相互作用（Tang and Zhang， 2021），模拟裂隙面微凸体的弹性变形和弹塑性变形破坏过程. 接触区域的变形来源于两部分，即微凸体自身的变形和微凸体之间相互作用引起的变形，如式（10）：

{\overline{u}}_{\mathrm{z}}(x,y)=\frac{\mathrm{1}}{\mathrm{\pi }{E}^{\mathrm{*}}}{\displaystyle \underset{\mathrm{\Omega }}{\overset{}{\iint }}}\frac{{\sigma }_{\mathrm{n}}(x\text{'},y\text{'})\mathrm{d}x\text{'}\mathrm{d}y\text{'}}{\sqrt[]{{(x-x\text{'})}^{\mathrm{2}}+{(y-y\text{'})}^{\mathrm{2}}}}.（10）

式（10）离散化后得到式（11）：

({\overline{u}}_{\mathrm{z}}{)}_l={\displaystyle \sum _{l=\mathrm{1}}^M}{C}_{kl}{p}_k，（11）

式中：k和l分别为接触应力位置和表面位移位置两个指标；M是发生几何干涉的初始接触点微凸体的总数；C_kl 是表征应力与裂隙面位移关系的影响系数矩阵.

将弹性接触模型推广到弹塑性接触分析，纳入考虑塑性变形引起的能量耗散，如式（12）：

V^{\mathrm{*}}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \underset{\mathrm{\Omega }}{\overset{}{\iint }}}{\sigma }_{\mathrm{n}}{\overline{u}}_z(x,y)\mathrm{d}x\mathrm{d}y-{\displaystyle \underset{\Omega }{\overset{}{\iint }}}{\sigma }_n[{\overline{u}}_z^{\mathrm{*}}(x,y)-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{u}_z^p(x,y)]\mathrm{d}x\mathrm{d}y，（12）

式中：Δu_z^p 为局部接触应力p达到岩石基质硬度H时；塑性变形区域的复合增量表面位移. 其中，局部接触应力p大于0的区域为接触区域.

在外部恒定应力水平下，对于单个接触的微凸体，随着时间推移，接触区域的弹性模量随时间递减（式8），接触应力和接触面积相互平衡，取其不同时刻的平均接触应力{\overline{\sigma }}_{\mathrm{n}}（t， x， y），得到不同时刻的法向变形{\overline{u}}_{\mathrm{z}}（t， x， y），如式（13）：

{\overline{u}}_{\mathrm{z}}(t,x,y)=\frac{\mathrm{1}}{\mathrm{\pi }{E}^{\mathrm{*}}(t,x\text{'},y\text{'})}{\displaystyle \underset{\mathrm{\Omega }}{\overset{}{\iint }}}\frac{{\overline{\sigma }}_{\mathrm{n}}(t,x\text{'},y\text{'})\mathrm{d}x\text{'}\mathrm{d}y\text{'}}{\sqrt[]{{(x-x\text{'})}^{\mathrm{2}}+{(y-y\text{'})}^{\mathrm{2}}}}.（13）

计算流程如图9所示：

3 模拟结果

分别模拟了红砂岩和石灰岩两种不同岩石裂隙在1、3、5、7、9 MPa下的时变闭合过程，并对比了试验前后损伤区域的差异，如图10和11所示. 在受压初始阶段，只有少数微凸体相互接触；随着时间的推移，接触的微凸体数量越来越多，接触面积不断扩大；主要接触微凸体被压溃，应力重分布并从高径比大的微凸体向高径比低的微凸体传递；裂隙法向刚度增大，闭合速率逐渐减慢. 对比图10和图11，红砂岩裂隙微凸体分布较均匀，面积扩大明显，应力分布更均衡；而石灰岩微凸体硬度高，体积大，抗压强度大，塑性变形能力较弱，快速压碎后面积呈现小幅增加，法向应力主要集中在几个关键微凸体上.

红砂岩和石灰岩裂隙面在不同法向应力下损伤面积增长速率的演化趋势如图12所示. 图中r₊为不同应力水平下72 h后损伤面积占整个裂隙面的比例，r_-为同一应力水平下初始损伤面积占72 h后损伤面积的比例. 如式（14）：

\left\{\begin{array}{l}{r}_{+}=\frac{A_{\mathrm{72}\mathrm{h}}}{A_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\times \mathrm{100}\mathrm{\%}\\ r_{-}=\frac{A_{\mathrm{0}}}{A_{\mathrm{72}\mathrm{h}}}\times \mathrm{100}\mathrm{\%}\end{array}\right.，（14）

式中：A_{72 h}为某一级应力水平下蠕变保持72 h后的损伤面积（mm²）；A_total为裂隙总面积（mm²）；A₀为施加完某一级应力后，未经蠕变时的损伤面积（mm²）.

随着法向应力增大，r₊增长速率逐渐变缓，而r_-增长速率逐渐减小，这是因为接触的微凸体逐渐被压平，法向接触刚度增大所致. 最终，在法向应力为9 MPa条件下，红砂岩和石灰岩的损伤面积分别占整个裂隙面的3.37%和1.56%，分别约306 mm²和145 mm². 扫描试验后裂隙面的三维轮廓，得到红砂岩和石灰岩的损伤面积分别为310 mm²和151 mm²，两者误差分别为1.31%和4.14%，试验前后损伤区域和法向变形基本吻合，如图13和图14所示.

选取裂隙面上高径比差异的典型微凸体并作重点分析，如图15和图16，展示了两种岩石裂隙面微凸体的变形和受力特性随时间的演化规律. 在相同时间内，裂隙表面的不同微凸体可能处于衰减蠕变阶段、稳定蠕变阶段或加速蠕变阶段. 对于红砂岩裂隙面，微凸体体积小、粗糙度低. 高径比较小的微凸体较扁平［如：No.31（λ₀=0.100 1）；No.32（λ₀=0.104 9）］，起伏度较小. 因此，这类微凸体接触时间通常晚于高径比较大的微凸体（No.31和No.32分别在第36小时和第14小时发生接触）. 但它们的稳定性更好，分担应力的水平更高，抵抗变形的能力更强，在相同时间内应变更小；反之，高径比较大的微凸体较尖锐［如：No.13（λ₀=0.203 4），No.36（λ₀=0.194 7）］，体积大、起伏度大，更容易在受压初始阶段就发生接触. 对于石灰岩裂隙面，微凸体强度高、硬度大，高径比大的微凸体在受压初始时就发生接触［No.23（λ₀=0.235 0）和No.70（λ₀=0.221 5）］，大部分微凸体的法向变形和接触面积变化并不明显. 实际上，岩石裂隙的时间依赖性闭合仅取决于少数相互接触的主要微凸体，这些微凸体承担了大部分应力，对裂隙的时间依赖性闭合发挥了决定性作用（Li et al.， 2015）.

图17和图18对比了试验后两种岩石裂隙面微凸体高度和高径比的分布情况. 受压时变闭合后的微凸体数量显著减少，相比裂隙面初始状态（图6），可识别的微凸体总数量分别为121个和104个，分别减少了26.2%和28.8%. 红砂岩和石灰岩裂隙面微凸体最大高度分别为0.48 mm和0.37 mm，高度和高径比均明显降低，最大高径比分别为0.23和0.28. 一方面，高度和高径比较小的微凸体发生塑性破坏被压平，另一方面，高度和高径比较大的微凸体被压碎形成更细小的微凸体. 说明试验后裂隙面粗糙度变小，整体趋于平缓. 相比来说，红砂岩裂隙面微凸体高度呈现出均匀下降趋势，微凸体通过渐进式压平实现粗糙度降低，而石灰岩裂隙面高径比0.24~0.26的微凸体出现“间断”，微凸体主要以脆性破裂为主导机制，导致局部高度分布非连续.

4 讨论

在裂隙法向受压闭合方面，学者们已经归纳总结出了对数模型，统一指数型和双曲模型等大量经验模型；在岩石流变特性方面，以往对岩石流变特性的研究主要分为两类：完整岩石的单轴/三轴压缩蠕变试验，以及含裂隙试样在渗流耦合作用下的流变试验（王振等， 2017；Zhang et al.， 2019），对法向应力作用下裂隙的时间依赖性闭合规律关注较少. 本文在前人的基础上，考虑了时间因素和形貌因素对法向应力作用下裂隙闭合规律的影响，从细观角度更真实地揭示了岩石裂隙在法向应力作用下的受压时变闭合机理.

为了对比两种岩石裂隙中不同形貌微凸体在裂隙受压时变闭合过程中的差异性，选取红砂岩下表面S⁃No.5（λ₀=0.125 9）微凸体和石灰岩下表面L⁃No.23（λ₀=0.235 0）两个微凸体，分析两者应变、蠕变变形、接触面积和法向应力的演化规律，如图19所示. 微凸体形貌受岩石类型、矿物组成和胶结方式等因素的影响显著，红砂岩中石英颗粒与黏土矿物的混合导致微凸体硬度小、高径比小（Qu et al.， 2025；Huang et al.， 2023），而石灰岩中的方解石成分导致微凸体硬度大、高径比大. L⁃No.23的高径比结构导致应变集中于局部压碎区，总体表现为应变量大但法向变形小，应变增长速度和幅度均大于S⁃No.5，表明接触区域会形成更显著的应力集中，而S⁃No.5的应变则会均匀分布于整个微凸体（图19a，19b）. L⁃No.23的接触面积增长速率大于S⁃No.5，但其法向应力下降程度更显著，而S⁃No.5的应变均匀分布于整个微凸体，颗粒滑移和错动导致整体变形缓慢增加，法向接触应力衰减程度较小（图19c，19d）.

试验前后，红砂岩和石灰岩裂隙面微凸体形貌的差异性为预测裂隙面的整体闭合速率和密封性提供了细观依据. 具体来说，红砂岩微凸体的高度均匀下降趋势反映了其矿物颗粒的均质弹-塑性变形特性，说明微凸体通过渐进式压平实现粗糙度降低；而石灰岩微凸体的高径比区间更容易出现“间断”，说明其表面微凸体主要受脆性破裂机制控制，即硬度较高、体积较大的微凸体为主要承载体，它们在应力集中下发生压碎，生成次级微凸体，导致局部高度分布非连续.

S⁃No.5和L⁃No.23在裂隙面上的位置如图20所示. 微凸体形貌是影响深地工程岩石裂隙力学行为的关键因素，其高径比特性虽不直接对应体积大小，但会显著影响接触刚度和应力分布. 当高径比大的微凸体被压碎后，应力会向更低的微凸体传递. 对比来看，红砂岩裂隙的渐进压平特性更适合长时间尺度密封；石灰岩裂隙的脆性破裂特征应采取更主动的裂隙调控和监测措施，重点关注局部密封失效风险，充分利用其应力重分布特性优化储层改造.

5 结论

本研究开展了不同形貌微凸体的受压时效性试验和不同类型岩石裂隙的受压时变闭合试验，得到了不同微凸体弹性模量随时间的衰减规律，提出了一种高效快速识别微凸体形貌的方法，并将上述创新应用于岩石裂隙受压时变闭合的边界元数值计算. 主要得到以下结论：

（1）通过开展不同高径比微凸体的受压时效性试验，发现微凸体具有明显的阶段性破坏特征，高径比较小的微凸体分担应力的能力更强；高径比大的微凸体尖端容易发生应力集中，在蠕变过程中易发生细微的脆性断裂. 微凸体尖端曲率半径越小，应力水平越高，对颗粒重排和应力重分布的驱动力越强，导致弹性模量的衰减速率越快.

（2）通过小波分析法、区域生长算法和参考面法，提出了一种裂隙面细观尺度微凸体形貌的识别方法，可简单高效地识别裂隙面上的关键微凸体，减小了数据噪声的干扰，弥补了分水岭等传统方法过分割和欠分割的缺陷，为从细观单个微凸体尺度分析裂隙面的受压时间依赖性闭合问题提供了有效手段. 基于此，通过边界元开展了真实岩石裂隙的受压时变闭合计算，实现了对单个微凸体力学参数随时间演化的定量分析.

（3）随着时间推移和法向应力增大，相互接触的微凸体数量逐渐增长，不同微凸体受细观形貌特征的影响处于不同的蠕变状态，但裂隙的时变闭合特性主要受高径比较大的几个关键微凸体控制. 相比而言，高径比较大的微凸体应变速率变化大，更容易发生加速蠕变而进入塑性阶段；而高径比较小的微凸体较稳定，分担应力的能力更强，在较长的时间内应变速率变化不大. 高径比大的微凸体一旦失稳后，应力向高径比低的微凸体传递，从而发生应力重分布，这些微凸体相互作用共同导致了裂隙的时效性闭合.

（4）红砂岩微凸体的高度均匀下降趋势反映了其矿物颗粒的均质弹-塑性变形特性，说明微凸体通过渐进式压平实现粗糙度降低；而石灰岩微凸体的高径比区间更容易出现“间断”，说明其表面微凸体主要受脆性破裂机制控制，即硬度较高、体积较大的微凸体为主要承载体，它们在应力集中下发生压碎，生成次级微凸体，导致局部高度分布非连续.

参考文献

[1] Alamos, F. J., Philo, M., Go, D. B., et al., 2022. Rough Surface Contact under Creep Conditions. Tribology International, 176: 107916. https://doi.org/10.1016/j.triboint.2022.107916

[2] Brot, C. C., Etsion, I., Kligerman, Y., 2008. A Contact Model for a Creeping Sphere and a Rigid Flat. Wear, 265(5/6): 598-605. https://doi.org/10.1016/j.wear.2007.12.003

[3] Brown, S. R., Scholz, C. H., 1985. Closure of Random Elastic Surfaces in Contact. Journal of Geophysical Research: Solid Earth, 90(B7): 5531-5545. https://doi.org/10.1029/JB090iB07p05531

[4] Chen, L., Zhao, X. G., Liu, J., et al., 2023. Progress on Rock Mechanics Research of Beishan Granite for Geological Disposal of High⁃Level Radioactive Waste in China. Rock Mechanics Bulletin, 2(3): 100046. https://doi.org/10.1016/j.rockmb.2023.100046

[5] Chung, J. C., 2010. Elastic⁃Plastic Contact Analysis of an Ellipsoid and a Rigid Flat. Tribology International, 43(1/2): 491-502. https://doi.org/10.1016/j.triboint.2009.08.005

[6] Ciavarella, M., Delfine, V., Demelio, G., 2006. A “Re⁃Vitalized” Greenwood and Williamson Model of Elastic Contact between Fractal Surfaces. Journal of the Mechanics and Physics of Solids, 54(12): 2569-2591. https://doi.org/10.1016/j.jmps.2006.05.006

[7] Goedecke, A., Jackson, R. L., Mock, R., 2010. Asperity Creep under Constant Force Boundary Conditions. Wear, 268(11/12): 1285-1294. https://doi.org/10.1016/j.wear.2010.01.025

[8] Greenwood, J. A., 2006. A Simplified Elliptic Model of Rough Surface Contact. Wear, 261(2): 191-200. https://doi.org/10.1016/j.wear.2005.09.031

[9] Greenwood, J. A., Williamson, J. B. P., 1966. Contact of Nominally Flat Surfaces. Proceedings of the Royal Society of London Series A, Mathematical and Physical Sciences, 295(1442): 300-319

[10] Huang, K., Yu, F., Zhang, W., et al., 2023. Experimental and Numerical Simulation Study on the Influence of Gaseous Water on the Mechanical Properties of Red⁃Layer Mudstone in Central Sichuan. Rock Mechanics and Rock Engineering, 56(4): 3159-3178. https://doi.org/10.1007/s00603⁃023⁃03228⁃z

[11] Kang, F. C., Li, Y. C., Tang, C. A., et al., 2022. Competition between Cooling Contraction and Fluid Overpressure on Aperture Evolution in a Geothermal System. Renewable Energy, 186: 704-716. https://doi.org/10.1016/j.renene.2022.01.033

[12] Kang, H., Einstein, H., Brown, S., et al., 2020. Numerical Simulation for Rock Fracture Viscoelastic Creep under Dry Conditions. Geofluids, 2020(1): 8879890. https://doi.org/10.1155/2020/8879890

[13] Kling, T., Vogler, D., Pastewka, L., et al., 2018. Numerical Simulations and Validation of Contact Mechanics in a Granodiorite Fracture. Rock Mechanics and Rock Engineering, 51(9): 2805-2824. https://doi.org/10.1007/s00603⁃018⁃1498⁃x

[14] Kumamoto, K. M., Thom, C. A., Wallis, D., et al., 2017. Size Effects Resolve Discrepancies in 40 Years of Work on Low⁃Temperature Plasticity in Olivine. Science Advances, 3(9): e1701338. https://doi.org/10.1126/sciadv.1701338

[15] Li, B., Cui, X. F., Mo, Y. Y., et al., 2021. Deformation Behavior of Dislocated Sandstone Fractures Subject to Normal Stresses. Rock and Soil Mechanics, 42(7): 1850-1860 (in Chinese with English abstract).

[16] Li, B., Mo, Y. Y., Zou, L. C., et al., 2022. An Extended Hyperbolic Closure Model for Unmated Granite Fractures Subject to Normal Loading. Rock Mechanics and Rock Engineering, 55(7): 4139-4158. https://doi.org/10.1007/s00603⁃022⁃02862⁃3

[17] Li, B., Zhao, Z. H., Jiang, Y. J., et al., 2015. Contact Mechanism of a Rock Fracture Subjected to Normal Loading and Its Impact on Fast Closure Behavior during Initial Stage of Fluid Flow Experiment. International Journal for Numerical and Analytical Methods in Geomechanics, 39(13): 1431-1449. https://doi.org/10.1002/nag.2365

[18] Ma, H. C., Cao, Y., Qian, J. Z., et al., 2023. Theoretical Study of the Mesoscopic Mechanism of Rock Fractures during Normal Deformation. Rock Mechanics and Rock Engineering, 56(8): 5719-5733. https://doi.org/10.1007/s00603⁃023⁃03372⁃6

[19] Malamut, S., Kligerman, Y., Etsion, I., 2009. The Effect of Dwell Time on the Static Friction in Creeping Elastic⁃Plastic Polymer Spherical Contact. Tribology Letters, 35(3): 159-170. https://doi.org/10.1007/s11249⁃009⁃9445⁃3

[20] Matsuki, K., Wang, E. Q., Sakaguchi, K., et al., 2001. Time⁃Dependent Closure of a Fracture with Rough Surfaces under Constant Normal Stress. International Journal of Rock Mechanics and Mining Sciences, 38(5): 607-619. https://doi.org/10.1016/S1365⁃1609(01)00022⁃3

[21] Ovcharenko, A., Halperin, G., Etsion, I., 2009. Experimental Study of a Creeping Polymer Sphere in Contact with a Rigid Flat. Journal of Tribology, 131: 011404. https://doi.org/10.1115/1.3002330

[22] Pyrak⁃Nolte, L. J., Nolte, D. D., 2016. Approaching a Universal Scaling Relationship between Fracture Stiffness and Fluid Flow. Nature Communications, 7: 10663. https://doi.org/10.1038/ncomms10663

[23] Qu, J. K., Xue, Y. G., Kong, F. M., et al., 2025. Multi⁃Scale Analysis of the Influence of Red⁃Bed Lithological Interface on Tunnel Deformation and Instability. Engineering Geology, 357: 108314. https://doi.org/10.1016/j.enggeo.2025.108314

[24] Rohmer, J., Pluymakers, A., Renard, F., 2016. Mechano⁃Chemical Interactions in Sedimentary Rocks in the Context of CO₂ Storage: Weak Acid, Weak Effects? Earth⁃Science Reviews, 157: 86-110. https://doi.org/10.1016/j.earscirev.2016.03.009

[25] Rostami, A., Goedecke, A., Mock, R., et al., 2014. Three⁃Dimensional Modeling of Elasto⁃Plastic Sinusoidal Contact under Time Dependent Deformation Due to Stress Relaxation. Tribology International, 73: 25-35. https://doi.org/10.1016/j.triboint.2013.12.020

[26] Serati, M., Alehossein, H., Williams, D. J., 2015. Estimating the Tensile Strength of Super Hard Brittle Materials Using Truncated Spheroidal Specimens. Journal of the Mechanics and Physics of Solids, 78: 123-140. https://doi.org/10.1016/j.jmps.2015.02.011

[27] Song Z. Y., Li B., Yan B. M., et al., 2024. Time⁃Dependent Closure Calculation Method of Rock Fractures Considering the Attenuation of Relaxation Modulus. In: China Rock Mechanics and Engineering Society, International Association for Geohazards and Risk Reduction. Proceedings of the 21st Chinese Rock Mechanics and Engineering Academic Conference, CHINA ROCK 2024. Department of Underground Architecture and Engineering, Tongji University; Department of Environmental Science and Engineering, Royal Institute of Technology, Sweden, 2024: 42-48 (in Chinese with English abstract).

[28] Song, Z. Y., Li, B., Chen, X. W., et al., 2025. Time⁃Dependent Deformation of Fracture Asperities with Different Height⁃to⁃Radius Ratios Subject to Normal Loading. Rock Mechanics and Rock Engineering, 2025: 1-18. https://doi.org/10.1007/s00603⁃025⁃04858⁃1

[29] Sun, J., 2007. Rock Rheological Mechanics and Its Advance in Engineering Applications. Chinese Journal of Rock Mechanics and Engineering, 26(6): 1081-1106 (in Chinese with English abstract).

[30] Tang, Z. C., Zhang, Q. Z., 2021. Elliptical Hertz⁃Based General Closure Model for Rock Joints. Rock Mechanics and Rock Engineering, 54(1): 477-486. https://doi.org/10.1007/s00603⁃020⁃02275⁃0

[31] Tian, X. F., Bhushan, B., 1996. A Numerical Three⁃Dimensional Model for the Contact of Rough Surfacesby Variational Principle. Journal of Tribology, 118(1): 33-42. https://doi.org/10.1115/1.2837089

[32] Viswanathan, H. S., Ajo⁃Franklin, J., Birkholzer, J. T., et al., 2022. From Fluid Flow to Coupled Processes in Fractured Rock: Recent Advances and New Frontiers. Reviews of Geophysics, 60(1): e2021RG000744. https://doi.org/10.1029/2021RG000744

[33] Wang, Z., Shen, M. R., Tian, G. H., et al., 2017. Characteristics of Aging Strength of Structural Planes with Different Roughness. Chinese Journal of Rock Mechanics and Engineering, 36(S1): 3287-3296(in Chinese with English abstract).

[34] Wen, Y. Q., Tang, J. Y., Zhou, W., et al., 2022. New Analytical Model of Elastic⁃Plastic Contact for Three⁃Dimensional Rough Surfaces Considering Interaction of Asperities. Friction, 10(2): 217-231. https://doi.org/10.1007/s40544⁃020⁃0419⁃7

[35] Xue, Y. C., Xu, T., Heap, M. J., et al., 2023. Time⁃Dependent Cracking and Brittle Creep in Macrofractured Sandstone. International Journal of Rock Mechanics and Mining Sciences, 162: 105305. https://doi.org/10.1016/j.ijrmms.2022.105305

[36] Zhang, L. L., Wang, X. J., 2020. Viscoelastic⁃Plastic Damage Creep Model for Rock. Chinese Journal of Geotechnical Engineering, 42(6): 1085-1092 (in Chinese with English abstract).

[37] Zhang, Q. Z., Wu, C. Z., Fei, X. C., et al., 2019. Time⁃Dependent Behavior of Rock Joints Considering Asperity Degradation. Journal of Structural Geology, 121: 1-9. https://doi.org/10.1016/j.jsg.2019.01.004

[38] Zhang, W. G., Lin, S. C., Wang, L. Q., et al., 2024. A Novel Creep Contact Model for Rock and Its Implement in Discrete Element Simulation. Computers and Geotechnics, 167: 106054. https://doi.org/10.1016/j.compgeo.2 023. 106054

[39] Zou, L. C., Li, B., Mo, Y. Y., et al., 2020. A High⁃Resolution Contact Analysis of Rough⁃Walled Crystalline Rock Fractures Subject to Normal Stress. Rock Mechanics and Rock Engineering, 53(5): 2141-2155. https://doi.org/10.1007/s00603⁃019⁃02034⁃w

基金资助

国家自然科学基金资助项目(42077252)

国家自然科学基金资助项目(42377162)