共振冲击钻井技术是一项以适当的冲击振幅及转速，通过调节破岩工具的冲击频率使其达到正在钻进的地层岩石的固有频率，从而使钻头下岩石产生共振来实现高速破岩的钻井提速技术^[1-2]。该技术可以有效解决深部地层硬度高、研磨性强造成的机械钻速缓慢、钻具磨损严重等技术难题。作为冲击钻井技术之一，由于轴向简谐动载的引入使钻进系统的动力学特性变得更加复杂，为延长钻具寿命、提高钻进效率，开展共振钻井系统动力学特性研究十分重要。国外有关钻进系统动力学特性的研究已经有数十年的历史。一方面是针对轴向-扭转耦合运动模型展开研究，Richard等^[3]首次引入了状态依赖时滞变量，对自激振荡系统的动力学特性及其黏滑现象进行了分析，Kamel等^[4]在此基础上通过参数优选实现了机械钻速提高、成本降低的目的；另一方面是针对轴向动静载复合冲击展开研究，Krivtsov等^[5]研究动静载对材料去除率的影响，Cao等^[6]通过数值仿真方法研究动静载作用下的动力学系统局部分岔问题。国内目前研究主要集中于扭转方向黏滑振动方面^[7-8]，吕苗荣等^[9]研发了多自由度钻柱振动仿真软件；许帅等^[10]通过分数阶控制器，主动抑制钻柱的黏滑振动等。然而综合现有研究成果发现已开展的研究或是针对钻进系统轴向或是扭转方向单一方向上的动力学特性研究，即使针对钻进系统轴向-周向耦合运动展开研究，轴向方向上也仅考虑了单一静载作用。因此为研究共振钻井系统的动力学特性，笔者在已有冲击钻井动力学模型基础上，综合考虑轴向简谐动载和静载复合冲击作用以及轴向-扭转方向耦合运动，分别基于定切深和状态依赖时滞方法建立共振冲击系统的钻进模型，开展简谐动载关键参数对冲击系统动力学特性的影响规律分析。

1 共振冲击钻井系统动力学模型

1.1 共振冲击系统物理模型

共振冲击钻井系统包括轴向简谐动载-静载复合作用以及轴向钻进-周向扭转耦合运动，因此用轴向振动模型模拟轴向钻进运动，用扭转振动模型模拟扭转切削运动，构建轴向和扭转方向二自由度模型。

共振冲击钻井系统力学模型如图1所示，K_a为冲击系统的轴向刚度系数；C_a为冲击系统的轴向阻尼系数；v₀为大钩恒定加载速度；M为钻杆和底部钻具组合的集中质量；F₀和f分别为施加在钻头上的简谐动载幅值和频率；W₀为钻头所受稳态钻压；W为钻头所受反作用力；X为冲击系统相对于静平衡位置的轴向位移；K_t为冲击系统的扭转刚度系数；C_t为冲击系统的扭转阻尼系数；Ω₀为转盘的恒定角速度；J为钻杆和底部钻具组合的集中转动惯量；T为钻头所受扭矩；Φ为冲击系统相对于静平衡位置的角位移。

为了便于分析共振冲击钻井系统动力学特性，假设：①钻柱系统组成部分为一个整体；②钻柱系统没有水平运动；③井眼与钻柱系统均为垂直状态，钻柱系统各元件轴线始终与井眼轴线重合；④不考虑钻头等元件的磨损情况；⑤岩石的固有属性不随钻进深度变化而变化。

1.2 基于定切深方法的共振冲击系统钻进模型

定切深模型是指整个钻头的垂直切削深度不随时间的变化而变化，始终是一个定值，钻头所受反作用力与扭矩也始终为一个定值。根据牛顿定律，其动力学方程表达式为

式中，J为钻杆和底部钻具组合的集中转动惯量，kg·m²；$\ddot{\Phi }$为冲击系统角加速度，rad/s²；$\dot{\Phi }$为冲击系统角速度，rad/s；Ω₀为转盘的恒定角速度，rad/s；T为钻头所受扭矩，N·m；$\ddot{X}$为冲击系统轴向加速度，m/s²；$\dot{X}$为冲击系统轴向速度，m/s；X为冲击系统相对于静平衡位置的轴向位移，m。

整理方程（1）、（2）可得

实际上冲击系统受到的外力作用形式多样，力函数可能是周期性的，如简谐激励，也可能是非周期性的，如阶跃力、斜坡力等^[11]。根据振动力学相关理论，可以使用杜哈梅积分对任意外加激励的线性系统进行求解^[12]。因此对运动微分方程（3）、（4）进行求解，可得

其中

$\begin{array}{c}{U}_1\left(t\right)=-T+K_{\mathrm{t}}{\mathrm{\Omega }}_{0}t,{\omega }_{\mathrm{n}1}=\sqrt{K_{\mathrm{t}}/J},{\omega }_{\mathrm{d}1}=\sqrt{1-{\zeta }_1^2}{\omega }_{\mathrm{n}1},\\ U_2\left(t\right)=F_{0}cos\left(2\pi ft\right)+W_0-W+K_{\mathrm{a}}{v}_{0}t,\\ {\omega }_{\mathrm{n}2}=\sqrt{K_{\mathrm{a}}/M},{\omega }_{\mathrm{d}2}=\sqrt{1-{\zeta }_2^2}{\omega }_{\mathrm{n}2}.\end{array}$

式中，ω_d1和ω_n1分别为扭转方向有阻尼和无阻尼固有角频率，rad/s；ω_d2和ω_n2分别为轴向方向有阻尼和无阻尼固有角频率，rad/s；ζ₁和ζ₂分别为扭转方向和轴向方向阻尼比。

结合转盘转速和大钩恒定加载速度等初始条件，共振冲击动力学系统微分方程在轴向和扭转方向上的位移响应和速度响应^[13]分别为

$\Phi =-{\zeta }_1{\omega }_{\mathrm{n}1}\Phi +\exp \left(-{\zeta }_1{\omega }_{\mathrm{n}1}t\right)\left\{\left[{\mathrm{\Omega }}_0+\right.\right.$

$\left.\left.{\zeta }_1{\omega }_{\mathrm{n}1}\Phi \left(0\right)\right]\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{\mathrm{d}1}t\right)\right\}-$

1.3 基于状态依赖时滞的共振冲击系统钻进模型

在实际钻井过程中钻头的垂直切削深度不是定值，而是由每个切削齿的垂直切削深度累加得到。由于每个切削齿垂直切削深度不断变化，钻头所受钻压W与扭矩T也不断变化，在考虑顺应性的基础上，钻压与扭矩的表达式^[14-15]分别为

其中

$\begin{array}{c}{T}_{\mathrm{c}}=\frac{1}{2}\epsilon R_{\mathrm{b}}n{d}_n,T_{\mathrm{f}}=\frac{1}{2}{R}_{\mathrm{b}}^2\gamma \mu l\sigma ,\\ W_{\mathrm{c}}=\zeta \epsilon R_{\mathrm{b}}n{d}_n,W_{\mathrm{f}}=R_{\mathrm{b}}l\sigma ,\end{array}$

$H\left(y\right)=\left\{\begin{array}{c}1,y>0;\\ \frac{1}{2},y=0;\\ 0,y<0.\end{array}\right.\mathrm{s}\mathrm{g}\mathrm{n}\left(y\right)=\left\{\begin{array}{c}1,y>0;\\ 0,y=0;\\ -1,y<0.\end{array}\right.$

式中，T_c为扭矩切削分量，N·m；T_f为扭矩摩擦分量，N·m；W_c为钻压切削分量，N；W_f为钻压摩擦分量，N；H（y）为单位阶跃函数; sgn（y）为符号函数；ε为岩石的破碎比能，MPa；R_b为钻头半径，m；n为切削齿数量；d_n为每个切削齿的垂直切削深度，m；ζ为与切削力方向有关的参数；γ为大于1的钻头几何特征参数；μ为摩擦系数；l为钻头磨损平面长度，m；σ为接触应力，MPa。

由于钻头每个切削齿之间的夹角相同，为了更接近钻进系统实际情况，引入状态依赖时滞变量t_n。在不考虑钻头磨损的条件下，每个切削齿的垂直切削深度d_n相同，可根据轴向位移X计算得出：

时滞系统的运动规律不仅与当前运动状态有关，还依赖于过去的运动状态，并且极大影响着系统的动力学性质^[16]:

由于状态依赖时滞变量t_n难以显性表达，只能通过泰勒展开式求出近似解，并且由于非自治微分方程的存在，上述方程的解析解很难求出，故采用数值计算软件求出方程的数值解^[17]。令

$\begin{array}{c}{Z}_1=\Phi ,Z_2=\Phi ,Z_3=X,Z_4=\dot{X},U_1\left(t\right)=-T+K_{\mathrm{t}}{\mathrm{\Omega }}_{0}t,\\ U_2=F_{0}cos\left(2\pi ft\right)+W_0-W+K_{\mathrm{a}}{v}_{0}t.\end{array}$

并将其带入式（3）、（4）中，记作矩阵形式为

其中

$\begin{array}{c}A=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{K_{\mathrm{t}}}{J}& -\frac{C_{\mathrm{t}}}{J}& 0& 0\\ 0& 0& 0& 1\\ 0& 0& -\frac{K_{\mathrm{a}}}{M}& -\frac{C_{\mathrm{a}}}{M}\end{array}\right],Z=\left[\begin{array}{c}{Z}_1\\ Z_2\\ Z_3\\ Z_4\end{array}\right],\\ B=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& \frac{1}{J}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{M}\end{array}\right],U=\left[\begin{array}{c}{U}_1\\ U_1\\ U_2\\ U_2\end{array}\right].\end{array}$

传统求解动力学系统模型主要采用欧拉法，但是欧拉法只估计了时间段t_n与t_n+1的左端点的导数，当迭代次数过多时其解的精确度下降，故本文中采用精度更高的龙格-库塔方法进行求解^[18]。龙格-库塔方法是通过已知方程式和其区间中可能存在的解作为初值，推导出新的方程求出近似解的方法，求解的阶数越高，结果精度越高^[19]。实际工程问题中也常采用四阶龙格-库塔方法来求解微分方程^[20]，表达式为

式中，h为步长，一般取10^-5~10^-4 s；k_i为龙格-库塔系数矩阵。

2 模型结果分析

应用数值计算软件，分别采用杜哈梅积分和四阶龙格-库塔方法对共振冲击钻井系统定切深模型和状态依赖时滞模型求解，分析冲击系统的动力学特性。模型计算参数均结合实验数据或现场工况确定，其中初始条件为：大钩加载速度v₀=0.015 m/s，转盘转速Ω₀=4π rad/s。

2.1 定切深模型计算结果

分析中简谐动载频率f为100 Hz，幅值F₀为60 kN，求解共振冲击系统的速度和位移响应。这里冲击系统非齐次微分方程的通解代表系统的全响应，可表示为齐次解和特解之和，齐次解代表系统的暂态响应，特解代表系统的稳态响应。如图2（a）所示，在轴向静载和简谐动载复合作用下共振冲击系统轴向速度全响应曲线呈简谐形式变化，开始时振荡剧烈，且最大反向速度出现了较大的负值。同时由其包络线形态可知，速度响应每个周期的峰值以负指数函数形式衰减。随着时间的增加，轴向速度全响应曲线逐渐衰减到与稳态响应曲线重合。由图2（b）可以看出，轴向振动冲击系统稳态响应位移与时间t成正比，其斜率与轴向刚度系数K_a和大钩速度v₀有关。在考虑初始条件、常值力W₀-W、斜坡力K_av₀t、简谐动载F₀cos（2πft）作用下，轴向位移全响应曲线在开始时呈周期性变化，且曲线上升、下降交替出现，说明此时钻头出现跳钻现象。随着振荡幅值下降，钻头跳钻现象逐渐减弱，轴向位移全响应曲线开始向稳态响应曲线贴近。

如图3（a）所示，扭转冲击系统的角速度全响应曲线与轴向速度全响应曲线类似，随着时间推移，最终也将与其稳态响应曲线重合。角速度稳态响应曲线为一条水平线，不随时间变化，只与转盘转速有关。如图3（b）所示，扭转冲击系统的稳态响应角位移与时间t成正比，其斜率与扭转刚度系数K_t和转盘转速Ω₀有关。在常值扭矩T和斜坡扭矩K_tΩ₀t的作用下，角位移全响应曲线出现周期性变化，且曲线上升、停滞交替出现，说明此时钻头发生黏滑现象。随着振荡幅值下降，钻头黏滑现象逐渐减弱，角位移全响应曲线开始向稳态响应曲线转变。

2.2 状态依赖时间迟滞模型计算结果

为分析简谐动载下状态依赖时滞模型动力学特性，对不同频率、不同幅值下共振钻井系统最大反向速度进行计算。如图4所示，无论是在冲击系统轴向还是扭转方向，钻进速度和角速度均随动载频率的增大而下降，随着动载幅值的增大而增加。根据计算结果，选取3组代表性的幅频组合进行分析，结果如图5所示。

由图5（a）可以看出：前20 s内，角速度、轴向速度每一周期的最大值均不断增大；在20~100 s内，其最大值几乎不变，均趋于稳定。轴向最大反向速度在前20 s内不断增加，出现了负值，在20~100 s内，轴向最大反向速度几乎保持不变，说明钻头此时产生了跳钻现象。对于角速度，在0~100 s内每个周期的角速度最小值均接近0，说明钻头一直存在黏滑问题。

如图5（b）所示，共振冲击系统的轴向速度和角速度均按负指数形式衰减，冲击系统轴向钻进和周向切削均进入稳定运动状态，轴向方向的跳钻现象和扭转方向的黏滑现象消失。随着时间的增加，轴向速度和角速度振动幅值越来越小，趋于平衡解，共振冲击系统几乎与大钩、转盘保持同步运动。

如图5（c）所示，当简谐动载频率f为25 Hz、幅值F₀为30 kN时，冲击系统扭转方向再次出现黏滑现象。从第3个周期开始，每个周期的角速度最大值未发生明显变化，此时角速度达到转盘转速的2.2倍左右。需要注意的是，钻头角速度过大会产生较大应力，加剧底部钻具组合失效。在轴向方向上，跳钻现象也再次出现，从第四个周期开始轴向最大反向速度达到最大且几乎保持不变。

基于对共振冲击系统动力学特性分析可知，高频低幅的简谐动载对钻井效率有利。高频的范围主要取决于共振模块的工作原理，通过机械方式驱动的其频率可达几十赫兹，通过电驱动的其高频超过200 Hz，低幅的范围一般使其不超过钻压为优。在实际工程应用中，可基于共振冲击系统的动力学模型，结合井下动力钻具类型、岩石力学特性等因素进行作业参数和施工方案优选。

3 共振冲击系统动力学模型对比

3.1 动力学特性响应对比

在频率f=100 Hz，幅值F₀=60 kN且相同初始条件下对共振冲击系统定切深和状态依赖时滞两种动力学模型进行对比分析。图6、7分别为基于两种模型条件下冲击系统轴向和扭转方向振动响应结果。

如图6所示：基于定切深模型的轴向速度响应曲线仅在初始阶段出现负值，随着时间的增加，轴向速度振动幅值逐渐减弱，也正因如此，其相应的轴向位移响应曲线在初始的剧烈波动后，逐渐稳定成正比于大钩加载速度的直线；而基于状态依赖时滞的轴向速度响应曲线在初始阶段没有出现负值，从18 s开始曲线呈现稳定的周期性变化，最大反向速度接近0.04 m/s，相应地其轴向位移响应曲线在初始时波动不明显，随着时间的增加，开始周期性的阶梯式变化。

如图7所示：基于定切深模型的冲击系统角速度响应曲线仅在初始阶段出现停滞现象，因此其相应的角位移响应曲线也仅在初始阶段出现阶梯式跃进变化；而基于状态依赖时滞模型的角速度响应曲线在经历了初始递增阶段后幅值趋于稳定，其最大反向速度接近0，相应角位移响应曲线始终出现周期性的阶梯式变化。

此外，基于定切深模型得到的轴向与扭转方向动力学特性响应频率与系统轴向、扭转方向无阻尼固有频率较为接近，而采用状态依赖时滞模型得到的两个方向动力学特性响应频率均小于系统的无阻尼固有频率。在相同参数条件下基于状态依赖时滞模型的系统响应滞后于基于定切深模型的系统响应。

在相同初始条件和简谐动载激励下基于两种模型得到的动力学特性响应曲线无论是形状还是频率均存在差异。这是由于定切深模型是将轴向和扭转方向运动视作两个互相独立的单自由度系统，其切削齿的垂直切削深度、系统所受的钻压和扭矩均是定值。而采用状态依赖时滞模型进行求解时，冲击系统轴向和扭转方向相互耦合，切削齿垂直切削深度以及系统所受钻压和扭矩基于钻进状态不断发生改变。因此状态依赖时滞模型更接近现场实际钻进状态，可以研究系统更复杂的动力学行为，模型的适应性更强，应用性更广，但是求解也更复杂。定切深模型虽然假设前提更多，但其求解结果也基本符合客观实际，求解也相对简单，在钻进深度较小、系统不产生复杂动力学特性时可作为替代方案。

3.2 结果分布情况对比

采用统计学箱线图方法^[21]对两种模型结果的分布情况进行对比分析。该统计方法利用上边缘、下边缘、中位数、上四分位数、下四分位数5个统计量来描绘数据分布和性质^[22]，其结构如图8所示。采用箱线图分析数据时通常设置上边缘、下边缘进行检测。当数据大于上边缘或小于下边缘时被认为是离群点。

共振冲击系统响应速度箱线图如图9所示。在500001个采样点中，基于定切深模型计算得到的共振冲击系统轴向速度数据中离群点数为110877，基于状态依赖时滞模型得到的离群点数为1405，说明采用状态依赖时滞模型得到的轴向速度数据点更为集中；基于定切深模型得到的角速度数据中偏离数据主要分布区域的采样点远多于状态依赖时间迟滞模型，两者分别为66744与0。

离群点数量过多说明模型波动较大，稳定性较差。在实际钻井过程中，钻头会出现黏滞、跳钻等现象，使钻进参数产生偏离，但其波动程度仍在一定范围内，系统仍相对稳定。因此通过箱线图对比可知，采用状态依赖时滞模型得到的数据更贴近真实钻柱在井下运动的状况

4 结论

（1）在一定钻压、转速、大钩速度条件下，提高轴向简谐动载的幅值或降低频率会加剧扭转方向的黏滑效应和轴向方向的跳钻现象。在工程实际中，采用高频低幅的简谐动载有利于提高钻井效率，可结合井下动力钻具类型、岩石力学特性等因素进行作业参数和施工方案优选。

（2）基于定切深模型的冲击系统跳钻和黏滞现象会随时间增加而减弱，而基于状态依赖时滞模型得到的系统响应则取决于钻进状态，状态依赖时滞模型更接近于钻井实际工况。

（3）无论是轴向还是扭转方向上，基于定切深模型得到的离群点数都远多于状态依赖时滞模型离群点数，采用状态依赖时滞模型得到的数据更集中，更贴近真实钻柱在井下运动的状况。

参考文献