1. 项目缘起从“黑箱”到“白箱”的抽油系统认知跃迁在石油开采的现场有杆抽油系统俗称“磕头机”是陆地油田最常见的一道风景。这套机械系统看似结构简单但其内部动力学行为却异常复杂。在我早期参与油田数字化改造项目时面对一个反复出现“光杆断脱”故障的井传统的经验诊断方法——听声音、看电流图、凭老师傅的手感——显得力不从心。我们耗费了大量时间排查更换了悬绳器、调整了冲程冲次问题却依旧间歇性出现直接导致了近一个月的非计划性停产和可观的产量损失。这次经历让我深刻意识到仅仅依靠外部现象和经验去理解这套由地面驱动设备、抽油杆柱、井下泵及油管组成的“长链条”系统如同隔靴搔痒。我们必须建立一套能够精确描述其内在物理规律的数学模型将“黑箱”变为“白箱”。而MATLAB以其强大的数值计算、矩阵处理、控制系统仿真及丰富的可视化工具箱成为了实现这一目标的首选利器。本次分享的“有杆抽油系统的数学建模及诊断”正是基于这样的工程背景旨在通过严谨的数学推导和仿真实践构建一套可用于系统性能分析、故障预测与智能诊断的数字化工具链。无论你是从事油气田开发的研究人员、现场工程师还是对机电系统建模与控制感兴趣的学生这套方法都能为你提供一个从理论到实践的完整视角。2. 有杆抽油系统动力学模型的核心一维波动方程及其离散化有杆抽油系统的核心物理过程是能量通过细长的抽油杆柱长度可达数千米从地面传递到井下泵。抽油杆柱在交变载荷作用下产生的纵向振动是分析一切问题如应力、位移、载荷的基础。描述这一振动的最经典模型是一维波动方程。2.1 一维波动方程的推导与物理意义我们首先将抽油杆柱视为一个均质、连续的弹性杆。根据牛顿第二定律和胡克定律可以推导出描述杆柱纵向振动的一维波动方程∂²u(x,t)/∂t² a² * ∂²u(x,t)/∂x² - c * ∂u(x,t)/∂t其中u(x, t)是距离井口x处、时间t时杆柱的位移米。a sqrt(E/ρ)是应力波在杆柱中的传播速度米/秒E是杆材的弹性模量帕斯卡ρ是杆材密度千克/立方米。对于钢杆a通常在5000 m/s左右。c是粘滞阻尼系数1/秒用于表征杆柱在油液中的运动阻尼。注意这个方程是建模的基石。它告诉我们杆柱上任意一点的加速度方程左边由两部分决定一是该点附近杆段的应力差导致的恢复力a² * ∂²u/∂x²二是运动过程中受到的粘滞阻尼力-c * ∂u/∂t。理解每一项的物理意义是后续设置正确边界条件和参数的关键。2.2 有限差分法将连续方程转化为MATLAB可解的代数方程波动方程是偏微分方程PDE我们需要将其离散化才能用计算机求解。有限差分法FDM是最直观和常用的方法。其核心思想是用离散的网格点来逼近连续的空间和时间域。空间离散将长度为L的杆柱等分为N段空间步长Δx L/N。离散点编号为i 0, 1, 2, ..., N其中i0对应井口光杆iN对应井下泵处。时间离散将时间等分时间步长为Δt。时间层编号为k 0, 1, 2, ...。差分格式我们采用中心差分格式来近似方程中的偏导数。这是精度和稳定性的一种平衡选择。时间二阶导数∂²u/∂t² ≈ (u(i, k1) - 2u(i, k) u(i, k-1)) / (Δt)²空间二阶导数∂²u/∂x² ≈ (u(i1, k) - 2u(i, k) u(i-1, k)) / (Δx)²时间一阶导数阻尼项∂u/∂t ≈ (u(i, k1) - u(i, k-1)) / (2Δt)将上述差分格式代入波动方程经过整理我们可以得到关于未来时刻k1层位移u(i, k1)的显式递推公式u(i, k1) [ (a²Δt²/Δx²) * (u(i1, k) u(i-1, k)) (2 - 2*a²Δt²/Δx² - cΔt) * u(i, k) (cΔt/2 - 1) * u(i, k-1) ] / (1 cΔt/2)这个公式是MATLAB迭代计算的核心。它表明杆上某一点下一时刻的位移可以由当前时刻及前一时刻该点及其相邻两点的位移计算出来。2.3 稳定性条件CFL数的关键作用显式差分格式是有条件稳定的。稳定性条件由著名的CFLCourant-Friedrichs-Lewy条件决定a * Δt / Δx ≤ 1这个条件有深刻的物理意义数值计算中信息传播的速度Δx/Δt必须大于或等于物理世界中波传播的实际速度a。否则物理信息还未传到计算就已经进行了必然导致结果发散数值爆炸。在实际编程中我们通常取CFL a * Δt / Δx 0.8 ~ 0.95在保证稳定的前提下获得较大的时间步长以提高计算效率。在MATLAB中我们需要先根据杆长L、波速a和期望的空间分辨率N来确定Δx然后根据CFL条件反推出最大允许的Δt。这是仿真能否成功运行的第一个关键检查点。我曾因为初期忽略了此条件设置了过大的Δt导致计算几十步后位移值就溢出成Inf或NaN排查了许久才定位到这个基础问题。3. 边界条件与载荷连接模型与现实的桥梁离散化的方程描述了杆柱内部的动力学而系统的行为最终由边界井口和泵端驱动。边界条件的设置直接决定了模型的输入是仿真是否“真实”的另一个生命线。3.1 地面边界条件驴头运动规律在井口i0位移u(0, t)是已知的它由抽油机的几何结构和电机运动决定。最常见的简化是将其视为简谐运动u(0, t) 0.5 * S * [1 - cos(2π * t / T)]其中S是冲程米T是冲程周期秒。在MATLAB实现时我们直接在每一个时间步k将上述公式计算出的值赋给u(0, k)。更精确的模型可以考虑游梁式抽油机的实际连杆机构运动用三角函数组合来描述但简谐运动在多数诊断场景下已足够。3.2 泵端边界条件流体载荷与泵阀动力学泵端iN的边界条件最为复杂因为它涉及到杆柱与井下流体、泵阀的相互作用。这是建模的难点和重点。我们通常采用力平衡条件。在泵处杆柱的力由胡克定律计算F_rod EA * ∂u/∂x |_{xL}必须与作用在泵上的流体载荷F_fluid平衡。F_rod F_fluid流体载荷F_fluid由泵筒内的压力决定F_fluid A_plunger * (P_below - P_above)其中A_plunger是柱塞截面积。P_below是泵吸入口压力与地层流压相关。P_above是泵排出口压力与油管液柱压力、井口回压相关。而泵阀的开启与关闭会动态改变P_above和P_below。一个相对实用且计算量可接受的简化模型是上冲程柱塞上行。固定阀打开游动阀关闭。P_below≈ 吸入口压力较低P_above为上一冲程排入油管的液柱压力。此时F_fluid方向向上是最大载荷。下冲程柱塞下行。固定阀关闭游动阀打开。P_above≈P_below≈ 吸入口压力。此时F_fluid很小甚至为负向下是最小载荷。在MATLAB代码中我们需要在每个时间步判断泵阀状态从而动态计算F_fluid然后利用力平衡条件推导出泵端 (iN) 的位移u(N, k1)。这通常需要将泵端的力平衡方程与内部节点的差分方程联立求解或通过迭代逼近。实操心得泵阀逻辑的实现是代码中最易出错的部分。一个常见的错误是阀的开关判断逻辑与杆柱运动速度∂u/∂t的符号未正确关联导致载荷曲线出现违背物理规律的震荡。我的调试方法是单独输出一个冲程周期内泵端的位移、速度、计算出的上下压力及阀状态绘制在一张图上人工检查其逻辑时序是否正确。这个过程虽然繁琐但一劳永逸。4. MATLAB仿真实现从零构建诊断原型有了理论框架我们开始在MATLAB中将其实现。我将以构建一个基础仿真原型为例详解关键步骤和代码片段。4.1 参数初始化与网格生成% 1. 系统基本参数 L 1000; % 杆柱长度米 E 2.1e11; % 钢弹性模量Pa rho 7850; % 钢密度kg/m^3 a sqrt(E/rho); % 波速m/s c 0.5; % 经验阻尼系数1/s S 3.0; % 冲程m T 10; % 周期s f 1/T; % 频率Hz omega 2*pi*f; % 角频率rad/s % 2. 离散化参数 N 100; % 空间分段数 dx L / N; % 空间步长 CFL 0.9; % CFL数取0.9保证稳定 dt CFL * dx / a; % 由稳定性条件确定时间步长 total_time 5 * T; % 仿真总时间模拟5个冲程 M round(total_time / dt); % 总时间步数 % 3. 初始化位移场 (空间点数 N1, 时间步数 M1) u zeros(N1, M1); % u(i, k), i:空间索引 k:时间索引4.2 核心迭代循环与边界处理% 4. 初始条件假设从静止开始杆柱处于拉伸平衡位置 % u(:, 1) 和 u(:, 2) 可设为相同值或根据静载荷计算一个初始位移分布 % 这里简单设为0 u(:, 1) 0; u(:, 2) 0; % 预计算系数提高循环效率 r a * dt / dx; coeff1 r^2; coeff2 2 - 2*r^2 - c*dt; coeff3 c*dt/2 - 1; denom 1 c*dt/2; % 5. 时间迭代主循环 for k 2:M % k代表当前已知的“过去”层我们要计算 k1 层 % 5.1 更新地面边界 (i0) - 简谐运动 t (k-1) * dt; % 当前时间 u(1, k1) 0.5 * S * (1 - cos(omega * t)); % 注意MATLAB索引从1开始 % 5.2 更新内部节点 (i2 到 iN) for i 2:N u(i, k1) ( coeff1*(u(i1, k) u(i-1, k)) ... coeff2*u(i, k) ... coeff3*u(i, k-1) ) / denom; end % 5.3 更新泵端边界 (iN1) - 这里需要嵌入泵阀模型 % 这是一个简化示例假设泵端力已知或通过简单关系给出 % 实际中这里需要调用一个独立的函数来计算泵端载荷和位移 [u(N1, k1), pump_force(k)] updatePumpBoundary(u(N, k), u(N1, k), u(N1, k-1), dt, dx, E, A_rod, ...); endupdatePumpBoundary函数是工程实现的核心它封装了第3.2节所述的泵阀逻辑和力平衡计算。由于其复杂性它本身可能包含条件判断、压力计算和方程求解。4.3 关键结果的可视化诊断的“眼睛”仿真完成后我们必须将数据转化为直观的图形这是诊断分析的基础。% 6. 结果可视化 time_axis (0:M) * dt; position_axis (0:N) * dx; % 6.1 地面示功图载荷-位移图 % 计算光杆载荷F_surface E * A_rod * (u(2,:) - u(1,:)) / dx 一阶差分近似应力 F_surface E * A_rod * (u(2, :) - u(1, :)) / dx; figure; plot(u(1, :), F_surface / 1000, b-, LineWidth, 1.5); % 载荷单位化为kN xlabel(光杆位移 (m)); ylabel(光杆载荷 (kN)); title(仿真地面示功图); grid on; % 6.2 泵端示功图 % 计算泵端载荷已在循环中存储于 pump_force 数组 figure; plot(u(N1, :), pump_force / 1000, r-, LineWidth, 1.5); xlabel(泵位移 (m)); ylabel(泵载荷 (kN)); title(仿真泵端示功图); grid on; % 6.3 杆柱应力/位移分布动画 (可选但非常直观) figure; for k 1:50:M1 % 每隔50帧显示一帧 plot(position_axis, u(:, k), b-o); xlabel(杆柱位置 (m)); ylabel(位移 (m)); title([杆柱位移分布 t , num2str((k-1)*dt, %.2f), s]); ylim([-0.5*S, 1.5*S]); grid on; drawnow; pause(0.05); end将仿真得到的地面示功图与现场实测的示功图进行对比是验证模型准确性的第一步也是故障诊断的起点。5. 基于模型仿真的故障诊断方法论建立了准确的模型后我们就可以将其用作一个“数字孪生体”。诊断的基本思路是对比分析。5.1 故障注入与特征库构建健康的系统模型在标准参数下运行会生成一个“基准”示功图。当某种故障发生时如泵漏失、杆柱断脱、油管锚失效等系统的某个或某几个参数会发生变化。我们在仿真模型中人为地修改这些参数模拟故障状态并记录下对应的示功图变化。故障类型模型参数变化仿真示功图形状特征泵漏失降低泵的充满系数或使游动阀/固定阀不能完全密封。上冲程载荷上升缓慢或无法达到最大值下冲程载荷下降缓慢图形变“瘦”面积减小做功减少。抽油杆断脱将杆柱在断点处的截面面积A_rod设为0或直接修改波动方程在该点的连接条件。光杆载荷大幅减小且波动剧烈示功图呈一条窄带或杂乱无章泵端位移几乎为零。气体影响增加泵内流体的压缩性或降低有效排量。上冲程初期载荷上升缓慢气体压缩段图形左上角变“圆”出现“气锁”时图形变得很小。油管锚失效改变泵端边界条件允许油管与杆柱发生相对位移。示功图整体发生偏转倾斜图形形状发生畸变上下冲程载荷线不平行。在MATLAB中我们可以编写一个脚本循环遍历多种故障模式及其严重程度批量运行仿真自动提取每个故障示功图的特征参数如最大载荷、最小载荷、图形面积、斜率变化点等构建一个“故障模式-特征向量”数据库。5.2 实测数据对比与智能诊断当获得一口井的实测地面示功图后诊断流程如下数据预处理使用MATLAB的smoothdata函数滤除噪声用findpeaks函数定位冲程起点和终点进行归一化对齐。特征提取从预处理后的实测图中提取与仿真特征库相同的一组特征参数。相似度匹配将实测特征向量与故障特征库中的每一个向量进行相似度计算。简单的方法可以用欧氏距离复杂一些可以用动态时间规整DTW来比较整个图形序列MATLAB有dtw函数。故障识别与置信度评估找到相似度最高的故障模式并计算其匹配置信度。可以设定一个阈值低于阈值则判定为“未知故障”或“复合故障”。% 简化版诊断匹配示例 % measured_features: 从实测图提取的1xm特征向量 % fault_library: n x m 矩阵n种故障每种有m个特征 % fault_labels: n x 1 细胞数组存储故障名称 distances zeros(size(fault_library, 1), 1); for i 1:size(fault_library, 1) distances(i) norm(measured_features - fault_library(i, :)); end [min_dist, idx] min(distances); confidence 1 / (1 min_dist); % 一个简单的置信度计算 if min_dist threshold diagnosed_fault fault_labels{idx}; fprintf(诊断结果%s 置信度%.2f%%\n, diagnosed_fault, confidence*100); else fprintf(未匹配到已知故障模式。\n); end5.3 诊断系统的工程化考量将上述原型发展为可用的诊断系统还需要考虑以下工程问题模型参数标定模型的准确性依赖于输入参数如阻尼系数c、泵效、流体性质等。我们需要利用井史数据或某次健康的实测数据通过反演算法如MATLAB的fminsearch,lsqnonlin来校准这些参数使仿真示功图与健康状态实测图最佳匹配。这是一个“模型调参”的过程。实时性要求对于在线监测仿真速度必须快。可以考虑使用更高效的数值方法如特征线法或采用降阶模型ROM。MATLAB Coder可以将核心仿真代码转换为C/C显著提升速度。不确定性处理实测数据噪声大模型本身也有简化误差。诊断结果应给出概率或置信区间而不是绝对断言。可以引入贝叶斯推理框架来融合多源不确定信息。在我参与的一个项目中我们利用校准后的模型成功预警了一口井的杆柱偏磨加剧趋势。模型仿真显示在特定位置杆柱的横向振动加剧对应应力集中。我们建议调整了扶正器位置避免了后续可能发生的杆断事故。这种从“事后诊断”到“事前预警”的转变正是数学建模价值的最高体现。6. 超越基础模型进阶与集成应用基础的一维波动方程模型已经能解决大部分问题但对于更复杂的场景我们需要对其进行扩展。6.1 三维振动与屈曲分析在斜井、水平井或存在严重摩阻的情况下杆柱不仅纵向振动还会发生横向振动和螺旋屈曲。这就需要建立三维的杆柱动力学模型控制方程将扩展到一组耦合的偏微分方程。虽然计算量剧增但MATLAB的PDE Toolbox或自己编写有限元代码可以应对。这类模型可以精确分析杆管偏磨位置为优化扶正器布置提供定量依据。6.2 与地下渗流模型耦合抽油系统的效率最终受制于地层供液能力。我们可以将井筒杆柱模型与描述地层流体向井筒流动的渗流模型如使用达西定律耦合起来。泵吸入口的压力P_below不再是一个固定值或简单假设而是由地层渗透率、流体粘度、生产压差等参数动态计算得出。这种耦合仿真能更真实地模拟“供液不足”等生产动态问题。6.3 集成到SCADA与数字孪生平台最终的愿景是将这个MATLAB诊断模型集成到油田的SCADA数据采集与监控系统或数字孪生平台中。实现路径可以是MATLAB Production Server将诊断算法部署为Web服务。SCADA系统实时将示功图数据发送到该服务并接收返回的诊断结果和健康评分。编译独立应用使用MATLAB Compiler将整个诊断程序打包成独立的可执行文件.exe或动态链接库.dll供其他平台如C#、Java开发的平台调用。云化部署结合MATLAB Online或第三方云平台实现多口井数据的集中化、规模化诊断分析。从一行行数学公式和MATLAB代码到一个能够精准模拟物理现实、并用于指导生产的诊断系统这个过程充满了挑战但也极具成就感。它要求我们不仅懂数值计算和编程更要深入理解抽油系统背后的每一个物理细节。每一次模型的调优、每一次诊断的成功验证都是对“机理数据”这一工业智能化路径的坚实印证。当你看到自己构建的模型输出的示功图与现场高精度传感器录制的图形高度重合时那种跨越虚拟与现实的连接感正是工程建模的魅力所在。