% 修正:轴向力从井口(最大)向井底(最小)递减
% 基于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