简介针对钻井作业中的钻孔顺序优化问题这套MATLAB代码提供了一组基于蚁群算法的改进实现主要面向石油工程、工业工程或运筹优化方向的开发者和学生。核心思路是将旅行商问题TSP映射到井位遍历路径规划通过最小化总行程成本提升作业效率压缩包共六个文件全部为.m脚本包含算法主程序、成本计算、信息素更新、初始放置等模块整体体积仅3KB结构精简。代码覆盖蚁群算法关键环节并带有对应原始版本的修改标记可能针对钻井约束做了适应性调整便于读者研究算法变体与工程落地方案。对于想了解智能优化算法实际应用的读者这是一份轻量且完整的参考案例已有111人学习适合作为算法研发或课程设计的起点。1. 为什么钻井工程需要一套自己的 MATLAB 代码钻井工程师做钻具组合校核或水力参数设计时往往面临两种尴尬商业软件能出漂亮报告但改一个地层系数就要重新填表单Excel 表格够用但一涉及井深网格、迭代求解或蒙特卡洛模拟公式拉扯到几十列就会把人逼疯。drilling-code 这种自建 MATLAB 计算包正是为了在“可追溯”和“够灵活”之间找平衡。它把井眼轨迹、钻具组合、流变参数、机械钻速回归这些钻井工程常算的东西拆成独立函数输入结构体输出结构体主脚本只负责组织数据。这样无论是算一口定向井的摩阻扭矩还是批跑 200 组钻压转速组合都能在一段可读的 MATLAB 代码里完成。这套方案的适用人群很明确石油工程专业的研究生、钻井工艺工程师、以及做井眼轨迹控制或钻柱动力学的仿真人员。你不会在一个 drilling-code 文件里找到花哨的图形界面但你能找到所有计算步骤的中间量、单位换算和注释。本文不依赖任何商业工具箱来打通主干计算只有流变拟合和回归部分用到 MATLAB 优化工具箱用来替代手写梯度迭代。如果你已经见过大型钻井软件的后台逻辑你会更容易理解这套代码为什么把“数据结构”放在第一步。2. drilling-code 的第一层井身轨迹与钻具组合数据结构钻井计算最常出错的点不是公式本身而是数据放不对地方。井深、井斜、方位、钻具外径、线重、地层压力这些参数来自不同报表单位也不一致。drilling-code 的常见做法是先用结构体struct把原始数据一次性装载再用函数把它变成后续计算需要的节点数组或表格。这样每个计算模块都不必关心原始文件长什么样只认结构体字段。2.1 用结构体和表格组织钻井数据先定义一个最小可用的井眼轨迹结构体字段名尽量和钻井报告中的英文缩写一致% 定义井眼轨迹原始数据MD为测量井深Incl为井斜Azi为方位 MD [0 300 600 900 1200 1500 1800 2100]; % 单位m Incl [0 2 5 8 12 20 30 40]; % 单位deg Azi [0 0 10 15 25 30 35 40]; % 单位deg traj struct(MD, MD, Incl, Incl, Azi, Azi);这种装载方式的好处是后续增加测斜计算方法时不需要改函数签名只要在traj里多放一个method字段。实际项目中我一般会在主脚本开头加一段assert检查确保MD、Incl、Azi三个列向量长度一致且MD严格递增避免数据拼接时出现错位。从结构体到计算坐标需要用测斜计算方法换算垂深TVD和北东坐标。最简单的平均角法公式是% 平均角法用相邻两个测点的平均井斜和平均方位递推 TVD zeros(size(MD)); N zeros(size(MD)); E zeros(size(MD)); for i 2:length(MD) dMD MD(i) - MD(i-1); inc1 deg2rad(Incl(i-1)); inc2 deg2rad(Incl(i)); azi1 deg2rad(Azi(i-1)); azi2 deg2rad(Azi(i)); inc_avg (inc1 inc2) / 2; azi_avg (azi1 azi2) / 2; TVD(i) TVD(i-1) dMD * cos(inc_avg); N(i) N(i-1) dMD * sin(inc_avg) * cos(azi_avg); E(i) E(i-1) dMD * sin(inc_avg) * sin(azi_avg); end这里每个测段的垂深增量是测量斜深增量乘以平均井斜余弦北向增量和东向增量再按平均方位分解。参数inc_avg和azi_avg用弧度制因为 MATLAB 的三角函数默认输入弧度这一点最容易和 Excel 公式混用。平均角法在井眼曲率不大时精度足够但造斜段或扭方位段通常需要最小曲率法。2.2 井眼轨迹的最小曲率修正最小曲率法假设相邻测点之间的井眼轴线是一段圆弧比平均角法更接近真实钻进曲线。核心是先把狗腿角beta算出来再求一个比例因子RF% 最小曲率法计算承接上面的变量 beta zeros(size(MD)); RF ones(size(MD)); TVD2 zeros(size(MD)); N2 zeros(size(MD)); E2 zeros(size(MD)); for i 2:length(MD) dMD MD(i) - MD(i-1); alpha1 deg2rad(Incl(i-1)); alpha2 deg2rad(Incl(i)); phi1 deg2rad(Azi(i-1)); phi2 deg2rad(Azi(i)); % 狗腿角余弦 cosBeta cos(alpha1)*cos(alpha2) ... sin(alpha1)*sin(alpha2)*cos(phi2-phi1); cosBeta max(min(cosBeta, 1), -1); % 防止数值溢出 beta(i) acos(cosBeta); % 比例因子beta很小时直接取1 if beta(i) 1e-9 RF(i) 1; else RF(i) 2 / beta(i) * tan(beta(i)/2); end TVD2(i) TVD2(i-1) dMD * (cos(alpha1) cos(alpha2)) / 2 * RF(i); N2(i) N2(i-1) dMD * (sin(alpha1)*cos(phi1) sin(alpha2)*cos(phi2)) / 2 * RF(i); E2(i) E2(i-1) dMD * (sin(alpha1)*sin(phi1) sin(alpha2)*sin(phi2)) / 2 * RF(i); end这里cosBeta被强制限制在[-1,1]区间因为浮点运算可能让acos的输入超出定义域这是测斜计算里最常见的 NaN 来源。RF在狗腿角接近零时趋近于1这时最小曲率法退化为正切法所以我们用阈值1e-9避免除零。实际用 drilling-code 时建议把两种方法都实现在报告里对比垂深差异用来评估仪器测斜误差对井身质量的影响。2.3 钻具组合离散化从一段管柱到有限元节点钻柱力学计算不能把整个组合看作一个刚体需要把 BHA 和钻杆按实际长度离散成若干单元。下段代码展示如何用结构体描述一段钻铤% 钻铤段参数外径/内径/线重/弹性模量 bha struct(OD, 0.165, ... % 外径 165.1mm ID, 0.071, ... % 内径 71.1mm w, 1.96e3, ... % 线重 196 kg/m E, 2.07e11); % 弹性模量 207 GPa % 将钻铤分成20个梁单元离散节点坐标 L 100; % 钻铤长度m Nele 20; % 单元数 z linspace(0, L, Nele1); % 节点位置从井底算起 % 计算截面积和惯性矩 A pi/4 * (bha.OD^2 - bha.ID^2); I_sec pi/64 * (bha.OD^4 - bha.ID^4);为什么把 BHA 离散成节点因为后面的摩阻扭矩递推需要知道每个节点上的单位长度重量和侧向力。z向量从 0 到 100 表示井底到井口方向注意后面钻压的符号约定要以井底为起点。A和I_sec是后续计算轴向应力、弯曲应力和屈曲特征值的公共参数把它们放在一个结构体里随时可以传给子函数。在这一层数据结构和基本几何已经就绪。很多新手在drilling-code上栽跟头是因为没有区分“测量井深”和“垂深”所以写代码前先把轨迹换算出 TVD 和 N/E 坐标后面做地层分层对比就方便多了。下一章会在这套离散节点上建立钻柱力学计算模块。3. 核心计算模块轴向力、扭矩与屈曲临界载荷钻柱力学是钻井工程里公认最难算准的部分因为井眼不规则、钻具变形、地层摩阻都是现场真实数据不可能用一个纯理论公式包打天下。drilling-code 的价值在于把“模型选择”和“计算过程”分开你可以在同一个数据结构里切换软绳模型、梁模型或屈曲判别式然后对比结果差异。3.1 浮力、摩阻与轴向载荷递推从井底向上逐段计算轴向力是大多数扭矩摩阻软件的底层做法。先算浮力系数再对每个节点计算重力分量和侧向力function [F, T, N] axial_force_torque(BHA, traj, WOB, rpm, mu) % BHA : 钻具组合结构体包含节点坐标z、外径OD、内径ID、线重w % traj : 井眼轨迹结构体包含MD、Incl弧度、Azi弧度 % WOB : 井底钻压N拉伸为正压缩为负 % rpm : 转速RPM用于扭矩计算此处仅作记录 % mu : 井眼摩阻系数 % rho_m : 钻井液密度kg/m^3 % rho_s : 钢材密度kg/m^3 rho_m BHA.rho_mud; rho_s 7850; BF 1 - rho_m / rho_s; % 浮力系数 n length(BHA.z); F zeros(n, 1); T zeros(n, 1); N zeros(n, 1); F(end) WOB; % 井底位置压缩钻压取负 for i n-1:-1:1 dL BHA.z(i1) - BHA.z(i); alpha deg2rad(traj.Incl(i)); % 该节点井斜 phi deg2rad(traj.Azi(i)); % 重力分量和侧向力简化计算 dW BHA.w * BF * dL; % 浮重单位N F_axial F(i1); % 当前轴向力 % 侧向力来自管柱重力的径向分量 轴向力方向变化简化模型 incl_change deg2rad(traj.Incl(i1)) - alpha; azi_change deg2rad(traj.Azi(i1)) - phi; F_side sqrt((dW * sin(alpha) F_axial * incl_change)^2 ... (F_axial * sin(alpha) * azi_change)^2); % 轴向力递推上提 mu*F_side下放 - mu*F_side % 这里按正常钻进下放计算 F(i) F(i1) dW * cos(alpha) - mu * F_side; N(i) F_side; % 扭矩递推 r_contact BHA.OD / 2; % 接触半径假设为外径一半 T(i) T(i1) mu * F_side * r_contact; end end这个函数的递推方向是“从井底到井口”因为井底钻压已知向上逐段加上浮重和摩阻就得到井口大钩载荷。注意这里的WOB的符号约定压缩取负所以井底位置的F(end)是负值。F_side假设侧向力由重力径向分量和井眼曲率造成没有考虑钻柱弯曲刚度引起的额外接触力如果想要更精确的结果下一步可以用梁单元求解横向位移再把接触力迭代回来。但作为工程快速评估这个简化模型已是很多批量敏感性分析的基础。3.2 用参数化函数替换屈曲公式钻柱屈曲判别式有很多版本不同钻井手册的系数略有差别。drilling-code 里最好的做法是把临界载荷计算写成独立函数参数化井斜角、径向间隙和抗弯刚度这样公式改动只动一处function [Fcr_sin, Fcr_hel] buckling_load(EI, w_b, alpha, r_c) % 计算钻柱在斜井眼中的屈曲临界载荷 % EI : 抗弯刚度 (N·m^2) % w_b : 单位长度浮重 (N/m) % alpha : 井斜角 (rad) % r_c : 管柱与井壁的径向间隙 (m) % 参考工程近似正弦屈曲临界载荷和螺旋屈曲临界载荷 % 系数2和4来自 Lubinski 经典公式在斜直井眼中的推广 Fcr_sin 2 * sqrt(EI * w_b * sin(alpha)) / r_c; Fcr_hel 4 * sqrt(EI * w_b * sin(alpha)) / r_c; if Fcr_hel Fcr_sin Fcr_hel Fcr_sin; % 防止奇异输入 end end调用时传入每个节点的局部井斜角alpha、该节点的浮重线重和外径半径差。径向间隙r_c的计算方式如果是裸眼段等于井眼半径减管柱外径半径如果在套管内则用套管内径半径减管柱外径半径。注意单位保持国际制否则临界载荷会差好几个数量级。这个函数返回两个值后面画安全窗口时把实际轴向力画成一条曲线再把正弦屈曲临界载荷和螺旋屈曲临界载荷画成包络线。3.3 输出井眼安全操作窗口有了轴向力曲线和屈曲临界载荷就能画出一张钻井工程常用的“钻压安全窗口”图% 计算并绘制安全窗口 alpha_deg traj.Incl; % 井斜角度 alpha_rad deg2rad(alpha_deg); Fcr_s zeros(size(z)); Fcr_h zeros(size(z)); for i 1:length(z) [Fcr_s(i), Fcr_h(i)] buckling_load(bha.E * I_sec, ... bha.w * 9.81 * BF, alpha_rad(i), r_c); end figure; plot(F / 1000, z, b, LineWidth, 1.5); hold on; plot(Fcr_s / 1000, z, r--, LineWidth, 1.2); plot(Fcr_h / 1000, z, g--, LineWidth, 1.2); legend(实际轴向力, 正弦屈曲临界, 螺旋屈曲临界); xlabel(轴向力 (kN)); ylabel(节点位置 (m)); set(gca, YDir, reverse); % 井深从上往下 title(钻柱轴向力与屈曲安全窗口);绘制时注意把纵轴方向反转让井口在顶部、井底在底部。实际轴向力曲线若位于临界载荷包络线左侧说明该深度段没有屈曲风险若穿越临界线就说明该段已经进入正弦或螺旋屈曲需要调整钻压或改用更重的下部钻具组合。这里的r_c需要根据实际井眼直径计算我在代码里用全局变量但更正式的做法是把它放进BHA结构体里一并传递。4. 水力参数与机械钻速回归让 drilling-code 真正能“算”钻柱力学负责“这段钻具受不受力”水力参数负责“能不能把岩屑带出来”机械钻速回归负责“钻得快不快”。这三个问题在钻井排班里往往互相耦合所以 drilling-code 的常见组织结构里会把水力计算和 ROP 回归放在同一层最终输出一份“操作对比表”。4.1 流变模型拟合从实测数据到 n、K钻井液的流变参数通常用旋转粘度计测出剪切速率和剪切应力。幂律模型的数学形式是 τ K * γ^n两边取对数就能用线性回归拟合% 实测数据剪切速率γ (1/s) 和剪切应力τ (Pa) gamma [5 10 20 50 100 200 300 600]; tau [3.2 4.1 5.6 8.9 14.0 23.0 31.0 49.0]; % 对数线性化log(τ) log(K) n*log(γ) log_gamma log(gamma); log_tau log(tau); p polyfit(log_gamma, log_tau, 1); n p(1); K exp(p(2)); fprintf(n %.3f, K %.3f Pa·s^n\n, n, K);polyfit做最小二乘拟合返回的p(1)是斜率对应幂律指数np(2)是截距取指数后得到稠度系数K。这个方法的优点是速度快、可解释性强但缺点是拟合对实测数据的测量误差比较敏感而且当体系存在屈服应力时线性拟合会失真。此时我一般会改用 MATLAB 优化工具箱里的lsqnonlin直接以τ - (K*γ^n τ0)为目标函数拟合三参数 Herschel-Bulkley 模型% 三参数流变模型τ τ0 K*γ^n model (x) x(1) x(2) * gamma.^x(3) - tau; x0 [1, 0.5, 0.8]; % 初值 [τ0, K, n] x lsqnonlin(model, x0, [0 0 0.1], [10 10 1.5]); tau0 x(1); K2 x(2); n2 x(3);lsqnonlin支持边界约束初值的选择对收敛影响很大。工程上一般把屈服应力初值设为几个 PaK 初值设为 0.5 左右n 初值设为 0.8 左右。如果拟合结果出现负屈服应力说明该钻井液可能不需要三参数模型这时可以退回幂律模型。4.2 钻头喷嘴压降与水力功率计算水力参数里最常反复计算的是钻头压降和钻头水力功率。压降公式为ΔP_bit (ρ * Q²) / (2 * Cd² * A_nozzle²)其中 ρ 是钻井液密度单位 kg/m³Q 是排量m³/sCd 是喷嘴流量系数通常取 0.95A_nozzle 是喷嘴总过流面积m²。用 MATLAB 写成一个函数输入排量 L/s 和喷嘴当量直径 mmfunction [dP_bit, HHP] bit_hydraulics(rho_kg_m3, q_lps, d_nozzle_mm) % rho_kg_m3 : 钻井液密度 (kg/m^3) % q_lps : 排量 (L/s) % d_nozzle_mm : 喷嘴当量直径 (mm) % 输出钻头压降 (MPa) 和钻头水功率 (kW) Q q_lps / 1000; % 转换为 m^3/s A pi/4 * (d_nozzle_mm/1000)^2; % 单喷嘴面积m^2 Cd 0.95; dP_bit (rho_kg_m3 * Q^2) / (2 * Cd^2 * A^2); % Pa dP_bit dP_bit / 1e6; % 转换为 MPa HHP dP_bit * 1e6 * Q / 1000; % 水功率 ΔP * Q (W) kW end注意喷嘴当量直径和喷嘴个数的关系。如果现场用的是三喷嘴当量直径是三个喷嘴直径的平方和开根号实际上当量流量面积是 π/4 * sum(d_i²)所以传入d_nozzle_mm时要么传总当量面积对应的直径要么传单个直径但内部乘喷嘴数量。为了避免歧义我的代码里只计算一个喷嘴的压降实际使用时会写一个包装函数把多喷嘴的面积合并后传入。这个细节很容易在单位换算里出错建议在函数内写assert(q_lps 0)防呆。4.3 ROP 回归与钻压-转速交互分析机械钻速回归常用的工程模型是 ROP a * W^b * N^c其中 W 是钻压单位通常为 kN 或 tN 是转速RPMa、b、c 是待拟合参数。取对数后是线性回归问题% 现场实测数据钻压(kN), 转速(RPM), 实际ROP(m/h) W [40 50 60 70 80 90]; N [60 60 80 80 100 100]; ROP_obs [8.2 10.5 13.1 14.8 17.2 19.0]; % 对数化后用 regress 做多元线性回归 X [ones(size(W)), log(W), log(N)]; y log(ROP_obs); coef regress(y, X); a exp(coef(1)); b coef(2); c coef(3); fprintf(ROP %.3f * W^%.3f * N^%.3f\n, a, b, c);回归后可以用nlpredci计算预测区间但更常见的是把拟合出的模型画成三维曲面再用contour等高线找出目标 ROP 下的钻压-转速组合。这个过程很容易忽略一点X里log(W)和log(N)的尺度差异会让回归矩阵病态所以工程上往往先把 W 和 N 归一化到 [0,1] 区间再回归得到系数后再反推原始参数。drilling-code 中我一般会保留归一化前后的系数方便现场直接代入。5. 用批量脚本做敏感性分析并验证代码当单个计算模块都跑通以后drilling-code 的优势才真正体现出来你可以不点一次鼠标就把几十组钻压、转速、排量和喷嘴组合全部算完并输出成表格。这一章讲批量执行的标准化写法以及我用来保障结果可信的一个回归断言技巧。5.1 参数扫描的标准化脚本敏感性分析不需要写复杂循环嵌套。把要扫描的参数定义为向量再用meshgrid生成网格然后对每个网格点调用前面的水力或 ROP 函数% 扫描对象钻压 W (kN) 和转速 N (RPM) W_vec [50 60 70 80]; N_vec [60 80 100]; [W_grid, N_grid] meshgrid(W_vec, N_vec); ROP_surf zeros(size(W_grid)); dP_surf zeros(size(W_grid)); for i 1:numel(W_grid) [ROP_surf(i), dP_surf(i)] simulate_drilling(W_grid(i), N_grid(i), ... rho_mud, d_nozzle, regr_coef); endsimulate_drilling可以是前面几个模块的封装输入钻压转速先算钻柱屈曲风险再算水力压降最后用回归模型算 ROP。封装时注意输入输出参数名要一致代码维护会更省心。meshgrid生成的网格是矩阵但numel遍历时顺序是按列来的如果你之后要surf绘图记得保持矩阵形状。5.2 并行与表格输出如果扫描点超过几千个可以用parfor代替for但要注意循环体内的函数都必须能序列化而且不要声明全局变量。常见做法是写一个独立函数把参数作为结构体传入% 用并行循环批量计算 results cell(numel(W_grid), 1); parfor i 1:numel(W_grid) res simulate_drilling(W_grid(i), N_grid(i), params); results{i} res; end % 汇总成表 T_out table(W_grid(:), N_grid(:), cellfun((x) x.ROP, results), ... cellfun((x) x.HHP, results), ... VariableNames, {W_kN, N_RPM, ROP_mh, HHP_kW}); writetable(T_out, sensitivity_results.csv);parfor的随机性问题如果simulate_drilling里有随机数生成需要用RandStream保证每个 worker 的随机序列可重复。上述代码中cellfun提取结果里的字段最终生成的table可以直接写 CSV也可以再用gather做可视化。实际工作中我建议第一次先不用parfor跑一遍小规模验证再开并行否则一旦远端parfor报错调试成本远高于计算省下的时间。5.3 一个检查中间结果的技巧批量计算最容易出现的问题是参数组合里某组数据让迭代不收敛、结果变成 NaN但主脚本不一定报错。我习惯在每个循环迭代结束后立即做一次断言检查把可疑结果当场拦下来if ~isfinite(res.ROP) || res.ROP 0 warning(异常结果: W%d kN, N%d RPM, W_grid(i), N_grid(i)); continue; % 跳过不让异常值污染统计结果 end更进一步的验证是能量或量纲检查。比如钻头水功率一定小于泵压乘以排量再比如井底钻压绝对值大于井口钩载时说明摩阻递推方向反了。把这些断言写进simulate_drilling函数内部就能在批量运行的第一轮发现模型或参数的错误而不是等到画图时才发现曲线形状完全不对。把这个断言放进你的主脚本以后改参数、改公式一眼就能看出哪个结果已经不再可信。本文还有配套的精品资源点击获取