International Open Access Journal Platform

人工裂缝导流能力是判断压裂施工后增产效果优劣的重要指标。随着裂缝导流能力损伤研究的发展，由岩石蠕变引起的导流能力损伤渐渐受到了较多关注。国内外学者对岩石的黏弹性力学行为也进行了深入的研究。当前裂缝导流能力衰减研究多聚焦支撑剂变形与嵌入的影响，未考虑岩石蠕变这种黏弹性力学行为的影响［1-2］。而大量岩石力学实验证实，在含水环境下岩石会普遍表现出较强的蠕变行为[3]。实验研究表明，蠕变行为对裂缝长期导流能力影响显著，实验测试20小时后裂缝导流能力下降20%[4]。同时，水-岩作用会显著促进岩石的蠕变特征，蠕变速率可提升46%～162%，杨氏模量与硬度明显降低［5-6］。

结合数值模拟与室内实验，众多学者建立了关于岩石蠕变的数学模型并在工程问题中开展应用。Cao等
人[7]采用实验测试和数值模拟的手段，研究了砂岩的蠕变与井眼失稳防治。Firme等人[8]通过软件，模拟了基于单参数、双参数和多参数本构模型的井眼和岩心蠕变。之后，Orlic等人[9]进行了盐膏盖层单井废弃过程中井眼蠕变闭合的数值模拟。2015年，Almasoodi等人[10]研究了油田工作液对Eagle Ford区块岩石粘弹性蠕变的影响。2016年，Rassouli和Zoback[11]通过实验对比了非常规储层中长期和短期蠕变的差异。同年，王萍[12]建立了能够考虑储层岩石水化和蠕变损伤的粘弹性本构方程。

在裂缝导流能力表征方面，现有数学模型多考虑支撑剂变形与嵌入，忽略了岩石蠕变的影响[13]。由于蠕变会降低导流能力并加剧支撑剂嵌入，将此类模型用于导流能力描述时会产生较大误差。为提升模型精度则需将蠕变模型与裂缝导流能力模型耦合。目前基于实验与模拟已建立多种经典组分模型描述页岩蠕变，包括Maxwell模型、Kelvin模型、Kelvin-Voigt模型等。但上述模型的精度目前已无法满足数值模拟的需求。近年来，众多学者围绕岩石蠕变行为的反演模拟开展了研究。Sheng等人[14]对比了几种拉普拉斯变换的数值反演方法。2016年，Meng等人[15]对时间依赖力学行为进行了分数阶描述。同年，周仁战和姚兆明[16]采用分数阶Burgers模型和模拟退火算法，分析了温度、应力水平对人工冻土黏弹性行为的影响。2017年，唐皓和王东坡[17]建立了基于分数阶导数的三元蠕变模型，并用盐膏的蠕变曲线验证了其有效性。同年，Xu和Jiang[18]拟合了不同应力水平下高密度聚乙烯、聚醚醚酮的蠕变，得到了不同分数阶模型的拟合结果。这些研究充分说明了分数阶模型在描述岩石蠕变过程中的有效性和可行性。

本文利用分数阶微积分原理建立了分数阶蠕变本构方程，并结合岩石蠕变实验数据对本构方程的准确性进行了验证；而后基于分数阶粘弹对应性原理建立了考虑支撑剂嵌入、变形和岩石蠕变的导流能力损伤模型；最终运用嵌入式离散裂缝方法建立了考虑导流能力损伤的气-水两相产能模型。通过模拟分析，讨论了不同导流能力损伤模式下气井产能的变化规律，确定了可用于表征页岩蠕变行为强弱的指标参数。

1 分数阶蠕变本构方程

黏弹性元件主要有：弹簧、黏壶、弹簧壶（图1）。其中，弹簧壶元件对于低速蠕变的拟合效果较好，对于岩石这一类明显偏向于弹性的黏弹性材料需要采用弹簧壶元件来进行模拟。

（a）弹簧元件 （b）黏壶元件 （c）弹簧壶元件

图 1 黏弹性元件

Figure 1 Viscoelastic element

分数阶Kelvin模型是指用弹簧壶元件替代经典Kelvin模型中的粘壶元件，分数阶Kelvin模型的本构方程为：

（1）

式中，——黏性模量，MPa；

——应力，MPa；

——应变，无量纲；

——分数阶稠度系数，MPa·sα；

——表示积分下限为0，上限为t的α阶分数阶微分算子。

结合Mittag-Leffler函数的拉普拉斯变换性质，可以对拉普拉斯空间中分数阶Kelvin模型的蠕变柔量进行反演得到实空间中的蠕变柔量：

（2）

式中，——分数阶Kelvin模型的蠕变柔量，MPa-1；

——分数阶弛豫时间，sα；

——求导阶数；

——特殊函数。

为验证分数阶Kelvin模型在描述岩石蠕变中的适应性，本文分别运用经典Kelvin蠕变模型和分数阶Kelvin蠕变模型与岩石蠕变数据进行拟合。拟合效果如图2所示。

图 2 完全饱和岩样的拟合对比结果

Figure 2 The fitting comparison results of fully saturated rock samples

由以上结果可以分析得出，分数阶Kelvin模型的拟合精度远高于经典Kelvin模型的拟合精度。因此，运用本文提出的分数阶Kelvin模型与岩石蠕变实验数据进行拟合，可以得到符合目标储层的岩石蠕变本构方程，以及分数阶弛豫时间、求导阶数、黏性模量等参数的具体数值。

2 考虑支撑剂变形、嵌入和岩石蠕变的导流能力损伤模型

（1）考虑支撑剂变形、嵌入的导流能力损伤模型

假设支撑剂多层铺置情况如图3所示，中间的球体假设为支撑剂，两侧挤压的平面假设为裂缝壁面。

图 3 支撑剂多层铺置情况示意图

Figure 3 Schematic diagram of multi-layer placement of proppant

有支撑剂变形引起的裂缝宽度变化计算模型为：

（3）

式中，——支撑剂球体的变形量，mm；

——闭合压力，MPa；

——支撑剂粒径，mm；

——支撑剂泊松比；

——支撑剂体积模量，MPa；

——支撑剂杨氏模量，MPa。

本文采用支撑剂单层铺置情况下的嵌入深度计算模型，嵌入模型如下：

（4）

式中，——裂缝壁面储层厚度，mm；

——储层杨氏模量，MPa。

支撑缝的宽度变化包括了支撑剂的变形和支撑剂嵌入两个部分，为了使计算结果更接近于实际，需要同时考虑支撑剂变形和嵌入带来的影响。因此，本文采用了如下模型进行计算：

（5）

式中，——缝宽变化量，mm。

在上述缝宽计算模型的基于上，可进一步得到考虑支撑剂变形、嵌入的导流能力损伤模型：

（6）

式中，——裂缝导流能力，D·cm；

——人工裂缝的渗透率，D；

——闭合应力下的缝宽，cm；

——闭合压裂为0时的孔隙度；

——闭合压力为0时的孔隙喉道的半径，；

——闭合压力为0时的孔隙迂曲度。

（2）考虑支撑剂变形、嵌入和岩石蠕变的导流能力损伤模型

将上述导流能力衰减模型和岩石蠕变本构方程结合，根据分数阶粘弹对应性原理，可以得到考虑支撑剂变形、嵌入和岩石蠕变的导流能力损伤模型：

（7）

式中，——拉氏空间中的储层岩石弹性模量，MPa；

——拉氏空间中的压力响应，MPa；

——储层岩石的弹性模量，MPa；

——储层岩石的黏性模量，MPa；

——分数阶弛豫时间，sα；

s——拉氏变量，s-1。

该裂缝导流能力损伤模型中的分数阶弛豫时间、求导阶数、黏性模量参数均可通过拟合得到的岩石蠕变本构方程获得。

3 考虑导流能力损伤的产能模型

本文结合了嵌入式离散裂缝模型和双重介质模型，运用已建立的导流能力损伤模型获取裂缝的导流能力，从而得到考虑导流能力损伤的气-水两相产能模型。

（1）主裂缝的渗流模型

在忽略毛管力的基础上，根据连续性方程和运动方程可以得到裂缝中液相的渗流微分方程：

（8）

式中，——裂缝液相相对渗透率，无量纲；

——水相的粘度，mPa·s；

——水相体积系数，m3/m3；

——主裂缝中水相的源汇项，m3/s；

——主裂缝与微裂缝网络网格之间的水相窜流量，m3/s；

——主裂缝孔隙度；

——主裂缝中水相饱和度，无量纲;

——主裂缝导流能力，通过考虑支撑剂变形、嵌入和岩石蠕变的导流能力损伤模型计算获得，mD。

气相的渗流微分方程：

（9）

式中，——裂缝气相相对渗透率，无量纲；

——裂缝中气相的源汇项，m3/s；

——裂缝中气相饱和度，无量纲；

——主裂缝中的气相饱和度；

——气体的体积系数。

（2）微裂缝的渗流模型

微裂缝中同时存在着主裂缝-微裂缝和微裂缝-基质两种渗流［19，20］。微裂缝中的水相渗流方程为：

（10）

式中，——微裂缝水相相对渗透率；

——基质水相相对渗透率；

——微裂缝水相饱和度；

——方程系数，当微裂缝网络网格中没有主裂缝嵌入时，取0，反之取1；

——微裂缝孔隙度。

微裂缝中气相渗流方程为：

（11）

式中。——微裂缝中气相相对渗透率；

——基质中气相相对渗透率；

——微裂缝中的气相饱和度。

（3）基质中的渗流模型

基质水相渗流微分方程：

（12）

式中，——基质中水相饱和度；

——基质孔隙度；

——方程系数，当基质网格中没有主裂缝嵌入时，取0，反之取1。

基质气相渗流微分方程：

（13）

式中，——基质中气相饱和度。

（4）边界条件

内边界条件为井底流压：

（14）

式中，，——井的坐标；

——井底流压，MPa。

模型中外边界为非渗透边界，外边界条件为：

（15）

（5）初始条件

初始压力取原始地层压力，且每个位置取相同的初始液相饱和度。

4 模型求解

先对裂缝和油藏基质中的渗流微分方程进行有限差分离散，再联立压力和饱和度求解差分方程组。

裂缝系统的有限差分方程用矩阵形式表示为：

（16）

式中，——为裂缝系统传导率矩阵；

——为裂缝系统累积项矩阵；

——为相交主裂缝之间的流量交换矩阵；

——为基质和微裂缝网络之间的窜流项矩阵；

——为微裂缝网络和主裂缝之间的窜流项矩阵；

——为源汇向量；

——为裂缝系统差分方程处理中的常数向量。

基质系统的有限差分方程用矩阵形式表示为：

（17）

式中：——为基质系统传导率矩阵；

——为基质系统累积项矩阵；

——为基质系统差分方程处理中的常数向量。

本文模型求解过程中，先求解裂缝系统中的压力分布和饱和度分布，然后求解基质系统中的压力分布和饱和度分布。

5 模型分析

5.1 模拟验证

为了验证模型的准确性，选取目标区块的储层物性参数进行模拟。模型验证参数如表1所示。

表 1 模型所需参数

Table 1 The model parameters

储层压力（MPa）	20	储层温度（K）	328
储层大小（m×m）	700×340	SRV大小（m×m）	540×220
基质孔隙度	0.08	主裂缝孔隙度	0.8
微裂缝网络孔隙度	0.00022	微裂缝网络渗透率（D）	2.08×10-6
主裂缝渗透率（D）	8.33	基质渗透率（D）	2.6×10-4
基质密度（kg/m3）	2500	地层厚度（m）	10
岩石压缩系数（MPa-1）	0.0001	井径（m）	0.05
主裂缝半长（m）	90	主裂缝宽度（m）	0.001

设定井底流压为10MPa，得到了本文模型的累积产气量和文献中累积产气量的对比曲线（图4），从图中可以看出，两条累积产气曲线相对接近，累积曲线之间的差异会越来越大，这也符合一般规律，从而验证了本文模型的正确性。

图 4 模型验证

Figure 4 Model validation

5.2 模拟结果分析

利用本文建立的产能模型分别对不同导流能力损伤情况下的产气情况进行模拟分析，得到了生产30天后的日产量曲线（图5）、累积产量曲线（图6）。

图 5 不同导流能力衰减模式的日产气量对比

Figure 5 Comparison of daily gas production under different conductivity decline modes

图 6 不同导流能力衰减模式的累积产气量对比

Figure 6 Comparison of cumulative gas production for different conductivity decline patterns

由图5和图6可知，只考虑支撑剂嵌入和变形影响时的产能模拟结果与未考虑导流能力损伤时的模拟结果差异较小；而当模型考虑了岩石蠕变的影响后，产能的模拟结果明显低于其他两种情况下的模拟结果。因此，由岩石蠕变引起的导流能力损伤会强于支撑剂嵌入和变形带来的影响。同时，考虑了支撑剂嵌入、变形和岩石蠕变影响的日产量模拟结果与其他两种情况的模拟结果差异有着先增大后减小的趋势。

6 蠕变损伤评价

对于不同的储层，其岩石蠕变特性会存在一定差异，这使得岩石蠕变引起的导流能力损伤程度强弱也会不同［21-23］。由于岩石蠕变本构方程中的分数阶弛豫时间、求导阶数、黏性模量三个参数是由分数阶Kelvin模型与岩石蠕变实验数据拟合得到的，因此，这三个参数的数值能够在一定程度上体现储层岩石自身的蠕变性。为此，本文分析对比了不同分数阶弛豫时间、不同求导阶数、不同黏性模量情况下的产能。

图 7 不同分数阶弛豫时间的累积产气量对比

Figure 7 Comparison of cumulative gas production with different fractional-order relaxation times

图 8 不同求导阶数的累积产气量对比

Figure 8 Comparison of cumulative gas production at different derivative orders

图 9 不同黏性模量的累积产气量对比

Figure 9 Comparison of cumulative gas production with different viscosity modulus

通过对上述模拟结果的对比分析可知，岩石蠕变本构方程中的分数阶弛豫时间和求导阶数两个参数对气井的产能影响较小；而黏性模量却对产能有着较为明显的影响。当储层黏性模量越高时，岩石蠕变越不明显，导流能力损伤程度也越弱，这种情况下便能获得较高的产能。若储层黏性模量增加到一个较大的数值时，由岩石蠕变引起的导流能力损伤便可忽略。因此，黏性模量能够用于表征储层蠕变特征，通过岩心实验获取黏性模量便可明确目标储层的蠕变强弱。

7 结论

（1）本文通过分数阶微积分原理和粘弹性理建立了分数阶Kelvin模型，运用该模型与岩石蠕变实验数据进行拟合，验证了本文提出的分数阶Kelvin模型的准确性。同时建立了综合考虑支撑剂变形、嵌入和岩石蠕变的气-水两相产能模型，并验证了该模型的准确性。

（2）通过对不同导流能力损伤情况下的气井产能进行模拟分析得出，由岩石蠕变引起的导流能力损伤强于支撑剂嵌入和变形带来的影响；同时，考虑了支撑剂嵌入、变形和岩石蠕变影响的日产量模拟结果与只考虑支撑剂嵌入变形和不考虑导流能力衰减两种情况下的模拟结果的差异有着先增大后减小的趋势。

（3）分析对比了不同分数阶弛豫时间、不同求导阶数、不同黏性模量对气井产能的影响得出，岩石蠕变本构方程中的分数阶弛豫时间和求导阶数两个参数对气井的产能影响较小；而黏性模量却对产能有着较为明显的影响，可作为储层蠕变行为的表征参数。

参考文献

[1] Meng Y，Wang Z，Zhang L，et al．Research Article Experimental Evaluation on the Conductivity of Branch Fracture with Low Sand Laying Concentration and Its Influencing Factors in Shale Oil Reservoirs［J］．Lithosphere，2021，1：4500227.

[2] Wei J，Zhou X，Fu X，et al．Experimental investigation of long-term fracture conductivity filled with quartz sand：Mixing proppants and closing pressure［J］．International Journal of Hydrogen Energy，2021，46（64）：32394-32402.

[3] Zhang J，Ouyang L，Zhu D，et al．Experimental and numerical studies of reduced fracture conductivity due to proppant embedment in the shale reservoir［J］．Journal of Petroleum Science and Engineering，2015，130：37-45.

[4] Zhang J，Kamenov A，Zhu D，et al．Laboratory Measurement of Hydraulic-Fracture Conductivities in the Barnett Shale［J］．Spe Production & Operations，2014，29（3）：1930-1855.

[5] Liu Y，Yang C，Wan J，et al．New insights into hydration-induced creep behavior of shale：A comparison study of brittle black shale and clayey oil shale at micro-scale［J］．Marine and Petroleum Geology，2022，138：105554.

[6] Song J，Xiang D，Hu D，et al．Creep characteristics of a fracturing fluid-softened shale investigated by microindentation［J］．International Journal of Rock Mechanics and Mining Sciences，2022，152：105067.

[7] Cao Y，Deng J，Yu B，et al．Analysis of sandstone creep and wellbore instability prevention［J］．Journal of Natural Gas Science and Engineering，2014，19：237-243.

[8] Firme P，Roehl D M，Romanel C，et al．Creep constitutive modeling applied to the stability of pre-salt wellbores through salt layers［C］．48th US Rock Mechanics/Geomechanics Symposium，American Rock Mechanics Association，2014.

[9] Orlic B，Buijze L．Numerical modeling of wellbore closure by the creep of rock salt caprocks［C］．48th US Rock Mechanics / Geomechanics Symposium，American Rock Mechanics Association，2014.

[10] Almasoodi M M，Abousleiman Y N，Hoang S K．Viscoelastic creep of eagle ford shale：investigating fluid/shale interaction［J］．Journal of Canadian Petroleum Technology，2015，54 （3）：142-143.

[11] Rassouli F S，Zoback M D．A comparison of short-term and long-term creep experiments in unconventional reservoir formations［C］．50th US Rock Mechanics / Geomechanics Symposium，American Rock Mechanics Association，2016.

[12] 王萍．泥页岩井壁水化损伤的蠕变失稳力学研究［D］．西北工业大学，2015.

[13] Li K，Gao Y，Lyu Y，et al．New mathematical models for calculating the proppant embedment and fracture conductivity［J］．Spe Journal，2015，20（3）：496-507.

[14] Sheng H，Li Y，Chen Y．Application of numerical inverse laplace transform algorithms in fractional calculus［J］．Journal of the Franklin Institute，2011，348 （2）：315-330.

[15] Meng R，Yin D，Zhou C，et al．Fractional description of time-dependent mechanical property evolution in materials with strain softening behavior［J］．Applied Mathematical Modelling，2016，40 （1）：398-406.

[16] 周仁战，姚兆明．基于模拟退火分数阶导数伯格斯模型的人工冻土蠕变研究［J］．安徽理工大学学报（自科版），2016，36 （6）：30-38.

[17] 唐皓，王东坡，裴向军，等．基于分段模拟的岩石蠕变模型［J］．水文地质工程地质，2017，44 （1）：41-47.

[18] Xu H，Jiang X．Creep constitutive models for viscoelastic materials based on fractional derivatives［J］．Computers & Mathematics with Applications，2017，73 （6）：1377-1384.

[19] Luo A，He J，Li J，et al．Experimental Study on the Influence of Water-rock Interaction on the Mechanical Characteristics and Creep Behavior of Shale［J］．Journal of Geo-Energy and Environment，2025，1（2）：61-69.

[20] Gong X，He J，Li J，et al．Professional Evaluation and Distribution Patterns of Shale Gas Reservoirs in the Wufeng Formation-Long 11 Sub-member of Well Block Z205，Sichuan Basin［J］．Journal of Geo-Energy and Environment，2025，1（2）：96-105.

[21] He S，Tan W C，Li Hu，et al．Mineralogical and lithofacies controls on gas storage mechanisms in organic-rich marine shales［J］．Energy & Fuels，2025，39（8）：3846-3858.

[22] Li H．Quantitative prediction of complex tectonic fractures in the tight sandstone reservoirs：a fractal method［J］．Arabian Journal of Geosciences，2021，14（19）：1986.

[23] Ma G，Huang R，Dong Y，et al．AHigh-Accuracy Cost Prediction Model for Shale Gas Drilling in Southern Sichuan Using PCA and BP Neural Network［J］．Journal of Geo-Energy and Environment，2026，2（1），46-55.