币型模型最初由Geertsma等［7］在1969年提出。由于当时计算机技术尚不成熟，无法对该模型的控制方程进行数值分析，所以仅给出了近似解。随着计算机技术的不断发展，美国工程院院士Detournay等［1， 3， 5］结合数值计算和解析求解对币型模型进行分析，从能量消耗角度阐释了币型裂缝发育的物理机理。币型裂缝发育受到“存储效应和滤失效应”、“压裂液黏度耗散和储集层岩石韧性抗裂”效应的影响。其中，存储效应指压裂液通过高压泵进入裂缝时，其携带的能量转换为裂缝扩展所需的新表面能，促使裂缝发育。因此，可以将存储效应理解为裂缝扩展所需的有效能量消耗。而滤失效应指压裂液进入裂缝后未参与裂缝发育，而是以渗透的方式经裂缝内表面进入储集层，导致裂缝不能利用该部分压裂液的能量进一步发育。所以，可以将滤失效应理解为裂缝扩展过程中的无效能量损失。另一方面，压裂液黏度耗散指的是当压裂液黏度较大而岩石韧性较小时，黏稠的压裂液在裂缝内部不易流动，岩石较易破碎起裂。这时，裂缝发育所需的有效能耗主要用来克服压裂液粘度，使得黏稠的压裂液在裂缝内部能够顺利流动，从而将岩石破碎，促进裂缝发育。而韧性抗裂指的是当压裂液黏度较小而岩石韧性较大时，压裂液在裂缝内部流动较容易，但岩石不易破碎起裂。这时裂缝发育所需的有效能耗主要用来克服岩石韧性，直接将岩石破碎，促进裂缝发育。
Detournay等［1， 2， 5］对上述4种物理机制进行量化分析，发现在4组参数条件下，可以得到币型模型的半解析解。它们分别是：（1）存储和黏度耗散占主导时的半解析解，记为M解；（2）滤失和黏度耗散占主导时的半解析解，记为[M]解；（3）存储和韧性抗裂占主导时的半解析解，记为K解；（4）滤失和韧性抗裂占主导时的半解析解，记为[K]解。
然而，韧性抗裂占主导会导致无量纲韧性κ取值变大，很难保证币型模型数值求解的精确性和稳定性［8］。因此，使用本课题组和外籍学者徐冠水共同开发的FrackOptima真三维水力压裂模拟软件，在上述4组特定参数条件下，开展币型模型的数值计算研究，重点研究κ取值较大时的计算结果。比较数值计算结果和半解析解运算结果，论证该软件的稳定性和精确性，为进一步使用FrackOptima优化实际压裂方案，尤其是优化低黏度滑溜水压裂方案，奠定了理论基础。
1 币型模型基本理论
币型裂缝模型如图1所示。由图1可知，币型裂缝关于油井轴对称。
<G:\武汉工程大学\2024\第1期\刘长昊-1.tif>[裂缝尖端][流动方向][w][p][弹性介质][油井井口][r][Q0]
图1 币型裂缝示意图
Fig. 1 Schematic diagram of penny-shaped fracture
该模型的几何特征和物理特性总结如下［1-3］：
（1） 岩石介质是无限大均匀且各向同性的弹性介质。其杨氏模量、泊松比、韧性分别记为E、v、KIC；
（2） 压裂液为不计重力的牛顿流体，仅沿r方向做层流运动。其井口注入流量为常数Q0，黏度为常数μ。压裂液流动前端和裂缝尖端重合。
该模型的流体流动控制方程是［1-2］：
[?w?t?+CLt-t0r=1μr??rrw3?p?r] （1）
固体线弹性断裂力学控制方程是［1-2］：
[w=8πER0RGrR,xRpx,txdx] （2）
记x=r/R、x′=x′/R，其中积分核心G（r/R，x′/R）可表示为G（x，x′），其表达式为：
[Gξ,ξ=1ξFarcsin1-ξ21-ξ2,ξ2ξ2,???????ξ>ξ1ξFarcsin1-ξ21-ξ2,ξ2ξ2,?????????ξ>ξ] （3）
式中E′=E/（1-v2）；μ′=12μ；x′是虚拟变量；t是压裂持续时间；R是裂缝半径，它是关于时间t的函数；p是裂缝内的流体压力，w是裂缝缝宽，它们都是关于坐标r和时间t的函数；CL是卡特模型滤失系数；t0是压裂液流至裂缝内某点所用时间，它是该点坐标r的函数t0（r）；F是第一类不完全椭圆积分。
上述控制方程的初始条件满足w（r，0）=0，同时边界条件满足［1-2］：
w（r，t）=0， r≥R （4）
?p/?x =0， r=R （5）
式（1）和式（2）都是无法通过传统数学分析方法进行解析求解的复杂方程。然而，当相关参数取值在若干特定范围内时，能够使用渐进式近似求解方法对其半解析求解［1-3］。引入无量纲滤失常数?与无量纲时间常数τ ［2］：
?=[μ3E11C4Q0K14]， τ = [K?9μ2.5E6.5Q1.50t] （6）
式中C′=2CL，K′=4（π/2）1/2KIC。
在给定压裂算例相关参数取值后，可以通过式（6）计算lgτ和lg?的取值。如果以（lgτ， lg?）为坐标的点落在图2所示的M、[M]、K、[K] 4个半解析解所对应的区域中的任意1个之内，则表明在该算例参数条件下，存在该区间所对应的币型模型半解析解。因此，可以在图2所示的参数条件图中确定M解、[M]解、K解、[K]解4组半解析解分别适用的参数条件［2］。图2中的叉号和圆点坐标分别代表了本文中满足[M]解和[K]解的2个算例参数所对应的lgτ和lg?的值。
<G:\武汉工程大学\2024\第1期\刘长昊-2.tif>[-10 0 10 20
lgτ][20
10
0
-10
-20
-30][lg?]
图2 币型模型半解析解参数条件图
Fig. 2 Parametric condition of semi-analytical solutions to penny-shaped fracture model
这4组半解析解的具体解析表达式如下：
M解 ［1-3］：
[R(t)=0.694 4EQ30μ19t49w0(t)=1.190 1μ2Q30E219t19pw(t)=2.410 9E2μ13Φt-13] （7）
[M]解［2， 9-10］：
[R(t)=0.450 2Q0C12t14w0(t)=1.057 4μ2Q30E2C18t116pw(t)=3.093 1μ2E6C3Q3018Φt-316] （8）
式中，Φ是包含了第一类和第二类完全椭圆积分及其微分的复杂数学表达式，具体形式请见文献［2］。其值需要通过数值计算来确定。本文采用MATLAB 2021a进行计算，可得Φ≈0.501 4。而w0是井口处（r=0）的裂缝最大缝宽，pw是井口处（r=0）处的井底压力。
另2组半解析解表达式如下：
K解 ［1-3］：
[R(t)=0.854 6EQ0K25t25w0(t)=0.653 7K4Q0E415t15pw(t)=0.300 4K6EQ015t-15] （9）
[K]解 ［2， 9-10］：
[R(t)=0.450 2Q0C12t14w0(t)=0.474 4K4Q0E4C14t18pw(t)=0.413 9K4CQ014t-18] （10）
此外，当不考虑滤失（CL=0）时，可以通过裂缝尖端无量纲化渐进式分析［3， 11］，取得更精确的无量纲半径γ的解［3］：
[γ=0.697 6,????????????????????????????????????????????????????????????????????????????????? κ≤13π225κ-25-54475π2κ-4,?????????κ≥3.5] （11）
式中，κ是无量纲韧性，其表达式为：
[κ=Kt2μ5Q30E13118] （12）
式（11）中的κ≤1和κ≥3.5分别对应不考虑滤失时M解和K解存在所需满足的参数条件。由式（12）可知，当黏度取值较小时，κ取值较大。在设计实际压裂方案时，如果选取低黏度的滑溜水作为压裂液，无量纲韧性κ的取值一般大于4［8］。
另一方面，若通过式（6）计算得到的lgτ和lg?取值使得以（lgτ， lg?）为坐标的点未进入图2中的4个区域，则表明币型裂缝处于发育过渡阶段［9］。在其对应的参数条件下，不存在半解析解。
2 币型模型的数值模拟计算
采用真三维水力压裂模拟软件FrackOptima，在图2和式（6）所示4组半解析解存在的参数取值范围内建立币型模型算例，开展数值仿真研究，重点关注κ≥4时的模拟研究。FrackOptima已在中国石油勘探开发研究院、康菲、壳牌等众多科研机构和油气公司得到了油田现场应用［8， 12-13］。相关技术细节请详见文献［12］。
在算例参数取值方面，币型模型理论忽略压裂液重力，即将压裂液密度ρ视为0。但Frack-Optima是针对现实压裂方案的工程软件，所以在算例中取压裂液密度为较小数值，即ρ=10 kg/m3，从而近似满足“压裂液密度为0”这一理论条件。
基于压裂现场工况［14-16］，设置岩石泊松比v=0.25。根据图2和式（6）所示的参数取值范围，设置分别满足M解、[M]解、K解、[K]解这4组半解析解存在所需的4组参数。其取值如表1所示。
表1 4组算例参数取值
Tab. 1 Parameters of four numerical examples
[参数 E /
GPa μ /
(Pa·s) KIC /
(MPa·m0.5) CL /
(m/s0.5) Q0 /
(m3/s) t / 103s M解 12 0.100 0.2 0 0.025 1.0 M?解 10 0.100 0.2 0.002 0.200 108.0 K解 10 0.001 8.0 0 0.060 1.8 K?解 20 0.001 8.0 0.002 0.200 108.0 ]
对于表1中的[M]解和[K]解的参数值，其对应的lgτ和lg?值分别如图2中的叉号和圆点坐标所示，所以它们分别满足[M]解和[K]解参数条件。而表1中M解参数取值对应的lgτ =-8.246，?=0，所以lg?→-∞。以（-8.246， -∞）为坐标的点落在图2中的M解所对应的区域内。类似地，表1中K解参数取值对应的lgτ =11.372，?=0，所以lg?→-∞。以（11.372， -∞）为坐标的点落在图2中的K解所对应的区域内。
另一方面，针对式（11）所示的无量纲半径解，如表2所示，设置11组参数条件，满足无量纲韧性κ的取值范围，进而开展数值计算。
根据水力压裂力学理论［1］和前期研究经验［13-14］，在任意维度上布置至少10个元素，就有望得到较为可靠的仿真结果。根据式（7）~式（10）中的半径计算式，表1~表2所选参数生成的裂缝半径都在60 m左右。因此，本文选取5 m（水平方向）×5 m（竖直方向）矩形网格进行数值计算，保证径向上有约12个元素，确保得到较可靠的计算结果。
表2 无量纲算例参数取值
Tab. 2 Parameters of dimensionless numerical examples
[参数 E /
GPa μ /
(Pa·s) KIC /
(MPa·m0.5) Q0 / (m3/s) t / s κ #1 9.446 0.250 0.2 0.14 4 876 0.1 #2 10.990 0.250 0.5 0.20 2 994 0.2 #3 9.107 0.100 0.8 0.16 3 541 0.5 #4 5.694 0.080 0.9 0.20 3 186 0.8 #5 5.284 0.050 1.0 0.25 2 441 1.0 #6 4.377 0.004 1.6 0.24 1 872 3.5 #7 7.269 0.002 2.0 0.15 1 972 4.0 #8 10.990 0.001 3.0 0.18 980 5.0 #9 10.850 0.001 5.0 0.22 844 8.0 #10 12.250 0.001 6.0 0.13 1 215 10.0 #11 11.340 0.001 8.0 0.25 759 12.0 ]
3 计算结果及分析
3.1 4组半解析解参数条件
图3（a-c）、图4（a-c）、图5（a-c）和图6（a-c）分别比较了M解、[M]解、K解、[K]解这4组半解析解存在的参数条件下，软件数值计算得到的币型裂缝半径R、井口最大缝宽w0、井底压力pw及其对应的半解析解运算结果。而图3（d）、图4（d）、图5（d）和图6（d）则展示了数值计算得到的币型水力裂缝形状。它们的裂缝平面形状和图1所示的币型裂缝平面理论形状相同，均为圆币型。
另一方面，在图3（a-c）、图4（a-c）、图5（a-c）和图6（a-c）中，软件数值计算结果与半解析解运算结果较好地吻合。所以，5 m×5 m可以作为币型模型数值模拟的网格尺寸。
在这4组参数条件下，数值模拟计算结果和半解析解相差小于10%，满足实际工程的误差要求。图3（c）中，存储效应和黏度耗散占主导（M解）条件下，2种方法计算得到的井底压力在开始阶段相差较大，但是相互收敛得很快。这是因为币型裂缝的井底压力计算式是关于时间t的奇异方程。对比式（7）~式（10）中的井底压力计算式可知，M解条件下t的奇异性最强，因此在压裂初始阶段，半解析解得到的压力和数值解相差较大。这一现象也在图4（c）、图5（c）和图6（c）中有所体现。此外，根据相关研究［1-2］，由于奇异性的存在，井底压力半解析解在压裂初始阶段（即t较小时）的精度并不具备较好的实际工程对比价值。
图3（a）、图4（a）、图5（a-c）和图6（a-c）的数值计算结果曲线不是光滑曲线，而是带有阶梯型变化，造成该现象的原因有2点。首先，在滤失效应占主导时（图4和图6），大部分压裂液渗入岩石，需要较多时间积累足够的压裂液来促进裂缝发育。这一现象在数据的体现上就是时间跨度较大的阶梯型半径增长。其次，根据水力压裂计算力学基本理论，计算得出时间步长Δt和流体黏度μ成正比［4］。所以，在韧性抗裂占主导时（图5和图6），流体黏度μ很小，这就要求Δt也必须很小，从而保证数值计算的稳定性。例如在图6对应的[K]解参数条件下，Δt量级维持在10-3 s。这使得计算难以保持长时间的连续稳定性，导致在若干连续时间步长内会得到一两个和其他数据相差较大的瞬时数据。因此，FrackOptima软件在图形用户界面采用“图像平滑”算法来修正这些瞬时数据，这也会导致数据的阶梯型变化。
根据式（12）可得，表1中K解相关参数所对应的无量纲韧性κ =18.53。图5表明，在该条件下FrackOptima软件能够获得较精确和稳定的计算结果，其数值计算结果和半解析解误差绝对值不超过5%。而在其他3组参数条件下，数值结果和半解析解的误差绝对值也不超过10%。另一方面，文献［8］表明，大多数压裂软件在κ≥3，即黏度μ较小时，就失去计算稳定性，无法获得较精确的结果。
3.2 无量纲半解析解参数条件
表2中的11组算例数值计算得出的数值解的无量纲半径[γn]和方程（11）的无量纲半径γ的对比如图7所示。其中，数值解的无量纲半径γn计算式是：
[γn=μQ30Et419Rn] （13）
式中，Rn是软件数值求解得到的裂缝半径。
由图7可知，数值解和无滤失条件下的无量纲半解析解误差绝对值小于5%，吻合度较理想。
4 结 论
使用FrackOptima压裂软件开展币型水力裂缝模型数值计算研究，通过比较数值模拟结果和半解析解运算结果可以得到如下结论：
（1） 数值模拟结果和半解析解计算结果误差绝对值小于10%，和无滤失条件下的无量纲半解析解误差绝对值小于5%，吻合度较好。
（2） 当压裂液黏度很小，无量纲韧性κ=18.53时，数值计算仍能获得较稳定、精确的结果。为进一步采用FrackOptima软件数值模拟优化实际压裂方案，尤其是优化滑溜水（黏度很小的压裂液）压裂方案，奠定了理论基础。