堵漏材料与天然裂缝宽度的匹配对于裂缝性地层堵漏作业及裂缝性储层屏蔽暂堵作业有着重要的作用^[1-3]，采用成像测井、声波测井和核磁共振等方法虽然能识别井下天然裂缝宽度^[4-5]，但在井漏时进行测井作业存在风险。若能利用实钻井漏录井数据对天然裂缝宽度进行预测，不仅可以快速为裂缝性地层堵漏和裂缝性储层屏蔽暂堵提供决策依据，而且能避免因测井等附加作业带来的风险。众多学者对钻井液漏失进行了理论计算和方法研究^[6-13]，然而国内利用实钻井漏数据预测天然裂缝宽度的研究开展得很少^[14-16]。国外O.Lietard等人^[17]利用实钻井漏数据建立了预测天然裂缝宽度的模型(以下简称为Lietard模型)，并在北海中部2口井中进行了应用，取得了较好的效果；F.Verga等人^[18]运用Lietard模型预测了3口井的天然裂缝宽度，预测结果与成像测井解释结果有较好的一致性。Lietard模型采用图版法求解，需要进行漏失曲线图版绘制、实钻井漏数据描点、理论漏失特征曲线和实钻漏失特征曲线人工匹配等一系列求解过程，存在求解过程繁琐、使用不便和易产生人为误差等问题。为此，笔者在Lietard模型的基础上，运用数值解法，并结合自适应收索法提出一种新解法。该解法可以实现无盲区自动匹配，提高Lietard模型的求解速度与精度，便于现场应用。

O.Lietard等人^[17]运用宾汉流体模型描述了发生井漏时钻井液在天然裂缝中的流动，其流动压力梯度为^[19]：

式中：p为压力，Pa；r为钻井液侵入半径，m；μ_p为钻井液塑性黏度，Pa·s；τ_y为动切力，Pa；v为钻井液漏失速度，m³/s；w为裂缝宽度，m。

漏失速度的表达式为：

式中：Q(t)为钻井液漏失速率，m³/s；V为钻井液漏失量，m³。

将式(2)代入式(1)，结合几何关系对式(1)积分，可得：

式中：Δp为井底压差，Pa；r_w为井眼半径，m。

为使式(3)便于求解，引入无因次时间及无因次半径：

其中

式中：t_D为无因次时间；r_D为无因次半径；t为时间，s；β为时间利用系数，s^-1。

将式(4)、式(5)代入式(3)，可得：

其中

式中：α为无因次有限侵入值。

式(7)的数值解为：

式(7)的初始条件为：r_D=1，t_D=0。

根据式(9)及初始条件，可求得一系列r_D，t_D值，再分别以lg t_D及lg(r²_D－1)为纵、横坐标，即可按一定步长的无因次有限侵入值α绘制理论漏失特征曲线，如图 1所示。

图1 理论漏失特征曲线及实钻漏失特征曲线平移示意 Fig.1 Shift sketch of theoretical and actual leakage characteristics curves

图选项

在钻井过程中，钻井液漏失体积计算公式为：

令纵、横坐标lg y，lg x分别为：

分析可知，不同无因次有限侵入值α有着与之对应的唯一理论漏失特征曲线。以及为纵横坐标的曲线称之为实钻漏失特征曲线，其与相对应的理论漏失特征曲线具有相同的形态及大小，只是余项lgw和－2lgw分别使实钻漏失特征曲线较与之对应的理论漏失特征曲线发生了垂直及水平平移(见图 1)。

图版法求解Lietard模型的步骤为：

1) 以及为纵横坐标，利用实钻井漏数据描点，绘制漏失曲线图版；

2) 采用垂直、水平平移方法寻找与实钻漏失特征曲线匹配的理论漏失特征曲线，得到与之对应的理论漏失特征曲线的无因次有限侵入值α，并反演出天然裂缝宽度。

在漏失曲线图版的绘制过程中，只能按照一定的步长进行，使2条相邻的理论漏失特征曲线间产生曲线盲区，在实钻漏失特征曲线匹配理论漏失特征曲线过程中，可能使找到的理论漏失特征曲线位于曲线盲区。从图 1可以看出，若实钻漏失特征曲线与α为0.06理论漏失特征曲线匹配，但该图版未能绘制α为0.06理论漏失特征曲线，则平移后实钻漏失特征曲线将位于α为0.1和0.01的2条理论漏失特征曲线之间，需人工估计相应的理论漏失特征曲线，这可能会产生人为误差。为避免或减小人为误差，可缩小α步长，加密漏失曲线图版的理论漏失特征曲线。但这不能从根本上消除人为误差，因为漏失曲线图版不能无限加密，总会存在实钻漏失特征曲线不能找到与之对应的理论漏失特征曲线而产生人为误差的可能性。

为解决Lietard模型求解速度慢且易产生人为误差的问题，分析可知，不同α值对应的理论漏失特征曲线均有一段近似直线段，为分析近似直线段的线性相关度^[20-21]，分别计算了1,0.1,0.001和0.000 1等代表不同数量级α值的线性相关系数，均达到0.999，可对近似直线段进行线性处理，近似直线段斜率分别为2.165 26，1.633 76，1.382 06，1.168 45和1.116 96(见图 1)。通过自适应搜索法自动设置α步长，求取不同α值所对应的理论漏失特征曲线直线段斜率。然后，取得实钻漏失特征曲线直线段数据，运用最小二乘法计算其斜率，比对二者斜率，直至其相对误差在合理范围内，确定对应斜率的理论漏失特征曲线，取其α值，反演出天然裂缝宽度。具体计算步骤为：

1) 输入基本参数Δp，r_w，μ_p，τ_y和允许匹配相对误差δ；

2) 求得实钻漏失特征曲线纵横坐标序列l和；

3) 建立方程组

4) 输入数据，得到矩阵，B=A^T，C=B×A；

5) 解方程组，得到b和k₂；

6) 计算相关系数

7) 设定步长α，α_max及α_min；

8) 计算不同α值对应理论漏失特征曲线近似直线段斜率k₁；

9) 计算相对误差

10) 若δ_r≤δ，记录α值，反演出裂缝宽度；否则，重复步骤(7)—(9)，直至δ_r≤δ。

为验证新求解方法的准确性及计算速度，选取X区块大、小不同井漏级别的A、B 2口井进行计算分析。A井在钻至垂深4 338.00 m时发生井漏，漏失速度为46 m³/h，由地质资料及邻井资料可知为裂缝性漏失。此时，井眼直径为311.1 mm，井底压差Δp为7.38 MPa，钻井液动切力τ_y为15 Pa，塑性黏度μ_p为39 mPa·s。结合其井漏数据(见图 2)，用新求解方法求解裂缝宽度，输入基本参数及相对误差δ=0.001%，将漏失数据序列V(n)和t(n)(取时间间隔为4 s)代入计算步骤2)，求出纵横坐标序列Y(n)、X(n)，求得实钻漏失特征曲线近似直线段的斜率为1.228 26，通过自适应搜索法自动找到与之匹配的理论漏失特征曲线，最终计算出α=0.001 8。将α=0.001 8代入α=，反演出裂缝宽度为527 μm。为对计算结果α=0.001 8进行验证，将纵横坐标为Y(n)、X(n)的实钻漏失特征曲线平移至α=0.1和α=0.001 8理论漏失特征曲线处进行匹配(见图 3)，对比分析与两理论漏失特征曲线的重合度可知，实钻漏失特征曲线与α=0.001 8的理论漏失特征曲线匹配度高。

图2 A井井漏数据 Fig.2 The leakage data of Well A

图选项

图3 A井实钻漏失特征曲线匹配结果 Fig.3 The matching results of actual leakage characteristics curve of Well A

图选项

B井在钻至垂深4 302.00 m时发生井漏，漏失速度为10.8 m³/h。井眼直径为311.1 mm，井底压差Δp为7.32 MPa，钻井液动切力τ_y为15 Pa，塑性黏度μ_p为38 mPa·s。在相对误差δ=0.001%情况下，用新求解方法计算得α=0.004 9，反演出裂缝宽度为195 μm。将实钻漏失特征曲线平移至α=0.1和α=0.004 9理论漏失特征曲线处进行匹配(见图 4)，实钻漏失特征曲线与α=0.004 9的理论漏失特征曲线匹配度高。

图4 B井实钻漏失特征曲线匹配结果 Fig.4 The matching results of actual leakage characteristics curve of Well B

图选项

从图 3、图 4可以看出，针对大、小不同级别井漏，Lietard模型的求解精度与实钻漏失特征曲线近似直线段的线性相关度呈正比关系；在同一相关度下，新方法因实现了无盲区自动匹配，匹配误差控制在0.001%，比图版求解法更准确地求得实钻漏失特征曲线所对应理论漏失特征曲线的无因次有限侵入值，从而提高了天然裂缝宽度的预测精度。

1) 分析Lietard模型漏失特征曲线可知，实钻漏失特征曲线及理论漏失特征曲线均有一段线性相关度为0.999的近似直线段，可作线性化处理，便于模型求解。

2) Lietard模型裂缝宽度求解新方法实现了计算机无盲区自动匹配理论漏失特征曲线，消除了图版法可能产生的人为误差，将匹配误差控制在0.001%内，提高了求解精度。

3) Lietard模型裂缝宽度求解新方法简化了求解过程，可以在数秒内完成模型求解，提高了求解速度。