基于状态依赖时滞的钻柱动力学稳定性分析

2022-12-01狄勤丰王文昌段浩宇

振动与冲击 2022年22期

张鹤，狄勤丰,2，王文昌,2，陈锋，段浩宇

(1.上海大学上海市应用数学和力学研究所，上海 200072；2.上海市能源工程力学重点实验室，上海 200072；3.上海大学机电工程与自动化学院，上海 200072)

钻井工程是石油、天然气、页岩气勘探开发的重要环节之一，主要为油、气提供生产通道。钻柱作为钻井工程中的核心部件，由钻杆和底部钻具组合(bottom hole assembly,BHA)构成。BHA主要由钻头、钻铤、稳定器以及井下测量和辅助钻井工具等组成，为一长度约为100～300 m的管柱结构。钻井过程中，钻柱在地面设备驱动下，在充满液体或气体的狭长井眼中转动，并带动钻头旋转破碎井底岩石，形成井筒(见图1)。钻头通过切削齿挤压或剪切的方式破坏岩石，与岩石的相互作用会激发钻柱的振动，其主要有三种形式[1]：轴向振动[2]、扭转振动[3]和横向振动[4-5]，三种振动在井下相互耦合，其最剧烈、最具破坏性的极端表现形式为跳钻、黏滑和涡动。

自二十世纪九十年代始，聚晶金刚石复合片(polycrystalline diamond compact，PDC)钻头得到广泛应用，但使用PDC钻头进行深井钻探时，钻柱容易发生扭转振动甚至恶化为黏滑振动，其主要表现为钻头处于停钻(钻头转速为0，称为黏滞相)和高速转动(钻头最大转速达到地面恒定转速的两倍及以上，称为滑脱相)的周期交替振荡状态。据统计[6-7]，在深井钻探中，钻柱发生黏滑振动的时长占钻进时间的一半以上。钻柱黏滑振动引起的交变应力容易造成钻具接头的疲劳失效，加剧钻头磨损，降低钻井效率。因此，本文主要以钻柱的扭转黏滑振动为研究对象，通过对描述钻柱黏滑振动的数学模型进行稳定性分析，并在此基础上研究钻井参数与钻井液阻尼对钻柱黏滑振动的影响。

图1 钻井和钻柱系统示意图(改自文献[8]) Fig.1 Schematic of drilling and drillstring systems (Modified from [8])

早期对钻柱黏滑振动的研究主要是对扭转振动单独建模求解，不考虑扭转振动与轴向或横向振动的相互耦合[9-10]，并认为钻头反作用扭矩随钻头转速增加而减小的速度弱化效应是造成钻柱发生黏滑振动的主要原因，并提出不同的函数关系来表示钻头反扭矩与转速之间的速度弱化规律[11-19]，在一定程度上重现了钻柱的黏滑振动。然而，随着对振动机理认识的不断深入，人们发现钻柱不同振动模式之间存在相互耦合，需要考虑耦合效应才能正确地理解钻柱的黏滑振动。同时，单个钻头切削齿切削岩石的试验结果与速度弱化效应并不相符，表明了钻头和岩石之间的作用力与钻头运动速度无关[20]。

Detournay等[21]在单个钻头切削齿切削岩石试验结果的基础上，提出了与钻头转速无关的钻头-岩石相互作用模型，该模型将钻头与岩石相互作用分为齿切削面对岩石的切削过程以及齿磨损面与岩石的接触摩擦过程。Richard等[22]在上述钻头-岩石相互作用模型的基础上，考虑了BHA的轴向和扭转运动，建立了低维钻柱动力学离散模型(文献中称为RGD模型)。RGD模型将钻头的切削深度表示为钻头在当前时刻与之前某一时刻轴向位置的差值，进而在模型中引入了取决于钻头的运动状态的时滞变量，使得钻头的轴向和扭转振动通过该状态依赖时滞相互耦合，其运动控制方程为状态依赖时滞微分方程。数值模拟结果表明，钻头切削岩石过程的状态依赖时滞是造成钻柱发生自激黏滑振动的根本原因，而钻头扭矩表现出的速度弱化效应仅为钻头发生黏滑振动后的宏观表现。由于RGD模型中的钻柱模型和钻头的几何形状较为简化，因此近年来很多国外学者对其进行了改进，其中主要包括：①在钻柱模型中考虑钻杆的轴向刚度和钻井液阻尼[23-25]；②建立钻柱的多自由度有限元模型或连续模型[26-29]；③研究几何形状更为复杂的钻头，其与岩石的相互作用涉及两个甚至多个状态依赖时滞[30-31]；④通过定义二元函数来表示时滞变量，避免向钻柱模型中引入状态依赖时滞[32-34]；但国内在这方面的研究鲜有报道。

对非线性动力学系统进行线性稳定性分析，可以定性地研究系统参数对系统稳态运动的影响。文献中对RGD及其改进模型的稳定性分析主要分为两种：一种是直接对状态依赖时滞微分方程进行线性化，得到线性系统的特征方程，进而得到特征方程根轨迹的解析解[35]；另一种是将状态依赖时滞微分方程转化为非线性耦合的偏微分-常微分方程组，并对其进行线性化求解[36]。针对RGD模型的线性稳定性分析指出，由于忽略阻尼的影响，因此其在任何参数下均不稳定。然而，以上解析方法仅适用于低维钻柱模型，难以应用到多自由度钻柱动力学系统的稳定性分析。鉴于此，本文考虑阻尼的影响，利用Zhang等研究中的钻头轨迹函数，将RGD模型的状态依赖时滞微分方程转化为非线性耦合的偏微分和常微分方程。在此基础上，利用谱方法和伽辽金方法将线性化后的偏微分和常微分方程离散为耦合的常微分方程组，以实现对RGD模型的稳定性分析。此外，本文提出的稳定性数值分析方法还容易推广到带状态依赖时滞的多自由度钻柱动力学系统的稳定性分析。

1 数学模型

1.1 钻柱模型

RGD模型在轴向上，仅考虑了BHA的集中质量Mb，忽略了钻杆的轴向刚度和轴向阻尼，这是因为实际钻井中钻柱在轴向上、下边界均为力边界条件(上边界为大钩载荷H0，下边界为钻头-岩石相互作用的钻压Wb)，因此钻杆轴向刚度和轴向阻尼对BHA的轴向运动几乎没有影响(见图2)；而钻柱扭转方向的上、下边界分别为位移边界和力边界(上边界为地面恒定转速Ω0，下边界为钻头和岩石相互作用的扭矩Tb)，因此在扭转方向将BHA简化为转动惯量为Jb的刚体，将钻杆简化为扭转刚度为Kp的扭转弹簧，忽略了扭转阻尼的影响，从而使得RGD模型在任何参数下均不稳定。

基于此，本文在RGD模型的基础上考虑了扭转方向的阻尼以研究其对RGD模型稳定性的影响。由图2可知，RGD模型考虑了BHA的轴向和扭转两个自由度，其轴向和扭转位移(速度)分别用Ub(Vb)和Φb(Ωb)表示。根据钻柱的上、下边界条件，可以推导出控制BHA轴向和扭转运动的常微分方程(ordinary differential equation,ODE)分别为

(1)

(2)

图2 钻柱动力学模型示意图 Fig.2 Schematic diagram of drill string dynamics model

1.2 钻头-岩石相互作用模型

钻头处的钻压Wb和扭矩Tb可以分别由两个分量构成[37]，即

Wb=Wc+Wf,Tb=Tc+Tf

(3)

式中：Wc和Wf分别为钻压的切削和接触摩擦分量；Tc和Tf分别为扭矩的切削和接触摩擦分量。

切削分量Wc和Tc的大小均与钻头的切削深度d成正比，即

(4)

式中：a为钻头半径；ε为岩石的本征比能；ζ描述了切削力在切削平面上的方位。钻头的准螺旋运动在井底形成一定的岩石轮廓，而切削深度取决于相邻刀翼间的井底岩石轮廓，如图3所示。

图3 相邻刀翼间井底岩石轮廓及切削深度示意图(改自Richard等的研究) Fig.3 Section of the bottom hole profile between two adjacent blades (Modified from Richard,et al)

Richard 等研究中通过引入轴向再生效应来表示钻头的切削深度，其大小为钻头轴向位移在当前时刻与此前某一时刻的差值

d=ndn,dn=Ub(t)-Ub(t-tn)

(5)

式中：n为钻头的刀翼个数；dn为单个刀翼的切削深度，由于钻头做刚体运动且刀翼沿钻头轴向对称均匀分布，因此每个刀翼的切削深度均相同；ub(t)为钻头当前时刻的轴向位移，Ub(t-tn)为钻头在t-tn时刻的轴向位移，时滞时间tn为钻头绕其轴线转过2π/n(相邻刀翼间的夹角)所需的时间，即

(6)

式中，Φb(t)和Φb(t-tn)分别为钻头在当前时刻和t-tn时刻的扭转角位移。由于钻头存在扭转振动，因此时滞tn不为常数，其大小取决于钻头的运动状态，也即tn为状态依赖时滞。

钻压和扭矩的接触摩擦分量Wf和Tf之间满足以下约束关系

(7)

式中：μ为钻头切削齿磨损面与岩石间的摩擦因数；γ为与钻头几何形状相关的常量，对于RGD模型中所用钻头，其取值为1。Wf的取值与钻头的轴向速度Vb相关

Wf=anlσ[1-g(Vb)]

(8)

式中：σ为钻头与岩石平均接触应力的最大值；g(Vb)为集值函数，其定义为

(9)

由于钻头切削过程中引入了状态依赖时滞tn，因此将式(3)～式(9)描述的钻头-岩石相互作用代入微分方程式(1)～式(2)中，可以得到控制钻头轴向和扭转运动的状态依赖时滞微分方程(state-dependent delay differential equation,SDDDE)；另一方面，钻头与岩石的接触摩擦过程向方程中引入了集值函数g(Vb)，从而使钻头处的边界条件呈现出极强的非线性。

1.3 钻头轨迹函数

为了描述钻头的空间运动轨迹，需要建立整体和局部两个坐标系，如图4所示。整体坐标系ORXΘ的原点O固定在井口；R为原点指向井壁的径向坐标；X沿井眼轴线方向，用以表示钻头的轴向位置；Θ正方向与钻头扭转运动正方向一致，用以表示钻头的扭转角位移。局部坐标系orxθ的原点固定在钻头上随钻头一起转动；r沿钻头半径方向，用以表示钻头切削齿在半径方向上的分布；x指向钻头钻进方向，表示钻头切削齿在钻头上的轴向位置；θ正方向与钻头转动正方向一致，表示钻头相对于整体坐标系的转动。RGD模型假设刀翼沿钻头轴线对称均匀分布，且钻头切削齿沿半径方向紧密排列，因此0≤r≤a；且钻头切削齿均在同一水平面上。若令参考刀翼的坐标为x=0,0≤r≤a,θ=0，则沿θ正向第i个刀翼的坐标可表示为x=0,0≤r≤a,θ=2π(i-1)/n。

图4 整体和局部坐标系示意图 Fig.4 Schematic of global and local cylindrical coordinate systems

由于钻头做刚体运动且刀翼沿钻头轴线对称均匀分布，因此钻头轨迹在相邻刀翼间形成的井底轮廓曲线均相同，如图5所示(4刀翼钻头运动轨迹)，其中刀翼的轨迹函数由方程X=H[Φb(t)-θ]表示。在此基础上定义钻头轨迹函数

h(θ,t)=-H[Φb(t)-θ],θ∈[0,2π/n],t>0

(10)

根据图5和式(5)可知，单个刀翼的切削深度dn可用钻头轨迹函数h(θ,t)表示为

(11)

此外，在整体坐标系ORXΘ中，刀翼的轨迹函数X=H(Θ),0<Θ<Φb(t)与时间t无关，因此在给定坐标Θ=Φb(t)-θ处

(12)

(13)

由式(11)可知，通过引入钻头轨迹函数表示钻头切削深度，避免了向方程中引入状态依赖时滞。此外，钻头轨迹函数的变化受式(13)控制，将钻头切削深度式(11)代入钻压和扭矩切削分量表达式(3)～式(4)并进一步代入到式(1)～式(2)后，即可将SDDDE转化为耦合的PDE-ODE。

图5 井底岩石轮廓的平面展开示意图 Fig.5 Planar schematic of bottom hole profile

1.4 无量纲化

首先定义系统的特征时间和特征长度

(14)

(15)

式中，U0和Φ0分别为钻头稳态钻进时的轴向和扭转位移。无量纲形式的切削深度δ0,δn和时间τ可分别表示为

(16)

由此，式(1)和式(2)的无量纲形式为

(17)

其中，

(18)

(19)

式中，v0为钻头稳态钻进时的轴向速度的无量纲形式。

(20)

(21)

ћ(0,θ)=-ub(τ)

(22)

式中，ω0为钻头稳态运动时扭转速度的无量纲形式。在此基础上，钻头切削深度表达式(11)可以用无量纲形式的钻头轨迹函数表示为

(23)

2 线性化及稳定性分析

2.1 线性化

在非线性耦合PDE-ODE稳态解的基础上施加小扰动，可以实现对该方程组的线性化处理。由式(15)及偏微分方程初始条件式(21)可知PDE-ODE方程的稳态解为

(24)

(25)

式中的扰动均为小量。

当钻头稳态运动时，g(vb)=0，因此将扰动表达式(25)代入式(17)～式(18)、式(20)～式(23)，忽略扰动小量的乘积并整理后，可得线性化后控制扰动变量的耦合PDE-ODE为

(26)

(27)

其中偏微分方程式(27)的初始和边界条件为

(28)

2.2 稳定性分析

利用Gupta等研究中的方法，可以推导出上述线性PDE-ODE方程的解析解，然而该解析法只适用于低维的钻柱模型(仅含有两个自由度)，无法推广到多自由度钻柱模型的稳定性分析。为此，本文提出了分析线性PDE-ODE稳定性的数值方法，该方法不受系统自由度的限制，可以容易地拓展到多自由度系统的稳定性分析。

(29)

(30)

(31)

(32)

(33)

将式(30)中的第二式代入常微分方程式(26)后可得

(34)

(35)

则上述耦合的常微分方程组写成矩阵矢量的表示形式为

(36)

式中，系数矩阵A和B中的元素均是由无量纲参数组成的已知量。由此，利用谱方法和伽辽金方法可以将线性PDE-ODE的稳定性分析转化为分析线性常微分方程式(36)的稳定性。线性常微分方程式(36)的稳定性可以通过其特征值实部的正负来判断，当所有特征值的实部均为负时，线性常微分方程式(36)的零解渐进稳定，则原非线性PDE-ODE方程的稳态解也渐进稳定；相反，若至少有一个特征值的实部为正，则线性常微分方程式(36)的零解不稳定，同样原非线性PDE-ODE方程的稳态解也不稳定；若除含有零实部特征值外，其余特征值的实部均为负，则线性方程组的零解及原非线性PDE-ODE方程的稳态解均为临界稳定。

3 稳定性结果验证及阻尼影响

3.1 稳定性结果验证

RGD模型存在两种不稳定性机制：快不稳定性机制和慢不稳定性机制。当地面转速ω0大于临界转速ωc时，系统的稳定性由扭转振动的不稳定极点主导，表现为慢不稳定性，钻头轴向和扭转速度为准周期运动，需要经过很长时间才能发展为自激的黏滑振动；而当ω0<ωc时，系统的不稳定性由轴向振动的不稳定极点主导，表现为快不稳定性，即钻头在短时间内即出现自激的黏滑振动。临界转速ωc的近似计算公式为

(37)