基于颗粒接触理论的岩石压密阶段本构模型-有色金属在线

段非线性变形与岩石内部颗粒之间的相互接触有关，引入基于凹凸体正态分布的G-W接触模型并建立接触体元件模型；将该模型与弹簧元件串联，建立压密阶段的本构模型。利用元件模型对实验结果进行模拟，同时对模型中参数进行分析。研究结果表明：岩石颗粒间的张开度及颗粒表面的粗糙度均对岩石压密阶段本构关系有一定影响。颗粒间的张开度越大，压密过程越长；在相同应变条件下，粗糙度越大，应力增长越缓慢。本文模型在模拟4种岩石压密阶段的本构关系过程中，拟合优度均大于0.95，说明模型能够准确描述岩石的应力-应变关系。同时，模型还能够反映围压对本构关系的影响，证明利用接触模型来描述压密阶段本构关系是合理的。

关键词:岩石力学；压密阶段；本构模型；接触理论；正态分布

中图分类号:TU45 文献标志码:A 文章编号:1672-7207（2019）08-1941-08

Constitutive model of rock compaction stage based on contact theory

MA　Qiufeng, QIN　Yueping, ZHOU　Tianbai, YANG　Xiaobin

(School of Emergency Management and Safety Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China)

Abstract: Considering that the nonlinear deformation in the compaction stage of rock is related to the contact between the particles in the rock,a G-W contact model based on the normal distribution of concave convex body was introduced and then the contact element model was established. The constitutive model of the compaction stage was established by connecting the model with the spring element in series. The element model was used to simulate the experimental results, and the parameters of the model were analyzed. The results show that both the aperture of rock particles and the roughness of the particle surface have certain effects on the constitutive relationship of rock compaction stage. The larger the aperture degree between particles is, the longer time the compaction process takes. Under the same strain condition, the greater the roughness is, the slower the stress growth is. In the process of simulating the constitutive relationship of four rock compaction stages, the goodness of fit of the model is greater than 0.95. The proposed model can accurately describe the stress-strain relationship of rock. At the same time, the model reflects the influence of confining pressure on the constitutive relationship, which verifies that it is reasonable to use contact model to describe the constitutive relationship in the compaction stage .

Key words: rock mechanics; compaction stage; constitutive model; contact theory; normal distribution

砂岩是最重要的油气储集层，世界上大部分的油气资源都储存在砂岩中。此外，砂岩中常常有煤和各种金属矿床。因此，对于砂岩力学性质的研究具有重要的意义。由于砂岩内部结构的特殊性，传统Hooke模型不能准确地描述其应力-应变关系，其原因如下：砂岩存在压密阶段，岩石表现出非线性特征；砂岩具有塑性硬化和软化现象，同样表现出非线性特征。人们对砂岩的力学性质进行了大量研究，发现砂岩的强度受到孔隙率的影响，孔隙率越高，砂岩的强度越低^[1^-3]。同时，砂岩屈服面、硬化函数、塑性变形同样受到孔隙率的影响^[4^-8]。BAUD等^[6]分析了不同孔隙率下岩石的单轴抗压强度，提出了影响岩石强度的含孔隙翼型裂纹模型。KULHAWY^[7]认为随着孔隙的压密，弹性模量迅速增大。AL-HARTHI等^[8]分析了孔隙率对岩石弹性模量的影响，并建立了考虑孔隙影响的弹性模量模型。ZHAO等^[9^-10]提出了各向异性应力条件下多孔岩石的三维模型，在该模型中，一部分孔隙内的变形用基于自然应变的胡克定律描述，另一部分孔隙内的变形用基于工程应变的胡克定律描述。李连崇等^[11]完善了上述模型，并利用该模型实现了FLAC^3D开发。王青元等^[12]将岩石峰前变形阶段分为2个部分，分别建立了本构模型。曹文贵等^[13^-14]在统计损伤模型的基础上，采用宏观与微观相结合的方法，将空隙岩石分为岩石骨架和空隙2个部分，建立空隙岩石变形分析模型；该模型不仅反映了空隙岩石的应变软化特征，而且还反映了空隙岩石在初始压密阶段的变形非线性特点。近年来，人们从颗粒接触入手，对砂岩的力学性质进行分析。LESZCZYNSKI^[15]利用颗粒的接触黏聚力，模拟颗粒材料的动载过程。NASSAUER等^[16]提出了多面体颗粒接触计算方法。FENG^[17]提出了非圆球形接触模型的计算方法，为离散元计算提供了接触模型。袁康等^[18]对岩石压缩荷载作用下内部颗粒组分的力学响应进行了研究，从细观层面揭示了岩石的破裂机制。目前，有关砂岩的压密阶段以及砂岩颗粒接触的研究较多，而通过接触理论来描述岩石压密阶段的本构关系的研究较少。本文作者以颗粒接触为研究基础，认为岩石压密阶段中碎屑颗粒之间的孔隙也被压密，并且岩石压密阶段非线性特征与颗粒间的接触有关。通过引入G-W接触模型，建立接触元件模型用于表征岩石中碎屑颗粒的接触，同时采用弹簧元件表征颗粒本身的线弹性特征，将二者串联，用于描述岩石压密过程中的应力-应变关系。利用建立的串联模型计算岩石压密阶段的本构关系，并与实验结果进行对比以验证模型的合理性。

1 接触模型的建立

FX_GRP_ID800008DA

图1　西原正夫模型示意图

Fig. 1　Diagram of Nishihara model

在沉积学中，根据砂岩内部结构，砂岩主要由碎屑颗粒、杂基、胶结物和孔隙这4个部分组成。其中杂基和胶结物又合称为填隙物。根据砂岩内部碎屑颗粒数量的不同，砂岩的承载形式主要分为2类：颗粒承载和杂基承载。本文主要研究碎屑颗粒数较多的颗粒承载。CHEN等^[19]认为，在成岩过程中，岩石存在压溶作用，岩石中碎屑颗粒之间呈现凹凸状的接触。由于颗粒间凹凸体的存在，岩石在初始加载过程中呈现出非线性特征。

在描述压密阶段的本构关系时，传统本构模型^[7^-10]往往采用孔隙率作为内变量来描述岩石的应力-应变关系。尽管模型在一定程度上能够描述压密过程，然而此类模型只能反映出压密阶段一部分力学行为，并没有从岩石压密过程的机理上反映出岩石的非线性变形规律。例如，模型没有考虑压密过程中渗透方向等问题。

本文作者通过分析沉积学中关于砂岩压密阶段的承载规律，提出利用颗粒间不平整接触面作为本构模型的建模基础。首先建立能够反映颗粒间不平整接触的元件模型，然后利用该元件模型建立岩石压密阶段本构模型。从模型结构来看，本文提出的模型更加简洁，同时物理意义更加明确。

传统的元件模型在岩石力学行为研究中得到了广泛应用。其中，经典西原正夫模型(Nishihara model)被广泛应用于描述岩石蠕变特性。西原正夫模型示意图如图1所示。该模型主要包括弹簧元件、黏壶元件以及塑性元件。人们在元件模型的基础上，提出修正原正夫模型或者新的模型^[20^-23]，但用于描述岩石压密阶段的模型还未见报道。本文作者基于上述颗粒接触理论，在G-W接触模型的基础上，提出岩石压密阶段的元件模型，用于描述岩石的本构关系，如图2所示。

FX_GRP_ID800028AA

图2 压密阶段本构关系的元件模型

Fig. 2　Element model in compaction stage

由图2可知：该元件存在一个倾斜面，在该倾斜面上存在满足赫兹接触理论的凸起，且凸起高度满足正态分布。2个接触面之间的平均距离为(即接触面之间的张开度，简称张开度)。接触面与垂直方向的夹角为(简称接触面倾角)。接触面法向方向的力-位移关系可以采用G-W模型进行计算。同时假设该元件宽度和厚度均为1 m。由于颗粒之间接触面的角度是随机的，并没有固定的方向，因此，接触面倾角用可于描述大量接触面接触方向的整体效果。

接触元件模型用于描述压密阶段颗粒间的接触效应。此外，假设颗粒本身满足线弹性特征。为了表征颗粒本身的变形规律，在接触元件上串联1个弹簧元件(见图2(b))。

由于压密阶段往往应力较小，本文假设碎屑颗粒不发生塑性变形，即颗粒一直处于弹性状态；由于碎屑颗粒之间存在一定的黏聚力，因此，假设接触面不发生相互滑动，仅考虑接触面的法向作用力。

2 G-W模型简介

经典G-W模型用于描述颗粒间接触面的相互作用。G-W模型中，粗糙表面被描述为一系列凸起的组合，凸起之间通过赫兹定律进行计算。

G-W模型假设如下：1) 粗糙面剖面凸起的高度服从正态分布；2) 凸起顶端具有相同的曲率；3) 凸起变形相互独立，互不影响；4) 忽略凸起下部的体积变化。

2.1　G-W模型中的粗糙度与概率密度

G-W模型中的粗糙表面示意图如图3所示。图中，为在长度l_c范围内横坐标为x处的高度。剖面粗糙度均方根(简称粗糙度)表示粗糙表面顶点高度与粗糙表面平均高度差值的标准差，其表达式为

(1)

G-W模型中粗糙面的概率密度示意图如图4所示。图中，黑色区域表示高度大于z的凸起的概率。概率密度函数用于描述剖面某点处高度的概率分布。若以粗糙面平均高度为基准，在高度为(z,z+dz)范围内出现凸起顶点的概率为p(z，z+dz)，则概率密度函数为

(2)

在G-W模型中，认为粗糙面凸起高度满足正态分布，其概率密度为

(3)

FX_GRP_ID80002C82

图3 粗糙表面示意图

Fig. 3　Diagram of surface roughness

FX_GRP_ID80002E05

图4　粗糙面概率密度的求解示意图

Fig. 4　Schematic map for probability density solution of rough surface

2.2　粗糙面之间的接触分析

2个不同粗糙度的平面相互接触问题可以等效为刚性平面与可变形的粗糙平面相互接触的问题。等效后平面的粗糙度_e可用等效前2个平面的粗糙度表示：

(4)

式中：和分别为2个接触面的粗糙度。

将研究对象取为高度为的凸起(为高于粗糙面平均高度)。当刚性平面与粗糙面平均高度平面的距离为d时：1) 若d>，则凸起与刚性平面不接触；若d<，则凸起与刚性平面相互接触，两者相互重叠的量为。根据赫兹定律，凸起与刚性平面的接触力f (z_s)为

(5)

式中：为凸起的弹性模量，为凸起的曲率。与刚性平面接触的所有凸起的概率p_rob(z_s>d)为

(6)

因此，单位面积上产生的接触力P(d)为

(7)

式中：N为接触凸起的数量。为了便于计算，假设当接触面位移为0 mm时，接触产生的压力近似为0 N。转化坐标后的接触面示意图如图5所示。以接触面零位移为坐标原点，凸起的平均高度距离起始点的距离视为张开度d₀。根据式(7)转化坐标后可得位移为时接触面的法向力P()为

(8)

由于颗粒间存在黏结物，忽略剪应力对接触面的影响。当一个斜向力作用在接触面时，假设只考虑法向分量对截面的影响。由于本文重点研究压密阶段岩石的本构关系，故暂未考虑当剪应力较大时接触面滑移以及岩石破坏的问题。

FX_GRP_ID80003323

图5　转化坐标后的接触面示意图

Fig. 5　Schematic diagram of contact surface after coordinate transformation

3 压密阶段本构关系

3.1　组合元件的本构关系

当元件两端产生压缩量后，令弹簧体的竖直方向位移为，接触体的竖直方向位移为，则有

(9)

弹簧中的力P₁为

(10)

式中：为弹性系数。根据几何关系，接触体的法向位移为

(11)

接触体的接触力P()计算公式见式(8)。

对接触力进行分解，得到接触元件在竖直方向的力P₂为

(12)

根据串联结构，令，联立式(10)和(12)，计算得到压缩量后组合元件中的力。已知模型尺寸后，便能够计算得到应力-应变关系。

通过上述关系，可计算得到在任意压缩量下组合单元中的力。为了求解组合单元中的应力-应变关系，需要定义每个元件的宽度与高度。已知岩石颗粒的弹性模量为E_s，则弹簧的弹性系数为

(13)

式中：为元件所代表的岩石水平方向的截面积，此处为单位面积；为元件中弹性部分的高度。

对于任意元件，

(14)

式中：h₀为元件高度，此处取为1 m。

3.2　二分法在计算本构关系中的应用

式(8)是反映接触元件力-位移关系的重要表达式，但无法对其进行积分求解。本文采用数值计算的方法，计算得到近似解。由组合元件的力-位移关系可得：

(15)

假设总压缩量已知，求解组合元件的受力。采用二分法进行数值计算，如图6所示。

FX_GRP_ID80003962

图6　二分法在计算过程中的应用

Fig. 6　Application of dichotomy in calculation process

首先将l分为两等分(见图6(a))，令，假设此时弹簧元件中的力大于接触元件中的力，这说明弹簧元件的压缩量过大，需进行下一步修正。将图中AB段进行等分(见图6(b))，此时令，。若此时>，则等分左侧AD段，令，；若<，则等分DB段，令，；若此时=，则说明原假设，就是此时压缩量分配的计算结果。

依照上述方法，不断进行划分，直到达到计算的精度要求，即(其中为精度判断参数)。

3.3　侧向压力对压密过程的影响

本文第3.1节中的受力分析过程仅考虑在单轴压缩条件下的本构关系，下面将分析侧向压力条件下的本构关系。

首先对元件施加侧向压应力，然后再施加竖直方向的位移。通过侧向压力计算得到接触面上的侧向力，然后计算得到侧向压力在接触面产生的法向初始位移。在接触面法向存在预压力：

(16)

此后，施加竖直方向的位移。假设接触元件竖直位移为，此时接触面的法向位移为

(17)

通过式(8)计算得到法向接触力P()，此时因竖直位移产生的力P₂为

(18)

联立，当已知时即可计算得到单个组合元件的力。在已知模型尺寸后，便能够计算得到应力-应变关系。

4 模型的验证

下面验证模拟孔隙压密阶段岩石本构关系的元件模型的合理性。本文模型的参数主要包括颗粒的弹性模量、单位面积凸起的数量N、凸起的曲率、接触面的张开度、接触面粗糙度和接触模型中的倾角。

本文选取4种不同孔隙率砂岩的实验结果^[24^-25]进行对比。由于文献[24-25]中没有给出具体的模型参数，本文参考文献[26-27]中微观砂岩颗粒的观察结果，推测出4种岩石的参数，如表1所示。由于孔隙率是确定颗粒间的张开度的重要参数，因此也列于表1中。由于对称性，颗粒之间接触面倾角的范围为[0°，90°]。对于颗粒接触面分布较为均匀的岩石而言，接触模型中的倾角可以通过求解倾角的数学期望得到，。

表1　模型参数

Table 1　Model parameters

利用MATLAB数值计算软件对模型进行编程后，模拟得到模型压密阶段的应力-应变关系。不同砂岩本构关系的实验结果与数值模拟结果对比如图7所示。图7中，R²为模型对实验结果的拟合优度。

FX_GRP_ID8000461E

图7　不同砂岩本构关系的实验结果与数值模拟结果对比

Fig. 7　Comparison of experimental results and numerical simulation results for different sandstone constitutive relations

从图7可以看出：在岩石压密阶段本构关系表现出明显的非线性特征。在初始加载过程中，岩石应力增长较为缓慢，随着轴向位移的增大，应力增长速率逐渐趋于稳定。对比图7(a)和图7(c)可知：绿砂岩的压密阶段对应的应变范围为0~0.1%，Boise砂岩的压密阶段对应的应变范围为0~0.23%，这说明对于不同岩石材料，压密阶段对应的应变范围不同。

本文模型对绿砂岩的拟合优度为0.984，对Berea砂岩的拟合优度为0.991，对Boise砂岩的拟合优度为0.972，对Fontainebleau砂岩的拟合优度为0.978。本文模型对实验结果模拟过程中，拟合优度均大于0.95，能够较准确地描述岩石在压密阶段的非线性变形规律，证明了模型的合理性。

利用Fontainebleau砂岩的模型参数，对不同侧向压力条件下进行压缩模拟，模拟过程中分别采用0，5，10，20和30 MPa的侧向压力。

模拟的加载过程如下：首先施加侧向压力，垂直方向位移清零后施加垂直压缩位移。图8所示为位移过程中竖向的应力和应变的关系。

从图8可以看出：在不同侧向压力的情况下，随着侧向压力的增大，初始压密阶段增长的速率越快。当侧压力接近30 MPa时，应力-应变曲线趋近于线性。由此可见，数值模拟结果反映了岩石压密阶段的变形特点，与实验规律基本一致，这说明模型具有合理性。

5 模型参数分析

FX_GRP_ID80005C46

图8　不同侧向压力条件下压密阶段的应力-应变关系

Fig. 8　Stress-strain relationship in compaction stage under different lateral pressure conditions

5.1　接触面粗糙度对压密阶段的影响

首先对接触面粗糙度进行分析。在本模型中，粗糙度表示接触表面凸起高度距离平均高度的离散程度。不同粗糙度的表面如图9所示。由图9可见：表面1的离散度更高，因此表面1的粗糙度更大。

取不同粗糙度在无侧向压力的条件下进行数值计算，计算结果如图10所示。

FX_GRP_ID8000498A

图9　不同粗糙度的表面

Fig. 9　Surfaces with different roughness

FX_GRP_ID80004ABD

图10　不同粗糙度条件下的应力-应变关系

Fig. 10　Stress-strain relationship under different roughness conditions

从图10可以看出：当粗糙度较小时，随着应变的增大，应力增长较快；当粗糙度较大时，应力增长的速率较慢。这与接触面中凸起的接触密切相关，当较小时，说明发生接触的凸起数量迅速增加，导致应力快速增长。

5.2　接触面的张开度对压密阶段的影响

张开度与岩石的孔隙率紧密相关。通过表1与图7的计算结果可以看出：当岩石孔隙率较大时，需采用较大的张开度才能得到较好的数值计算结果；当岩石孔隙率较小时，需采用较小的张开度才能得到较好的数值计算结果。

取不同张开度在无侧向压力的条件下进行数值计算，计算结果如图11所示。从图11可以看出：当张开度较小时，在经历较小的应变后，应力迅速增大；当张开度较大时，在经历较大的应变后，应力才开始逐渐增长。张开度综合反映了压密过程的长短，当张开度较大时，需要经历较长的压密过程，而小的张开度只需要经历较少的压密过程。

6 结论

1) 本文建立的接触元件模型在模拟4种岩石弹性阶段本构关系过程中，拟合优度均大于0.95，模型的模拟结果与实验结果一致，说明模型能够准确地描述岩石弹性阶段的本构关系。

2) 本文提出的接触元件模型能够反映围压对岩石轴向本构关系的影响规律，证实了利用接触模型用于描述压密阶段本构关系的合理性。

3) 岩石颗粒间的张开度对岩石压密阶段具有一定影响。颗粒间的张开度越大，压密过程经历的应变越多；相反，张开度越小，岩石颗粒经历较少的应变后能够进入线弹性阶段。

4) 岩石接触表面的粗糙度对岩石压密阶段具有一定影响。在相同应变条件下，粗糙度越大，应力增大越缓慢。

FX_GRP_ID800048EC

图11　不同张开度条件下的应力-应变关系

Fig. 11　Stress-strain relationship under different opening conditions

参考文献：

[1] JAMSHIDI A, ZAMANIAN H, ZAREI SAHAMIEH R. The effect of density and porosity on the correlation between uniaxial compressive strength and P-wave velocity[J]. Rock Mechanics and Rock Engineering, 2018, 51(4): 1279-1286.

[2] VAJDOVA V, ZHU W, NATALIE CHEN T M, et al. Micromechanics of brittle faulting and cataclastic flow in Tavel limestone[J]. Journal of Structural Geology, 2010, 32(8): 1158-1169.

[3] MENNDEZ B, ZHU Wenlu, WONG T F. Micromechanics of brittle faulting and cataclastic flow in Berea sandstone[J]. Journal of Structural Geology, 1996, 18(1): 1-16.

[4] SUN Lizhi, HUANG Zhuping. Dynamic void growth in rate-sensitive plastic solids[J]. International Journal of Plasticity, 1992, 8(8): 903-924.

[5] BESSON J, ABOUAF M. Rheology of porous alumina and simulation of hot isostatic pressing[J]. Journal of the American Ceramic Society, 1992, 75(8): 2165-2172.

[6] BAUD P, WONG T F, ZHU W. Effects of porosity and crack density on the compressive strength of rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2014, 67: 202-211.

[7] KULHAWY F H. Stress deformation properties of rock and rock discontinuities[J]. Engineering Geology, 1975, 9(4): 327-350.

[8] AL-HARTHI A A, AL-AMRI R M, SHEHATA W M. The porosity and engineering properties of vesicular basalt in Saudi Arabia[J]. Engineering Geology, 1999, 54(3/4): 313-320

[9] ZHAO Yu, LIU Huihai. An elastic stress-strain relationship for porous rock under anisotropic stress conditions[J]. Rock Mechanics and Rock Engineering, 2012, 45(3): 389-399.

[10] LIU Huihai, RUTQVIST J, BERRYMAN J G. On the relationship between stress and elastic strain for porous and fractured rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2009, 46(2): 289-296.

[11] 李连崇, LIU Huihai, 赵瑜. 基于双应变胡克模型的岩石非线性弹性行为分析[J]. 岩石力学与工程学报, 2012, 31(10): 2119-2126.

LI Lianchong, LIU Huihai, ZHAO Yu. Investigation on nonlinear elastic behaviour of rocks based on a two-part Hooke's model[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(10): 2119-2126.

[12] 王青元, 朱万成, 徐涛, 等. 考虑孔隙压密的岩石非线性变形行为计算分析[J]. 中国科学:技术科学, 2018, 48(5): 565-574.

WANG Qingyuan, ZHU Wancheng, XU Tao, et al. Numerical analysis on nonlinear deformation behavior of rock considering compaction of pores[J]. Scientia Sinica(Technologica), 2018, 48(5): 565-574.

[13] 曹文贵, 杨尚, 张超. 考虑弹性模量变化的岩石统计损伤本构模型[J]. 水文地质工程地质, 2017, 44(3): 42-48.

CAO Wengui, YANG Shang, ZHANG Chao. A statistical damage constitutive model of rocks considering the variation of the elastic modulus[J]. Hydrogeology & Engineering Geology, 2017, 44(3): 42-48.

[14] 曹文贵, 张超, 贺敏, 等. 考虑空隙压密阶段特征的岩石应变软化统计损伤模拟方法[J]. 岩土工程学报, 2016, 38(10): 1754-1761.

CAO Wengui, ZHANG Chao, HE Min, et al. Statistical damage simulation method of strain softening deformation process for rocks considering characteristics of void compaction stage[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(10): 1754-1761.

[15] LESZCZYNSKI J S. A discrete model of a two-particle contact applied to cohesive granular materials[J]. Granular Matter, 2003, 5(2): 91-98.

[16] NASSAUER B, KUNA M. Contact forces of polyhedral particles in discrete element method[J]. Granular Matter, 2013, 15(3): 349-355.

[17] FENG Y T. A general contact theory for non-spherical particles[C]// Proceedings of the 7th International Conference on Discrete Element Methods. Singapore: Springer, 2016, 181: 29-35.

[18] 袁康, 蒋宇静, 李亿民, 等. 基于颗粒离散元法岩石压缩过程破裂机制宏细观研究[J]. 中南大学学报(自然科学版), 2016, 47(3): 913-922.

YUAN Kang, JIANG Yujing, LI Yimin, et al. Macro-micro mechanical research on failure mechanism of rock subjected to compression loading based on DEM[J]. Journal of Central South University(Science and Technology), 2016, 47(3): 913-922.

[19] CHEN Zhihe, FANG Hongwei. Analysis of the complex morphology of sediment particle surface based on electron microscope images[J]. Science China Technological Sciences, 2013, 56(2): 280-285.

[20] PENG Yu, ZHAO Jinzhou, LI Yongming. A wellbore creep model based on the fractional viscoelastic constitutive equation[J]. Petroleum Exploration and Development, 2017, 44(6): 1038-1044.

[21] HOU Feng, LAI Yuanming, LIU Enlong, et al. A creep constitutive model for frozen soils with different contents of coarse grains[J]. Cold Regions Science and Technology, 2018, 145: 119-126.

[22] HAN Bing, XIE Huibing, ZHU Li, et al. Nonlinear model for early age creep of concrete under compression strains[J]. Construction and Building Materials, 2017, 147: 203-211.

[23] YANG Liu, LI Zhida. Nonlinear variation parameters creep model of rock and parametric inversion[J]. Geotechnical and Geological Engineering, 2018, 36(5): 2985-2993.

[24] WANG Q Y, ZHU W C, XU T, et al. Numerical simulation of rock creep behavior with a damage-based constitutive law[J]. International Journal of Geomechanics, 2017, 17(1): 04016044.

[25] BAUD P, WONG T F, ZHU W. Effects of porosity and crack density on the compressive strength of rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2014, 67: 202-211.

[26] 刘惟庆, 吴伟, 秦成岗, 等. 基于扫描电镜的砂岩储层分析: 以白云凹陷珠江组三角洲前缘砂体为例[J]. 电子显微学报, 2014, 33(2): 117-122.

LIU Weiqing, WU Wei, QIN Chenggang, et al. Analysis on the sandstone reservoir by SEM: A case study of the delta front sand body of zhujiang formation, Baiyun sag[J]. Journal of Chinese Electron Microscopy Society, 2014, 33(2): 117-122.

[27] 许江, 田傲雪, 程立朝, 等. 砂岩双面剪切细观损伤特性试验研究[J]. 岩石力学与工程学报, 2015, 34(S2): 3642-3651.

XU Jiang, TIAN Aoxue, CHENG Lichao, et al. Experimental research on meso-damage characteristic of double sheared sandstone[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(S2): 3642-3651.

(编辑伍锦花)

收稿日期： 2018 -09 -13; 修回日期： 2018 -12 -13

基金项目:国家自然科学基金资助项目(50904071, 51274207)

Foundation:(Projects (50904071,51274207) supported by the National Natural Science Foundation of China)

通信作者：马秋峰，博士研究生，从事岩石力学与煤岩动力研究；E-mail：maqiufeng666@sina.com

摘要:鉴于岩石压密阶段非线性变形与岩石内部颗粒之间的相互接触有关，引入基于凹凸体正态分布的G-W接触模型并建立接触体元件模型；将该模型与弹簧元件串联，建立压密阶段的本构模型。利用元件模型对实验结果进行模拟，同时对模型中参数进行分析。研究结果表明：岩石颗粒间的张开度及颗粒表面的粗糙度均对岩石压密阶段本构关系有一定影响。颗粒间的张开度越大，压密过程越长；在相同应变条件下，粗糙度越大，应力增长越缓慢。本文模型在模拟4种岩石压密阶段的本构关系过程中，拟合优度均大于0.95，说明模型能够准确描述岩石的应力-应变关系。同时，模型还能够反映围压对本构关系的影响，证明利用接触模型来描述压密阶段本构关系是合理的。