“把 if 往上提,for 往下放!”

本文经授权转自公众号优快云(ID:优快云news)

作者 | Alex Kladov,翻译 | 苏宓

很多程序员写代码的时候都会遇到这种情况:一个判断条件到底该放在函数里面还是外面?循环里是不是可以加个 if?这些看起来无关紧要的小选择,实际会影响代码的清晰度、性能。最近,一位热爱简单代码和编程语言的程序员 Alex Kladov 在博客上分享了两条简单但特别实用的经验法则,引发不少开发者的注意。其表示:“把 if 往上提,把 for 往下放(Push Ifs Up, Push Fors Down)”。这会对提升代码质量有哪些作用,我们不妨通过本文一起来看看。

关于两条相关经验法则的简短说明。

1、把 If 条件往上推

如果函数内部有 if 条件判断,可以考虑是否能把它移到调用者那里:

// 好的写法fn frobnicate(walrus: Walrus) {    ...}// 不好的写法fn frobnicate(walrus: Option<Walrus>) {  let walrus = match walrus {    Some(it) => it,    None => return,  };  ...}

就像上面举的例子一样,函数里经常会有“前置条件”的检查——比如你传进来的参数是不是有效、是不是空之类的。很多时候,程序员会在函数内部去检查这些情况,如果不符合就什么都不做。但更好的做法,其实是把这些判断提前,在调用这个函数之前就处理好。这样一来,函数本身就可以更专注地完成它的核心工作,不用管太多额外的事。

为什么要这么做?有两个原因:

首先,如果每个函数都自己检查一次条件,那整个程序里可能会有很多重复判断,既费事又容易出错。把判断提前处理好,可以让代码整体更简洁高效。

其次,if 语句(也就是“如果……就……”)本身就是让程序变复杂的东西。判断多了,程序流程就绕了,Bug 也容易藏在里面。所以与其到处乱插 if,不如集中放在一个地方处理——比如放在主函数里统一决定怎么走,然后把具体的执行任务交给其他小函数去做,这样结构更清晰,出错的概率也更低。

而且,如果所有的判断都集中在一个地方,程序员在看代码的时候也能更容易发现有没有重复判断,或者根本不会被执行的“死代码”。比如:

fn f() {  if foo && bar {    if foo {    } else {    }  }}fn g() {  if foo && bar {    h()  }}fn h() {  if foo {  } else {  }}

对于 f() 来说,更容易看出无用分支,而分散在 g() 和 h() 的组合里不容易发现。

这里还有一个相关的模式,叫做“dissolving enum”重构。有时候代码会长这样:

enum E {  Foo(i32),  Bar(String),}fn main() {  let e = f();  g(e)}fn f() -> E {  if condition {    E::Foo(x)  } else {    E::Bar(y)  }}fn g(e: E) {  match e {    E::Foo(x) => foo(x),    E::Bar(y) => bar(y)  }}

这里有两处分支,通过往上提,可以发现它们其实是同一个条件被重复了三次(第三次变成了数据结构):

fn main() {  if condition {    foo(x)  } else {    bar(y)  }}
2、把 For 循环往下推

这条建议来自“面向数据”的编程思想。简单说就是:少就是少,多就是多。现实中的程序,大多数时候都不是在处理一个东西,而是一大批东西。特别是在那些对性能要求高的关键部分(也叫“热点路径”),几乎总是在同时处理很多个对象——正是因为数量多,这部分代码才会“热”。

所以,一个常见的好习惯是:从一开始就把“批量处理”当成主要方式,而不是只想着一个个地处理。如果真的需要处理单个对象,那就把它当成“批量里的一个特例”来对待。这样做,不光能让程序结构更清晰,还能为后面提升性能打下基础。比如:

// 好的写法frobnicate_batch(walruses)// 不好的写法for walrus in walruses {  frobnicate(walrus)}

这样做最主要的好处是性能,尤其在极端情况下非常明显。

如果有一整批数据可以一起处理,就能摊销启动成本,还能灵活调整处理顺序。实际上,你甚至不必按特定顺序处理实体,可以先对所有实体的某个字段统一处理,再处理其他字段,利用向量化或结构化数组技巧。

一个有趣的例子是基于 FFT 的多项式乘法:在一堆点上同时计算多项式,比一个点一个点计算要快得多!

如果 if 往上推,for 循环往下推,两条规则还能组合使用:

// 好的写法if condition {  for walrus in walruses {    walrus.frobnicate()  }} else {  for walrus in walruses {    walrus.transmogrify()  }}// 不好的写法for walrus in walruses {  if condition {    walrus.frobnicate()  } else {    walrus.transmogrify()  }}

好的写法避免了在循环内部反复判断条件,减少了循环里的分支,提高了性能,甚至可能解锁向量化优化。

这个模式既适用于微观层面,也适用于宏观架构层面,比如 TigerBeetle 的数据平面就是对批量对象操作,以摊销控制平面决策的成本。

虽然性能是“for 循环往下推”的主要动机,但它也能提升表达力。比如早期很成功的 jQuery 就是针对元素集合的操作,抽象向量空间的语言思维往往比逐个坐标方程更清晰。

最后用一句话总结:把 if 条件往上推,把 for 循环往下推!

本文转自公众号“优快云”,ID:优快云news

---END---

% 修正:轴向力从井口(最大)向井底(最小)递减 % 基于Cerberus/Maxwell等先进软件的电缆力学模型 % 作者:Advanced Cable Mechanics Analysis System % 版本:3.4 Professional (Axial Force Corrected) tic; clc; clear; close all; fprintf('========================================================================================================\n'); fprintf(' 测井电缆高精度力学分析系统 v3.4 - 轴向力修正版\n'); fprintf(' 基于Cerberus/Maxwell算法核心 + 最新工业标准API RP 9B\n'); fprintf('========================================================================================================\n'); fprintf('模型特性:\n'); fprintf(' • 压缩点自动识别与标记 (Compression Point Detection)\n'); fprintf(' • 弱点受力精确分析 (Weak Point Analysis)\n'); fprintf(' • 多层电缆结构建模 (Multi-layer Cable Model)\n'); fprintf(' • 实时屈曲状态监测 (Real-time Buckling Monitor)\n'); fprintf(' • 轴向力正确计算 (Axial Force Corrected)\n'); fprintf(' • Dawson-Paslay模型屈曲判断 (Dawson-Paslay Buckling Model)\n'); fprintf('========================================================================================================\n\n'); %% ===================== 1. 数据准备与参数定义 ===================== fprintf('[1/8] 初始化系统参数...\n'); % 检查并加载井眼轨迹数据 if ~exist('井眼轨迹.xlsx', 'file') fprintf(' → 创建示例井眼轨迹...\n'); depth_example = (0:50:5000)'; inclination_example = [zeros(20,1); linspace(0, 60, 30)'; 60*ones(51,1)]; azimuth_example = [zeros(20,1); linspace(0, 135, 30)'; 135*ones(51,1)]; well_trajectory = [depth_example, inclination_example, azimuth_example]; writematrix(well_trajectory, '井眼轨迹.xlsx'); else well_trajectory = readmatrix('井眼轨迹.xlsx'); end % ===== 工况选择 ===== WORKING_CONDITION = 2; % 1=上提, 2=下放, 3=静止 condition_names = {'上提作业', '下放作业', '静止悬挂'}; fprintf(' → 当前工况: %s\n', condition_names{WORKING_CONDITION}); % ===== 电缆结构参数(多层模型)===== cable = struct(); cable.OD = 0.0118; % 外径 11.8mm (7/16") cable.length = 5619; % 长度 m cable.breaking_strength = 8500; % 破断强度 kg % 电缆多层结构定义(从内到外) cable.layers = [ % 直径(m), 密度(kg/m³), 弹性模量(GPa), 泊松比, 屈服强度(MPa) 0.002, 8900, 110, 0.34, 250; % 铜芯导体 0.004, 1200, 3.5, 0.40, 50; % 绝缘层 0.006, 2700, 70, 0.33, 180; % 铝护套 0.008, 1400, 2.0, 0.42, 40; % 外护套 0.0118, 7850, 200, 0.30, 450; % 钢丝铠装 ]; % 计算等效参数 [cable.density, cable.E, cable.nu, cable.yield_strength] = calculate_composite_properties(cable); cable.weight_in_air = cable.density * pi * (cable.OD/2)^2 * 9.81; % N/m % ===== 弱点参数定义 ===== weak_points = struct(); weak_points.locations = [500, 1500, 2500, 3500, 4500]; % 弱点位置 m weak_points.types = {'接头', '磨损', '腐蚀', '接头', '疲劳'}; % 弱点类型 weak_points.strength_factors = [0.85, 0.90, 0.88, 0.85, 0.82]; % 强度折减系数 weak_points.stiffness_factors = [0.95, 1.00, 0.98, 0.95, 0.92]; % 刚度折减系数 % ===== 工具串参数 ===== tool = struct(); tool.OD = 0.0432; % 外径 43.2mm (1-11/16") tool.length = 17.4; % 总长 m tool.weight = 2500; % 总重 N tool.sections = [ % 长度(m), 外径(m), 重量(N), 名称 5.0, 0.0432, 500, "密度工具"; 3.0, 0.0380, 300, "中子工具"; 2.5, 0.0432, 250, "伽马工具"; 4.0, 0.0350, 350, "电阻率"; 3.5, 0.0432, 300, "声波工具"; 2.0, 0.0300, 100, "张力短节"; ]; % ===== 井筒与流体参数 ===== wellbore = struct(); wellbore.casing_ID = 0.1778; % 7" 套管内径 m wellbore.tubing_ID = 0.0762; % 2-7/8" 油管内径 m wellbore.casing_shoe = 3000; % 套管鞋深度 m fluid = struct(); fluid.density = 1150; % kg/m³ fluid.viscosity = 35; % mPa·s 钻井液粘度 fluid.yield_point = 15; % Pa 钻井液动切力 fluid.compressibility = 4.5e-10; % 1/Pa % ===== 摩擦参数(Stribeck模型)===== friction = struct(); friction.static = 0.3; % 静摩擦系数 friction.kinetic = 0.2; % 动摩擦系数 friction.viscous = 0.001; % 粘性摩擦系数 friction.transition_velocity = 0.01; % m/s friction.cased_hole_factor = 0.7; % 套管段系数 friction.open_hole_factor = 1.2; % 裸眼段系数 % ===== 运动参数 ===== if WORKING_CONDITION == 1 cable.velocity = 0.25; % m/s 上提速度 cable.acceleration = 0.03; % m/s² elseif WORKING_CONDITION == 2 cable.velocity = -0.25; % m/s 下放速度 cable.acceleration = -0.02; % m/s² else cable.velocity = 0; cable.acceleration = 0; end % ===== 计算参数 ===== calc = struct(); calc.step = 5; % 计算步长 m calc.max_iterations = 50; % 最大迭代次数 calc.tolerance = 1e-6; % 收敛容差 calc.output_interval = 100; % 输出间隔 m %% ===================== 2. 井眼轨迹预处理 ===================== fprintf('[2/8] 处理井眼轨迹数据...\n'); % 插值到计算步长 depth = (0:calc.step:well_trajectory(end,1))'; n_points = length(depth); % 高精度插值 inc_deg = interp1(well_trajectory(:,1), well_trajectory(:,2), depth, 'pchip'); azi_deg = interp1(well_trajectory(:,1), well_trajectory(:,3), depth, 'pchip'); % 转换为弧度 inc = inc_deg * pi/180; azi = azi_deg * pi/180; % 计算3D轨迹参数 fprintf(' → 计算曲率和挠率...\n'); [curvature, torsion, dogleg] = calculate_advanced_trajectory(depth, inc, azi); % 计算3D坐标(用于轨迹显示) [north, east, tvd] = calculate_3d_coordinates(depth, inc, azi); % 计算井壁接触特性 contact_diameter = zeros(n_points, 1); for i = 1:n_points if depth(i) <= wellbore.casing_shoe contact_diameter(i) = wellbore.casing_ID; else contact_diameter(i) = wellbore.tubing_ID; end end %% ===================== 3. 初始化结果数组 ===================== fprintf('[3/8] 初始化计算数组...\n'); % 主要力学参数 results = struct(); results.tension = zeros(n_points, 1); % 张力 N results.axial_force = zeros(n_points, 1); % 轴向力 N results.apparent_weight = zeros(n_points, 1); % 表观悬重 N results.torque = zeros(n_points, 1); % 扭矩 N·m results.bending_moment = zeros(n_points, 1); % 弯矩 N·m results.normal_force = zeros(n_points, 1); % 法向力 N results.friction_force = zeros(n_points, 1); % 摩擦力 N results.side_force = zeros(n_points, 1); % 侧向力 N results.cumulative_friction = zeros(n_points, 1); % 累积摩阻 N % 状态参数 results.buckling_state = zeros(n_points, 1); % 0=直线, 1=正弦, 2=螺旋 results.compression_flag = zeros(n_points, 1); % 压缩标志 results.stress_state = zeros(n_points, 1); % von Mises应力 Pa () results.safety_factor = zeros(n_points, 1); % 安全系数 results.fatigue_damage = zeros(n_points, 1); % 疲劳损伤 % 压缩点和零张力点 results.compression_points = []; % 压缩点位置 results.neutral_point = NaN; % 中性点位置 results.zero_tension_point = NaN; % 零张力点 % 弱点分析 results.weak_point_stress = zeros(length(weak_points.locations), 1); results.weak_point_safety = zeros(length(weak_points.locations), 1); %% ===================== 4. 高级T&D力学计算(轴向力修正版)===================== fprintf('[4/8] 执行高级T&D力学分析(轴向力修正版)...\n'); % 创建进度条 prog_bar = waitbar(0, '正在计算...', 'Name', 'T&D Analysis'); % 初始化数组 T = zeros(n_points, 1); % 张力 W = zeros(n_points, 1); % 悬重(电缆重量累积) A = zeros(n_points, 1); % 轴向力 N = zeros(n_points, 1); % 法向力 F = zeros(n_points, 1); % 摩擦力 F_friction = zeros(n_points, 1); % 单位长度摩擦力 (N/m) % 井口边界条件 T(1) = 0; % 井口张力初始值为0 W(1) = 0; % 井口悬重初始值为0 % T&D迭代求解 - 从井口向下计算张力 fprintf(' → 执行T&D迭代计算...\n'); % 预计算工具重量的浮力系数 buoyancy_tool_avg = 1 - fluid.density / 7850; % 使用平均流体密度 for iter = 1:calc.max_iterations T_old = T; % 从井口向下计算张力 for i = 1:n_points-1 depth_i = depth(i+1); ds = calc.step; % 流体密度(考虑温压效应) temp = 20 + 0.03 * depth_i; pressure = 0.01 * depth_i; % MPa rho_fluid = fluid.density * (1 + fluid.compressibility * pressure * 1e6); % 电缆浮重 buoyancy = 1 - rho_fluid / cable.density; w_cable = cable.weight_in_air * buoyancy; % N/m % ===== 法向接触力计算 ===== if dogleg(i+1) > 0 && T(i) > 0 N_dogleg = T(i) * sin(dogleg(i+1) * ds); else N_dogleg = 0; end N_gravity = w_cable * ds * sin(inc(i+1)); if abs(cable.velocity) > 0 && curvature(i+1) > 0 N_centrifugal = cable.density * pi * (cable.OD/2)^2 * cable.velocity^2 * curvature(i+1) * ds; else N_centrifugal = 0; end N(i+1) = N_dogleg + N_gravity + N_centrifugal; % ===== 摩擦系数 ===== if abs(cable.velocity) > 0 v_ratio = abs(cable.velocity) / friction.transition_velocity; mu = friction.kinetic + (friction.static - friction.kinetic) * exp(-v_ratio^2); else mu = friction.static; end if depth_i <= wellbore.casing_shoe mu = mu * friction.cased_hole_factor; else mu = mu * friction.open_hole_factor; end % ===== 摩擦力计算 ===== F_friction(i+1) = mu * N(i+1) / ds; % N/m if WORKING_CONDITION == 1 % 上提 F(i+1) = F_friction(i+1); elseif WORKING_CONDITION == 2 % 下放 F(i+1) = -F_friction(i+1); else F(i+1) = 0; end % ===== 流体阻力 ===== if abs(cable.velocity) > 0 Re = abs(cable.velocity) * cable.OD * rho_fluid / (fluid.viscosity * 1e-3); if Re < 1 Cd = 24 / Re; elseif Re < 1000 Cd = 24 / Re * (1 + 0.15 * Re^0.687); else Cd = 0.44; end F_drag = 0.5 * Cd * rho_fluid * cable.velocity * abs(cable.velocity) * cable.OD; if WORKING_CONDITION == 1 F_drag = abs(F_drag); elseif WORKING_CONDITION == 2 F_drag = -abs(F_drag); else F_drag = 0; end else F_drag = 0; end % ===== 张力递推 ===== dT_weight = w_cable * cos(inc(i+1)); dT_friction = F(i+1); dT_drag = F_drag; W(i+1) = W(i) + w_cable * ds; T(i+1) = T(i) + (dT_weight + dT_friction + dT_drag) * ds; % 修改:工具重量在迭代中的影响 - 基于有效质量和井斜 % 计算该点到井底的距离比例 depth_ratio = (depth(end) - depth(i+1)) / depth(end); % 根据工况调整工具重量的影响方式 if WORKING_CONDITION == 1 % 上提时,工具重量对所有点都有影响 tool_effect = tool.weight * buoyancy_tool_avg * depth_ratio / n_points; elseif WORKING_CONDITION == 2 % 下放时,工具重量的影响需要考虑摩擦方向 % 计算工具重量对当前点的有效影响 effective_weight = tool.weight * buoyancy_tool_avg; % 下放时,摩擦方向向下,工具重量的影响被部分抵消 tool_effect = effective_weight * depth_ratio / n_points * 0.5; % 根据物理特性调整权重 else % 静止时 tool_effect = tool.weight * buoyancy_tool_avg * depth_ratio / n_points; end T(i+1) = T(i+1) + tool_effect; % 累积摩阻 results.cumulative_friction(i+1) = results.cumulative_friction(i) + abs(F_friction(i+1)) * ds; % 限制最大张力 weak_idx = find(abs(depth_i - weak_points.locations) < calc.step/2); if ~isempty(weak_idx) strength_factor = weak_points.strength_factors(weak_idx(1)); max_T = cable.breaking_strength * 9.81 * strength_factor; else max_T = cable.breaking_strength * 9.81; end if T(i+1) > max_T T(i+1) = max_T; fprintf(' 警告: %.0fm处张力达到极限\n', depth(i+1)); end if mod(i, 100) == 0 waitbar(i / n_points, prog_bar, sprintf('计算进度: %.1f%%', i / n_points * 100)); end end if max(abs(T - T_old)) < calc.tolerance fprintf(' → T&D迭代收敛,迭代次数: %d\n', iter); break; end end % ===== 轴向力计算(修正版:从井口最大到井底最小)===== fprintf(' → 计算轴向力分布(井口最大→井底最小)...\n'); % 方法:计算每点以下的所有载荷 for i = 1:n_points % 计算该点以下的所有载荷 remaining_weight = 0; % 剩余电缆重量 remaining_friction = 0; % 剩余摩阻 % 累积该点以下的电缆重量和摩阻 for j = i:n_points-1 depth_j = depth(j); temp = 20 + 0.03 * depth_j; pressure = 0.01 * depth_j; rho_fluid = fluid.density * (1 + fluid.compressibility * pressure * 1e6); buoyancy = 1 - rho_fluid / cable.density; w_cable_j = cable.weight_in_air * buoyancy; % 电缆重量贡献 remaining_weight = remaining_weight + w_cable_j * calc.step * cos(inc(j)); % 摩阻贡献 if WORKING_CONDITION == 1 % 上提时摩阻增加轴向力 remaining_friction = remaining_friction + abs(F_friction(j)) * calc.step; elseif WORKING_CONDITION == 2 % 下放时摩阻减小轴向力 remaining_friction = remaining_friction - abs(F_friction(j)) * calc.step; end end % 加上工具重量(在最底部) if i <= n_points depth_bottom = depth(end); temp = 20 + 0.03 * depth_bottom; pressure = 0.01 * depth_bottom; rho_fluid = fluid.density * (1 + fluid.compressibility * pressure * 1e6); buoyancy_tool = 1 - rho_fluid / 7850; tool_contribution = tool.weight * buoyancy_tool; else tool_contribution = 0; end % 轴向力 = 下部所有载荷的总和 A(i) = remaining_weight + tool_contribution + remaining_friction; end % 确保轴向力从井口到井底递减 % 井口轴向力应该最大(等于总悬重) A(1) = T(end); % 井口轴向力 = 井底张力(总悬重) % 从井口向下重新计算轴向力,确保递减 depth_i = depth(1); for i = 2:n_points temp = 20 + 0.03 * depth_i; pressure = 0.01 * depth_i; rho_fluid = fluid.density * (1 + fluid.compressibility * pressure * 1e6); buoyancy = 1 - rho_fluid / cable.density; w_cable_segment = cable.weight_in_air * buoyancy; % 轴向力递减量 dA = w_cable_segment * calc.step * cos(inc(i)); if WORKING_CONDITION == 1 % 上提 dA = dA + abs(F_friction(i)) * calc.step; elseif WORKING_CONDITION == 2 % 下放 dA = dA - abs(F_friction(i)) * calc.step; end A(i) = A(i-1) - dA; end % 存储结果 results.tension = T; results.axial_force = A; % 修改表观悬重计算,现在基于正确的张力分布 % 表观悬重直接等于考虑了工具重量和摩擦方向的张力 results.apparent_weight = T; results.cable_weight = W; results.normal_force = N; results.friction_force = abs(F_friction) * calc.step; results.friction_per_meter = abs(F_friction); results.surface_weight = T(end); % 压缩点判断 for i = 2:n_points-1 if A(i) < 0 if inc_deg(i) > 60 || curvature(i) > 0.015 || dogleg(i) > 0.008 if isempty(results.compression_points) || abs(depth(i) - results.compression_points(end)) > 100 results.compression_points = [results.compression_points; depth(i)]; results.compression_flag(i) = 1; fprintf(' → 压缩点: %.0f m (轴向力: %.0f N, 井斜: %.1f°)\n', depth(i), A(i), inc_deg(i)); [buckling_state, ~] = analyze_buckling(A(i), cable, contact_diameter(i), inc(i)); results.buckling_state(i) = buckling_state; end end end end close(prog_bar); % 输出关键结果 fprintf(' → 井口张力: %.2f kN (深度0m)\n', T(1)/1000); fprintf(' → 井底张力: %.2f kN (深度%.0fm)\n', T(end)/1000, depth(end)); fprintf(' → 井口轴向力: %.2f kN (最大)\n', A(1)/1000); fprintf(' → 井底轴向力: %.2f kN (最小)\n', A(end)/1000); fprintf(' → 累积摩阻: %.2f kN\n', results.cumulative_friction(end)/1000); fprintf(' → 表观悬重: %.2f kN\n', results.surface_weight/1000); %% ===================== 5. 应力分析与安全评估 ===================== fprintf('[5/8] 执行应力分析与安全评估...\n'); A_cross = pi * (cable.OD/2)^2; I = pi * cable.OD^4 / 64; J = pi * cable.OD^4 / 32; for i = 1:n_points sigma_axial = abs(results.tension(i)) / A_cross; if curvature(i) > 0 M_bend = cable.E * I * curvature(i); sigma_bend = M_bend * (cable.OD/2) / I; else sigma_bend = 0; end if torsion(i) > 0 && results.torque(i) > 0 tau_torsion = results.torque(i) * (cable.OD/2) / J; else tau_torsion = 0; end for i = 1:n_points % 1. 计算当前点的最大允许张力(考虑弱点折减) strength_factor = 1.0; weak_idx = find(abs(depth(i) - weak_points.locations) < calc.step/2, 1); if ~isempty(weak_idx) strength_factor = weak_points.strength_factors(weak_idx); end max_tension = cable.breaking_strength * 9.81 * strength_factor; % 2. 计算实际张力(绝对值) actual_tension = abs(results.tension(i)); % 3. 安全系数 = 最大允许张力 / 实际张力 if actual_tension > 0 results.safety_factor(i) = max_tension / actual_tension; else results.safety_factor(i) = 10; % 压缩状态设为安全值 end end if WORKING_CONDITION ~= 3 C = 1e20; m = 3; stress_range = abs(results.stress_state(i)); if stress_range > 0 N_cycles = C * (stress_range)^(-m); results.fatigue_damage(i) = 1 / N_cycles; else results.fatigue_damage(i) = 0; end end end fprintf(' → 分析弱点受力...\n'); for wp = 1:length(weak_points.locations) wp_idx = find(abs(depth - weak_points.locations(wp)) < calc.step/2); if ~isempty(wp_idx) idx = wp_idx(1); stress_concentration = 1.5; switch weak_points.types{wp} case '接头' stress_concentration = 1.8; case '磨损' stress_concentration = 1.4; case '腐蚀' stress_concentration = 1.6; case '疲劳' stress_concentration = 2.0; end results.weak_point_stress(wp) = results.stress_state(idx) * stress_concentration; weak_yield = cable.yield_strength * 1e6 * weak_points.strength_factors(wp); results.weak_point_safety(wp) = weak_yield / abs(results.weak_point_stress(wp)); end end %% ===================== 6. 累积计算 ===================== fprintf('[6/8] 计算累积参数...\n'); results.cumulative_torque = zeros(n_points, 1); for i = 2:n_points if torsion(i) > 0 dM = cable.E * J * torsion(i) * calc.step; results.cumulative_torque(i) = results.cumulative_torque(i-1) + dM; else results.cumulative_torque(i) = results.cumulative_torque(i-1); end end results.cumulative_fatigue = cumsum(results.fatigue_damage); if abs(cable.velocity) > 0 results.power_loss = results.cumulative_friction * abs(cable.velocity) / 1000; max_power = max(results.power_loss); fprintf(' → 最大功率损失: %.2f kW\n', max_power); else results.power_loss = zeros(n_points, 1); end [max_stress, max_stress_idx] = max(results.stress_state); fprintf(' → 最大von Mises应力: %.1f MPa @ %.0f m\n', max_stress/1e6, depth(max_stress_idx)); %% ===================== 7. 高级可视化 ===================== fprintf('[7/8] 生成分析图表...\n'); fig = figure('Position', [20, 20, 1900, 1000], 'Name', 'Cable Analysis - Axial Force Corrected', 'Color', 'w'); % --- 子图1:张力与轴向力分布(修正版)--- subplot(4,4,[1,5]); hold on; plot(results.tension/1000, depth, 'b-', 'LineWidth', 2.5, 'DisplayName', '张力'); plot(results.axial_force/1000, depth, 'g-', 'LineWidth', 2.5, 'DisplayName', '轴向力'); xline(0, 'k-', 'LineWidth', 1.5, 'DisplayName', '零力线'); % 初始化屈曲状态统计计数器 count_compression = 0; % 压缩未屈曲计数 count_sin_buckling = 0; % 正弦屈曲计数 count_hel_buckling = 0; % 螺旋屈曲计数 % 找出所有压缩区域 compression_zones = find(results.axial_force < 0); if ~isempty(compression_zones) % 提前计算所有压缩点的屈曲状态 buckling_states = zeros(size(compression_zones)); for i = 1:length(compression_zones) idx = compression_zones(i); [buckling_states(i), ~] = analyze_buckling(results.axial_force(idx), cable, contact_diameter(idx), inc(idx)); end % 按屈曲状态分组显示 % 1. 压缩未屈曲区域(黄色) unbuckled_idx = compression_zones(buckling_states == 0); count_compression = length(unbuckled_idx); if ~isempty(unbuckled_idx) plot(results.axial_force(unbuckled_idx)/1000, depth(unbuckled_idx), 'y.', ... 'MarkerSize', 10, 'DisplayName', '压缩未屈曲'); % 填充区域 for i = 1:length(unbuckled_idx) if i == 1 || unbuckled_idx(i) ~= unbuckled_idx(i-1) + 1 start_idx = unbuckled_idx(i); end_idx = unbuckled_idx(i); else end_idx = unbuckled_idx(i); end if i == length(unbuckled_idx) || unbuckled_idx(i) ~= unbuckled_idx(i+1) - 1 fill([-100, 0, 0, -100], [depth(start_idx), depth(start_idx), depth(end_idx), depth(end_idx)], ... 'y', 'FaceAlpha', 0.2, 'EdgeColor', 'none'); end end end % 2. 正弦屈曲区域(蓝色) sin_buckled_idx = compression_zones(buckling_states == 1); count_sin_buckling = length(sin_buckled_idx); if ~isempty(sin_buckled_idx) plot(results.axial_force(sin_buckled_idx)/1000, depth(sin_buckled_idx), 'b.', ... 'MarkerSize', 10, 'DisplayName', '正弦屈曲'); % 填充区域 for i = 1:length(sin_buckled_idx) if i == 1 || sin_buckled_idx(i) ~= sin_buckled_idx(i-1) + 1 start_idx = sin_buckled_idx(i); end_idx = sin_buckled_idx(i); else end_idx = sin_buckled_idx(i); end if i == length(sin_buckled_idx) || sin_buckled_idx(i) ~= sin_buckled_idx(i+1) - 1 fill([-100, 0, 0, -100], [depth(start_idx), depth(start_idx), depth(end_idx), depth(end_idx)], ... 'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none'); end end % 添加蓝色星形标记正弦屈曲中心 if ~isempty(sin_buckled_idx) mid_idx = floor(length(sin_buckled_idx)/2); plot(results.axial_force(sin_buckled_idx(mid_idx))/1000, depth(sin_buckled_idx(mid_idx)), 'b*', 'MarkerSize', 12); end end % 3. 螺旋屈曲区域(红色) hel_buckled_idx = compression_zones(buckling_states == 2); count_hel_buckling = length(hel_buckled_idx); if ~isempty(hel_buckled_idx) plot(results.axial_force(hel_buckled_idx)/1000, depth(hel_buckled_idx), 'r.', ... 'MarkerSize', 10, 'DisplayName', '螺旋屈曲'); % 填充区域 for i = 1:length(hel_buckled_idx) if i == 1 || hel_buckled_idx(i) ~= hel_buckled_idx(i-1) + 1 start_idx = hel_buckled_idx(i); end_idx = hel_buckled_idx(i); else end_idx = hel_buckled_idx(i); end if i == length(hel_buckled_idx) || hel_buckled_idx(i) ~= hel_buckled_idx(i+1) - 1 fill([-100, 0, 0, -100], [depth(start_idx), depth(start_idx), depth(end_idx), depth(end_idx)], ... 'r', 'FaceAlpha', 0.3, 'EdgeColor', 'none'); end end % 添加红色星形标记螺旋屈曲中心 if ~isempty(hel_buckled_idx) mid_idx = floor(length(hel_buckled_idx)/2); plot(results.axial_force(hel_buckled_idx(mid_idx))/1000, depth(hel_buckled_idx(mid_idx)), 'r*', 'MarkerSize', 12); end end end breaking_force = cable.breaking_strength * 9.81; xline(breaking_force/1000, 'r--', 'LineWidth', 1.5, 'DisplayName', '破断强度'); xline(breaking_force*0.8/1000, 'y--', 'LineWidth', 1.5, 'DisplayName', '80%破断'); xlabel('力 (kN)'); ylabel('深度 (m)'); % 添加屈曲状态统计信息到标题 title(sprintf('张力与轴向力分布 (压缩未屈曲: %d点, 正弦屈曲: %d点, 螺旋屈曲: %d点)', ... count_compression, count_sin_buckling, count_hel_buckling)); grid on; set(gca, 'YDir', 'reverse'); xlim([min([results.axial_force/1000; -5]), max([results.tension/1000; results.axial_force/1000])+5]); legend('Location', 'best'); % 保留原有文本标注代码 text(results.tension(1)/1000, 0, sprintf('井口张力: %.1f kN', results.tension(1)/1000), ... 'VerticalAlignment', 'bottom', 'FontSize', 8); text(results.tension(end)/1000, depth(end), sprintf('井底张力: %.1f kN', results.tension(end)/1000), ... 'VerticalAlignment', 'top', 'FontSize', 8); text(results.axial_force(1)/1000, 0, sprintf('井口轴向力: %.1f kN', results.axial_force(1)/1000), ... 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left', 'FontSize', 8, 'Color', 'g'); text(results.axial_force(end)/1000, depth(end), sprintf('井底轴向力: %.1f kN', results.axial_force(end)/1000), ... 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left', 'FontSize', 8, 'Color', 'g'); % --- 其他子图保持不变 --- % 子图2:悬重分析 subplot(4,4,[2,6]); hold on; plot(results.apparent_weight/1000, depth, 'b-', 'LineWidth', 2, 'DisplayName', '表观悬重'); plot(results.cable_weight/1000, depth, 'g:', 'LineWidth', 1.5, 'DisplayName', '电缆重量'); xlabel('悬重 (kN)'); ylabel('深度 (m)'); title('悬重分布分析'); grid on; set(gca, 'YDir', 'reverse'); legend('Location', 'best'); % 子图3:摩阻分析 subplot(4,4,[3,7]); hold on; plot(results.cumulative_friction/1000, depth, 'r-', 'LineWidth', 2.5); xlabel('累积摩阻 (kN)'); ylabel('深度 (m)'); title('累积摩阻分布'); grid on; set(gca, 'YDir', 'reverse'); % 子图4:弱点分析 subplot(4,4,[4,8]); hold on; plot(results.safety_factor, depth, 'g-', 'LineWidth', 2); for wp = 1:length(weak_points.locations) wp_depth = weak_points.locations(wp); wp_sf = results.weak_point_safety(wp); if wp_sf < 1.8 color = 'r'; elseif wp_sf < 2.5 color = 'y'; else color = 'g'; end plot(wp_sf, wp_depth, 'o', 'MarkerSize', 10, 'MarkerFaceColor', color, 'MarkerEdgeColor', 'k'); text(wp_sf+0.1, wp_depth, sprintf('%s\nSF=%.2f', weak_points.types{wp}, wp_sf), 'FontSize', 8); end xline(1.8, 'r--', 'LineWidth', 1.5); xline(2.5, 'y--', 'LineWidth', 1.5); xlabel('安全系数'); ylabel('深度 (m)'); title('弱点安全系数分析'); grid on; set(gca, 'YDir', 'reverse'); xlim([0 5]); % 子图5-12:其他分析图表(保持不变) subplot(4,4,9); hold on; colors = {'g', 'y', 'r'}; labels = {'直线', '正弦屈曲', '螺旋屈曲'}; for bs = 0:2 idx = find(results.buckling_state == bs); if ~isempty(idx) plot(ones(size(idx))*bs, depth(idx), '.', 'Color', colors{bs+1}, 'MarkerSize', 6); end end xlabel('屈曲状态'); ylabel('深度 (m)'); title('电缆屈曲状态分布'); set(gca, 'XTick', [0 1 2], 'XTickLabel', labels); set(gca, 'YDir', 'reverse'); grid on; xlim([-0.5 2.5]); subplot(4,4,10); hold on; plot(results.stress_state/1e6, depth, 'b-', 'LineWidth', 2); xline(cable.yield_strength, 'r--', 'LineWidth', 1.5); xlabel('von Mises应力 (MPa)'); ylabel('深度 (m)'); title('应力分布'); grid on; set(gca, 'YDir', 'reverse'); subplot(4,4,11); hold on; plot(results.normal_force/1000, depth, 'b-', 'LineWidth', 2); xlabel('法向接触力 (kN)'); ylabel('深度 (m)'); title('接触力分布'); grid on; set(gca, 'YDir', 'reverse'); subplot(4,4,12); hold on; yyaxis left; plot(inc_deg, depth, 'b-', 'LineWidth', 2); ylabel('井斜角 (°)'); yyaxis right; plot(azi_deg, depth, 'r-', 'LineWidth', 1.5); ylabel('方位角 (°)'); xlabel('角度'); title('井眼轨迹'); grid on; set(gca, 'YDir', 'reverse'); subplot(4,4,[13,14]); plot3(east, north, tvd, 'b-', 'LineWidth', 2); hold on; plot3(0, 0, 0, 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g'); plot3(east(end), north(end), tvd(end), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); for wp = 1:length(weak_points.locations) wp_idx = find(abs(depth - weak_points.locations(wp)) < calc.step); if ~isempty(wp_idx) plot3(east(wp_idx(1)), north(wp_idx(1)), tvd(wp_idx(1)), 'ks', ... 'MarkerSize', 8, 'MarkerFaceColor', 'y'); end end xlabel('东坐标 (m)'); ylabel('北坐标 (m)'); zlabel('垂深 (m)'); title('3D井眼轨迹'); grid on; set(gca, 'ZDir', 'reverse'); view(45, 30); axis equal; subplot(4,4,15); if WORKING_CONDITION ~= 3 semilogy(results.cumulative_fatigue, depth, 'r-', 'LineWidth', 2); xlabel('累积疲劳损伤'); ylabel('深度 (m)'); title('疲劳损伤累积'); grid on; set(gca, 'YDir', 'reverse'); else text(0.5, 0.5, '静止状态\n无疲劳累积', 'HorizontalAlignment', 'center', 'FontSize', 12); axis off; end subplot(4,4,16); axis off; monitor_text = { sprintf('\bf实时监控数据\rm'); ''; sprintf('\bf井口参数(0m):\rm'); sprintf(' 张力: %.2f kN', results.tension(1)/1000); sprintf(' 轴向力: %.2f kN', results.axial_force(1)/1000); ''; sprintf('\bf井底参数(%.0fm):\rm', depth(end)); sprintf(' 张力: %.2f kN', results.tension(end)/1000); sprintf(' 轴向力: %.2f kN', results.axial_force(end)/1000); ''; sprintf('\bf系统状态:\rm'); sprintf(' 表观悬重: %.2f kN', results.surface_weight/1000); sprintf(' 累积摩阻: %.2f kN', results.cumulative_friction(end)/1000); ''; sprintf('\bf安全状态:\rm'); sprintf(' 最小SF: %.2f', min(results.safety_factor(results.safety_factor>0))); sprintf(' 利用率: %.1f%%', max(results.tension)/breaking_force*100); sprintf(' 压缩点: %d个', length(results.compression_points)); }; text(0.05, 0.95, monitor_text, 'VerticalAlignment', 'top', 'FontSize', 9); sgtitle(sprintf('电缆力学分析(轴向力修正版) - %s | 悬重: %.1f kN | 井口轴向力: %.1f kN | 井底轴向力: %.1f kN', ... condition_names{WORKING_CONDITION}, results.surface_weight/1000, ... results.axial_force(1)/1000, results.axial_force(end)/1000), ... 'FontSize', 14, 'FontWeight', 'bold'); %% ===================== 8. 生成分析报告 ===================== fprintf('[8/8] 生成分析报告...\n'); fprintf('\n'); fprintf('========================================================================================================\n'); fprintf(' 电缆力学分析报告(轴向力修正版)\n'); fprintf('========================================================================================================\n'); fprintf('分析时间: %s\n', datestr(now)); fprintf('工况类型: %s\n', condition_names{WORKING_CONDITION}); fprintf('井深范围: 0 - %.0f m\n', depth(end)); fprintf('电缆规格: %.2f mm | 破断强度: %.0f kg\n', cable.OD*1000, cable.breaking_strength); fprintf('\n'); fprintf('【关键结果】\n'); fprintf('────────────────────────────────────────────────────────────────────────────────\n'); fprintf(' 参数 井口(0m) 井底(%.0fm) 变化\n', depth(end)); fprintf(' ────────── ────────── ────────── ──────────\n'); fprintf(' 张力 %8.2f kN %8.2f kN ↑ %.1f kN\n', ... results.tension(1)/1000, results.tension(end)/1000, (results.tension(end)-results.tension(1))/1000); fprintf(' 轴向力 %8.2f kN %8.2f kN ↓ %.1f kN\n', ... results.axial_force(1)/1000, results.axial_force(end)/1000, (results.axial_force(1)-results.axial_force(end))/1000); fprintf(' 累积摩阻 %8.2f kN %8.2f kN ↑ %.1f kN\n', ... results.cumulative_friction(1)/1000, results.cumulative_friction(end)/1000, results.cumulative_friction(end)/1000); fprintf('\n'); fprintf(' 表观悬重: %.2f kN\n', results.surface_weight/1000); fprintf(' 最小安全系数: %.2f\n', min(results.safety_factor(results.safety_factor>0))); if ~isempty(results.compression_points) fprintf(' 首个压缩点: %.0f m\n', results.compression_points(1)); end fprintf('\n【力学参数分布特征】\n'); fprintf('────────────────────────────────────────────────────────────────────────────────\n'); fprintf(' • 张力: 从井口向井底递增(0 → %.1f kN)\n', results.tension(end)/1000); fprintf(' • 轴向力: 从井口向井底递减(%.1f → %.1f kN)\n', results.axial_force(1)/1000, results.axial_force(end)/1000); fprintf(' • 井口轴向力代表整个系统的真实载荷\n'); fprintf(' • 井底轴向力<0表示电缆处于压缩状态\n'); fprintf('\n========================================================================================================\n'); fprintf('分析完成 | 总耗时: %.2f 秒\n', toc); fprintf('========================================================================================================\n\n'); %% ===================== 辅助函数定义 ===================== function [density_eq, E_eq, nu_eq, yield_eq] = calculate_composite_properties(cable) total_area = pi * (cable.OD/2)^2; density_sum = 0; E_sum = 0; nu_sum = 0; yield_sum = 0; for i = 1:size(cable.layers, 1) if i == 1 area = pi * (cable.layers(i,1)/2)^2; else area = pi * ((cable.layers(i,1)/2)^2 - (cable.layers(i-1,1)/2)^2); end volume_fraction = area / total_area; density_sum = density_sum + cable.layers(i,2) * volume_fraction; E_sum = E_sum + cable.layers(i,3) * 1e9 * volume_fraction; nu_sum = nu_sum + cable.layers(i,4) * volume_fraction; yield_sum = yield_sum + cable.layers(i,5) * volume_fraction; end density_eq = density_sum; E_eq = E_sum; nu_eq = nu_sum; yield_eq = yield_sum; end function [curvature, torsion, dogleg] = calculate_advanced_trajectory(depth, inc, azi) n = length(depth); curvature = zeros(n, 1); torsion = zeros(n, 1); dogleg = zeros(n, 1); for i = 2:n-1 ds = depth(i+1) - depth(i); u1 = [sin(inc(i-1))*sin(azi(i-1)); sin(inc(i-1))*cos(azi(i-1)); cos(inc(i-1))]; u2 = [sin(inc(i))*sin(azi(i)); sin(inc(i))*cos(azi(i)); cos(inc(i))]; u3 = [sin(inc(i+1))*sin(azi(i+1)); sin(inc(i+1))*cos(azi(i+1)); cos(inc(i+1))]; cos_dogleg = dot(u1, u2); cos_dogleg = max(-1, min(1, cos_dogleg)); dogleg(i) = acos(cos_dogleg) / ds; du = (u3 - u1) / (2*ds); curvature(i) = norm(du); if i > 2 && i < n-1 u0 = [sin(inc(i-2))*sin(azi(i-2)); sin(inc(i-2))*cos(azi(i-2)); cos(inc(i-2))]; u4 = [sin(inc(i+2))*sin(azi(i+2)); sin(inc(i+2))*cos(azi(i+2)); cos(inc(i+2))]; t1 = (u2 - u0) / (2*ds); t2 = (u3 - u1) / (2*ds); t3 = (u4 - u2) / (2*ds); dt = (t3 - t1) / (2*ds); if norm(cross(t2, dt)) > 1e-10 torsion(i) = dot(cross(t2, dt), t2) / norm(t2)^2; end end end curvature(1) = curvature(2); curvature(n) = curvature(n-1); torsion(1) = torsion(2); torsion(n) = torsion(n-1); dogleg(1) = dogleg(2); dogleg(n) = dogleg(n-1); curvature = min(curvature, 0.05); torsion = min(abs(torsion), 0.02) .* sign(torsion); dogleg = min(dogleg, 0.05); end function [north, east, tvd] = calculate_3d_coordinates(depth, inc, azi) n = length(depth); north = zeros(n, 1); east = zeros(n, 1); tvd = zeros(n, 1); for i = 2:n md = depth(i) - depth(i-1); inc_avg = (inc(i) + inc(i-1)) / 2; azi_avg = (azi(i) + azi(i-1)) / 2; dtvd = md * cos(inc_avg); dnorth = md * sin(inc_avg) * cos(azi_avg); deast = md * sin(inc_avg) * sin(azi_avg); tvd(i) = tvd(i-1) + dtvd; north(i) = north(i-1) + dnorth; east(i) = east(i-1) + deast; end end function [buckling_state, F_critical] = analyze_buckling(T, cable, diameter, inclination) if T >= -50 buckling_state = 0; F_critical = 0; return; end I = pi * cable.OD^4 / 64; clearance = (diameter - cable.OD) / 2; w_e = cable.weight_in_air * (1 - 1200/cable.density); % 有效浮重 if inclination > 0.1 % 斜井情况 r = clearance; % 径向间隙 EI = cable.E * I; % 抗弯刚度 % Dawson-Paslay模型的临界力计算公式 F_sin = 2 * sqrt(EI * w_e * sin(inclination) / r); % 正弦屈曲临界力 F_hel = 2*sqrt(2) * F_sin; % 螺旋屈曲临界力(约为正弦屈曲的2.8倍) else % 垂直井情况 % 使用垂直井的简化公式 F_sin = 3.0 * (cable.E * I * w_e^2 / clearance)^(1/3); F_hel = 5.0 * (cable.E * I * w_e^2 / clearance)^(1/3); end F = abs(T); % 当前轴向力绝对值 if F < F_sin buckling_state = 0; % 压缩未屈曲 F_critical = F_sin; elseif F < F_hel buckling_state = 1; % 正弦屈曲 F_critical = F_hel; else buckling_state = 2; % 螺旋屈曲 F_critical = F_hel; end end
11-26
function [TAi,MPi,FAi,Sk,fh,NAi,SKi,fs,Xs,Ys,Zs,Depth,wc,Locs,Loch,wg,DEP,depth,SF,asw,azhangli,amozu,ajiechu,aa]=abcfunc3(wc,a,b,Holedia,Strden,ml,qqlimit,vis,vs,tao,Tp,tractionForce) %% 基础参数输入 working_condition=wc; data=a; % 井眼轨迹数据 C=b; % 电缆的相关数据 VS=vs; % 上提下放速度 omega=4*pi; % 钻柱旋转角速度,对于测井电缆,旋转情况需根据实际调整 rhoi=0; % 内流体密度,测井电缆内部一般无流体,设为 0 miu=vis; % 钻井液塑性粘度,取值范围为 25 - 60 s Dw=Holedia; % 油管直径,单位为 m Den=Strden; % 电缆密度,单位为 kg/m3 taof=tao; % 钻井液动切力,即屈服值,取值范围为 5 - 25 pa miut=0.30; % 摩擦系数 miua=miut; E=ml; % 弹性模量 thegma=qqlimit; % 电缆屈服极限,单位为 MPa rhoo=rhoi+50; % 钻柱外流体密度 g=9.81; % 重力加速度 t=numel(C(:,1)); ltrans=zeros(t,1); dtrans=zeros(t,1); Dtrans=zeros(t,1); mtrans=zeros(t,1); for i=1:t ltrans(i)=round(C(i,3)); dtrans(i)=C(i,1); Dtrans(i)=C(i,2); mtrans(i)=Den*9.8*0.25*pi*(Dtrans(i)^2-dtrans(i)^2); end aa=12*10^(-6); %管材线膨胀系数,1/℃ TP=Tp; %油压 %% %钻具组合参数 Ntrans=numel(Dtrans); %所加载的钻具组合数量 nntrans=zeros(Ntrans,1); %即nntrans的前n项累加 ntrans=round(data(end,1)); Rt=zeros(ntrans,1); %各段钻具组合外半径/m rt=zeros(ntrans,1); %各段钻具组合内半径/m Aot=zeros(ntrans,1); %各段钻具组合外截面积/m^2 Ait=zeros(ntrans,1); %各段钻具组合内截面积/m^2 qt=zeros(ntrans,1); %各段钻具组合线重/N.m^-1 qmt=zeros(ntrans,1); %各段钻具组合浮重/N.m^-1 Kft=zeros(ntrans,1); %浮力系数 It=zeros(ntrans,1); %惯性矩 ht=zeros(ntrans,1); %壁厚 Pi=zeros(ntrans,1); %环空流体压力,MPa PI1=Pi; PI2=Pi; Tc=PI1; Po=zeros(ntrans,1); %管内流体压力,MPa PO1=Po; PO2=Pi; FH=zeros(ntrans,1); %正弦屈曲临界载荷,N FS=zeros(ntrans,1); %正弦屈曲临界载荷,N %计算nnntrans,便于后续计算 nntrans(1)=ltrans(1); for i=2:Ntrans nntrans(i)=nntrans(i-1)+ltrans(i); end %计算钻具组合各分段的半径、截面积、线重 for i=1:Ntrans if i==1 for j=1:nntrans(i) Rt(j)=Dtrans(i)/2; rt(j)=dtrans(i)/2; Aot(j)=pi*Rt(j)^2; Ait(j)=pi*rt(j)^2; qt(j)=mtrans(i); qmt(j)=qt(j)-(Aot(j)-Ait(j))*rhoi*g; Kft(j)=qmt(j)/qt(j); It(j)=pi*(Rt(j)^4-rt(j)^4)/4; ht(j)=Rt(j)-rt(j); end else for j=nntrans(i-1)+1:nntrans(i) Rt(j)=Dtrans(i)/2; rt(j)=dtrans(i)/2; Aot(j)=pi*Rt(j).^2; Ait(j)=pi*rt(j).^2; qt(j)=mtrans(i); qmt(j)=qt(j)-(Aot(j)-Ait(j))*rhoi*g; Kft(j)=qmt(j)/qt(j); It(j)=pi*(Rt(j)^4-rt(j)^4)/4; ht(j)=Rt(j)-rt(j); end end end for i=working_condition if i==1 sign1=-1;%-1为上提 sign2=0; v1=VS; %正为输出(下放),负为注入(上提),油管内流 v2=VS; %套管外流 Mium=15.7*10^(-6);%流体运动粘度 Rei=2*rhoi*abs(v1)*2*rt(100)/Mium; derta=0.05; if Rei<2320 nonmda=64/Rei; elseif Rei<=22*(2*rt(100)/derta)^(8/7) nonmda=0.316*Rei^(-0.25); elseif Rei<=597*(2*rt(100)/derta)^(9/8) nonmda=(1.14-2*log10(derta/2/rt(100)+21.25/Rei^0.9))^(-2); else nonmda=0.11*(2*rt(100)/derta)^0.25; end Fflow=sign(v1)*rhoi*nonmda*v1^2*pi*rt(100)/4; elseif i==2 sign1=1;% 1为下放 sign2=0; v1=VS; %正为输出(下放),负为注入(上提),油管内流 v2=VS; %套管外流 Mium=15.7*10^(-6);%流体运动粘度 Rei=2*rhoi*abs(v1)*2*rt(100)/Mium; derta=0.05; if Rei<2320 nonmda=64/Rei; elseif Rei<=22*(2*rt(100)/derta)^(8/7) nonmda=0.316*Rei^(-0.25); elseif Rei<=597*(2*rt(100)/derta)^(9/8) nonmda=(1.14-2*log10(derta/2/rt(100)+21.25/Rei^0.9))^(-2); else nonmda=0.11*(2*rt(100)/derta)^0.25; end Fflow=sign(v1)*rhoi*nonmda*v1^2*pi*rt(100)/4; end end %计算参数 Dw=Dw+0.001; ds=1; %步长 len=ntrans-1; %计算长度 sspan=0:ds:len; % 预分配 ft ca ba 变量 ft = zeros(ntrans,1); ca = zeros(ntrans - 1,1); ba = zeros(ntrans - 1,1); %% %处理井眼轨迹函数 [Mk,mk,Sk,alphak,phik]=deal_curve_data2(data); %初始化k和tao等参数 [alpha,phi,ks,dks,ddks,kphis,kalphas,taos]=prepare_data(sspan,Mk,mk,Sk,alphak,phik); DAL=zeros(ntrans,1); ALpha=flipud(alpha*180/pi); for i=1:ntrans-1 DAL(i)=abs((ALpha(i+1)-ALpha(i))*30)*10; end [~,Ys,Zs,Xs]=deal_trank(Sk,alpha,phi); %输出从底部到井口位置的自然长度以及各点坐标 for i=1:ntrans PI1(i)=rhoi*g*cos(alpha(ntrans-i+1)); PO1(i)=rhoo*g*cos(alpha(ntrans-i+1)); Tc(i)=4.5/100*cos(alpha(i)); ft(i)=aa*E*(Aot(i)-Ait(i))*Tc(i); end PI2(1)=PI1(1); PO2(1)=PO1(1); for i=2:ntrans PI2(i)=PI2(i-1)+PI1(i); PO2(i)=PO2(i-1)+PO1(i); end PI2=PI2-9.81*rhoi;%油管静液柱压力 PO2=PO2-9.81*rhoo;%环空静液柱压力 AiT=flipud(Ait);%正序外截面积 AoT=flipud(Aot);%正序内截面积 df=PI2; for i=1:ntrans df(i)=(PI2(i)+TP*10^6-TP/ntrans/2*i*10^6)*AiT(i)-PO2(i)*AoT(i); end Pi=flipud(PI2); Po=flipud(PO2); %% %识别井眼轨迹 smalpha=smooth(alpha,201); B=mean(smalpha(1:500)); B1=mean(smalpha(1:300)); Z=mean(smalpha(end-300:end)); Ls=0; Lh=0; if B<0.44 %若所有井斜角小于25°,则当成直井处理 Ls=ntrans; else for i=1:ntrans-1 ca(i)=abs(smalpha(end-i+1)-Z); Ls=i+1;%垂直井段 if ca(i)>0.11 && smalpha(end-i+1)>0.25 break end end end if B1>0.44 && B<1.23 %判断25-70°的定向井 for i=1:ntrans-1 if Ls==ntrans break else ba(i)=abs(smalpha(i)-B1); Lh=i;%水平井段 if ba(i)>0.11 && smalpha(i)<0.44 break end end end else for i=1:ntrans-1 if Ls==ntrans break else ba(i)=abs(smalpha(i)-B); Lh=i;%水平井段 if ba(i)>0.11 && smalpha(i)<1.35 break end end end end Lc=ntrans-Lh-Ls; %造斜段 for i=1:ntrans if i<Lh+1 FH(i)=-2.75*(E*It(i)*qt(i)*sin(alpha(i))/(0.5*Dw-Rt(i)))^(0.5);%斜直,水平 FS(i)=-2*(E*It(i)*qt(i)*sin(alpha(i))/(0.5*Dw-Rt(i)))^(0.5); elseif i<Lh+1+Lc FH(i)=-7.56*E*It(i)/(Rt(i)*(0.5*Dw-Rt(i)))/(10^3);%弯曲 FS(i)=-4*E*It(i)/(Rt(i)*(0.5*Dw-Rt(i)))/(10^3); else FH(i)=-5.55*(E*It(i)*qt(i)^2)^(1/3);%垂直 FS(i)=-2.55*(E*It(i)*qt(i)^2)^(1/3); end end Fh=FH; Fs=FS; Miud=zeros(ntrans,1); % 初始化摩擦系数数组 Miudi = zeros(i, 1); % 转换 alpha 为角度 alpha_deg = alpha * 180 / pi; % 根据工作条件和井斜角设置摩擦系数 for j = 1:i current_alpha = alpha_deg(ntrans - j + 1); % 注意 alpha 数组的索引逻辑 if wc == 1 % 上提工况 if current_alpha <= 5 Miudi(j) = 0.25; elseif current_alpha <= 45 Miudi(j) = 0.35; elseif current_alpha <= 85 Miudi(j) = 0.45; else Miudi(j) = 0.5; end else % 下放工况 if current_alpha <= 5 Miudi(j) = 0.2; elseif current_alpha <= 45 Miudi(j) = 0.3; elseif current_alpha <= 85 Miudi(j) = 0.4; else Miudi(j) = 0.5; end end end %% %开始计算 Ttemp=0;%钻压 Mtemp=0;%扭矩 Nbtemp=0; Nntemp=0; A(1)=0; B(1)=0; C(1)=0; shendu=round(data(end,1)); DI=50; if mod(shendu,50)==0 Sdep=50; else Sdep=mod(shendu,50); end DEP=zeros(round((ntrans-Sdep+DI)/DI),1); for i=1:round((ntrans-Sdep+DI)/DI) DEP(i)= Sdep+DI*(i-1); end TAi=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); MAi=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); FAi=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); fai=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); NAi=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); SKi=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); MPi=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); Ks=flipud(ks); Dks=flipud(dks); Ddks=flipud(ddks); Kphis=flipud(kphis); Kalphas=flipud(kalphas); Taos=flipud(taos); MMk=flipud(Mk); mmk=flipud(mk); Alphak=flipud(alphak); Phik=flipud(phik); RT=flipud(Rt); rT=flipud(rt); qmT=flipud(qmt); AiT=flipud(Ait); AoT=flipud(Aot); IT=flipud(It); fh=flipud(FH)/1000; fs=flipud(FS)/1000; for i=Sdep:DI:ntrans SKi(1:i,(i-Sdep+DI)/DI)=1:1:i; ksi=flipud(Ks(1:i)); dksi=flipud(Dks(1:i)); ddksi=flipud(Ddks(1:i)); kphisi=flipud(Kphis(1:i)); kalphasi=flipud(Kalphas(1:i)); taosi=flipud(Taos(1:i)); sspani=sspan(1:i); Rti=flipud(RT(1:i)); rti=flipud(rT(1:i)); qmti=flipud(qmT(1:i)); Aiti=flipud(AiT(1:i)); Aoti=flipud(AoT(1:i)); Iti=flipud(IT(1:i)); Mki=flipud(MMk(1:i)); mki=flipud(mmk(1:i)); Ski=Sk(1:i); alphaki=flipud(Alphak(1:i)); phiki=flipud(Phik(1:i)); fti=ft(1:i); Pii=Pi(1:i); Poi=Po(1:i); Fs1=flipud(Fs); Fh1=flipud(Fh); Fsi=flipud(Fs1(1:i)); Fhi=flipud(Fh1(1:i)); Miudi=Miud(1:i); opts=odeset('events',@eventfun); [~,y]=ode45(@(s,y)odefunc2(s,y,ksi,dksi,ddksi,kphisi,kalphasi,taosi,sspani,v1,v2,omega,taof,miu,Rti,rti,Dw,Miudi,miut,qmti... ,Aiti,Aoti,rhoi,rhoo,E,Iti,g,Mki,mki,Ski,alphaki,phiki,Ttemp,Mtemp,Nbtemp,Nntemp,sign1,sign2,Fsi,Fhi,fti,Pii,Poi,Fflow,tractionForce),0:ds:(i-1)*ds,[A(1);B(1);C(1)],opts); Depth=length(Ski); numl=min(y(:,1))/1000; if numl<=-200 break; end TAi(1:i,(i-Sdep+DI)/DI)=flipud(y(:,1))/1000; %轴力,KN MAi(1:i,(i-Sdep+DI)/DI)=flipud(y(:,2))/1000; FAi(1:i,(i-Sdep+DI)/DI)=flipud(y(:,3))*miua/1000; %摩阻,KN fai(1:i,(i-Sdep+DI)/DI)=flipud(y(:,3)); NAi(1:i-1,(i-Sdep+DI)/DI)=-diff(fai(1:i,(i-Sdep+DI)/DI)); PP=MAi(1:i,(i-Sdep+DI)/DI); for j=1:i PP(j)=abs(y(j,1))/(Aot(j)-Ait(j))/1000000; end MPi(1:i,(i-Sdep+DI)/DI)=flipud(PP)/thegma; Depth=length(Ski(:,1)); end thegma1=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); thegma2=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); thegma3=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); Thegma=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); Thegma1111=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); SF=zeros(ntrans,round((ntrans-Sdep+DI)/DI)); TAii=TAi(1,any(TAi~=0)); NUM=length(TAii(1,:)); for i=1:1:length(TAi(1:end,1)) for j=1:length(TAi(1,1:NUM)) thegma1(i,j)=((RT(i)^2+rT(i)^2)*PI2(i)-2*RT(i)^2*PO2(i))/(RT(i)^2-rT(i)^2); thegma2(i,j)=-PI2(i); thegma3(i,j)=TAi(i,j)*1000/(AoT(i)-AiT(i)); Thegma(i,j)=(((thegma2(i,j)-abs(thegma1(i,j)))^2+(thegma1(i,j)-abs(thegma3(i,j)))^2+(thegma3(i,j)-abs(thegma2(i,j)))^2)/2)^0.5; Thegma1111(i,j)=(((thegma2(i,j)-thegma1(i,j))^2+(thegma1(i,j)-thegma3(i,j))^2+(thegma3(i,j)-thegma2(i,j))^2)/2)^0.5; SF(i,j)=(Thegma(i,j)/1000000)/thegma; end end %% TAii=TAi(:,any(TAi~=0)); %提取非零列 k=length(TAii(1,:)); wg=TAi(1,:); %悬重 depth=length(fs(:,1)); KS=zeros(ntrans,1); KH=zeros(ntrans,1); for i=1:ntrans if TAi(i,k)<fs(i) && TAi(i,k)>fh(i) KS(i)=TAi(i,k); else KS(i)=0; end end for i=1:ntrans if TAi(i,k)<=fh(i) KH(i)=TAi(i,k); else KH(i)=0; end end Locs=find(KS<0); %取出正弦屈曲的位置 Loch=find(KH<0); %取出螺旋屈曲的位置 function [value,isterminal,direction]=eventfun(~,y) value=y(:,1)/1000+300; isterminal=1; direction=0; end % 提取并转置数据 asw = TAi(1, :)'; azhangli = TAi(:, end); amozu = FAi(end:-1:1, end); ajiechu = NAi(:, end); end
08-09
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值