CO₂ 注入引起的流固耦合作用对井筒损伤研究

Study on Wellbore Damage Caused by CO₂ Injection‒Induced Hydromechanical Coupling

摘要

在CO₂地质封存过程中，井筒周围受到多场耦合作用，易产生不同程度的损伤，进而成为潜在的CO₂泄漏通道。为了研究CO₂注入对井筒完整性的影响，建立了两相渗流流固耦合模型，分析了储层岩石初始孔隙度、渗透率、弹性模量，以及应力比和注入速率对CO₂运移和井筒、储层岩石的损伤影响。结果表明：初始孔渗大小和注入速率对储层岩石和水泥环损伤影响显著；较低的初始孔渗和较大的注入速率更容易造成井筒损伤。采用库仑失效准则评估储层岩石和水泥环的损伤演变。分析发现：储层岩石在与水泥环界面处容易损伤，而水泥环在与套管界面处更易发生损伤而失效，损伤可能由最初的点损伤演变为面损伤；对于现场注入CO₂，如果维持较高的注入速率，则需要改善低渗透性储层的孔渗条件，降低井筒失效风险。

Abstract

Objective Wellbores can be damaged due to multiphysical coupling during geological CO₂ storage in deep saline aquifers, which compromises wellbore integrity and creates potential leakage pathways for CO₂. This study proposes a fully coupled hydromechanical model to simulate the two-phase flow of brine and CO₂ near the wellbore and to investigate the mechanical responses of the reservoir rock and cement sheath. Methods The wellbore structure in the model consisted of casing, a cement sheath, and the surrounding reservoir rock. The reservoir was considered a saturated, homogeneous, isotropic elastic medium, with deformation characterized as small, linear elastic deformation. Pore pressure diffusion was governed by the Biot storage model, while the two-phase flow for brine and CO₂ was described by Darcy’s law and the Brooks‒Corey model. Since pore pressure increase led to wellbore damage during CO₂ injection, the indirect variation in reservoir porosity and permeability was also considered, showing a correlation with changes in average effective stress. The mechanical response and pore pressure variation were fully coupled, and the Coulomb failure criterion was utilized to evaluate the damage risk of both the reservoir rock and the cement sheath. It was found that damage to the cement sheath depended on its inherent properties, the surrounding stress state, and changes in pore pressure. Therefore, the effects of initial porosity, permeability, and elastic modulus of the reservoir rock, as well as the stress ratio and injection rate, on CO₂ migration and damage to the wellbore and reservoir rock were analyzed. The initial porosity and permeability of the reservoir rock were described using the Weibull distribution. Results and Discussions The influence of various factors on CO₂ migration and damage to the wellbore and reservoir rock was examined. Analyses of the benchmark model showed that CO₂ saturation and average pressure distributions were closely related, with the highest pore pressure occurring near the injection well. Pore pressure changed significantly from 0.5 to 4.0 days, with maximum average pressure changes initially reaching 5 MPa and increasing up to 15 MPa as the pore pressure diffused. This increase reflected the influence of boundary conditions on the migration of CO₂. Heterogeneity in reservoir porosity and permeability significantly affected distributions of both pore pressure and CO₂. Initial CO₂ saturation exhibited a distinctly nonuniform pattern. As CO₂ migrated, its distribution became more uniform but still depended on variations in porosity and permeability. Previous field simulations indicated that reservoir heterogeneity and the presence of fractures can significantly alter CO₂ migration patterns. For the wellbore model, despite boundary constraints limiting further pressure diffusion, heterogeneity impacted CO₂ migration and pore pressure distribution, influencing potential damage to the wellbore and surrounding areas. It was found that the injection rate, reservoir porosity, and permeability controlled CO₂ migration, with the stress ratio and elastic modulus being largely negligible, as indicated by analyzing CO₂ saturation and pore pressure at 4 days. In cases of low porosity and permeability, a higher injection pressure was required, resulting in greater CO₂ saturation and extended migration distances, though excessive pressure posed a risk of fracturing. Since CO₂ migration was linked to pore pressure diffusion, pressure-induced damage was mainly influenced by injection rate and changes in reservoir characteristics. Damage to the reservoir rock was attributed to increased porosity and permeability resulting from CO₂ injection, as well as mechanical damage at the cement sheath interface. Analysis showed that porosity changes induced by effective stress variation ranged between 0 and 0.001. Permeability changes reached up to 0.03 mD and were even higher locally, with local increments reaching up to 0.06 mD. The Coulomb failure criterion indicated varying degrees of damage under different conditions. Scenarios such as low injection rate and high permeability delayed damage onset beyond 1 day, while low permeability or high injection rates accelerated damage, potentially fracturing the rocks. Damage initially appeared as point damage, potentially evolving into surface damage with continuous injection, occurring at the reservoir-cement sheath interface. The integrity of the cement sheath depended on its properties, the surrounding stress conditions, and variations in pore pressure. For high CO₂ injection volumes, increasing the reservoir's porosity and permeability reduced pressure buildup. Visual analysis revealed damage at the interface between the cement sheath and the casing. Nonuniform changes in the Coulomb failure criterion at these interfaces, influenced by heterogeneity in reservoir rock properties, can result in radial cracks within the cement sheath. Conclusions This study employs a two-phase hydromechanical coupling model to assess the impact of multiple factors on wellbore and surrounding rock damage caused by CO₂ injection. The findings indicate that CO₂ migration is aligned with pore pressure diffusion, resulting in both compression and expansion of internal pores. Although reservoir heterogeneity influences CO₂ migration, the injection rate and initial permeability play a more critical role in determining CO₂ flow and pore pressure distribution. Damage typically occurs at the reservoir-cement sheath interface and the cement-casing interface. Reducing injection rates and increasing reservoir permeability can help mitigate damage to reservoir rocks and the cement sheath. The elastic modulus and stress ratio have limited influence on damage. As the damage evolves from point to surface, enhancing permeability in low-permeability reservoirs is recommended to sustain high CO₂ injection rates and reduce the risk of cement sheath failure.

Graphical abstract

关键词

流固耦合 / 两相渗流 / 井筒损伤 / CO₂注入 / 水泥环失效 / 非均质性

Key words

hydromechanical coupling / two-phase flow / wellbore damage / CO₂ injection / cement sheath failure / heterogeneity

本刊网刊

CO₂过度排放导致全球温室效应显著，由此引发的次生灾害，如海平面上升和极端天气效应，严重影响人类的生存环境。CO₂地质封存（GCS）旨在将多余CO₂注入到地质体中，通过不同的封存机制进行封存（如构造封存、残余封存、溶解封存和矿化封存），从而减少大气中CO₂含量^[1‒3]。超临界CO₂通过井筒被注入到地层中，将导致井筒周围的温度场、渗流场、应力场和化学场（THMC）发生改变。井筒的完整性受到多场耦合作用影响可能产生不同程度的损伤，因此，井筒也成为潜在的CO₂泄漏通道，需要对其泄漏风险进行评估分析^[4‒6]。

对井筒损伤的研究在地下能源开采方面具有重要意义，如页岩气开采^[7]和地热能提取^[8]，相关的研究理论和结果在CO₂地质封存中也同样适用。不同点在于，能源开采过程中需要通过水力压裂等手段增加储层岩石的渗透性，提高产量，由此可能引发井筒损伤。而CO₂地质封存过程中则需要保证井筒的完整性，确保CO₂不会沿井筒泄漏。

目前，许多研究关注耦合效应造成的井筒损伤。艾池等^[9]建立了完整的井筒模型，通过理论计算分析了温降和孔压变化对水泥环的影响，结果表明，水泥环和温度同时作用会增加水泥环损伤的风险，而单独作用时风险减小。Yan等^[10]则结合试验和模拟开展了水力压裂过程中射孔水泥环完整性的耦合分析，发现极低的套管压力有利于保持水泥环的完整性，并且只有同时优化水泥环的弹性模量和泊松比才能较好地维持水泥环的完整性。Li等^[11]基于模拟开展了相似的研究，主要分析了参数非均质性对水泥环的影响。而Samarakoon等^[12]则关注于温度和压力耦合作用下不同地层中固井水泥的热力学行为。李琦等^[13]基于耦合‒去耦合理论，建立了非均质岩石中注入CO₂的热孔弹耦合模型；结果表明，超临界CO₂注入会导致一定程度损伤，造成渗透率增加，且非均质岩石相比均质岩石更容易产生损伤；但是该研究中并未考虑套管、水泥环和岩石之间的相互作用。

井筒的强度是决定其是否发生损伤的重要因素。由于涉及流体注入，井筒周围岩石的孔渗条件和完整性同样会影响井筒的损伤程度。例如：Li等^[14]进一步地考虑了储层岩石中天然裂缝对井筒稳定性的影响；Peng等^[15]则考虑了CO₂注入过程中套管腐蚀对井筒完整性的影响，但这是一种长期效应；尹飞等^[16]从力学角度总结了储层压实可能导致的井筒失效机制，主要包括上覆岩层井筒拉伸破坏、套管失稳屈曲、套管剪切变形、套管径向变形和挤毁。CO₂注入的增压过程中可能遇到相似的问题。

通过上述研究发现，影响井筒完整性的因素众多。一方面，虽然目前许多研究基于耦合作用来分析井筒的损伤，然而，大多数研究中多场效应并未完全耦合，同时，缺乏对CO₂注入过程中井筒结构之间相互影响的研究。另一方面，虽然温度也会对损伤产生影响，但注入之后的力学损伤行为依然是流固耦合效应占主导^[17]。因此，本文结合两相渗流理论和流固耦合理论，建立了两相渗流流固耦合模型。模型考虑套管、水泥环和储层岩石之间的相互作用及其随CO₂注入的演变，对CO₂注入过程中的井筒损伤进行评估。进一步分析了不同条件下的孔渗变化和井筒损伤情况，为井筒完整性评价提供参考依据。

1 井筒模型建立与损伤评估

1.1 井筒模型描述

本文以宁夏盐池地区延长组底部地层为研究目标区域^[18]。延长组底部埋深约为3 km，地层温度约95 ℃。宁夏盐池县地处中国西北地区东部，鄂尔多斯盆地西缘。其地形复杂，总体呈南部高、北部低、中部高、东西两侧低的特点，局部地区起伏较大。该地层交错分布多层泥岩和砂岩，致密储层渗透率基本小于1×10^-16 m²，孔隙度小于0.1，弹性模量在10～45 GPa之间^[18]。

为了便于分析井筒损伤，建立2维井筒模型。井筒模型主要结构包含套管、水泥环和储层岩石，示意图如图1所示。图1中，σ_x 和σ_y 分别为x、y方向应力。套管选择N80套管，其内径为121.4 mm，外径为139.7 mm，厚度为9.15 mm。水泥环的厚度为54.35 mm，孔隙度和渗透率分别为0.01和1×10^-19 m²。套管、水泥环和储层岩石的物性参数如表1所示。完整模型为宽度10 m的正方形，选择此宽度主要是为了削减边界条件对CO₂运移的影响。

考虑到储层物性的非均质性，采用Weibull分布来生成储层岩石初始状态的孔隙度、渗透率和弹性模量。Weibull分布的概率密度函数f (β)表示如下^[19]：

f(\beta )=\frac{m}{a}{\left(\frac{\beta }{a}\right)}^{m-\mathrm{1}}·\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{\beta }{a}\right)}^m\right)（1）

式中：β为目标物性参数值；a为目标物性参数的平均值；m为非均匀指数，描述Weibull分布的形状，其值越小目标物性参数变化越大，这里取m=4^[20]。对于水泥环和套管，只考虑为均质的模型。

本研究主要分析储层岩石的非均质性（包括孔隙度θ、渗透率k和弹性模量E）、应力比R及注入速率Q对井筒及其周围损伤的影响，即储层和水泥环的损伤。在改变参数时，每次只变化1个参数，其余参数取基准值，并且孔隙度和渗透率同时发生变化。应力比为x方向应力σ_x 与y方向的应力σ_y 比值。考虑到模型埋深3 km，垂直应力约为80 MPa，所以假定水平应力σ_x 与垂直应力σ_y 的比值为0.7，σ_y 为56 MPa，因此，不同的应力比对应于不同的σ_x 变化。基准模型注入速率为0.025 kg/(m²∙s)，假定注入深度为1 km，则该速率对应的CO₂注入量约6×10⁵ t/a。模型#1、#2、#3中改变的参数如表2所示，参数值为目标物性参数的平均值，其中#2为基准模型。不同模型的初始孔隙度、渗透率和弹性模量分布如图2所示。图2中，▲表示最大值，▼表示最小值，下文同。

1.2 流固耦合模型

CO₂注入储层将导致流体与固体之间相互作用。Biot理论建立了固体变形和孔隙压力扩散之间的关系。该理论假定多孔介质为饱和均质、各向同性的弹性介质，变形为线弹性小变形，流体流动服从达西定律（不考虑惯性力的影响），并且单位时间内多孔介质单元体水量的变化量等于介质体积的变化量。因此，孔压的变化可由下述方程描述^[21]：

{\rho }_{\mathrm{f}}S\frac{\partial p}{\partial t}-\nabla \left(\frac{{\rho }_{\mathrm{f}}k}{{\mu }_{\mathrm{d}}}\nabla p\right)=Q_{\mathrm{m}}-{\rho }_{\mathrm{f}}\alpha \frac{\partial {\epsilon }_{\mathrm{v}\mathrm{o}\mathrm{l}}}{\partial t}（2）

S=\theta {\chi }_{\mathrm{f}}+\frac{(\mathrm{1}-\alpha )(\alpha -\theta )}{K_{\mathrm{d}}}（3）

K_{\mathrm{d}}=\frac{\mathrm{2}G(\mathrm{1}+\upsilon )}{\mathrm{3}(\mathrm{1}-\mathrm{2}\upsilon )}（4）

式（2）～（4）中，ρ_f为流体密度，S为储水系数，p为孔隙压力（当考虑两相渗流时，其大小为两相流体压力的饱和度平均值），t为时间，\nabla为Nabla算子，k为渗透率，μ_d为动力黏度，Q_m为源项，α为比奥系数，ε_vol为体积应变，θ为孔隙度，χ_f为流体压缩系数，K_d为体积模量，G为剪切模量，υ为泊松比。

对于固体变形，其力学平衡方程为：

\nabla \left(\mathit{\sigma }\text{'}-\alpha p\mathit{I}\right)+\mathit{F}=\mathrm{0}（5）

式中， σ'为有效应力张量， I 为单位矩阵， F 为体积力。

对于流固耦合模型，渗透率和孔隙度随流体注入发生改变，这里采用下述方程进行描述^[22]：

\theta =\alpha +({\theta }_{\mathrm{0}}-\alpha )\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{\Delta }{\sigma }_{\mathrm{e}\mathrm{f}\mathrm{f}}}{K_{\mathrm{d}}}\right)（6）

{\sigma }_{\mathrm{e}\mathrm{f}\mathrm{f}}=({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}})/\mathrm{3}+\alpha p（7）

k=k_{\mathrm{0}}(\theta /{\theta }_{\mathrm{0}}{)}^n（8）

式（6）～（8）中：θ₀和k₀分别为初始孔隙度和初始渗透率；∆σ_eff为平均有效应力σ_eff的变化量；σ₁、σ₂和σ₃分别为3个方向的主应力；n为依赖于岩石性质的无量纲常数，这里取n=3。力的方向设置为拉伸为正，压缩为负。

1.3 两相渗流模型

考虑超临界CO₂注入到咸水层的过程，由连续性方程可以得到^[23]：

\frac{\partial (\theta {\rho }_{\mathrm{f}i}{s}_i)}{\partial t}=\nabla \cdot \left({\rho }_{\mathrm{f}i}\frac{k{k}_{\mathrm{r}i}}{{\mu }_{\mathrm{d}i}}\nabla p_i\right)（9）

式中：下标i可以为w和nw，分别表示湿润相流体（咸水）和非湿润相流体（超临界CO₂）；s为相饱和度；k_r为相对渗透率。这里采用Brooks‒Corey模型来描述相对渗透率和毛细压力^[24]：

k_{\mathrm{r}\mathrm{w}}={\left(\frac{s_{\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}{\mathrm{1}-s_{\mathrm{r}\mathrm{n}\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}\right)}^{\mathrm{3}+\frac{\mathrm{2}}{\lambda }}（10）

k_{\mathrm{r}\mathrm{n}\mathrm{w}}={\left(\frac{s_{\mathrm{n}\mathrm{w}}-s_{\mathrm{r}\mathrm{n}\mathrm{w}}}{\mathrm{1}-s_{\mathrm{r}\mathrm{n}\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}\right)}^{\mathrm{2}}\left(\mathrm{1}-{\left(\frac{s_{\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}{\mathrm{1}-s_{\mathrm{r}\mathrm{n}\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}\right)}^{\mathrm{1}+\frac{\mathrm{2}}{\lambda }}\right)（11）

p_{\mathrm{c}}=p_{\mathrm{e}\mathrm{c}}{\left(\frac{s_{\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}{\mathrm{1}-s_{\mathrm{r}\mathrm{n}\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}\right)}^{-\frac{\mathrm{1}}{\lambda }}（12）

式（10）～（12）中，k_rw、s_w和s_rw分别为湿润相流体的相对渗透率、饱和度和残余饱和度，k_rnw、s_nw和s_rnw分别为非湿润相流体的相对渗透率、饱和度和残余饱和度，λ为孔隙分布指数，p_c为毛细压力，p_ec为毛细进入压力。CO₂初始相饱和度为0，s_rw=0.2，s_rnw=0，λ=2，p_ec=0.02 MPa^[2]。

在CO₂注入过程中，流体孔隙压力增加导致固体介质变形，产生的应力反作用于孔隙压力，而CO₂的运移又受到孔隙压力扩散的影响。基于上述关系，建立了两相渗流流固耦合模型。上述公式通过有限元软件COMSOL Multiphysics进行完全耦合求解。采用的物理场包括多孔介质相传递、达西定律和固体力学。多孔介质相传递和达西定律耦合计算CO₂运移，而达西定律和固体力学耦合计算CO₂注入引起的应力变化和变形。孔渗的变化则用函数定义。在计算前进行应力平衡，总计算时间t为4 d。

对于力学边界，左侧和底部为辊支撑，右侧和上部为荷载边界。对于流体边界，初始孔隙压力设置为30 MPa，模型四周均设置为无流动边界，主要考虑孔压积聚对井筒造成的损伤，对应于最不利的情况。模型的计算网格采用三角形网格划分（图1），网格数量约2.7万个。储层岩石与水泥环之间考虑为整体，并对该处的网格进行加密，保证计算精度。由于只考虑储层岩石的流固耦合作用，CO₂注入位置设置在储层岩石和水泥环交界处，不考虑流体进入水泥环。

1.4 井筒损伤评估

为了分析在注入过程中水泥环和储层岩石发生损伤的可能性，这里采用库仑失效准则进行评估^[25]，失效函数如下：

\begin{array}{c}\psi =({\sigma }_{\mathrm{3}}+\mathrm{\Delta }{p}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}})-\frac{\mathrm{1}+\mathrm{s}\mathrm{i}\mathrm{n}\varphi }{\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\varphi }({\sigma }_{\mathrm{1}}+\mathrm{\Delta }{p}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}})+\\ c\frac{\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\varphi }{\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\varphi }\left(\mathrm{1}+\frac{{\sigma }_{\mathrm{2}}-{\sigma }_{\mathrm{1}}}{{\sigma }_{\mathrm{3}}-{\sigma }_{\mathrm{1}}}\right)\end{array}（13）

式中，\psi为失效值，ϕ为库仑摩擦角，c为库仑黏聚力，Δp_well表示孔隙压力p_well的变化。失效值\psi=0为临界点，\psi>0表示稳定，而\psi<0表示有失效风险，\psi值越小损伤风险越大。

2 影响因素分析

2.1 基准模型CO₂ 运移的影响机制

在分析不同因素对CO₂运移规律影响前，首先分析基准模型在不同时刻平均压力和CO₂饱和度的分布规律，其饱和分布云图如图3所示。平均压力p_aver计算式如下：

p_{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}=s_{\mathrm{w}}{p}_{\mathrm{w}}+s_{\mathrm{n}\mathrm{w}}{p}_{\mathrm{n}\mathrm{w}}（14）

式中，p_aver为等于各相的饱和度与压力乘积之和，p_w为湿润相流体压力，p_nw为非湿润相流体压力。时间选择第0.5、2.0和4.0 d。由图3可以发现，CO₂饱和度分布与平均压力分布基本一致。由于CO₂注入是增压过程，最大的孔隙压力位于井口处，孔压向四周扩散并逐渐减小。当t=0.5 d时，最大平均压力变化为5 MPa，此时，边界处的压力已经开始增加；当t=4.0 d时，最大平均压力变化达到15 MPa，压力扩散至整个模型，边界处的压力发生明显增长。可以预见此时CO₂运移一定程度上受到边界条件的影响。此外，储层孔渗的非均质性也会对孔压和CO₂饱和度分布造成一定影响。在非均质介质中，CO₂饱和度在初始段呈现出十分明显的非均匀分布，井口位置CO₂饱和度较高。随着不断扩散，CO₂运移分布逐渐趋于均匀，但是饱和度分布依然是非均匀的，取决于孔隙度和渗透率分布。

以往的场地尺度模拟表明储层孔渗的非均质性可能会严重影响CO₂的运移，储层中是否存在裂隙也会导致不一样的运移路径^[26‒27]。就井筒周围情况而言，即使由于边界条件限制，压力未能向远处进一步扩散，储层岩石孔渗的非均质性也会显著影响CO₂的运移和孔压分布，进而影响井筒及其周围的损伤情况。

为了分析不同因素对CO₂运移的影响，选择截线（图1）处t=4.0 d时CO₂饱和度和孔压进行分析，结果如图4所示。图4中结果表明：t=4.0 d时的CO₂运移范围主要受到注入速率和储层孔渗大小的影响。图例表示改变的参数（Q、R、E、θ、k）和参数值对应的模型编号（#1、#2、#3），如Q‒#1表示Q取模型#1设定的值。应力比和弹性模量对CO₂运移的影响几乎可以忽略。由于考虑了流固耦合作用，应力比主要通过改变变形影响孔隙大小，而弹性模量的作用与此相似，也是通过影响变形量的大小来影响孔渗变化。在相同的质量流注入情况下，孔隙度和渗透率越低（图4中θ‒k‒#1），需要更大的注入压力才能注入CO₂，由此导致CO₂饱和度相对较大，运移距离更远，但同时过大的压力可能会导致压裂^[28]。而相同条件下注入量的提升同样会增加CO₂运移距离。由于CO₂运移与压力扩散相关，因此注入速率和储层岩石孔渗大小是影响压力诱发损伤的主要因素。

2.2 储层岩石在不同条件下的损伤分析

储层岩石的损伤主要来自两方面；一方面是CO₂注入引起的储层岩石孔渗增加，另一方面是与水泥环界面接触所导致的损伤。首先，选取截线处的孔隙度和渗透率变化进行分析，变化值为t= 4 d与初始时刻两个变量的差值，结果分别如图5和6所示。根据式（6），孔隙度的变化来源于有效应力的变化，即多孔介质骨架压力和孔隙压力。CO₂注入将导致孔隙扩张，而扩张的孔隙可能会导致相邻的孔隙压缩。所以即使整体的孔隙度增加，内部依然存在孔隙度减小的情况。整体的孔隙度变化范围在0~0.001内，且θ‒k‒#1情况下的孔隙度增量最大。减小注入速率（Q‒#1），增加岩石孔渗（θ‒k‒#3）以及较大的弹性模量（E‒#3）都导致相对较小的孔隙度增量。根据式（8），渗透率的变化取决于孔隙度的变化，因此，图6中渗透率和孔隙度的变化趋势一致。在不同条件下，渗透率变化范围基本位于0～0.03 mD，局部的最大渗透率增量达到0.06 mD。

储层岩石孔渗的变化预示不同的损伤情况，根据库仑失效准则对储层岩石的损伤进行进一步评估。徐海强^[18]给出了樊学地区长7段岩石力学参数测试结果，同时给出了测井解释剖面图。这里选取该地区的岩石参数，岩石的黏聚力取为10 MPa，内摩擦角为40°。

首先分析最小失效函数值（ψ越小表明越可能发生破坏）随注入时间的变化。图7展示了不同条件下储层岩石最小失效函数值随时间的变化规律。图7中，ψ>0表示水泥环处于安全状态，以阴影部分表示。图7的结果表明，在相同条件下，最终储层内的岩石都会发生一定程度的损伤，但是发生损伤的时间和损伤程度受到不同因素的影响。对于低注入速率（Q‒#1），高孔渗条件（θ‒k‒#3）及低弹性模量（E‒#1）的岩石，发生损伤的时间均在注入1 d以后，其中，减小注入速率会显著降低损伤程度。而在其他条件下，CO₂注入迅速导致损伤的发生，产生损伤的速率取决于不同参数。低孔渗条件（θ‒k‒#1）和高注入速率（Q‒#3）形成十分显著的损伤趋势，二者共同影响可能会压裂储层岩石^[29]。

由于上述为失效函数的最小值变化，所以损伤在初始时刻可能为点损伤，随着不断注入，点损伤可能转化为面损伤。进一步提取截线处t=4 d时的失效函数进行分析，且只选择3个不同的注入速率工况（Q‒#1～#3），结果如图8所示。由图8可知，损伤ψ<0主要发生在注入井位置，注入速率越大，井口处越容易发生损伤，且整体也更倾向于产生损伤。但是在当前的模型中，储层岩石整体都未发生破坏。图8中的云图展示了注入井位置附近区域的失效函数分布，结果同样表明在注入井井口位置更容易发生破坏。失效函数整体分布规律与截线处一致。

2.3 水泥环在不同条件下的损伤分析

上述内容主要分析了储层岩石在不同情况下的损伤情况，发现损伤主要位于储层岩石与水泥环的交界处。因此，交界处也可能是水泥环最容易发生破坏的位置。为了分析水泥环在不同条件下的损伤情况，同样借助库仑失效准则来进行判断。图9为不同条件下水泥环失效函数值ψ随时间的变化情况。这里依然选取失效函数的最小值。由图9可知：在所有条件下，水泥环在t=4.0 d时都有发生损伤的可能性，但是发生损伤的时间和损伤程度不同，与储层岩石损伤相似。改变应力比、增大注入速率及储层孔渗较低情况都会导致水泥环迅速发生损伤，而减小注入速率和增加孔渗则推迟了损伤起始的时间。与储层岩石的损伤情况不同，弹性模量的减小并不会缓解水泥环的损伤。由于并未考虑CO₂进入水泥环，因此水泥环的损伤主要来源于和储层岩石接触产生的力学损伤，即来自主应力的变化；而储层岩石的损伤来自流固耦合作用，孔压变化的同时导致有效主应力变化，二者损伤机制存在差异。

根据库仑失效函数，可以发现水泥环是否失效取决于水泥环本身的性质（c和ϕ）、周围的应力状态（σ₁、σ₂和σ₃）及孔隙压力的变化（Δp_well）。因此，为了防止水泥环破坏，最直接的方法是加强水泥环的强度，即增加c和ϕ，提升其完整性。如果需要较高的CO₂注入量，则需要增加储层的孔隙度和渗透率，减少压力在井口处积聚。

为了直接观察水泥环的损伤情况，图10展示了失效函数在水泥环中的变化规律。选择最小注入速率情况，时间选择第0.1、0.5、1.0、2.0、3.0和4.0 d。可以发现，随着CO₂不断注入，在水泥环与套管界面处失效函数逐渐减小，在第4.0 d时界面处ψ<0，预示着发生损伤的可能性。因此，水泥环最容易发生损伤的位置位于其与套管的交界面处。水泥环与储层岩石和套管的界面处都容易发生损伤，这可能导致水泥环径向裂纹的扩展^[30]。此外，储层岩石的孔渗和弹性模量是非均匀的，这导致失效函数在水泥环与储层岩石界面处呈现出非均匀变化，同样会影响水泥环内部的损伤情况。

3 结 论

基于建立的两相渗流流固耦合模型，分析了不同因素对CO₂注入导致井筒及其周围岩石损伤的影响。主要得到以下结论：

1）CO₂运移与压力扩散规律基本一致，CO₂注入同时导致岩石孔隙的压缩和扩张。储层岩石孔渗非均质性会影响CO₂运移，但相比于其他因素，注入速率和初始孔渗大小对CO₂运移和孔压分布影响更显著。

2）储层岩石在与水泥环的交界处更容易损伤，而水泥环在与套管界面处更容易形成损伤。减少注入速率和提升储层孔渗条件能够减弱储层岩石和水泥环损伤程度。弹性模量和应力比对损伤影响并不明显。

3）如果要维持较高的CO₂注入量，损伤可能由最初的点损伤演变为面损伤。因此，可以增加低渗透储层岩石的孔渗以减少水泥环失效风险。

参考文献

[1] Li Qi, Song Ranran, Kuang Dongqin,et al.Status and advances of abandoned process of wells for CO₂ geological storage[J].Advances in Earth Science,2016,31(3):225‒235.

[2] 李琦,宋然然,匡冬琴,等.二氧化碳地质封存与利用工程废弃井技术的现状与进展[J].地球科学进展,2016,31(3),225‒235.

[3] Yang Yanlin.The implementation technology of CO₂ geological storage in complex geological structure and applying it to Orods Basin[D].Changchun:Jilin University,2014.

[4] 杨艳林.CO₂地质储存中地质特征实现技术与应用——以鄂尔多斯CCS示范工程为例[D].长春:吉林大学,2014.

[5] Li Xiaochun, Li Qi, Bai Bing,et al.The geomechanics of Shenhua carbon dioxide capture and storage (CCS) demonstration project in Ordos Basin,China[J].Journal of Rock Mechanics and Geotechnical Engineering,2016,8(6):948‒966. doi:10.1016/j.jrmge.2016.07.002

[6] He Kai.Research on integrity evolution and leakage of CO₂ geological storage system[D].Daqing:Northeast Petroleum University,2019.

[7] 贺凯.CO₂地质封存系统完整性演化及其泄漏研究[D].大庆:东北石油大学,2019.

[8] Zhang Zhichao, Bai Mingxing, Gao Shuo,et al.Leakage risk assessment of geological CO₂ storage system[J].Petroleum Geology and Recovery Efficiency,2023,30(2):135‒143.

[9] 张志超,柏明星,高硕,等.CO₂地质封存系统泄漏风险评价[J].油气地质与采收率,2023, 30(2):135‒143.

[10] Song Li, Tan Yu.A two-dimensional model analysis of casing-cement-strata in geological storage of CO₂ [J].Chinese Journal of Applied Mechanics,2015,32(2):288‒293.

[11] 宋力,谭雨.二氧化碳地质埋存中套管‒水泥环‒地层二维模型分析[J].应用力学学报,2015,32(2):288‒293.

[12] Bai Bing, Chen Mian, Jin Yan,et al.The thermoporoelastic coupling analysis of wellbore stability in shale formation under supercritical CO₂ drilling conditions[J].Journal of Petroleum Science and Engineering,2022,220:111146. doi:10.1016/j.petrol.2022.111146

[13] Suo Yu, Dong Muyu, Fu Xiaofei,et al.Experimental and numerical simulation research on hot dry rock wellbore stability under different cooling methods[J].Geothermics,2024,119:102977. doi:10.1016/j.geothermics.2024.102977

[14] Ai Chi, Li Jing, Li Ziping,et al.Research on cement sheath stress integrity of CO₂ buried well in the process of injection[J].Science Technology and Engineering,2013,13(8):2057‒2061.

[15] 艾池,李静,李自平,等.CO₂埋存井注入过程中水泥环应力完整性研究[J].科学技术与工程,2013,13(8):2057‒2061.

[16] Yan Yan, Guan Zhichuan, Han Lihong,et al.Coupled thermo-hydro-mechanical analysis of perforated cement sheath integrity during hydraulic fracturing[J].Journal of Petroleum Science and Engineering,2022,218:110950. doi:10.1016/j.petrol.2022.110950

[17] Li Xiaorong, Zeng Yijin, Ding Zechen,et al.Cement sheath integrity in anisotropic shale formations:A numerical investigation[J].Geofluids,2022,2022:1‒16. doi:10.1155/2022/3553796

[18] Samarakoon M H,Ranjith P G Thermomechanical behaviour of well cement in different geological formations under the coupled effects of temperature and pressure[J].Geoenergy Science and Engineering,2023,223:211535. doi:10.1016/j.geoen.2023.211535

[19] Li Qi, Jing Miao.Thermo-poroelastic coupling analysis of rock damage around wellbore due to CO₂ injection[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(11):2205‒2213.

[20] 李琦,井淼.CO₂注入引起井筒周围岩石损伤的热孔弹耦合分析[J].岩石力学与工程学报,2013,32(11):2205‒2213.

[21] Li Xiaorong, Yang Henglin, Wei Yaoran,et al.Wellbore instability in naturally fractured formations:Experimental study and 3D numerical investigation[J].Gas Science and Engineering,2024,124:205265. doi:10.1016/j.jgsce.2024.205265

[22] Peng Haoyan, Xu Zhaodong, Xia Zhiheng,et al.Closed wellbore integrity failure induced by casing corrosion based on solid-chemical coupling model in CO₂ sequestration[J].Geoenergy Science and Engineering,2024,241:213140. doi:10.1016/j.geoen.2024.213140

[23] Yin Fei, Gao Deli, Zhao Jingfang,et al.Study on prediction of reservoir compaction and assessment of directional well integrity[J].Chinese Journal of Rock Mechanics and Engineering,2015,34(Supp2):4171‒4177.

[24] 尹飞,高德利,赵景芳,等.储层压实预测与定向井筒完整性评价研究[J].岩石力学与工程学报,2015,34():4171‒4177.

[25] Lei Hongwu, Li Jiaqi, Xu Tianfu,et al.Numercal simulation of coupled thermal‒hydrodynamic‒mechanical(THM) processes for CO₂ geological sequestration in deep saline aquifers at Ordos Basin,China[J].Journal of Jilin University (Earth Science Edition),2015,45(2):552‒563.

[26] 雷宏武,李佳琦,许天福,等.鄂尔多斯盆地深部咸水层二氧化碳地质储存热‒水动力‒力学(THM)耦合过程数值模拟[J].吉林大学学报(地球科学版),2015,45(2):552‒563.

[27] Xu Haiqiang.Study on rock mechanics parameters and In-Situ stress in Chang 7 member of Fan xue area in Ordos Basin[D].Beijing:China University of Geosciences (Beijing),2018.

[28] 徐海强.鄂尔多斯盆地樊学地区长7段储层岩石力学参数与地应力研究[D].北京:中国地质大学(北京),2018.

[29] Zhang Wei, Guo Tiankui, Qu Zhanqing,et al.Research of fracture initiation and propagation in HDR fracturing under thermal stress from meso-damage perspective[J].Energy,2019,178:508‒521. doi:10.1016/j.energy.2019.04.131

[30] Cao Wengui, Zhao Minghua, Liu Chengxue.Study on the model and its modifying method for rock softening and damage based on Weibull random distribution[J].Chinese Journal of Rock Mechanics and Engineering.2004,23(19):3226‒3231.

[31] 曹文贵,赵明华,刘成学.基于Weibull分布的岩石损伤软化模型及其修正方法研究[J].岩石力学与工程学报,2004,23(19):3226‒3231.

[32] Gao Dawei, Honn Kao, Wang Bei,et al.Complex 3D migration and delayed triggering of hydraulic fracturing‐induced seismicity:A case study near Fox Creek,Alberta[J].Geophysical Research Letters,2022,49(2):e2021GL093979. doi:10.1029/2021gl093979

[33] Alissar Yehya, Yang Zhuo, Rice James R.Effect of fault architecture and permeability evolution on response to fluid injection[J].Journal of Geophysical Research:Solid Earth,2018,123(11):9982‒9997. doi:10.1029/2018jb016550

[34] Kim M, Kim K Y, Han W S,et al.Density‐driven convection in a fractured porous media: Implications for geological CO₂ storage[J].Water Resources Research,2019,55(7):5852‒5870. doi:10.1029/2019wr024822

[35] Gershenzon N I, Ritzi J R W, Dominic D F,et al.Comparison of CO₂ trapping in highly heterogeneous reservoirs with Brooks‒Corey and van Genuchten type capillary pressure curves[J].Advances in Water Resources,2016,96:225‒236. doi:10.1016/j.advwatres.2016.07.022

[36] Hu Huaigang, Guan Zhichuan, Xu Yuqiang,et al.Bottom-hole stress analysis of ultra-deep wells based on theory of poroelastic mechanics[J].Journal of China University of Petroleum(Edition of Natural Sciences),2020,44(5):52‒61.

[37] 呼怀刚,管志川,许玉强,等.基于多孔弹性力学理论的深井井底应力场分析[J].中国石油大学学报(自然科学版),2020,44(5):52‒61.

[38] Shafabakhsh P, Ataie‒Ashtiani B, Simmons C T,et al.Convective-reactive transport of dissolved CO₂ in fractured-geological formations[J].International Journal of Greenhouse Gas Control,2021,109:103365. doi:10.1016/j.ijggc.2021.103365

[39] Qafoku N P, Lawter A R, Bacon D H,et al.Review of the impacts of leaking CO₂ gas and brine on groundwater quality[J].Earth‒Science Reviews,2017,169:69‒84. doi:10.1016/j.earscirev.2017.04.010

[40] Shao Changyue, Pan Pengzhi, Zhao Decai,et al.Effect of pumping rate on hydraulic fracturing breakdown pressure and pressurization rate[J].Rock and Soil Mechanics,2020,41(7):2411‒2421.

[41] 邵长跃,潘鹏志,赵德才,等.流量对水力压裂破裂压力和增压率的影响研究[J].岩土力学,2020,41(7):2411‒2421.

[42] Chen Liqiang, Tian Shouceng, Li Gensheng, et al.Initiation pressure models for supercritical CO₂ fracturing and sensitivity analysis[J].Rock and Soil Mechanics,2015,36(Supp2):125-131.

[43] 陈立强,田守嶒,李根生,等.超临界CO₂压裂起裂压力模型与参数敏感性研究[J].岩土力学,2015,36():125‒131.

[44] Gao Deli, Dou Haoyu, Dong Xuelin.Research progress in wellbore cement sheath integrity under conditions of CO₂ injection and storage[J].Journal of Yan'an University(Natural Science Edition),2022,41(3):1‒9.

[45] 高德利,窦浩宇,董雪林.二氧化碳注入条件下井筒水泥环完整性若干研究进展[J].延安大学学报(自然科学版),2022,41(3):1‒9.

基金资助

NSFC区域联合基金项目(U23A20671)

湖北省自然科学基金创新群体项目(2021CFA030)

岩土力学与工程国家重点实验室开放基金项目(SKLGME023018)