0 引言

近年来，(超)深层海相碳酸盐岩油气藏已成为四川盆地油气勘探开发的重要领域，如川中高石梯—磨溪地区龙王庙组、震旦系灯影组储层以及川西下二叠统栖霞组、茅口组等，这类油气藏储层埋藏深，具有高温、非均质性强等特征^[1]。

川西栖霞组埋藏深度超过7 200 m，储集空间以晶间孔、晶间溶孔和中小溶洞为主，天然裂缝较发育，属于低孔中低渗储层。栖二段下部储层厚度薄，邻层为水层，纵向上为白云岩和灰岩交替分布，层间应力差较小，在改造过程中易出现裂缝高度不受控、压开水层的风险。因此，该区块的改造难点为提高酸压裂缝导流能力的同时控制裂缝高度。

在早期的压裂/酸压设计数值计算时，一般将油气层的有效厚度视为裂缝高度。所以早期的PKN和CGD裂缝计算模型均将裂缝高度设置为常数^[2]。在早期的现场施工中，石油工程师们也普遍认为，水力裂缝不足以将盖层和底层压开，但在实际的增产作业中，常出现增产作业后产能不升反降的情况^[3]，学者们开始意识到裂缝高度影响产能的问题。在国内，李宾元^[4]基于椭圆裂缝模型，考虑层间物性参数差异对裂缝高度公式进行了修正，并以四川某井为例进行计算，结果显示该井裂缝高度大于储层厚度，解释了造成该井压裂减产的原因是压穿了薄盖层。此后专家学者们对裂缝高度问题展开了一系列的研究，胡阳明等^[5]、李年银等^[6]、李勇明等^[7]认为影响裂缝高度的因素包括地层参数、压裂液性能以及施工参数，并通过数学公式推导进一步分析了各参数对裂缝高度的影响规律。此外，裂缝模型也在进一步发展^[8]，根据文献[9]报到，Clifton等提出了代表性的三维裂缝模型，各石油公司也发展了全三维或拟三维的压裂、酸压模拟软件。

川西栖霞组储层情况复杂，现有的商业软件模拟结果与施工监测结果相差较大。使用有限元数值模拟方法，建立裂缝高度模拟的有限元模型，探索工程地质因素对裂缝高度发育的影响规律。在ABAQUS有限元模拟软件已有的数学模型基础上，通过用户子程序接口，自行编写孔渗参数耦合作用模块，更好地模拟非均质储层酸压过程中孔渗参数的耦合变化，并基于软件模拟结果，优化目标井X气井的施工参数，通过现场施工参数拟合、微地震测井监测结果^[10]和井温测井监测结果，确定实际裂缝高度，验证模型的合理性，以期更好地推广应用到实际生产中。

1 数学模型

基于有限元数值模拟方法，以ABAQUS有限元模拟软件流-固耦合离散模型作为基础，进行软件二次开发，并通过子程序接口^[11]编写用户自定义程序，考虑孔渗参数在模拟过程中的非均质演化，实现孔渗参数的耦合作用。

(1) 孔隙介质颗粒控制方程

流-固耦合模型采用有效应力原理，其表达式为

$ \bar \sigma = \sigma + m{p_{\rm{w}}} $

(1)

式中：$\bar \sigma $为有效应力矩阵，Pa；σ为总应力矩阵，Pa；m为[1，1，1，0，0，0]^T；p_w为孔隙压力，Pa。p_w符号表示为“压正拉负”，因流体仅能承受压力，故可认为p_w恒为正值，而$\bar \sigma $和σ的符号与之相反。

视川西下二叠统碳酸盐岩为孔隙介质，由虚功原理^[12]可推导出其岩体的固体平衡方程

$ \int_{\rm{V}} {\left( {\bar \sigma - m{p_{\rm{w}}}} \right)} \delta \varepsilon {\rm{d}}V = \int_{\rm{S}} t \delta v{\rm{d}}S + \int_{\rm{V}} f \delta v{\rm{d}}V $

(2)

式中：δ ε为单元应变率，s^-1；t为单元面力矢量；f为不考虑流体重力的单位体积力，N。

(2) 孔隙介质流体控制方程

使用达西定律描述流体渗流，推导连续性方程；基于质量守恒定律，由“入=出”原则，d t 时间内流入某一体积岩石内的液量等于其内部液量的增加，推导得连续性方程^[13]为

$ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( { \int\limits_{\Omega} \varphi {\rm{d}}\Omega } \right) = - \int\limits_{\rm{S}} {} \varphi n{v_{\rm{w}}}{\rm{d}}S $

(3)

式中：φ为岩石孔隙度，%；n为与积分外表面法线平行的方向矢量；v_w岩石孔隙间流体流动速度，m/s。

最终推导得到应力-渗流耦合方程表达式^[15]为

$ \left[\begin{array}{cc} K & C \\ E & G \end{array}\right] \frac{\mathrm{d}}{\mathrm{d} t}\left\{\begin{array}{l} \bar{u} \\ \bar{p}_{0} \end{array}\right\}+\left[\begin{array}{cc} 0 & 0 \\ 0 & F \end{array}\right]\left\{\begin{array}{c} \bar{u} \\ \bar{p} \end{array}\right\}=\left[\begin{array}{c} \frac{\mathrm{d} f}{\mathrm{d} t} \\ \hat{f} \end{array}\right] $

(4)

(3) 二次开发数学模型

在模拟计算过程中，设岩石不可压缩，则体积变化等于孔隙变化(体积)，定义$\varepsilon_{\mathrm{V}}=\varepsilon_{\mathrm{X}}+\varepsilon_{\mathrm{Y}}+\varepsilon_{\mathrm{Z}} $；同时

$ {\varepsilon _{\rm{v}}} = \frac{{\Delta V}}{{{V_0}}} = \frac{{\Delta {V_{\rm{p}}}}}{{{V_0}}} = \frac{{{V_{\rm{s}}}\Delta \varphi }}{{{V_{\rm{s}}}(1 + \Delta \varphi )}} = \frac{{\varphi - {\varphi _0}}}{{1 + {\varphi _0}}} $

(5)

式中：ε_V为岩石体积应变，m³；V₀为岩石初始体积，m³；φ₀为岩石初始孔隙度，%；V_s为固体体积，m³；V_p为孔隙体积，m³。

并由渗透率和孔隙度关系式，推导渗透率变化关系式

$ k=k_{0}\left[\left(\frac{1}{\varphi_{0}}\right)\left(1+\varepsilon_{\mathrm{V}}\right)^{3}-\left(\frac{1-\varphi_{0}}{\varphi_{0}}\right)\left(1+\varepsilon_{\mathrm{V}}\right)^{-\frac{1}{3}}\right]^{3} $

(6)

式中：k为渗透率，D；k₀为初始渗透率，D。

2 物理模型建立

选取川西X气井栖霞组储层作为建模依据，根据图 1中X气井的测井解释曲线，设计物理模型储隔层分布。X气井储层累计厚度为41.2 m，上层为11.0 m，下层为30.2 m，中间隔层为15.2 m。对储、隔层厚度数据进行归整处理，便于模型建立和网格划分，设计“五层式”模型，从上至下依次为上隔层、1#储层、中间隔层、2#储层和下隔层。1#储层厚度为10 m，中间隔层厚度为15 m，2#储层厚度为30 m。在1#储层上部和2#储层下部均设置一定厚度，防止裂缝压开整个模型，模型示意图如图 2所示。参考X气井栖霞组岩石物性参数，对参数进行处理，地层参数设置如表 1所列。因模型中无法设置水层含水，故在模拟计算时假设某层为水层，若裂缝压开该层，则视为压开水层。

预制裂缝采用Cohesive单元，预制人工裂缝沿最大水平主应力方向，设计变量为工作液黏度、注液排量和用液规模，探究不同施工参数下计算所得的缝高结果，并使用对节点进行参数提取，计算最终裂缝高度值，探究其变化规律。对模型进行离散化网格处理，采用Structured Hex网格，即结构正六面体网格，可以提高软件计算速度，精细化计算结果。模型网格化如图 3所示。

3 结果分析 3.1 工作液黏度

分析工作液黏度、工作液排量和工作液规模对缝高影响的单因素规律；设置排量为5 m³/min，工作液规模为250 m³，黏度分别为1 mPa·s，25 mPa·s，50 mPa·s和75 mPa·s，黏度变化范围涵盖目前现场施工大部分工作液黏度。计算结果如图 4所示。

提取裂缝边缘的节点x，y坐标，绘制裂缝边缘坐标图，可以表征裂缝的高度和长度。设置物理模型中间点为y坐标0点，裂缝起点为x坐标0点；各工作液黏度裂缝的高度和长度如图 5所示。

由图 5可知，单因素情况下，工作液黏度和裂缝高度呈正相关关系；工作液黏度高于25 mPa·s，裂缝会压开下隔层，按照建立物理模型时的假设，视为压开水层；裂缝长度与裂缝高度呈负相关关系；因为高黏工作液在缝内摩阻更高，滤失更少，酸压裂缝更易“憋压”，导致裂缝高度较大；在工作液规模相同的情况下，工作液更多作用在“缝高方向”，更少贡献给“缝长方向”。

3.2 施工排量

设置工作液黏度为25 mPa·s，规模为250 m³，排量为3~6 m³/min，梯度为0.5 m³/min，计算结果如图 6所示。

本次模拟所设计排量范围涵盖川西栖霞组气井增产施工设计的排量，满足现场酸压设计的要求，适用于X气井酸压改造。提取裂缝坐标点，绘制坐标图(图 7)。

由图 7可知，当排量大于5.0 m³/min时，裂缝延伸到下隔层，根据假设，视为压开水层；在X气井酸压施工时设计排量低于5.0 m³/min，可以满足控制裂缝高度的要求。

3.3 液体规模

设置排量为4.5 m³/min，工作液黏度为25 mPa· s，注液规模为100~350 m³，变化梯度为50 m³，绘制裂缝高度随工作液规模的变化图(图 8)。

由图 8可知，随工作液规模的扩大，裂缝高度呈增大趋势。当工作液规模大于250 m³后，中间隔层被完全压开。当工作液规模大于300 m³后下隔层被压开，视为出水。故在仅考虑单因素的情况下工作液规模尽量不超过300 m³，有利于在储层改造的同时不压开下隔层。

3.4 层间应力差

川西栖霞组层间应力差不大，这导致了施工过程中控制裂缝高度比较困难。根据李年银等^[15]研究，层间应力差越大，裂缝高度越小。故在酸压施工参数设计时，可考虑对X气井使用人工隔层，增大层间应力差，以达到更好地控制缝高的效果。

设置层间应力差为1~9 MPa，变化梯度为2 MPa，将计算出的裂缝高度值、单位厚度注入排量、注入工作液黏度进行插值绘图，可得到三者的变化规律(图 9)。图中排量为单位厚度注入排量，因为模型分为1#储层和2#储层，为了便于研究，将排量改为单位厚度注入排量，对应排量为3~6 m³/min。

图 9表示模型中所有储隔层均被压开。由图 9可知，层间应力差越大，裂缝高度越小。若X气井目标层位层间应力差增大，则在设计施工参数时可设计更大施工排量和更高黏度工作液(图 10)。根据X气井及其所在区块地层情况，推荐工作液黏度≤ 40 mPa·s；单位厚度排量≤ 0.12 m²/min，换算为排量≤ 4.8 m²/min；推荐用液规模不超过250 m³；使用这些参数进行模拟，结果如图 10所示。

由图 10可知，隔层并未被压开，上层裂缝仅在1#储层内发育，下层裂缝也仅在2#储层内发育；满足控制裂缝高度的施工要求；模拟结果中，1#储层内裂缝长度为30.35 m，2#储层内裂缝长度为91.32 m。

4 施工曲线拟合

X气井经参数优化后已完成酸压施工，工作液为胶凝酸，施工时间共计263 min，施工过程中最大排量为4.8 m³/min，一般排量为3.8~4.8 m³/min；最高油压为112.9 MPa，一般油压为104~112 MPa；最高环空压力为81.3 MPa，一般环空压力为35~ 47 MPa；X气井施工曲线如图 11所示。通过FracPT压裂模拟软件对酸压数据进行拟合，通过施工曲线拟合得到酸压裂缝几何参数，判断本次施工是否满足控制裂缝高度的酸压施工要求，证明本模型是否合理指导X气井施工。

X气井施工排量最高达4.8 m³/min，工作液黏度为17~18 mPa·s，用液规模为250 m³，按照本次模拟结果推荐的施工参数进行施工。再使用施工参数用FracPT软件进行净压力拟合，其结果如图 12所示。

由图 12可知，在X气井的施工参数下，经拟合形成2个酸压主缝；X气井成像测井显示，X气井井周鲜有天然裂缝分布，结合施工曲线2处油压降，可判断本次施工形成了2条主缝。拟合结果显示，1#储层裂缝高度为13.2 m，长度为41.7 m；2#储层裂缝高度为32.5 m，长度为98.0 m；FracPT拟合结果与本文模型计算结果高度吻合，本文模型计算裂缝几何形态与FracPT拟合裂缝形态基本一致。由X气井现场施工参数拟合结果显示，裂缝总高度为50.3 m。X气井生产中井温显示裂缝高度结果与现场施工参数拟合结果相当。本模型计算结果、净压力拟合结果与现场监测结果基本吻合。X气井最终测试产量为10.45万m³/d，本井施工效果好。

5 结论

(1) ABAQUS有限元模拟软件功能完善，内置有限元离散方程，便于二次开发，可较好模拟X气井储层情况，能较好模拟X气井裂缝高度的延伸，可用于川西下二叠统碳酸盐岩油气藏油气井控缝高酸压设计。

(2) 本模型模拟分析表明，影响酸压裂缝高度的主要因素为工作液黏度、施工排量、液体规模和层间应力差。

(3) 针对X气井情况，主要目标为控制裂缝不压开下隔层，在施工时加入上浮式和下沉式隔离剂，控制2#储层裂缝高度。本模型也应该针对此种工况作进一步修正。