从井底到井口的方式计数问题

358 篇文章 ¥29.90 ¥99.00
使用动态规划解决一个人从井底每次爬一步或三步到达井口的不同方式数问题。定义递推关系,初始化数组并进行循环计算,通过Python实现代码展示解题过程。

从井底到井口的方式计数问题

问题描述:
假设一个人位于井底,并希望通过爬行到达井口。这个人每次可以爬上一步或者三步。我们的任务是确定他有多少种不同的方式能够到达井口。

解决方案:
为了解决这个问题,我们可以使用动态规划的方法。我们将创建一个数组来存储到达每个位置的方式数量,并使用递推关系来计算每个位置的值。

我们假设井底为位置0,井口为位置n。我们可以定义一个长度为n+1的数组ways[],其中ways[i]表示从位置0到位置i的方式数量。

递推关系:
对于位置i,我们有两种方式可以到达它:

  1. 从位置i-1爬一步到达位置i;
  2. 从位置i-3爬三步到达位置i。

因此,我们可以使用以下递推关系计算ways[i]:
ways[i] = ways[i-1] + ways[i-3]

初始条件:
我们可以设置ways[0] = 1,表示从井底到井底的方式数量为1。并且当i < 0时,我们可以设置ways[i] = 0。

下面是使用Python编写的示例代码:

def count_ways(n):
    <
% 修正:轴向力从井口(最大)向井底(最小)递减 % 基于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 = 1; % 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.0132, 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 = 2490; % 总重 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 = 1180; % 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 results.stress_state(i) = sqrt(sigma_axial^2 + 3*tau_torsion^2 + sigma_bend^2); if results.stress_state(i) > 0 results.safety_factor(i) = cable.breaking_strength/(w_cable(i)+F_friction(i));%拉断力/(电缆浮重+工具串浮重+摩阻) else results.safety_factor(i) = 10; 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-18
源码地址: https://pan.quark.cn/s/d1f41682e390 miyoubiAuto 米游社每日米游币自动化Python脚本(务必使用Python3) 8更新:更换cookie的获取地址 注意:禁止在B站、贴吧、或各大论坛大肆传播! 作者已退游,项目不维护了。 如果有能力的可以pr修复。 小引一波 推荐关注几个非常可爱有趣的女孩! 欢迎B站搜索: @嘉然今天吃什么 @向晚大魔王 @乃琳Queen @贝拉kira 第三方库 食用方法 下载源码 在Global.py中设置米游社Cookie 运行myb.py 本地第一次运行时会自动生产一个文件储存cookie,请勿删除 当前仅支持单个账号! 获取Cookie方法 浏览器无痕模式打开 http://user.mihoyo.com/ ,登录账号 按,打开,找到并点击 按刷新页面,按下图复制 Cookie: How to get mys cookie 当触发时,可尝试按关闭,然后再次刷新页面,最后复制 Cookie。 也可以使用另一种方法: 复制代码 浏览器无痕模式打开 http://user.mihoyo.com/ ,登录账号 按,打开,找到并点击 控制台粘贴代码并运行,获得类似的输出信息 部分即为所需复制的 Cookie,点击确定复制 部署方法--腾讯云函数版(推荐! ) 下载项目源码和压缩包 进入项目文件夹打开命令行执行以下命令 xxxxxxx为通过上面方式或取得米游社cookie 一定要用双引号包裹!! 例如: png 复制返回内容(包括括号) 例如: QQ截图20210505031552.png 登录腾讯云函数官网 选择函数服务-新建-自定义创建 函数名称随意-地区随意-运行环境Python3....
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值