本文的研究对象是上述经典二维模型中的PKN模型。该模型最初由Perkins和Kern提出［9］。这两位学者忽略了裂缝扩展时裂缝缝宽变化所产生的的存储效应、以及压裂液通过裂缝表面进入储集层的滤失效应。他们提出的模型被称为PK模型。Nordgren［10］从能量消耗角度阐释了上述两种效应所代表的物理意义。存储效应指的是压裂液注入时，裂缝消耗压裂液所具有的能量形成新的裂缝表面，从而进一步扩展发育。也就是说，存储效应代表的是裂缝发育的有效能耗。而滤失效应指的是裂缝内部的压裂液通过裂缝表面进入储集层，带走了这部分压裂液所具有的能量，使得裂缝无法有效使用这部分能量形成新的裂缝表面，即裂缝无法使用这部分能量进行扩展发育。因此，滤失效应代表的是裂缝发育的无效能耗。这两种能耗模式在裂缝扩展过程中相互影响。通过数学推导和物理分析，Nordgren［10］提出了无量纲时间tD，对这两种能耗模式进行了量化。他通过计算发现，当tD<0.01时，存储效应占主导，这时可以得到模型的一组半解析解；当tD>1时，滤失效应占主导，可以得到模型的另一组半解析解；而当0.01≤ tD ≤ 1.00时，裂缝发育处于过渡阶段，无法获得模型的半解析解。Nordgren［10］提出了模型半解析解成立的参数判定条件，完善了模型。因此，业界以这3位学者姓氏首字母给模型命名，即PKN模型。Kemp［11］通过渐进式分析，改进了无滤失条件下的PKN模型裂缝尖端计算。Kovalyshen等［8］在存储效应占主导阶段、滤失效应占主导阶段、以及过渡阶段这3种参数条件下，进行了更普适的裂缝尖端渐进式分析，提出了更精确的半解析解成立参数判定条件，进一步改进了PKN模型。
采用本课题组和华裔学者徐冠水课题组联合研发的全三维水力压裂软件FrackOptima，在PKN 模型半解析解成立的参数条件下，建立相关算例进行数值模拟计算。通过对比数值模拟结果和半解析解计算结果，验证了该软件的精确性。为使用该软件进行实际复杂工况条件下的压裂优化设计提供了理论依据。
1 PKN模型理论
PKN型水力裂缝如图1所示。该裂缝几何与物理特征总结如下［8， 10］：
（1） 裂缝所处岩石介质是均匀、各向同性的无限大弹性介质，其中岩石韧性KIC忽略不计；
（2） 裂缝缝长L远大于缝高H，且H保持为常数。一般L/H大于3。因此，在每个和x轴垂直的平面上，可以假设平面应变条件成立；
（3） 牛顿流体以常数流量Q0/2进入裂缝，以层流方式仅沿x轴方向运动。流体重力忽略不计，同时不考虑裂缝尖端的流体滞后现象。
<G:\武汉工程大学\2023\第1期\龙恭博-1.tif>[缝宽w ][L(缝长)][最大缝宽
w0（x=0） ][H(缝高为常数)][y][x][z][L>>H>>w][压裂液流动方向]
图1 PKN型水力裂缝示意图（单位：m）
Fig. 1 Schematic view of PKN hydraulic fracture （unit： m）
该模型的控制方程是［8， 10］：
[E’128μH?2?x2(w40)=8CLπt-t0+?w0?t] （1）
式中E′是杨氏模量E和泊松比v的表达式：E′=E/（1-v2）；μ是流体黏度；t是压裂持续时间；w0是裂缝横截面的最大缝宽，它是时间t和坐标x的函数w0 （x， t）；CL是卡特滤失系数；t0是压裂液到达裂缝某点处所用时间，它是该点坐标x的函数t0 （x）。本文所有物理量的单位均为国际单位制单位。
其对应的初始条件是w0 （x， 0）=0；边界条件是［8， 10］：
w0（x，t）=0， x≥L（t） （2）
[?w40?xx=0=-512μQ0πE’] （3）
上述方程是一个高度非线性的偏微分方程。不借助基于计算机的数值方法，很难对其精确求解。但是在特定参数条件下，可以对其进行渐进式半解析求解。引入无量纲时间tD，其表达式为［8， 10］：
[tD=64E’C5LHt1.5π3Q2023] （4）
当tD<1.657 8×10-5时，裂缝处于存储效应占主导阶段，上述方程可得渐进式半解析解［8， 10］：
[L(t)=0.68×16-15E’Q30μH415t45w0(0, t)=2.5×2-15μQ20E’H15t15pw(t)=2.5×64-15E’4Q02μH615t15] （5）
当tD>206.31时，裂缝处于滤失效应占主导阶段，上述方程可得渐进式半解析解［8， 10］：
[L(t)=Q02πCLHt12w0(0, t)=4μQ20π3E’CLH14t18pw(t)=2E’3μQ20π3CLH514t18] （6）
方程（5）、（6）中的pw是裂缝入口（x=0）处的井底净压力。
而当1.657 8×10-5≤ tD ≤ 206.310 0时，裂缝处于过渡阶段［8］。在过渡阶段参数条件下，无法对控制方程进行较精确的解析求解，必须借助数值方法求解。
2 PKN模型的数值模拟计算
2.1 FrackOptima压裂软件简介
采用本课题组和华裔学者徐冠水课题组联合研发的全三维水力压裂软件FrackOptima，在PKN 模型半解析解成立的参数条件下，建立相关算例进行数值模拟计算。该软件通过基于全三维位错理论的“综合性边界元”方法进行离散处理，因此多个非平面裂缝之间的相互作用可以通过严格的弹性力学理论进行计算。裂缝壁面之间的流体运动通过经典润滑理论进行计算，压裂液的滤失则通过卡特模型进行建模。固流耦合、支撑剂的运移和沉降均通过时间步进算法进行计算。所以，该软件可以用来模拟计算地应力不均匀、岩石属性不均匀条件下非平面水力裂缝的扩展。该软件已在壳牌、中国石油勘探开发研究院等国内外大型石油公司得到了实际应用［12-13］。该软件的具体技术细节、相关方程见文献［14］。
2.2 算例参数与网格设置
PKN模型理论忽略流体重力和岩石韧性。而FrackOptima是面向实际复杂工况条件下的全三维压裂软件，因此将算例中的流体密度和岩石韧性设为较小数值，近似PKN模型的理论条件。设流体密度ρ=10 kg/m3、岩石韧性KIC=1.1 MPa?m0.5。同时设置裂缝所在区域的最小地应力比岩石其余区域最小地应力小20 MPa，从而限制裂缝的垂直缝高，满足缝高为常数的PKN模型理论要求。
根据实际工况［15-16］，设置裂缝扩展区域高度H=20 m、岩石杨氏模量E=30 GPa、泊松比v=0.2、压裂液黏度μ=0.05 Pa?s。根据PKN模型理论中的参数条件，如表1所示，设置下列参数分别满足存储效应占主导阶段、滤失效应占主导阶段的条件。
表1 2种阶段下的算例参数取值
Tab. 1 Parameters of two stages
[取值项 CL / （10-5m/s0.5） Q0 / （m3/s） t / s tD 存储主导 0 0.1 900 0 滤失主导 3.54 0.12 8 100 375.5 ]
FrackOptima压裂软件通过边界元方法进行离散处理。根据边界元方法的特点，只需对裂缝平面进行离散处理和网格划分。另一方面，根据水力压裂计算力学相关理论［14］，缝长和缝高的维度上至少需要设置5个元素，才能获得较为可靠的数值计算结果。同时，压裂软件的计算时间也是实际工程应用中需要考虑的重要因素之一，所以网格不宜划分过细。因此，本文在裂缝平面上选取4 m（垂直方向）×8 m（水平方向）和2 m（垂直方向）×4 m（水平方向）的两种矩形网格分别进行数值计算。
3 计算结果及分析
3.1 存储效应占主导阶段
在存储效应占主导阶段，FrackOptima数值模拟（4 m×8 m网格）得到的PKN型水力裂缝如图2所示。该图表明，在FrackOpitma的数值计算中，裂缝缝高保持常数，裂缝平面形状和图1中PKN型水力裂缝平面的理论形状一致。
FrackOptima数值计算结果（PKN模型裂缝缝长L、入口处最大缝宽w0、入口处流体井底净压力pw）和半解析解对比如图3所示。该图表明，FrackOptima的数值模拟计算结果和半解析解计算结果吻合程度较好。其中，4 m×8 m网格和更精细的2 m×4 m的数值计算结果的差别小于5%，但是二者的计算耗时分别是约20 min和约78 min。前者耗时约是后者的1/4。由于压裂软件的计算时间是实际工程应用中需要考虑的重要因素之一，因此，4 m×8 m可以作为存储效应占主导阶段条件下FrackOptima进行PKN模型数值计算的网格尺寸。
3.2 滤失效应占主导阶段
在滤失效应占主导阶段，FrackOptima数值模拟（4 m×8 m网格）得到的PKN型水力裂缝如图4所示。该图中的裂缝缝高保持常数，裂缝平面形状也和图1中PKN型水力裂缝平面的理论形状一致。
<G:\武汉工程大学\2023\第1期\龙恭博-4.tif>[裂缝宽度 /
mm][7
1]
图4 滤失效应占主导阶段条件下FrackOptima数值模拟得到的PKN型裂缝图
Fig. 4 PKN hydraulic fracture obtained by numerical
simulation of FrackOptima in leak-off-dominated stage
FrackOptima数值计算结果（PKN模型裂缝缝长L、最大缝宽w0、井底净压力pw）和半解析解对比如图5所示。该图表明，在滤失效应占主导阶段条件下，FrackOptima的数值模拟计算结果和半解析解计算结果吻合程度也较好。两种网格的数值计算结果的差别小于3%，但是二者的计算耗时分别是约88 min和约251 min。前者耗时约是后者的1/3。所以，4 m×8 m也可以作为滤失效应占主导阶段条件下FrackOptima进行PKN模型数值计算的网格尺寸。在上述两个算例中，本文采取两种网格尺寸进行计算，所得结果相差较小，完成了网格无关性检验。
3.3 结果分析
图3和图5表明，不管是存储效应占主导阶段，还是滤失效应占主导阶段，FrackOptima数值计算结果和半解析解结果都有较小的差别，而且数值计算结果比半解析解偏小。这是因为：
（1） PKN模型的建模基础包含了流体无重力、一维流动、岩石无韧性等不符合物理实际的理论假设。而FrackOptima是面向实际工况条件的压裂软件。在计算时对流体密度、岩石韧性取较小值来近似满足该条件，同时并未控制流体流动。因此，数值计算中流体的多维流动、岩石韧性对于裂缝扩展的限制都会使得数值计算中的流体相对于理论模型中的流体而言，多损耗一部分能量，使得数值结果比半解析解偏小。
（2） PKN模型假设裂缝所处岩石介质是完全均匀的，且裂缝高度保持常数不变。但是根据水力压裂力学基本理论，裂缝在完全均匀介质中会在裂缝高度方向上不断扩展［5］。为了满足缝高保持常数不变的理论假设，数值计算采取了设置最小地应力屏障的参数设定。这一设定能够保证缝高保持不变，但是也会使得一部分能量在裂缝对抗应力屏障的过程中损耗，导致数值计算结果比半解析解偏小。
4 结 论
通过对比基于FrackOptima的PKN型水力裂缝数值计算结果和相关半解析解计算结果，得出结论如下：
（1） 数值计算得到的PKN型水力裂缝结果和半解析解计算结果吻合度较好。只要在缝高维度上保证5个元素，就可以既获得较精确的数值计算结果，又能保证较合适的计算耗时。
（2） 数值计算所采取的相关参数设置措施能够近似满足PKN模型的建模理论条件。尽管这些措施会导致数值计算结果比半解析解偏小，但误差不超过10%，满足实际工程应用的需求，为进一步采用数值模拟优化水力压裂实际方案提供了理论基础。