阳春沟区块位于四川盆地南川页岩气田西部。综合评价认为，阳春沟五峰组−龙马溪组页岩气资源潜力大，已成为四川盆地下一步勘探开发的重点区域。该区块埋深大，应力高，导致现场施工压力较高，排量难以达到要求而无法完成加砂目标，亟需采用新技术提高储层改造效果[1-6]。
超临界CO2具有低黏度、低表面张力、高渗透的物理特性，在储层改造应用中表现出了巨大的潜力[7-9]。近年来，学者们针对超临界CO2压裂开展了大量研究工作。2015年，陈立强等[10]考虑超临界CO2黏度、可压缩性及增压速率等因素，建立了超临界CO2压裂起裂压力预测模型，但未考虑温度场和压力场对CO2性质的影响。Zhang和Liu等[11-12]基于有限元算法建立了超临界CO2压裂数值模型，研究了水、N2和超临界CO2对裂缝起裂的影响。结果表明，超临界CO2压裂具有较低的起裂压力。2019年，苏建政等[13]引入Pen-Robison方程，研究了超临界CO2和滑溜水压裂裂缝形态差异。结果表明，超临界CO2良好的扩散性及渗透能力可以增加孔隙压力，减少地应力对裂缝扩展的约束。因此，压后的裂缝形态更复杂。Mollaali等[14]应用相场模型，通过S-W方程和修正的达西定律，考虑CO2的密度变化以及可压缩性，模拟超临界CO2压裂裂缝起裂，但模型假设所有过程均为等温条件。2021年，Feng等[15]采用相场模型对比分析了水、CO2和N2压裂过程中裂缝扩展特征，证明CO2压裂可以降低破裂压力。
学者们基于数值模拟开展了大量关于超临界CO2压裂的研究，但多数模型未考虑CO2与地层之间的热交换及温压对CO2流体性质的影响，无法真实模拟超临界CO2压裂过程。因此，基于相场法建立了超临界CO2压裂热流固多场耦合模型，探究了相关参数对起裂压力和裂缝形态的影响，并在阳春沟区块开展现场应用，以期为超临界CO2压裂技术的现场应用提供借鉴作用。
根据Zhou等[16]的研究结果，多孔介质岩石总势能公式见式（1）。
式中：$ \psi $为多孔介质总势能，J；u为位移场，m；$ {\varOmega } $为计算域；Γ为裂缝域；$ \psi _{\text{k}}^{} $为动能密度，J/m3；$ \psi _\varepsilon ^{} $为弹性应变能密度，J/m3；$ \alpha $为Biot系数；$ p $为流体压力，Pa；$ {G_{\mathrm{c}}} $为临界能量释放率，J/m2；b为体积力，N；f为牵引力，N。
计算能量泛函L的一阶变分，设其值为0，即可得到式（2）。
式中：$ {\sigma _{ij}} $为应力张量分量，Pa；$ {\rho _{}} $为密度，kg/m3；$ {l_0} $为长度尺度参数，m；$ \phi $为标量相场；k为模型稳定性参数；$ \psi _\varepsilon ^ + $为拉伸弹性应变能密度，J/m3。
为保证裂缝扩展的不可逆特性，在此引入应变历史场H后，得到式（3）。
将式（2）中的ψε+用H(x,t)替换，相场方程可改写为式（4）。
χR和χF为两个线性函数。在基质区域：χR=1，χF=0；在裂缝区域：χR=0，χF=1。过渡区域表达式见式（5）~式（6）。
式中：c1和c2为过渡区域边界范围，当ϕ≤c1时，为裂缝区域；当ϕ≥c2时，为基质区域；当c1 < ϕ < c2时，为过渡区域。
假设流体在多孔介质内流动符合达西定律，流体流动控制方程见式（7）。
式中：密度$ \rho = {\rho _{\mathrm{R}}}{\chi _{\mathrm{R}}} + {\rho _{\mathrm{F}}}{\chi _{\mathrm{F}}} $，ρR和ρF分别为岩石和流体密度，kg/m3；S为储水系数，1/Pa；t为时间，s；渗透率$ K = {K_{\mathrm{R}}}{\chi _{\mathrm{R}}} + {K_{\mathrm{F}}}{\chi _{\mathrm{F}}} $，KR和KF分别为岩石和流体渗透率，m2；μ为流体黏度，Pa·s；g为重力加速度，m/s2；qm为质量源项，kg/(m3·s)；Biot系数$ \alpha = {\alpha _{\mathrm{R}}}{\chi _{\mathrm{R}}} + {\chi _{\mathrm{F}}} $，αR实际为岩石基质Biot系数，而$ \alpha $为模型考虑多孔介质（多孔介质包括流体+岩石基质）流体和基质共同作用下的Biot系数；εvol$=\nabla $·u，为ΩR（t）的体积应变，无因次。
CO2注入过程处于高压状态，应用Helmholtz自由能理论和S-W方程计算不同压力和温度条件下的CO2物性参数（见图1）。Helmholtz自由能表达式见式（8）[17]。
式中：Φ为Helmholtz自由能，无量纲；ΦO为气体部分Helmholtz自由能，无量纲；Φr为残余流体部分Helmholtz自由能，无量纲；δ为折算密度，无量纲；τ为折算温度，无量纲。
裂缝内超临界CO2与岩石发生热交换会产生热应力，影响裂缝的起裂和扩展。为了研究热应力对裂缝起裂和扩展的影响，用式（9）计算温度变化而产生的热应力。
式中：${\sigma _{ij}} $为热应力，Pa；λL为Lamé 第一常数，Pa；δij为克罗内克符号，无因次；G为Lamé 第二常数，Pa；αT为线性热膨胀系数，1/K；Kvol为体积模量，Pa；T为压裂后的地层温度，K；T0为地层初始温度，K。
多孔介质内流体与岩石的热量传递方程式见式（10）[18]。
式中：cp,eff为等效体积比定压热容，J/(m3·K)；λeff为等效导热系数，W/(m·K)；Q为热量源汇项，W/m3。
cp,eff、λeff可分别表示为式（11）~式（12）。
式中：cp,R、cp,F分别为岩石和流体比定压热容，J/(m3·K)；λR、λF分别为岩石导热系数和流体导热系数，W/(m·K)。
基于COMSOL建立的压裂模型如图2所示。前期地质分析表明，阳春沟区块天然裂缝发育，因此在模型中加入天然裂缝，天然裂缝长与宽分别设定为0.2 m与0.012 m，天然裂缝与初始压裂裂缝中心距离为0.12 m，为消除模型边界效应对裂缝扩展的影响[19-20]，模型尺寸的长、宽设置为1.5 m×1.5 m，初始流体压力为0，其余参数见表1[15]。
模型涉及5个耦合求解模块，分别为固体力学模块、相场模块、历史应变模块、达西流动模块和传热模块。通过对Helmholtz方程进行修正，建立相场模型，并基于ODE和DAE方程，确定历史应变模型。最后基于相场模型和历史应变模型对ϕ和H进行求解。
σh设定为1 MPa，水平应力差（p）分别设定为1 MPa、2 MPa、3 MPa和4 MPa。模拟结果如图3～图5所示。
从图3可知，水平应力差相同时，水力压裂裂缝偏转角度要小于超临界CO2压裂。裂缝偏转角度的大小体现了水平最大主应力对裂缝扩展控制作用的强弱，偏转角度越大，代表水平应力差对裂缝控制作用越强。结果表明，在水平应力差相同时，超临界CO2压裂裂缝受水平应力差影响更大。
从图4可知，在不同水平应力差下，与水力压裂相比，超临界CO2压裂时裂缝起裂压力较小，起裂时间更短。随着水平应力差增大，超临界CO2压裂与水力压裂的起裂压力与起裂时间均随之减小，但超临界CO2压裂的降幅较大。这表明，在实际施工过程中，对于高水平应力差地层，超临界CO2增能压裂的破岩效果会更明显。
从图5可看出：在p=1 MPa和2 MPa时（见图5（a）、图5（c）），超临界CO2压裂裂缝沿天然裂缝尖端形成4条裂缝；在p=3 MPa和4 MPa时（见图5（e）、图5（g）），超临界CO2压裂裂缝沿天然裂缝扩展分别形成3条和2条裂缝。这是因为超临界CO2的超低黏度和高渗透的物理特性，更易渗入并开启天然裂缝，即使在高水平应力差下也可形成复杂裂缝。根据Li等[21]的研究，裂缝网格数（相场ϕ>0.95）与模型总网格数的比值为裂缝面积占比，能够反映裂缝的复杂程度。计算发现，当p从1 MPa增加至4 MPa，超临界CO2压裂裂缝复杂度从0.034降低至0.031，水力压裂裂缝复杂度从0.032降低至0.028。结果表明，水平应力差增大，超临界CO2压裂裂缝和水力压裂裂缝复杂程度均有所减弱，但在相同水平应力差情况下，超临界CO2压裂裂缝复杂程度仍要高于水力压裂裂缝。
为研究天然裂缝逼近角（θ）对裂缝起裂的影响，设定σh为1 MPa，σH为4 MPa，θ分别为45°、60°、75°和90°，其他参数与2.1节相同时，实验模拟结果如图6所示。
从图6可看出：当θ分别为45°、60°和75°时，天然裂缝主要在左下和右上两个端点沿最大水平主应力开启扩展形成两条裂缝；当θ为90°时，裂缝在天然裂缝尖端扩展形成3条裂缝；随着θ的增加，裂缝复杂程度逐渐增加。
图7为不同天然裂缝逼近角下裂缝起裂时间与起裂压力的变化。从图7可看出，随着天然裂缝θ值增大，起裂时间和起裂压力均随之增大。这是因为天然裂缝θ值越小，天然裂缝面应力分量受最大主应力影响越少，水力裂缝越容易沿着天然裂缝扩展延伸，因此起裂压力较小；当θ值越大时，作用于天然裂缝面上的最大水平主应力就越大，导致起裂压力较大。
σH与σh分别设定为3 MPa和1 MPa，改变CO2注入温度，以20 K差值从278.15 K增加至338.15 K，其余参数与2.1节相同时，模拟结果如图8所示。
对比4种注入温度下裂缝扩展形态（见图8）可以看出，CO2与地层的温差越高，越有利于形成复杂的裂缝形态，CO2注入温度为278.15 K时，即与地层温差为87 K时，产生了地应力诱导裂缝。这是由于超临界CO2与地层之间的温差会在岩石内部产生热应力，这种热应力有助于降低岩石起裂压力，产生更多的裂缝。随着后期温差不断减小，热应力减小，裂缝的扩展模式相对简单。
图9为不同注入温度时CO2压裂裂缝起裂时间与起裂压力的变化情况。从图9可看出，CO2注入温度从338.15 K降低至278.15 K，起裂压力降低。这是因为随着注入温度降低，CO2与地层温差增加，热应力诱导作用增强，起裂压力降低，起裂时间缩短，这与卢义玉等[22]的研究结果一致。但当注入温度高于318.15 K后，CO2与地层温差减小，由热应力公式可知，温差减小，热应力减小，因此起裂压力和起裂时间降幅减小。
Y井位于阳春沟南斜坡，目的层为上奥陶统五峰组，水平段长1581 m，分为22段，水平段最小水平主应力为81.7 MPa，水平段最大水平应力为92.4 MPa，地层温度为115.9 ℃。2023年10月25日，采用常规水力压裂完成10段压裂施工。由于该井施工压力较高，从第11段开始，采用前置超临界CO2压裂施工。CO2在压力高于7.38 MPa、温度高于31.06 ℃时处于超临界状态。阳春沟Y井原始地层压力约为54.0 MPa，地层温度115.9 ℃，CO2通过套管注入到达地层时，必定处于超临界状态。施工前基于该井实际应力数据和地层温度进行模拟（见表1），模拟压力曲线及现场压裂曲线如图10所示。
由图10可知，现场施工排量8.0~17.5 m3，CO2注入排量为0.6 m3/min，最高砂比12%，胶液300 m3，酸30 m3。在注入前置超临界CO2 176.8 m3后，继续采用水基压裂液压裂，压力达到101.0 MPa时地层就破裂，相比第10段破裂压力（111.0 MPa，见图11）降低了10.0 MPa。同时，对比实际施工曲线和模拟压力曲线可知，破裂压力分别为101.0 MPa和97.0 MPa，误差为3.9%；而延伸压力分别为104.0 MPa和88.0 MPa，误差为15.3%。这主要是因为数值模拟并未模拟破裂后的加砂过程，因此地层破裂后压力并不会上升。结果表明，所采用的数值模拟方法可以在一定程度上指导现场超临界CO2增能压裂施工。
从图11可看出，在第1段~10段，地层破裂压力范围为109.8~116.2 MPa，而第11段~22段地层破裂压力为94.8~104.2 MPa，破裂压力显著下降。这是由于：①超临界CO2易扩散到水基压裂液不能到达的微孔裂隙，增大了孔隙压力；②超临界CO2与岩石发生热交换，在岩石内部产生了热应力。由于该井地层温度和杨氏模量均较高，热交换引起的热应力较大，在有效应力超过岩石破坏强度的热应力区域会产生诱导微裂缝。微裂缝的存在则会形成弱面，降低裂缝壁面岩石的破坏强度，最终表现为破裂压力降低。
在第1段~10段，停泵压力为73~81 MPa，而在第11段~22段，停泵压力为87.8~93.2 MPa。可以看出，加入前置超临界CO2不仅可以降低地层破裂压力，还可以提高后续常规水力压裂的停泵压力。停泵压力反映了裂缝缝内净压力的大小，缝内净压力的提高有利于增强压裂缝网复杂程度[23]。
1） 基于相场方法建立了热流固耦合超临界CO2压裂模型，通过S-W方程计算了不同温压下CO2流体性质的变化，基于能量变分和最小耗能原理求解相场断裂方程。模型考虑了压裂时温压对CO2流体性质的影响以及超临界CO2与岩石的热交换，可以更准确地模拟超临界CO2压裂裂缝起裂压力和裂缝形态。
2） 模拟结果表明，与常规水力压裂相比，超临界CO2压裂时裂缝起裂压力更小。在水平应力差相同时，超临界CO2可以获得更复杂的裂缝形态。天然裂缝逼近角与起裂压力和起裂时间成正相关关系。同时，天然裂缝逼近角越大，裂缝形态越复杂。CO2注入温度与地层温度差值越大，热应力诱导作用越明显，压裂裂缝形态越复杂。
3） 对模拟结果在阳春沟南斜坡Y井开展的超临界CO2增能压裂现场实践表明，采用前置超临界CO2增能压裂不仅可以降低地层破裂压力，还有利于提高后续水力压裂的停泵压力，增强储层改造效果。
2024, Vol. 53 Issue (5): 77-83 2024, Vol. 53 Issue (5): 77-83