A new calculation method of well trajectory of joint control

在滑动钻进时，钻具组合是给定的，钻具组合造斜能力和工具面角基本上保持不变，因此用恒工具面角模型来描述滑动钻进所形成的井眼轨迹比空间圆弧模型更恰当一些，因为在空间圆弧模型中，工具面角是随井深而变化的，而恒工具面角模型^[1-4]是满足“井眼曲率和工具面角恒定”的曲线。

井斜角变化规律为

方位角变化规律为

式中，$\omega $为工具面角，$\kappa $为井眼曲率，$\Delta L$为井深增量，$\alpha $和${\alpha _0}$分别为井眼轨迹上任意井深和初始井深处的井斜角，$\phi $和${\phi _0}$分别为井眼轨迹上任意井深和初始井深处的方位角。

定义

根据洛必达法则^[5]可知：$\mathop {\lim }\limits_{\alpha \to {\alpha _0}} f(\alpha ,{\alpha _0}) = \sin {\alpha _0}$，将式(1)代入式(2)，可将方位角变化规律简写为

记$\Delta N = N - {N_0}$、$\Delta E = E - {E_0}$、$\Delta H = H - {H_0}$分别为北坐标增量、东坐标增量和垂深增量，则当$\alpha = {\alpha _0}$时，

当$\alpha \ne {\alpha _0}$时，

式(8)、式(9)中的积分需要使用数值积分公式^[5]来计算，例如Simpson数值积分公式，但是计算量比较大。方敏等^[6]使用幂级数展开技巧给出了一种计算量较小的新方法。

假设钻具组合造斜率${\kappa _{\rm{h}}}$是已知的，单根长度为$\Delta {L_{\rm{g}}}$，复合钻进进尺为$\Delta {L_1}$，滑动钻进进尺为$\Delta {L_2}$、工具面角为${\omega _2}$，滑动钻进形成的井眼轨迹为恒工具面角曲线，则成立以下方程组

式中，${\alpha _{\rm{A}}}$和${\phi _{\rm{A}}}$分别为当前井底的井斜角和方位角，${\alpha _{\rm{B}}}$和${\phi _{\rm{B}}}$分别为预定钻达的井斜角和方位角。

在方程组(11)，有3个未知数$\Delta {L_1}$、$\Delta {L_2}$、${\omega _2}$，只要任给一个$\Delta {L_1}$值，就可以确定出$\Delta {L_2}$值，再从后2个方程确定出工具面角${\omega _2}$。从方程组(11)可得

工具面角的定义范围是$0 \leqslant {\omega _2} < 2{\rm{\pi }}$。在工程上常常使用−180°$ \leqslant {\omega _2} < $180°，故本文公式中约定$ - {\rm{\pi }} \leqslant {\omega _2} < {\rm{\pi }}$。

由于三角函数的周期性，单独从式(11)的第1式或第2式并不能唯一确定角度${\omega _2}$，参见图1，正确的计算公式是：

(1)当${s^2} + {c^2} \ne 1$时无解；

(2)当$\left| s \right| \leqslant 1$并且$\left| c \right| \leqslant 1$时：

需要注意的是，如果已知参数给定的不合适，可能求不出滑动钻进的工具面角${\omega _2}$。$\Delta {L_1}$必须满足的约束条件为

其中，${\rm{\delta }}L = $0.5 m。

Figure 1.  Variation range of tool face angle

由于适当给定复合进尺$\Delta {L_1}$就可以求出滑动钻进工具面角${\omega _2}$、使得单根结束时的井斜角和方位角达到预定值，这说明控制方案可以有很多种。这么多种设计方案中，必有一种设计方案是最佳的。

在做设计方案时，除了知道要控制的井斜角${\alpha _{\rm{B}}}$和方位角${\phi _{\rm{B}}}$之外，实际上还知道目标点(B)的NEH坐标$({N_{\rm{B}}},{E_{\rm{B}}},{H_{\rm{B}}})$。不同的设计方案结束点的NEH坐标是不同的，它与目标点(B)的距离也是不同的，认为这个距离达到最小时的设计方案是最佳的。

设井底(A)的NEH坐标为$({N_{\rm{A}}},{E_{\rm{A}}},{H_{\rm{A}}})$，复合钻进结束时的NEH坐标为$({N_{\rm{1}}},{E_{\rm{1}}},{H_{\rm{1}}})$，滑动钻进结束时的NEH坐标为$({N_{\rm{2}}},{E_{\rm{2}}},{H_{\rm{2}}})$，$\Delta {N_{\rm{2}}}$、$\Delta {E_{\rm{2}}}$、$\Delta {H_{\rm{2}}}$由式(5)~(7)或者式(8)~(10)来计算，则有：

滑动钻进结束点与目标点B的距离为

改变段长$\Delta {L_1}$，如果距离达到最小值，则是最佳方案，此时的$\Delta {L_1}$、$\Delta {L_2}$、${\omega _2}$即为最佳控制参数。图2给出了某算例距离随$\Delta {L_1}$变化的曲线，可见当$\Delta {L_1}$增加时，距离逐渐减小，达到最小值后又逐渐增大；距离曲线的几何特征是(单峰)下凸连续曲线，故可使用二分法求最小值点，下面给出两种数值算法。

Figure 2.  Variation of the distance from the target point with the compound drilling footage

给定复合进尺分辨率dL=0.01 m，在$\Delta {L_1}$的允许取值区间(式(14))中取等间距点

若$\Delta L_1^{(n)} > \Delta {L_{1\max }}$，则取$\Delta L_1^{(n)} = \Delta {L_{1\max }}$，取点结束。

对于所有的$\Delta L_1^{(n)}$，按第2节方法计算工具面角${\omega _2}$，再按式(16)~(18)计算滑动钻进结束点与目标点B的距离${d^{(n)}}$，令

则与$\Delta L_1^{(m)}$相对应的$\Delta {L_1}$、$\Delta {L_2}$、${\omega _2}$即为最佳控制参数的近似值。

如果单根长度$\Delta {L_{\rm{g}}} \approx $10 m，则算法一至多需要1 000次计算就可以求出最佳控制参数的近似值，误差至多为${10^{ - 2}}$，这对于工程应用来说，近似值已经足够准确，但是缺点是计算量比较大。

简记$Y = \{ \Delta {L_1},{\omega _2},d\} $，称之为一个计算点。给定复合进尺精度${\rm{d}}L = {10^{ - m}}$和距离精度$\varepsilon = {10^{ - p}}$，其中$m \geqslant 2$、$p \geqslant 4$,均为正整数。

(0)初始：令

分别计算滑动钻进结束点与目标点B的距离，相应的计算点记为${Y^{[1]}}$、${Y^{[3]}}$、${Y^{[3]}}$。

如果${d^{[m]}} = \min \{ {d^{[1]}},{d^{[2]}},{d^{[3]}}\} \leqslant \varepsilon $，则计算点${Y^{[m]}}$记为最小值点，结束。

(1)二分：再令$ \Delta {L}_{1}^{[a]}=0.5(\Delta {L}_{1}^{[1]}+\Delta {L}_{1}^{[2]})$为前2个段长的平均值，$ \Delta {L}_{1}^{[b]}=0.5(\Delta {L}_{1}^{[2]}+\Delta {L}_{1}^{[3]})$为后2个段长的平均值，分别计算滑动钻进结束点与目标点B的距离，相应的计算点记为${Y^{[a]}}$和${Y^{[b]}}$。

(2)判断：如果${d^{[m]}} = \min \{ {d^{[a]}},{d^{[2]}},{d^{[b]}}\} \leqslant \varepsilon $，则计算点${Y^{[m]}}$为最小值点，循环结束。

如果$\Delta L_1^{[3]} - \Delta L_1^{[1]} \leqslant {\rm{d}}L$，则在给定的精度条件下无解，但是可以将${d^{[m]}} = \min \{ {d^{[1]}},{d^{[a]}},{d^{[2]}},{d^{[b]}},{d^{[3]}}\} \leqslant \varepsilon $的计算点${Y^{[m]}}$作为近似最小值点，循环结束。

(3)选择：当${d^{[a]}} \leqslant {d^{[1]}}$且${d^{[a]}} \leqslant {d^{[2]}}$，称$1/a/2$是一个下凸构型；当${d^{[b]}} \leqslant {d^{[2]}}$且${d^{[b]}} \leqslant {d^{[3]}}$，称$2/b/3$是一个下凸构型。

如果$1/a/2$和$2/b/3$都是下凸构型，则：当${d^{[a]}} \leqslant $$ {d^{[b]}}$时，令${Y^{[3]}} = {Y^{[2]}}$、${Y^{[2]}} = {Y^{[a]}}$，转(1)；当${d^{[a]}} > {d^{[b]}}$时，令${Y^{[1]}} = {Y^{[2]}}$、${Y^{[2]}} = {Y^{[b]}}$，转(1)。如果只有$1/a/2$是下凸构型，则令${Y^{[3]}} = {Y^{[2]}}$、${Y^{[2]}} = {Y^{[a]}}$，转(1)。如果只有$2/b/3$是下凸构型，则令${Y^{[1]}} = {Y^{[2]}}$、${Y^{[2]}} = {Y^{[b]}}$，转(1)。否则令${Y^{[1]}} = {Y^{[a]}}$、${Y^{[3]}} = {Y^{[b]}}$，转(1)。

由于复合进尺区间每次循环后都缩小1/2，故至多计算M次就可以求出近似最小值点。

取${\rm{\delta }}L = $0.5 m，$\Delta {L_{1\max }} = $9.5 m，${\rm{d}}L = {10^{ - 6}}$ m，则M≈ 23，这个次数比算法一的1 000次减少50多倍。

为检验算法二的性能，制作一个算例。设计钻具组合造斜能力${\kappa _{\rm{h}}} = $7.87(°)/30 m，单根长度$\Delta {L_{\rm{g}}} = $10 m，复合钻进长度$\Delta {L_1} = $3 m，滑动钻进长度$\Delta {L_2} = $7 m，工具面角${\omega _2} = $60°，井底井斜角${\alpha _{\rm{A}}} = $30°、方位角${\phi _{\rm{A}}} = $120°，使用这些数据计算出来的单根结束时的井斜角、方位角、坐标增量见表1，其中实数至多保留小数点后8位有效数字是为了考察算法二的性能。算法二的精度参数取${\rm{d}}L = \varepsilon = {10^{ - 10}}$。其求解结果为：${\omega _2} = $60.000 000 012 7°，$\Delta {L_1} = $2.999 999 997 3，$\Delta {L_2} = $7.000 000 002 7，$d = $9.263×10⁻¹¹，其他计算结果见表2。求解结果表明，算法二非常正确地求出了最佳设计参数。

现场计算中，数据一般只能保证小数点后2位有效数字。将表1中数据保留小数点后2位有效数字，见表3。算法二的求解结果为：${\omega _2} = $60.277 1°，$\Delta {L_1} = $2.926 7，$\Delta {L_2} = $7.073 3，$d = $0.001 85，其他计算结果见表4。求解结果表明，算法二求出的最佳设计参数与制作算例时使用的参数非常接近。

待钻井眼轨设计数据的井身点的数据见表5，由于井眼轨道设计使用的是最小曲率法(空间圆弧模型)，工具面角在设计井段上不是常数，表中的11.97°是设计井段开始点的工具面角，11.37°是结束点的工具面角。使用“复合+滑动”钻进方式进行施工参数优化设计，钻具组合造斜能力${\kappa _{\rm{h}}} = $7.87(°)/30 m，使用算法二求得：${\omega _2} = $19.352 1°，$\Delta {L_1} = $2.970 0，$\Delta {L_2} = $7.029 97，$d = $0.029 9，井身参数计算结果见表6。

(1)对于现场单根控制问题，用“线段+恒工具面曲线”建立了“复合+滑动”钻进数学模型，给出了使用正弦值和余弦值计算工具面角的正确公式。

(2)对于单根控制设计的求解，使用最优化思想给出了两个数值算法，尤其是算法二的计算量很小，满足实时控制的要求。

(3)使用理论模型算例和实际钻井数据对算法进行了大量验证，结果表明，在理论模型算例的情形，算法二能够无限精确地反求出理论模型参数；在实际钻井数据的情形，算法二也能比较精确地求出控制参数。

(4)本文是针对现场定向井工程师研制的数学模型和求解算法，既可以适用于现场定向井工程师用于单根钻进时进行施工参数设计，也可以适用于电控系统的实时反馈控制。