钻柱粘滑振动模型：从机理到仿真与控制

简介面向钻井工程与钻柱动力学研究人员这是一份聚焦钻柱粘滑stick-slip现象建模与预防控制的MATLAB/Simulink资源包。内容包含线性与非线性钻柱模型、MPC控制实现及相关附录覆盖从模型搭建到控制策略验证的完整链条。压缩包共22个文件以.m脚本、.mat数据与.mdl模型为主整体仅90KB适合通过仿真快速复现粘滑特征并测试不同防滑控制方案。已有157人学习适用于掌握一定钻井力学基础和Simulink操作的研究生、工程师。借助该资源可系统理解钻柱非线性动力学特性对照模型代码与数据学习模型预测控制在粘滑抑制中的实际用法附录还提供了补充说明与参考便于进一步扩展研究。1. 拿到“Drillstringmodel_stick-slip”文件时先搞清楚是什么在振动第一次看到这个文件命名我第一反应是——这又是一个从油田现场传回来的“问题样本”。347977_ATTACHMENT01是内部编号Drillstringmodel指的是钻井井柱钻柱的动力学模型后面跟的两个几乎重复的词stick-slip和stickslip都指向同一个工程难题粘滑振动。如果你在钻井工程、井下工具研发、钻井仿真或者岩石力学相关的岗位上待过对这个词一定不陌生。它是钻井作业里最让人头疼的异常工况之一轻则磨损钻头、缩短钻具寿命重则导致井下工具疲劳断裂、螺杆钻具损坏甚至引发井控风险。这份文件的核心价值在于它把一个看起来只存在于现场现象层面的“粘滑振动”转化成了一套可以计算、可以复现、可以拿来做控制仿真的数学模型。打开压缩包后典型的内容应该包括钻柱集中质量动力学模型、扭转振动微分方程、摩擦扭矩模型、求解程序或仿真工程文件。对钻井工程师来说这份模型是分析井下动力学的钥匙对做控制算法的人来讲它是验证“如何消除粘滑振动”的试验床。这篇文章不打算写教科书而是以我实际打开和处理这类模型的操作为主线把模型背后的物理逻辑、代码实现、参数标定、仿真输出和现场印证这五件事讲透。不管你是刚接触钻井仿真的研究生还是想给自己的控制算法找一套工程对象的算法工程师读完应该能直接上手。2. 粘滑振动模型到底在建模什么2.1 一粘一滑到底是什么意思要把模型读懂先得在脑子里建立起物理画面。现场最常见的画面是转盘或者顶驱以恒定转速驱动钻柱顶端旋转但钻头接触岩石后受到巨大的摩擦阻力当阻力矩超过系统能够传递的扭矩时钻头会突然卡住不再转动。此时顶驱还在继续旋转钻柱像一个巨大的扭簧一样不断扭曲储能等到储存的弹性势能超过静摩擦阻力矩的临界值钻头瞬间释放、高速旋转把存储的扭转变形能一口气放完然后又因为阻力重新卡住。这个“卡住—储能—释放—再卡住”的循环就是粘滑振动。从动力学角度看粘滑振动本质上是自激振荡。系统没有外部周期激励而是由“静摩擦与动摩擦的切换”这个非线性环节自己激发了周期运动。这也是为什么模型里一定会包含非线性摩擦扭矩这个模块——没有它模型就是一根普通的线性扭转弹簧永远不可能自发产生振荡。理解这点是第一步也是最重要的一步因为之后你调模型、做控制、分析结果全都是在跟这个非线性环节打交道。2.2 模型的数学骨架扭转摆与集中质量工程上最常用的钻柱粘滑振动模型是图1所示的“二自由度集中质量扭转摆模型”。所谓集中质量是把整根钻柱等效成两个惯性元件和一个弹性元件顶端转动惯量代表顶驱和转盘旋转部分底端转动惯量代表钻铤、钻头和钻柱下部的大质量段中间的扭转刚度则代表整个钻柱的弹性变形能力。模型的核心控制方程可以写成下面两组二阶微分方程顶端转动惯量 $J_t$ 的动力学方程$$J_t \ddot{\theta}_t C_t(\dot{\theta}_t - \Omega_0) K(\theta_t - \theta_b) T_d$$底端转动惯量 $J_b$ 的动力学方程$$J_b \ddot{\theta}_b C_b \dot{\theta}b - K(\theta_t - \theta_b) -T{friction}(\dot{\theta}_b, F_w, \dots)$$其中 $\theta_t$ 是顶驱旋转角$\theta_b$ 是钻头旋转角$K$ 是整个钻柱的等效扭转刚度$C_t$ 和 $C_b$ 分别为顶驱和底部等效阻尼$\Omega_0$ 是顶驱设定转速$T_{friction}$ 就是之前强调过的非线性摩擦扭矩。物理图像上钻柱就是一根大型的“扭转弹簧”顶驱像一只手不断拧弹簧钻头侧则被岩石摩擦紧紧压住。用生活类比就是你用螺丝刀拧一颗严重生锈的螺丝螺丝不动螺丝刀杆越拧越紧积蓄弹性变形能等到力量够了螺丝突然“啪”地转了大半圈然后又被卡住这就是钻头粘滑振动的缩小版。这个方程体系看着简单但它已经把粘滑振动最关键的两个元素捕捉到了第一弹性储能过程由刚度项 $K(\theta_t - \theta_b)$ 描述第二摩擦与滑动的交替过程由函数 $T_{friction}(\dot{\theta}_b)$ 描述。后者的非线性特性直接决定系统是否会进入极限环振荡。2.3 摩擦扭矩模型选型库仑模型够用吗摩擦模型是粘滑模型里最敏感的一环。最简单的做法是采用库仑摩擦模型当钻头转速接近零时使用静摩擦系数 $T_s$当钻头转动时使用动摩擦系数 $T_c$。为了避免零速附近的数值奇异工程上通常把“近似为零”的速度区间做一个线性过渡。更精细的做法是采用斯托里贝克Stribeck摩擦曲线该曲线描述了摩擦扭矩随相对速度从静摩擦到动摩擦的下降过程并且通常带一个指数衰减项$$T_f(\dot{\theta}_b) \left[T_c (T_s - T_c) e^{-|\dot{\theta}_b / v_s|^\delta}\right] \cdot \operatorname{sgn}(\dot{\theta}_b)$$其中 $v_s$ 是过渡速度$\delta$ 是速度衰减指数通常取 1 或 2。从我的实践经验来看如果只是做机理分析库仑模型就够了但如果要做控制验证或者和现场录井曲线对比必须上斯托里贝克模型因为它的 $T_s$ 到 $T_c$ 下降过程会影响粘滑振动的极限环幅值和频率你拿库仑模型算出来的振荡周期可能和实测对不上。实操心得文件包里的模型如果只给了一组固定摩擦参数我建议你别直接用先拿现场历史钻压、扭矩和转速数据做一次简单反演标定。常见做法是把零转速附近的摩擦扭矩设为静摩擦段取起钻、下钻时“倒刹”工况的扭矩极值做参考把正常钻进时的平均扭矩作为动摩擦段参考。整个标定过程不需要什么高深算法Excel 也能做但参数不准的话后面仿真和控制设计都会跟着跑偏。3. 仿真实操怎么把模型跑起来3.1 工程参数录入模型能不能跑出现场那种“锯齿形转速曲线”绝大部分取决于初始参数给得对不对。以下是典型的某 215.9mm 井眼、3000 米井深钻柱参数从这一模型的工程实例中总结参数符号数值单位说明顶驱转动惯量$J_t$1580kg·m²包含顶驱转子、减速机构底部转动惯量$J_b$285kg·m²钻铤加钻头组合的等效值钻柱扭转刚度$K$420N·m/rad全井钻柱等效扭转刚度顶驱阻尼$C_t$90N·m·s/rad电机与传动系统的等效阻尼底部阻尼$C_b$30N·m·s/rad钻头与地层接触的等效粘性阻尼静摩擦扭矩$T_s$8.5kN·m低转速时岩石—钻头摩擦动摩擦扭矩$T_c$4.2kN·m正常滑动钻进时的平均扭矩顶驱设定转速$\Omega_0$90r/min约 9.42 rad/s要是文件里没有给全这些参数可以用估算公式补齐。转动惯量部分钻铤和钻柱按空心圆杆计算$J \frac{1}{2} \rho L \pi (R_o^4 - R_i^4)$其中 $\rho$ 是钢材密度 7850 kg/m³$L$ 是段长$R_o$ 和 $R_i$ 分别为外径和内径。刚度部分按 $K \frac{G I_p}{L}$ 折算其中 $G$ 是钢材剪切模量约 79.3 GPa$I_p$ 是极惯性矩 $\frac{\pi}{32}(D_o^4 - D_i^4)$。这些公式手册里都有实际工作中75%的模型跑不出发散振荡原因都不是“模型错了”而是参数给得离谱。3.2 程序实现四阶龙格库塔时的注意点模型方程建立后接下来要用数值方法求解。我习惯直接拿 MATLAB 或 Python 写一个四阶龙格库塔求解器。状态变量可以设置为 $[\theta_t, \omega_t, \theta_b, \omega_b]^T$其中 $\omega$ 为角速度。这样方程组就可以改写为状态空间形式$$\frac{d}{dt} \begin{bmatrix} \theta_t \ \omega_t \ \theta_b \ \omega_b \end{bmatrix} \begin{bmatrix} \omega_t \ \left(T_d - C_t(\omega_t - \Omega_0) - K(\theta_t - \theta_b)\right)/J_t \ \omega_b \ \left(-C_b \omega_b K(\theta_t - \theta_b) - T_f(\omega_b)\right)/J_b \end{bmatrix}$$注意摩擦项的处理必须放在每次求解器内部、每个子步都更新一次而不是在循环外面算好。如果你把 $T_f$ 当成常数算完整个时间序列那模型必然失真。另一个容易犯的坑是 $\operatorname{sgn}$ 函数在零速处的不连续建议用 tanh 做平滑否则高频抖动会让龙格库塔法的步长被迫降得很小严重影响求解效率。求解时间建议至少仿真 60 秒步长取 0.00050.001 秒。粘滑振动的一个周期通常在 25 秒量级60 秒足够抓到 10 个以上完整周期便于统计分析振动幅值。如果取大步长你可能会看到阻尼反而变大的假象因为数值耗散把高频成分吞掉了。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 模型参数 J_t 1580.0 # 顶驱转动惯量 kg·m^2 J_b 285.0 # 底部转动惯量 kg·m^2 K 420.0 # 扭转刚度 N·m/rad C_t 90.0 # 顶驱阻尼 N·m·s/rad C_b 30.0 # 底部阻尼 N·m·s/rad T_s 8500.0 # 静摩擦力矩 N·m T_c 4200.0 # 动摩擦力矩 N·m v_s 0.05 # Stribeck 过渡速度 rad/s Omega0 90 * np.pi / 30 # 顶驱设定转速 rad/s def T_friction(wb): # 平滑库仑 Stribeck 衰减 return (T_c (T_s - T_c) * np.exp(-abs(wb) / v_s)) * np.tanh(wb / 0.01) def rhs(t, y): th_t, om_t, th_b, om_b y dth_t om_t dom_t ( -C_t*(om_t - Omega0) - K*(th_t - th_b) ) / J_t dth_b om_b dom_b ( -C_b*om_b K*(th_t - th_b) - T_friction(om_b) ) / J_b return [dth_t, dom_t, dth_b, dom_b] y0 [0.0, Omega0, 0.0, Omega0] sol solve_ivp(rhs, [0, 60], y0, methodRK45, max_step0.001, dense_outputTrue) t sol.t theta_t sol.y[0] * 180 / np.pi # 转为角度显示 theta_b sol.y[2] * 180 / np.pi omega_t sol.y[1] * 30 / np.pi # 转为 rpm 显示 omega_b sol.y[3] * 30 / np.pi plt.figure(figsize(10, 4)) plt.plot(t, omega_b, lw0.8, label钻头转速 (rpm)) plt.plot(t, omega_t, lw0.8, label顶驱转速 (rpm)) plt.xlabel(时间/s); plt.ylabel(转速/rpm) plt.legend(); plt.grid(True) plt.show()3.3 结果输出与现象复现跑完 60 秒仿真正常的粘滑振动模型输出应该是这样钻头转速曲线在 15~50 秒区间出现周期性剧烈波动转速在零附近停留相当长的时间段这就是“粘”随后瞬间冲到 200 rpm 以上这就是“滑”顶驱转速则相对平稳在设定值 90 rpm 附近小幅波动。二者的“阶梯状”差异是判断模型是否成功的直观标志。我在第一次复现时只调了 20 秒结果曲线只出现了一次低频振荡当时差点以为模型错了。排查后发现是摩擦模型里过渡速度 $v_s$ 设得太大导致干摩擦切换变成了近似线性阻尼系统直接回到稳定状态。把 $v_s$ 从 0.5 降到 0.05 后粘滑振荡立刻出现了。这类问题用仿真软件的参数扫描功能很快就能定位关键在于你要先知道“该扫哪个参数”。4. 参数计算与调试细节4.1 扭转刚度、阻尼等效值的意义钻柱模型里扭转刚度 $K$ 的取值直接决定了振荡频率。实际钻柱是连续体刚性分布很复杂简化成单根等效弹簧时$K$ 可以通过现场实测的扭转自由振动频率来反标定。具体做法是在顶驱上施加一个阶跃扭矩测量井口扭矩响应的自然频率再代入公式 $f \frac{1}{2\pi}\sqrt{K(1/J_t 1/J_b)}$ 反算 $K$。如果手头没有现场实验数据就按 $K GI_p/L$ 计算注意这里 $L$ 要用整根钻柱长度而不是某个局部段长否则会差一到两个数量级。底部阻尼 $C_b$ 则比 $K$ 更难标定通常只能给一个经验范围20100 N·m·s/rad。它代表钻头与地层切削、钻柱与井壁碰撞的等效能量耗散。调参时你会发现$C_b$ 太小则模型轻微扰动就发散振荡太大则粘滑振荡被完全抑制所以 $C_b$ 本质上控制着粘滑振荡出现的“强度门槛”。仿真时建议作为扫描参数一次性做 50200 N·m·s/rad 区间的批处理观察极限环幅值变化。4.2 摩擦模型参数与钻压的耦合关系摩擦扭矩不是独立于工况的常量。现场钻压WOB是控制摩擦扭矩的直接因素钻压越大钻头压入岩石越深接触面积越大$T_s$ 和 $T_c$ 也越高。近似线性关系可以用 $T_s \approx \mu_s R_b W_f$、$T_c \approx \mu_c R_b W_f$ 描述其中 $R_b$ 是钻头等效半径$W_f$ 是钻压$\mu_s$ 与 $\mu_c$ 分别为静、动摩擦系数。实际工程文件的摩擦扭矩如果是一个常数那么它对应的其实是某一固定钻压工况做敏感性分析时不能只改钻压而不同步改摩擦。实操心得我在一次项目里连续换了三种不同型号的PDC钻头结果粘滑振动的剧烈程度差异巨大。问题根源就在牙齿形状改变了对地层的摩擦特性。所以在模型里做“改变钻头型号”仿真时最有效的改法不是去改 $J_b$其实底部转动惯量变化不大而是去改摩擦扭矩和它的速度衰减指数 $\delta$。这个细节很少有人写进论文里但对工程复现特别重要。4.3 边界条件与仿真时长求解时还有一个容易犯错的地方初始条件。如果让初始钻头转速等于顶驱转速即系统从完全同步状态启动有些参数组合下模型需要很长时间才能“自发”进入粘滑振荡因为非线性系统需要扰动才能偏离同步状态。我的习惯是初始钻头转速加一个 5% 的小扰动或者直接给初始时刻一个非零钻头角位移差这样能缩短仿真“预热”时间更快看到极限环。另外仿真时长建议至少覆盖 5 个粘滑周期瞬态阶段最开始 2~5 秒要舍弃不要纳入统计分析否则会低估粘滑振荡的幅值。5. 模型跑通后能做什么抑制粘滑的控制方向5.1 从模型到控制斜率就是方向模型的价值不只是复现故障更在于设计和验证抑制方案。工程上最常见的措施有如下几条提高顶驱设定转速 $\Omega_0$增加顶驱转速会减小钻头端的静、动摩擦扭矩切换造成的相对速度变化范围在一定窗口内可以弱化粘滑。代价是机械钻速可能下降、井斜控制压力增大。降低钻压WOB降低摩擦扭矩水平直接缩小“储能—释放”的能量差。这是现场最直观的应急操作但会牺牲机械钻速。施加主动阻尼控制在顶驱电机输出上叠加与钻头速度偏差成比例的附加扭矩虚拟增大 $C_t$。控制律通常很简单相当于一个 PID 控制器的 D 项但参数整定必须基于模型仿真做预验证。从仿真结果看第三种方案是效果最好且最值得投入的方向。控制律的基本思路是$T_d T_{set} - K_d(\omega_t - \omega_b)$这里 $K_d$ 是主动阻尼系数。在原本的模型中相当于把顶驱阻尼从 $C_t$ 提升到 $C_t K_d$。但要注意的是这个反馈量用的是钻头与顶驱的转速差实际现场中钻头转速很难直接获得通常需要用观测器从井口扭矩和顶驱转速中在线估计。5.2 把控制律加进仿真里验证以 Python 代码为例在原有 filter 函数里把控制扭矩从固定阀值改为反馈形式def rhs_control(t, y): th_t, om_t, th_b, om_b y K_d 300.0 # 主动阻尼系数 T_ctrl K_d * (om_t - om_b) # 反馈控制扭矩 dth_t om_t dom_t ( -C_t*(om_t - Omega0) T_ctrl - K*(th_t - th_b) ) / J_t dth_b om_b dom_b ( -C_b*om_b K*(th_t - th_b) - T_friction(om_b) ) / J_b return [dth_t, dom_t, dth_b, dom_b]从波形上看加入反馈控制后钻头转速曲线的周期性剧烈波动会明显削弱转速波动幅值从 0~250 rpm 缩小到 70~110 rpm 区间系统从极限环振荡回归稳定。实际现场中类似的软扭矩控制系统Soft Torque Rotary System正是通过修改顶驱电机的扭矩输出将粘滑振动引起的转速波动降低 80% 以上。这套控制策略已经在很多钻井平台上得到验证几乎成了高端顶驱的标配功能。5.3 不要只盯仿真曲线控制鲁棒性仿真好用不等于现场好用因为现场参数存在严重不确定性。比如井深变化导致 $K$ 改变钻头磨损导致摩擦特性变化泥浆密度改变导致阻尼变化。设计控制方案时至少要扫一组参数做鲁棒性分析比如把 $K$ 在 300~550 N·m/rad 范围内改变或把 $T_s$ 在 7~10 kN·m 范围内改变看控制效果是否仍然稳定。如果只在一个工作点调好参数到现场极有可能失效。6. 常见问题与排查技巧实录6.1 常见问题速查表现象可能原因排查方法仿真转速无振荡永远平稳摩擦模型过渡速度 $v_s$ 过大变成近似线性阻尼将 $v_s$ 调到 0.01~0.1 区间重新试振荡有但频率过高和现场不符扭转刚度 $K$ 虚高复核钻柱长度和外径确认折算公式里的 $L$ 是否取整长钻头转速长时间卡零无法滑脱$T_s$ 大于系统能传递的最大扭矩增大顶驱主动扭矩或降低静摩擦系数曲线初期杂乱后期勉强稳定初始条件离极限环太远瞬态干扰查看 5 秒以后的数据或从某个已收敛状态作为初值粘滑振荡幅度远大于现场录井$C_b$ 设置偏小能量耗散不足增加底部阻尼到 80~150 区间测试摩擦切换处数值发散$\operatorname{sgn}$ 函数不连续导致改用 tanh 平滑或采用低阶隐式求解器6.2 几条亲测有效的避坑心得第一不要迷信原始文件里的参数。这类共享文件里的参数常常来自某一口特定井的工况直接套用到别的井大概率对不上。最稳妥的方式是拿到文件后先运行基线工况然后从现场历史曲线中挑一段“正常钻进”和一段“粘滑振荡”的数据做简单的均值和幅值对比把参数按比例缩放修正再继续后面的控制设计。第二记录所有参数版本。粘滑模型参数多、版本迭代快前一天的参数可能只是摩擦系数微调了 0.1结果却完全改变系统是否收敛。我的习惯是给每次仿真配置一个 JSON 文件文件名带日期和井号。别提什么高大上的模型管理系统就这个土办法帮我省了不知道多少小时的重复排查时间。第三先用简化模型验证逻辑再上完整模型。如果你拿到的是很大很完整的模型文件比如几十个模块的 Simulink 工程不要一开始就全部启用。先把底部钻具组合、钻头、井壁接触等子系统全部屏蔽只留顶驱、钻柱刚度、底部惯量三个基本模块跑通基础版粘滑振荡再逐步恢复复杂环节。这样一旦波形异常你能立刻判断是哪个模块引入的问题而不是在一堆嵌套子系统里大海捞针。最后的实操建议至于这个文件到底该怎么进一步利用我的建议是从“反向验证”开始拿现场实测的钻头转速或扭矩数据跟模型输出做同一时间尺度的对齐先别管幅值是否完全一致重点看振荡频率、相位关系和波形形态是否匹配。匹配度达到六七成就说明模型的物理架构是对的可以大胆基于它做控制设计如果完全对不上先别动控制器回过去检查参数标定。这个文件是一个很好的起点但“能用”和“好用”之间的距离全靠参数标定和反复仿真迭代来填。本文还有配套的精品资源点击获取