随着油气田开发进程的不断发展，包括页岩油气在内的非常规油气藏在油田勘探开发中的地位日益重要。水力压裂技术作为开采这些非常规油气资源的主流方式之一，起到了至关重要的作用^[1]。利用数值模拟、物理实验及现场观测等手段揭示水力压裂过程，对促进油气增产具有重要意义。

水力压裂过程十分复杂，涉及缝内流体压力引起的岩体变形及裂缝扩展、缝内流体的流动及滤失^[2-3]等力学问题。基于不同假设，众多学者提出了水力压裂解析和半解析模型。Khristianovic等^[4]在假设裂缝高度方向形态完全一致的基础上，提出了KGD模型，推导出缝高远大于缝长时水力裂缝的解析解。Geertsma等^[5]考虑了裂缝边界的流体滤失，修正了KGD模型。在同一时期，Perkins等^[6]提出了PKN模型，Nordgren^[7]对PKN模型做了进一步修正。随后Simonson等^[8]建立了水力压裂的拟三维模型。Spence等^[9]在证明KGD模型中裂缝扩展具有自相似特性的基础上，提出了模型的半解析解。该半解析解不能显式地体现能量耗散的影响，同时不能考虑有滤失的情形。Bunger等^[10]系统性地研究了KGD模型的半解析解, 并利用表征能量耗散比和表征压裂液在岩体中赋存模式建立了KGD模型半解析解的完整体系。Shel等^[11]提出了轻微滤失的情况下拟三维模型的半解析解。

水力压裂常用的数值方法包括边界元法和有限元法。边界元法^[12]作为处理裂缝起裂和扩展的不连续界面问题的有效方法，利用微分算子的解析解作为边界积分方程的核函数，具有解析解与数值解相结合的特点，通常具有较高的计算精度。位移不连续法^[13-14]作为一种特殊的边界元法，其核心思想是对裂缝面进行离散，以不连续位移为基本解，在裂缝扩展时只需在新产生的裂缝面和裂缝尖端增加少许单元，不需要重新划分网格。文[15-17]以边界元方法为基础建立了改进的拟三维模型。Chen等^[18]针对边界元法求解拟三维水力压裂模型效率不高的问题，提出了一种基于Runge-Kutta-Legendre方法的显式时间步长算法。Adachi等^[19]利用其提出的拟三维模型，研究在2个对称应力边界上的水力裂缝的扩展高度。侯冰等^[20]讨论了页岩储层多级水力裂缝的扩展规律。

传统有限元法^[21]及其衍生的扩展有限元法^[22]在模拟非均质岩石中裂缝的扩展方面具有极大优势，适用于各种复杂受力状态和各种工程实际问题，目前已成为水力压裂数值计算模拟的强大工具。薛炳等^[23]针对油层岩石在水力载荷作用下的裂缝扩展及渗流行为，采用渗流应力耦合模型对水力压裂过程进行了三维有限元方法数值模拟。潘林华等^[24]在损伤力学理论的基础上，结合损伤单元模拟裂缝起裂和扩展，建立了水力压裂裂缝扩展的三维有限元模型。王涛等^[25]利用扩展有限元方法模拟了裂缝在平面内的扩展。曾青冬等^[26]针对页岩气藏水力压裂复杂缝网形成机理，建立了基于扩展有限元方法的多裂缝同步扩展数学模型。Gupta等^[27]在假设裂缝内流体压力完全一致的基础上，采用广义有限元方法建立水力压裂数值模拟的真三维模型，并讨论了井筒附近裂缝的曲面扩展。

水力压裂数值模型研究方面取得了长足的进步，从拟三维模型发展到全三维模型甚至真三维模型，从过去边界元占主导地位的情形发展到现今边界元方法和有限元方法共同主导的情形。总体而言，数值模型在物理实验验证方面仍然亟待进一步的研究，近些年这方面的主要研究成果包括Lecampion等^[28]基于边界元方法的轴对称模型与经典轴对称物理实验^[29]的对比以及Settgast等^[30]基于有限元方法的全三维模型与经典平面应变物理实验^[31]的对比。本文在建立表述水力压裂裂缝扩展、变形的固体有限元方程和表述裂缝内流体流动、滤失的流体有限元方程的基础上，通过耦合同步求解2套有限元方程，建立了全三维全耦合有限元水力压裂数值模型。为了验证数值模型的可靠性，数值模型与前述2个经典的水力压裂实验进行全面对比。数值模型完全遵从水力压裂的核心理论，对其物理过程未做任何简化，同时对其边界条件没有做任何增删，通过数值模型与物理实验的对比可以讨论水力压裂理论的适用性。

1 物理模型

如图 1所示，水力压裂是在远场最小围压σ₀和内水压力p_f的作用下裂缝在宽度、长度和高度上不断演化的过程。根据线弹性理论，在σ₀和p_f的作用下，岩石和裂缝力学响应的相关方程描述如下。

忽略重力的作用，在水力压裂准静态过程中岩石的应力平衡方程为

$ \nabla \cdot \boldsymbol{\sigma}=0. $

(1)

其中：▽∙为散度算子，σ为岩石的应力张量，它与应变张量ε的关系为

$ \boldsymbol{\sigma}=\boldsymbol{D}: \boldsymbol{\varepsilon}. $

(2)

其中：D为弹性系数张量，∶为张量的双点积。ε为位移梯度的对称部分，即

$ \boldsymbol{\varepsilon}=\left[\nabla \boldsymbol{u}+(\nabla \boldsymbol{u})^{\mathrm{T}}\right] / 2. $

(3)

其中：▽为梯度算子，u为位移。

记分析域的外表面为Γ_σ，裂缝表面为Γ_p，则应力边界条件为

$ \boldsymbol{\sigma} \cdot \boldsymbol{n}=\boldsymbol{\sigma}_{0}, \quad \text { 在 } \varGamma_{\sigma} \text { 上 } ; $

(4)

$ (\boldsymbol{\sigma} \cdot \boldsymbol{n}) \cdot \boldsymbol{n}=-p_{\mathrm{f}}, \quad \text { 在 } \varGamma_{\mathrm{p}} \text { 上. } $

(5)

其中：n为图 1中分析域表面单位外法线矢量。裂缝的宽度为裂缝两表面的位移差在其法线方向上的投影，因此其表达式为

$ w=[\boldsymbol{u}] \cdot \boldsymbol{n}. $

(6)

根据线弹性断裂力学理论，当裂缝尖端应力强度因子达到岩石断裂韧性时裂缝扩展，因此裂缝的扩展准则为

$ K_{\mathrm{I}}=K_{\mathrm{IC}}. $

(7)

其中：K_I为裂缝端部的应力强度因子，K_IC为岩石断裂韧性。

水力压裂中流体在裂缝中主要包括流动和滤失2个物理过程。忽略流体的压缩性，压裂液在地层中的质量守恒方程为^[2]

$ \nabla \cdot q+\frac{\partial w}{\partial t}+g-\delta(x, y) Q=0. $

(8)

其中：q为流量，g为流体滤失，δ为Dirac函数，Q为排量。当压裂液为Newton流体，根据边界层理论, q满足立方体定律^[32]，相应的方程为

$ q=\frac{w^{3}}{12 \mu} \nabla p_{\mathrm{f}}. $

(9)

其中：μ为动力黏度。采用Carter滤失模型^[33]，g的表达式为

$ g=\frac{2 C_{\mathrm{L}}}{\sqrt{t-t_{0}}}. $

(10)

其中：C_L为流体滤失系数，t₀为流体前缘达到时间。

2 数值方法

根据虚功原理，联立式(1)—(5)可得到任意时刻分析域Ω内可得到固体平衡方程的积分弱形式：

$ \begin{gathered} \int_{\varOmega}(\delta \boldsymbol{\varepsilon}): \boldsymbol{\sigma} \mathrm{d} \varOmega-\int_{\varGamma_{\sigma}}(\delta \boldsymbol{u})^{\mathrm{T}} \cdot \boldsymbol{\sigma}_{0} \mathrm{~d} \varGamma- \\ \int_{\varGamma_{\mathrm{p}}}(\delta \boldsymbol{u} \cdot \boldsymbol{n}) \cdot p_{\mathrm{f}} \mathrm{d} \varGamma=0. \end{gathered} $

(11)

其中：δu为任意容许虚位移，δε为其相应的虚应变。用四面体单元将分析域Ω离散为有限个连续单元。对于任一单元区域Ω_e，有

$ \boldsymbol{u}=\varSigma \boldsymbol{N}_{i} \boldsymbol{u}_{i}, \delta \boldsymbol{u}=\boldsymbol{N}_{i} \delta \boldsymbol{u}_{i} \quad(i=1,2,3,4). $

(12)

其中：N_i为单元内部节点的形函数，u_i和δu_i为相应的节点位移和虚位移。

根据式(1)—(5)和式(11)，有

$ \begin{gathered} \sum\limits_{\varOmega_{e}}\left(\delta \boldsymbol{u}_{i}\right)^{\mathrm{T}} \int_{\varOmega_{e}} \boldsymbol{N}_{i}^{\mathrm{T}} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B} \boldsymbol{N}_{j} \boldsymbol{u}_{j} \mathrm{~d} \varOmega_{e}= \\ \sum\limits_{\varOmega_{e}}\left(\delta \boldsymbol{u}_{i}\right)^{\mathrm{T}} \int_{\varGamma_{\sigma}^{e}} \boldsymbol{N}_{i}^{\mathrm{T}} \boldsymbol{\sigma}_{0} \mathrm{~d} \varGamma+ \\ \sum\limits_{\varOmega_{e}}\left(\delta \boldsymbol{u}_{i} \cdot \boldsymbol{n}\right)^{\mathrm{T}} \int_{\varGamma_{p}^{e}} \boldsymbol{N}_{i}^{\mathrm{T}} p_{\mathrm{f}} \mathrm{d} \varGamma. \end{gathered} $

(13)

其中：B为单元应变矩阵。式(13)中δu_i为任意容许虚位移，因此有

$ \begin{gathered} \sum\limits_{\varOmega_{e}} \int_{\varOmega_{e}} \boldsymbol{N}_{i}^{\mathrm{T}} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B} \boldsymbol{N}_{j} \boldsymbol{u}_{j} \mathrm{~d} \varOmega_{e}= \\ \sum\limits_{\varOmega_{e}} \int_{\varGamma_{\sigma}^{e}} \boldsymbol{N}_{i}^{\mathrm{T}} \boldsymbol{\sigma}_{0} \mathrm{~d} \varGamma+\sum\limits_{\varOmega_{e}} \int_{\varGamma_{\mathrm{p}}^{e}} \boldsymbol{N}_{i}^{\mathrm{T}} p_{\mathrm{f}} \mathrm{d} \varGamma. \end{gathered} $

(14)

式(14)可以简写成

$ \boldsymbol{K}_{s} \boldsymbol{U}=\boldsymbol{F}_{1}+\boldsymbol{F}_{2}. $

(15)

其中：K_s为Ω的整体刚度矩阵，U为节点位移向量，F₁、F₂分别为σ₀和p_f的等效节点力。式(15)即为岩石弹性响应的有限元方程。

采用隐式积分，根据式(8)—(10)，建立表述压裂液流动和滤失的有限元方程，即

$ \boldsymbol{K}_{\mathrm{f}} \boldsymbol{P} \Delta t+\boldsymbol{L}\left(\boldsymbol{W}-\boldsymbol{W}_{\mathrm{e}}\right)+\boldsymbol{G} \Delta t+\boldsymbol{Q} \Delta t=0. $

(16)

其中：K_f为流体的整体刚度矩阵，P和W为待求的流体压力和裂缝宽度矢量，L为裂缝面积矩阵，G为等效节点滤失矢量，Q为等效节点排量，W_e为上一平衡步中的裂缝宽度矢量，Δt为时间步长。

采用刚度凝聚技术^[34]，根据式(15)建立裂缝宽度和流体压力的关系方程为

$ \boldsymbol{W}=\boldsymbol{K} \boldsymbol{P}+\boldsymbol{C}. $

(17)

其中：K为反映裂缝宽度对流体压力全局依赖性的稠密矩阵，C为σ₀引起的裂缝闭合度。

求解以W和P为未知量的耦合方程式(16)和(17)即可得到新的变量。式(16)中K_f是W的非线性函数。采用Newton-Raphson法以及与文[34]相同的迭代收敛准则同步求解式(16)和(17)，从而实现数值模型的全耦合模拟。从式(16)和(17)可以看出，耦合方程的规模仅与裂缝面的节点数目有关。刚度凝聚技术有效地降低了耦合方程的规模。模型采用与文[34]相同的迭代收敛准则和应力强度因子计算方法。当裂缝前缘节点上的应力强度因子达到临界值时，裂缝在该点继续扩展并采用文[35]的方法更新有限元网格。这种方法将网格更新引起的计算开销降至最低。

根据式(16)和(17)求得裂缝面上各点流体压力后，由式(14)计算由σ₀和p_f引起的岩体位移和应力，并以此为基础采用M积分方法^[34]计算应力强度因子。M积分方法的优点是在裂缝前缘采用常规的有限单元同样可以获得良好的计算精度。因此，完全自主开发的全耦合数值模型采用常规四面体单元模拟岩体的变形，采用四面体单元在裂缝面内的三角形网格模拟裂缝的扩展，并将它们作为模拟流体流动和滤失的有限单元。当裂缝前缘节点上的应力强度因子达到临界值时，裂缝在该点继续扩展并采用文[35]中节点分裂的方法更新迭代有限元网格。这种方法将网格更新引起的计算开销降至最低。

从所采用的数值方法可以看出，水力压裂模型完全遵从水力压裂理论，完整地反映了水力压裂中4个核心的物理过程^[10]，即岩体的变形、裂缝的扩展、压裂液在裂缝中的流动和滤失。与常规的解耦方法不同^[30]，数值方法采用耦合方法而不需要假设任何流体压力边界条件，这符合水力压裂理论和实际工程中裂缝任意点流体压力待解的情形。数值方法对水力压裂理论中的边界条件未做任何增删。

3 实验验证

数值模型与2个经典物理实验进行对比以验证其可靠性。第1个物理实验为有滤失的轴对称实验^[29]，第2个物理实验为没有滤失的平面应变实验^[31]。2个物理实验包括裂缝宽度、流体压力和裂缝几何尺寸在内详实的观测数据，因此被用来验证各种数值模型^{[28, 30]}。

3.1 轴对称实验

文[29]中C12物理实验试件是边长为0.3 m的立方体，实验中压裂液通过半径为1 cm的孔眼注入到试件中心。在侧向围压的作用下裂缝沿水平面扩展。试件由混凝土和岩石的混合体组成，其弹性模量为24 GPa，Possion比为0.25，断裂韧性为0.6 MPa·m^1/2，压裂液的动力黏度为0.13 Pa·s，排量为2.952 mm³/s。模型的液滤失系数为3.0×10^-5 m·s^1/2[28]。该物理实验中裂缝近似圆形扩展，因此将此物理实验称为轴对称实验。根据对称性，数值方法采用图 2a所示的1/4试件模型模拟压裂过程。数值实验共有151 066个四面体单元和26 757个初始网格节点。如图 2b所示，裂缝面内有3 557个初始网格节点和13 856个三角形单元。裂缝面内有限元网格的特征尺寸为0.002 5 m。

图 3为轴对称物理实验和数值实验得到的裂缝半径、孔眼附近流体净压力及其裂缝宽度。从图 3a可以看出，物理实验中左半部分的半径和右半部分的半径差别始终不大。物理实验中裂缝近似圆形扩展。从图 3可以看出，数值实验与物理实验的主要差异体现在裂缝扩展初期流体净压力值上。随着裂缝的扩展，数值实验中裂缝半径、孔眼附近流体净压力及其裂缝宽度和物理实验结果具有很好的吻合度。

图 4为轴对称数值实验得到的4个不同时刻流体净压力和裂缝宽度分布图。与文[28]的轴对称模型不同，数值实验采用全三维模型模拟裂缝的扩展。尽管如此，从图 4可以看出，数值模拟结果中裂缝始终近似圆形扩展，这与物理实验具有相同的裂缝扩展模态。此外，当裂缝近似圆形扩展时，其前缘的净压力梯度逐渐减小，这与理论分析的结果一致^[36]。

3.2 平面应变实验

如图 5a所示，文[31]的无滤失平面应变物理实验是由三层透明有机玻璃组成，层间通过氯仿进行粘接。有机玻璃的弹性模量为3.28 GPa，Possion比为0.367，断裂韧性为1.2 MPa·m^1/2。压裂液的动力黏度为97.7 Pa·s，排量为73.2 mm³/s。根据对称性，数值实验采用图 5b所示的1/2试件模型模拟压裂过程，数值实验有103 976个四面体单元，19 101个初始网格节点。如图 5c所示，裂缝面内有5 160个三角形单元，1 365个初始网格节点，孔眼在裂缝面内的坐标为(0, 0)。物理实验在图 5c中(0, -0.028 5)坐标处有一裂缝宽度观测点。

图 6为平面应变实验物理实验和数值实验得到的裂缝半径、孔眼附近流体净压力和观测点裂缝宽度。虽然物理实验中左右初始裂缝半长有所差别(左边半长为22 mm，右边半长为26 mm)，但是从图 6a可以看出，实验仍然具有很好的对称性。因此数值实验采取对称分析是合理的。数值实验中初始裂缝半长为2.5 mm(详见图 7a)，流体净压力为1 MPa。初始条件不同是图 6中数值实验结果与物理实验在裂缝扩展初期有差异的一个重要原因。但是从图 6可以看出，随着裂缝的扩展，数值实验结果与物理实验逐渐趋于一致，初始条件的影响逐渐减弱。

图 7为平面应变数值实验得到的6个不同时刻流体净压力和裂缝宽度分布图。图 7a为压裂开始时裂缝的情形。从图 7b可以看出，刚开始裂缝由于远离上、下遮挡层的约束，其扩展模态为圆形；从图 7c可以看出，当裂缝前缘达到上、下遮挡层后，其扩展模态逐渐发生变化；从图 7d—7e可以看出，裂缝最终扩展模态为其前缘逐渐依次推进并且符合平面应变模型，这与物理实验观察的结果一致。

从与2个经典物理实验的对比中可以看出，数值实验结果与物理实验结果的差异主要是裂缝扩展初期的流体净压力。这种差异同样见诸于文[30]。随着裂缝的扩展，数值实验的结果与物理实验的结果渐渐趋于一致，其裂缝扩展模态也相同。数值模型得到经典实验的验证。裂缝中流体的流态受控于裂缝宽度。实际工程中裂缝宽度通常为毫米级别，而从图 3c和6c可以看出，2个物理实验中裂缝宽度为微米级别。数值模型完全遵从水力压裂理论，并采用流体流动的立方体定律(式(9))模拟实验中压裂液的流动。数值模拟结果与物理实验结果具有很好的一致性，说明裂缝宽度在微米级别立方体定律在水力压裂理论中仍然适用可行。

4 结论

在建立两套耦合有限元方程描述水力压裂核心物理过程岩体变形、裂缝扩展、流体流动和流体滤失的基础上，通过耦合同步求解2套方程，建立了基于有限元方法的水力压裂全三维全耦合数值模型。采用刚度凝聚技术使得耦合方程的规模仅仅依赖于裂缝面上的节点数目。为了验证数值模型的可靠性，与2个经典的水力压裂物理实验对比。在定量分析方面，数值实验得到的裂缝长度、流体净压力和裂缝宽度与物理实验结果具有非常好的一致性；在定性分析方面，2个数值实验与相应的物理实验具有相同的裂缝扩展模态。数值模型得到经典物理实验的验证。数值模型完全遵从水力压裂理论，对其物理过程未做任何简化，同时对其边界条件没有做任何的增删。因此数值模型可以用于讨论水力压裂理论的适用性，并且通过与物理实验的对比说明裂缝宽度在微米级别立方体定律在水力压裂理论中仍然适用可行。