%% ------------- 用户设置 -------------
xlsx_path = '附件1 放热能力数据.xlsx'; % Excel文件路径
sheetIndex = 1; % 工作表索引
scanRate_C_per_min = 5; % 扫描速率 °C/min
%% 参数初始化
h_body = 1.70; w = 60; % 身高(m), 体重(kg)
S_body = 0.0061*h_body + 0.0128*w - 0.1529; % 人体表面积(m²)
S = 1.25 * S_body; % 防护服面积(m²)
% 材料参数
h1 = 0.7e-3; % 内层厚度(m)
k1 = 0.068; c1 = 4803.8; rho1 = 208; P1_per_kg = 1000; % 500元/500g -> 1000元/kg
h2_options = [0.4e-3, 0.45e-3]; % 中间层厚度选项(m)
k2 = 0.06; c2 = 2400; rho2 = 552.3; P2_per_sqm = 10; % 中间层价格: 10元/m²
h3_options = (0.3:0.3:1.2)*1e-3; % 外层厚度选项(m) [0.3,0.6,0.9,1.2]mm
k3 = 0.0527; c3 = 5463.2; rho3 = 300; P3_per_kg = 300; % 150元/500g -> 300元/kg
% 环境参数
h_c = 3.0; % 对流换热系数(W/m²·K) 静止状态
T0 = 37; % 人体温度(°C)
T_env = -40; % 环境温度(°C)
% 原始资金计算 (h1=0.7mm, h2=0.4mm, h3=0.3mm)
inner_weight_orig = h1 * S * rho1;
outer_weight_orig = 0.3e-3 * S * rho3;
P0 = (inner_weight_orig * P1_per_kg) + ...
(outer_weight_orig * P3_per_kg) + ...
(S * P2_per_sqm);
max_P = 1.5 * P0; % 最大允许资金
%% ------------- 读取 DSC 数据并构造 dQ/dT(T) -------------
% 检查文件是否存在
if ~exist(xlsx_path, 'file')
error('找不到Excel文件: %s', xlsx_path);
end
[~, sheets] = xlsfinfo(xlsx_path);
tbl = readtable(xlsx_path, 'Sheet', sheets{sheetIndex});
T_dsc = tbl{:,1}; % 温度数据 (°C)
DSC_mW_per_mg = tbl{:,2}; % DSC数据 (mW/mg)
% 单位转换和扫描速率校正
scanRate = scanRate_C_per_min / 60; % 转换为 °C/s
DSC_W_per_kg = DSC_mW_per_mg * 1e3; % mW/mg -> W/kg
dQdT_vals = DSC_W_per_kg / scanRate; % J/(kg·K)
% 创建插值函数
dQdT_fun = @(Tq) interp1(T_dsc, dQdT_vals, Tq, 'pchip', 0);
% 为了兼容原代码，也创建DSC插值函数
DSC_interp = @(T) interp1(T_dsc, DSC_mW_per_mg, T, 'linear', 0);
%% 声明全局变量
global R1 R2 R3 R_conv C1 C2 C3 T0 T_env DSC_interp rho2 h2 S
% 遍历所有方案
max_time = 0; best_h2 = 0; best_h3 = 0;
results = [];
idx = 1;
fprintf('开始计算各种厚度组合...\n');
for i = 1:length(h2_options)
h2 = h2_options(i);
for j = 1:length(h3_options)
h3 = h3_options(j);
% 检查资金约束
outer_weight = h3 * S * rho3;
P_new = (inner_weight_orig * P1_per_kg) + ...
(outer_weight * P3_per_kg) + ...
(S * P2_per_sqm);
if P_new > max_P
fprintf('方案 h2=%.2fmm, h3=%.1fmm 超出资金限制，跳过\n', h2*1000, h3*1000);
continue;
end
% 计算重量约束时间
M_c = (h1*rho1 + h2*rho2 + h3*rho3) * S;
t_weight = (100 - M_c) / 0.05; % 重量限制时间(s)
% 计算热阻和热容
R1 = h1 / (k1 * S);
R2 = h2 / (k2 * S);
R3 = h3 / (k3 * S);
R_conv = 1 / (h_c * S);
C1 = c1 * rho1 * h1 * S;
C2 = c2 * rho2 * h2 * S;
C3 = c3 * rho3 * h3 * S;
% 设置ODE求解器
tspan = [0, 20000]; % 时间范围: 0~20000秒
y0 = [37, 37, 37]; % 初始温度: [T1, T2, T3]
options = odeset('Events', @tempEvent, 'RelTol', 1e-6, 'AbsTol', 1e-8);
try
% 求解ODE
[t, y, te] = ode15s(@(t,y) heatEqns(t, y), tspan, y0, options);
% 获取温度降至15℃的时间
if ~isempty(te)
t_temp = te(1);
else
t_temp = tspan(end);
end
% 有效时间
t_eff = min(t_temp, t_weight);
% 保存结果
results(idx,:) = [h2*1000, h3*1000, t_eff, P_new, M_c];
idx = idx + 1;
% 更新最优解
if t_eff > max_time
max_time = t_eff;
best_h2 = h2;
best_h3 = h3;
end
fprintf('h2=%.2fmm, h3=%.1fmm: t=%.1fs (%.2fmin), 成本=%.1f元, 重量=%.1fkg\n', ...
h2*1000, h3*1000, t_eff, t_eff/60, P_new, M_c);
catch ME
fprintf('方案 h2=%.2fmm, h3=%.1fmm 计算失败: %s\n', h2*1000, h3*1000, ME.message);
end
end
end
%% 输出结果
fprintf('\n========== 优化结果 ==========\n');
fprintf('最优方案: h2=%.2fmm, h3=%.1fmm, t=%.1fs (约%.2f分钟)\n', ...
best_h2*1000, best_h3*1000, max_time, max_time/60);
if ~isempty(results)
fprintf('\n所有方案结果:\n');
result_table = array2table(results, 'VariableNames', ...
{'h2_mm', 'h3_mm', 't_eff_s', 'cost_yuan', 'weight_kg'});
disp(result_table);
% 按有效时间排序
[~, sort_idx] = sort(results(:,3), 'descend');
fprintf('\n按有效时间排序的前5个方案:\n');
top5_results = results(sort_idx(1:min(5,size(results,1))),:);
top5_table = array2table(top5_results, 'VariableNames', ...
{'h2_mm', 'h3_mm', 't_eff_s', 'cost_yuan', 'weight_kg'});
disp(top5_table);
else
fprintf('警告: 没有找到满足约束条件的方案！\n');
end
%% ========== 辅助函数 ==========
function dydt = heatEqns(~, y)
% 解包变量
T1 = y(1); T2 = y(2); T3 = y(3);
% 全局参数
global R1 R2 R3 R_conv C1 C2 C3 T0 T_env DSC_interp rho2 h2 S
% 相变放热计算 (仅在14.7°C < T2 < 25°C时)
if T2 > 14.7 && T2 < 25
DSC_val = DSC_interp(T2); % mW/mg
phase_power_per_kg = -DSC_val * 1000; % W/kg (1 mW/mg = 1000 W/kg)
Q_phase = phase_power_per_kg * (rho2 * h2 * S); % W
else
Q_phase = 0;
end
% 热流计算
Q01 = (T0 - T1) / (0.5*R1);
Q12 = (T1 - T2) / (0.5*R1 + 0.5*R2);
Q23 = (T2 - T3) / (0.5*R2 + 0.5*R3);
Q34 = (T3 - T_env) / (0.5*R3 + R_conv);
% 温度变化率
dT1dt = (Q01 - Q12) / C1;
dT2dt = (Q12 - Q23 + Q_phase) / C2;
dT3dt = (Q23 - Q34) / C3;
dydt = [dT1dt; dT2dt; dT3dt];
end
function [value, isterminal, direction] = tempEvent(~, y)
value = y(1) - 15; % T1 > 15°C时继续
isterminal = 1; % 触发事件时停止
direction = -1; % 温度下降时触发
end我想获得最优的方案，但最后为啥中间层和最外层和最初的一样呢
最新发布
扫一扫
635
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
