致密油多级压裂水平井流-固全耦合产能数值模拟

张东旭，张烈辉，唐慧莹，赵玉龙

（西南石油大学油气藏地质及开发工程国家重点实验室，成都 610500）

0 引言

致密油资源丰富、开发潜力大，在中国油气资源中占据十分重要的地位[1-4]。与常规油气储集层相比，致密油储集层表现为多尺度性、强非均质性，储集层有效应力状态改变时，具有极强的流-固耦合效应[5-6]。传统渗流理论一般将多孔介质骨架视作刚性模型，但致密油储集层多孔介质变形与孔隙流体流动的耦合作用不可忽略，建立在传统渗流理论基础上的油藏数值模拟技术难以准确预测致密油储集层开发动态[7-8]。

Terzaghi有效应力原理和 Biot三维固结理论共同奠定了流-固耦合研究的基础[9]，许多学者在此基础上建立了多相饱和流体渗流和多孔介质变形耦合作用模型[10-12]。致密油储集层作为一种典型的双重介质，包含基质孔隙系统和天然裂缝系统，且二者通过窜流的方式进行流体交换，具有孔喉狭小、渗透率低的特点[5,13]。通常采用水平井水力压裂的方式改造致密油储集层，除原有小尺度的天然裂缝外，还出现了大尺度的导流裂缝[14-16]。随着致密油的采出，岩石骨架有效应力重新分布，骨架发生变形，储集层物性发生动态变化，这些变化又会影响流体在储集空间的流动和压力分布[17]。因此，致密油开发过程属于多尺度、多物理过程，且存在强非线性的力学问题，如何建立考虑多尺度流动、多物理场耦合的模型并开展压裂井生产动态预测，成为致密油研究的难点。本文提出一种有效的数学模型来模拟致密油储集层基质孔隙、天然裂缝及人工裂缝中的流-固耦合，在一定程度上完善了致密油储集层多级压裂水平井多尺度多场耦合渗流理论。

油藏数值模拟技术是预测压裂水平井生产动态及其增产效果的有效方法，目前采用的数值模型主要包括连续介质模型、离散裂缝模型和连续介质-离散裂缝混合模型，其中连续介质-离散裂缝混合模型具有双重介质和离散裂缝模型的优点，在裂缝储集层模拟中具有广阔的应用前景[18-19]。流-固耦合问题的系统方程本身具有不可独立求解和无法显式消去两类独立变量（压力、位移）的特征，根据求解方法的不同，可分为间接耦合方法（显式迭代耦合、隐式迭代耦合）和直接求解的全耦合方法。间接耦合方法的方程积分交错进行，耦合参数在两场间交叉迭代。Liu等[20]基于嵌入式离散裂缝模型和有限元方法建立的适用于致密油储集层多级压裂水平井流-固耦合的数值模拟，采用迭代耦合求解渗流力学和地质力学控制方程；Zhang等[8]基于离散裂缝模型和有限元方法建立的适用于裂缝性油藏注水开发过程的流-固耦合数值模拟，采用显式耦合求解。此类方法运算速度快，但精度较低，且对时间步的要求较高，一般只能以较短的时步进行求解，同时可能存在不收敛的情况。全耦合法可实现渗流场和应力场的同步求解，不存在时间滞后现象，因其计算结果稳定、误差较小，计算结果与实际物理过程更为一致，逐渐被国内外学者采用[21-22]。流-固耦合问题常用的数值模拟方法包括有限差分方法、有限体积法、有限元方法等，其中基于非结构网格的有限元方法因其能灵活处理任意边界条件和复杂的几何形状，并且相对于其他数值方法在多物理场耦合分析方面具有优势，被广泛应用于致密储集层生产动态模拟[5,8,23-24]。

以往的研究大多集中于应力对储集层渗流能力的影响，但对于致密油储集层压裂水平井开采的动态耦合及产能研究相对较少。因此，本文针对致密油储集层开采的多尺度流体流动机理、多物理场耦合作用机理和多级压裂水平井渗流理论等关键问题，建立了表征双重介质致密油储集层流体流动和储集层变形的数学模型，构建了致密油储集层物性演化模型；采用有限元方法全耦合求解模型，分析多重介质变形规律以及储集层流体流动的复杂机制，并对不同参数条件下的产能进行了研究。

1 致密油储集层流-固耦合模型

本文应用多孔介质弹性理论与流-固耦合作用机理，描述致密油储集层人工裂缝、天然裂缝、基质系统的变形，以及流体与固体之间的相互耦合作用。致密油储集层多尺度空间如图1所示，模型的基本假设条件如下：①储集层由基质系统与天然裂缝系统组成，二者之间存在拟稳态窜流；②多孔介质视为完全饱和且各向同性的线弹性体；③储集层为弹性变形，满足小变形假设，遵从Terzaghi有效应力原理；④单相流体等温渗流，其流动遵循达西定律，忽略重力影响；⑤原油从基质系统经天然裂缝流向压裂裂缝和井筒。

图1 致密油储集层多尺度空间示意图（根据文献[16]修改）

1.1 储集层岩体变形模型

1.2 流体渗流模型

1.3 流-固动态耦合模型

流-固耦合的过程包括：流体流动引起的压力变化会影响作用在岩石骨架上的有效应力，导致岩体发生形变；岩体变形使岩石孔隙空间、孔喉半径、裂缝开度等发生改变，从而改变岩体的渗透率、孔隙度和裂缝导流能力，同时又会影响流体的流动（见图2）。

图2 致密油储集层多物理场耦合关系示意图

1.4 初始和边界条件

2 模型求解与验证

2.1 离散与求解

（6）式、（12）式和（13）式构成了描述岩体变形和流体流动的控制方程，其中每个方程均是复杂的偏微分方程，而且多物理场耦合问题又是高度非线性化的，只能求解其数值解。本文采用有限元法对控制方程进行离散，将控制方程中位移和压力表示为结点变量的插值函数形式：

基于有限元方法与离散裂缝模型，考虑人工裂缝内渗流以及汇-源项影响，系统力学平衡控制方程（6）式和系统内质量守恒方程（12）式、（13）式可在空间上分别离散成如下矩阵形式[32]：

流-固全耦合有限元平衡方程在时域进行离散，控制方程的 3个主变量（储集层节点位移、基质系统压力、天然裂缝系统压力）均以隐式方式在同一时间步内联立解出；与主变量相关的次要变量及其耦合关系在每个时间步进行迭代更新，即可求得系统位移和压力的值。流-固耦合具体求解流程如图3所示：①模型初始化，包括流动初始化和地质力学初始化，得到t=0时刻空间不同位置流压场和岩体应力状态；②流-固全耦合，将求解过程划分成若干个时间步增量，采用时间步推进的方法联立耦合求解压力场和位移场，并检查其收敛性，若不收敛，返回重复迭代；③利用第②步得到的位移、压力（主要变量）计算该时间步下孔隙度、渗透率、裂缝导流能力，更新参数（次要变量）；④重复②—③步骤，直至模拟时间结束。

图3 流-固全耦合求解流程图

2.2 模型验证

为验证上述流-固全耦合模型的正确性，建立了致密油储集层多级压裂水平井模型，模型基本参数如表1所示，数值模拟参数取值主要参考文献[5]、[8]、[33]。模型采用三角形网格进行网格剖分，其中基质系统与天然裂缝系统共用一套网格，而压裂裂缝在网格划分过程中进行降维处理；由于裂缝周围压力梯度较大，所以在压裂裂缝周边进行网格加密，由此得到整个模型储集层区域网格划分图（见图4）。由数值计算程序得到全耦合刚性模型的模拟结果，并与商业软件Eclipse采用有限差分模型得到的日产油量和累计产油量进行对比（见图5），本文全耦合刚性模型数值模拟结果与Eclipse的模拟结果吻合性较好。同时，在二维均质、线性弹性、不渗透的岩石中，均匀受压裂缝沿裂缝长度方向的开度分布满足解析解，则将本文模型和 Jiang等[34]模型计算的裂缝开度与解析解进行对比（见图6）。可以看出，3个模型计算的裂缝开度相近，而本文模型与解析解更接近。从以上两个方面验证了本文流-固耦合模型的合理性与正确性。

表1 模型基本参数表

图4 储集层区域网格剖分示意图

图5 全耦合刚性模型与Eclipse模型模拟结果对比

图6 不同模型计算的裂缝开度分布与解析解对比

2.3 模型对比

本文对比了全耦合模型（模型1）、基质+天然裂缝刚性模型（模型2）、人工裂缝刚性模型（模型3）、全刚性模型（模型 4）模拟计算的多级压裂水平井产量（见图7），并以全耦合模型累计产油量为基数，计算其他模型累计产油量与基数的相对误差，结果见表2。可以看出，全耦合与全刚性模型产量预测差异大；开发前期，全耦合模型计算日产油量递减快，生产3 000 d后全耦合模型计算累计产量最小，全刚性模型计算累计产量最高。忽略基质与天然裂缝系统的流-固耦合作用，压裂水平井生产3 000 d产量预测误差为0.92%，主要原因在于实际开采过程中基质与天然裂缝系统物性逐渐变差，其向水平井筒的供给能力逐步减弱，而模型2不考虑该因素，因此模型2计算累计产量较模型1偏大。忽略人工裂缝系统的流-固耦合作用，压裂水平井生产3 000 d产量预测误差为36.46%，主要由于实际开采过程中人工裂缝导流能力强，裂缝内压力逐渐下降、开度不断减小，而模型 3不考虑该因素，因此模型3计算累计产量较模型1更大。对比模型2—模型4产量预测相对误差值，可以看出模型 4产量预测相对误差（38.30%）并非模型2和模型3产量预测相对误差简单的线性叠加（37.38%），进一步说明基质、天然裂缝、人工裂缝之间的影响具有相互耦合特性。

图7 不同模型计算的多级压裂水平井产量

表2 不同模型计算的累计产油量对比

综上所述，流-固耦合作用在致密油产能预测中的影响不可忽略，未耦合模型预测的多级压裂水平井产能偏高，流-固全耦合模型预测的产能更为精确；人工裂缝物性是决定压裂水平井产能的重要因素，有必要优化人工裂缝参数设计。

3 储集层物性演化规律

采用模型验证部分所用的数值模型，计算致密油储集层多级压裂水平井投入生产后不同时刻基质系统的渗透率分布和不同时刻天然裂缝系统的渗透率分布（见图8—图11）。流体经基质系统和天然裂缝系统流入人工裂缝及井筒时，由于人工裂缝内导流能力大、流动阻力小，近压裂裂缝区域存在较大的流体压力梯度，储集层渗透率发生显著变化。水平井生产2 000 d后，基质系统渗透率由 0.120×10-3μm2最小减至约0.117×10-3μm2，天然裂缝系统渗透率由 1.500×10-3μm2最小减至约1.070×10-3μm2，由于基质系统更为致密，其渗透率变化幅度小于天然裂缝系统。储集层渗透率的最小值分布于人工裂缝周围，随着时间的延续，渗透率降低的范围从人工裂缝附近逐渐向储集层四周扩大。

图8 生产100 d后基质系统渗透率分布

图9 生产2 000 d后基质系统渗透率分布

图10 生产100 d后天然裂缝渗透率分布

图11 生产2 000 d后天然裂缝渗透率分布

致密油储集层多级压裂水平井生产2 000 d后，人工裂缝开度由2.884×10-4m最小减至1.440×10-4m左右（见图12、图13）。图14为人工裂缝开度损失幅度和导流能力损失幅度曲线。开发前期，由于人工裂缝内流体迅速采出，其附近压力降低较快，导致人工裂缝开度与导流能力损失幅度大；开发500 d后，人工裂缝开度损失幅度为 40.37%，导流能力损失幅度为78.80%，随后人工裂缝物性损失幅度缓慢上升；3 000 d后，人工裂缝开度损失幅度达52.12%，导流能力损失幅度达89.02%。

图12 生产100 d后压裂裂缝开度分布

图13 生产2 000 d后压裂裂缝开度分布

图14 不同时刻人工裂缝开度和导流能力损失幅度曲线

4 压裂参数优化

4.1 参数敏感性

设计案例 1—案例 9，人工裂缝压缩系数分别为0.060，0.062，0.064 MPa-1，人工裂缝开度分别为2.884×10-4，2.978×10-4，3.065×10-4m，天然裂缝压缩系数分别为0.15，0.25，0.35 MPa-1，启动压力梯度分别为0，0.05，0.10 MPa/m，模拟计算9个案例多级压裂水平井产量，并以案例 1的累计产油量为基数，计算其他 8个案例累计产油量与基数的相对误差，结果见表3。图15—图18对比了不同参数下生产3 000 d的累计产油量，其中案例 5累计产油量最大，案例 9累计产油量最小，可见水平井产量对启动压力梯度最敏感，人工裂缝开度次之。启动压力梯度减缓了储集层压力波的传播，对致密油井的产量有较大的抑制作用，提高人工裂缝初始导流能力将有助于降低产量损失，同时水平井累计产量随人工裂缝和天然裂缝压缩系数的增加而逐渐降低。

表3 不同敏感性参数对产量的影响

图15 不同人工裂缝压缩系数下的累计产油量

图16 不同人工裂缝开度下的累计产油量

图17 不同天然裂缝压缩系数下的累计产油量

图18 不同启动压力梯度下的累计产油量

4.2 裂缝参数优化

人工裂缝几何分布是影响致密油储集层产能的重要因素之一，本文在保持压裂裂缝数量与压裂裂缝长度乘积不变的情况下，设计了11个案例（案例10—案例20），分别计算了不同压裂裂缝参数下的多级压裂水平井产量，并以案例10的累计产油量为基数，计算其他案例累计产油量与基数的相对误差，案例设计与模拟结果见表4，图19对比了不同案例生产3 000 d的累计产油量。可以看出，案例20累计产油量最大，案例13累计产油量最小，随着人工裂缝数量的增多，水平井产量不断增加；各案例累计产油量曲线在生产后期接近平行，因此水平井人工裂缝加密对产量的影响主要集中于生产早期。图20展示了不同人工裂缝数量与产量预测相对误差的关系。分别取产量预测相对误差曲线前5个点和后6个点，做拟合曲线，形成交会图，裂缝条数为8～9条时，产量预测相对误差曲线出现了拐点；当裂缝数量超过10条之后，水平井产量增幅减弱，裂缝条数对产量的影响程度也逐渐减小。因此压裂施工优化设计需要考虑人工裂缝间距、数量、长度的综合影响，只增加裂缝条数无法取得预期的增产效果；在本文模拟条件下，最佳裂缝条数为9条。

表4 不同裂缝参数对产量的影响

图19 不同案例多级压裂水平井累计产油量对比

5 结论

建立了致密油储集层多级压裂水平井多重孔隙介质变形与流体流动的全耦合数值模型，模型综合考虑了储集层基质系统、天然裂缝和人工裂缝的变形特征与分区流体渗流规律，能够准确模拟预测致密油储集层多级压裂水平井产能。

模型计算结果表明，致密油储集层水平井生产过程中近人工裂缝区域储集层物性变差，特别是人工裂缝开度和裂缝导流能力损失幅度分别达到 52.12%和89.02%；流-固耦合作用在致密油储集层产能预测中的影响不可忽略，模拟致密油储集层水平井生产3 000 d后，全耦合模型与未耦合模型产能预测的误差达38.30%；参数敏感性分析表明，致密油储集层水平井产能对启动压力梯度最敏感，人工裂缝开度次之，提高人工裂缝初始导流能力有助于提高致密油井产能；致密油储集层压裂施工设计需考虑人工裂缝导流能力、间距、数量、长度的综合影响，片面追求增加裂缝条数无法取得预期的增产效果。

符号注释：

b——比奥校正系数；B——应变矩阵；c——压缩系数，Pa-1；C——柔度矩阵；Cklrt,Ⅰ——基质系统柔度张量，Pa-1；Cklrt,Ⅱ——天然裂缝系统柔度张量，Pa-1；D——弹性矩阵；Dijkl,Ⅰ,Ⅱ——双重介质系统弹性张量，Pa；fi——重力项矢量，N/m3；——边界处已知应力，Pa；G——启动压力梯度，Pa/m；h——人工裂缝面高度，m；H——流体通量，kg/（m2·s）；k——迭代步；K——渗透率，m2；K0——初始渗透率，m2；m——单位张量；Mb——比奥模量，Pa；Ms,I——基质系统岩石骨架体积模量，Pa；n——边界法线方向；n——边界单位法向量；Np——压力场的形函数；Nu——位移场的形函数；p——孔隙压力，Pa；p0——初始孔隙压力，Pa；p——节点压力，Pa；Q——源汇项，kg/（m3·s）；Q——流动方程；sF——人工裂缝存储系数，Pa-1；t——时间，s；ts——时间步；u0——系统初始位移，m；u——位移矢量，m；ui,j——ui对j方向求偏导数；u——节点位移，m；us——边界处已知位移，m；V——流体速度，m/s；wF——人工裂缝宽度，m；wF,z——z单元压裂裂缝宽度，m；wF,0——初始人工裂缝宽度，m；W——单位体积流体质量，kg/m3；x，y——储集层坐标位置，m；z——裂缝单元；α——窜流系数，m-2；∇——哈密顿算子；μ——流体黏度，Pa·s；υ——岩石泊松比；δ——Kronecker符号；ε——应变张量；εij——应变分量；εv——岩石体应变；ρ——流体密度，kg/m3；σβn——作用于裂缝表面的法向应力，Pa；σ——总应力张量；σ0——系统初始应力，Pa；σij——应力分量，Pa；σij,j——σij对j方向求偏导数，Pa/m；σi′j——有效应力张量的分量，Pa；Γm——基质系统边界；Γf——天然裂缝系统边界；ΓF——人工裂缝系统边界；Ωe——整个储集层系统；ΩF——人工裂缝系统；Ωs——储集层边界；φ——孔隙度，%；φ0——初始孔隙度，%。下标：F——人工裂缝；i，j，k，l，r，s——Voigt 标记方法中的指标；η——基质系统或天然裂缝系统；τ——窜流；Ⅰ——基质系统；Ⅱ——天然裂缝系统。