如图1所示，边长为5cm的冰块，初始温度为-2℃，放在25℃的环境中自然冷却，对流换热系数为10W/m²K，本文将通过matlab编程求解冰块融化的过程，计算其温度场。

图1 案例示意图

02

温度场计算

本文通过matlab分别计算t=1min、5min、10min和20min后的温度云图及瞬态温度变化过程（计算总时间30min），计算结果展示如下表：

1min	5min
10min	20min

 

% 参数设置
Lx = 0.05; Ly = 0.05;   % 冰块尺寸（5cm x 5cm）
nx = 31; ny = 31;         % 网格数量
dx = Lx/(nx-1); dy = Ly/(ny-1);
T_initial = -2;           % 初始温度（℃）
T_env = 25;               % 环境温度（℃）
h = 10;                   % 对流换热系数（W/m²K）% 材料属性（冰）
rho_ice = 920;           % 密度（kg/m³）
k_ice = 2.18;            % 导热系数（W/mK）
c_ice = 2100;            % 比热容（J/kgK）
L = 334000;              % 潜热（J/kg）% 材料属性（水）
k_water = 0.6;           % 导热系数（W/mK）
c_water = 4200;          % 比热容（J/kgK）% 时间参数
alpha_ice = k_ice/(rho_ice*c_ice);
dt = 1 * dx^2/(4*alpha_ice); % 时间步长
total_time = 605;             % 总模拟时间（秒）
n_steps = round(total_time/dt);% 初始化
T = T_initial * ones(nx, ny);
f = zeros(nx, ny);       % 液相分数
k = k_ice * ones(nx, ny);
c = c_ice * ones(nx, ny);% 创建图形窗口
figure;
h_plot = pcolor(T);
shading interp;          % 双线性插值
colormap(jet(1024));     % 1024级颜色渐变
colorbar;
axis equal tight;
title('Temperature (℃)');% 预计算常数
mass = rho_ice * dx^2;  % 每个节点的质量（假设厚度为1m）% 主循环
for step = 1:n_stepsT_old = T;f_old = f;k_old = k;c_old = c;% ========== 向量化热流计算 ========== %% 初始化各方向热流矩阵[Q_east, Q_west, Q_north, Q_south] = deal(zeros(nx, ny));% 东向热流 (i < nx)if nx > 1k_east = 2 * k_old(1:end-1,:) .* k_old(2:end,:) ./ (k_old(1:end-1,:) + k_old(2:end,:) + eps);Q_east(1:end-1,:) = k_east .* (T_old(2:end,:) - T_old(1:end-1,:));end% 西向热流 (i > 1)if nx > 1k_west = 2 * k_old(2:end,:) .* k_old(1:end-1,:) ./ (k_old(2:end,:) + k_old(1:end-1,:) + eps);Q_west(2:end,:) = k_west .* (T_old(1:end-1,:) - T_old(2:end,:));end% 北向热流 (j < ny)if ny > 1k_north = 2 * k_old(:,1:end-1) .* k_old(:,2:end) ./ (k_old(:,1:end-1) + k_old(:,2:end) + eps);Q_north(:,1:end-1) = k_north .* (T_old(:,2:end) - T_old(:,1:end-1));end% 南向热流 (j > 1)if ny > 1k_south = 2 * k_old(:,2:end) .* k_old(:,1:end-1) ./ (k_old(:,2:end) + k_old(:,1:end-1) + eps);Q_south(:,2:end) = k_south .* (T_old(:,1:end-1) - T_old(:,2:end));end% 边界对流条件Q_boundary = zeros(nx, ny);Q_boundary(1,:)   = h*(T_env - T_old(1,:))*dx;   % 西边界Q_boundary(end,:) = h*(T_env - T_old(end,:))*dx; % 东边界 Q_boundary(:,1)   = h*(T_env - T_old(:,1))*dx;   % 南边界Q_boundary(:,end) = h*(T_env - T_old(:,end))*dx; % 北边界% 总热流Q_total = Q_east + Q_west + Q_north + Q_south + Q_boundary;delta_Q = Q_total * dt;% ========== 向量化相变计算 ========== %T_new = T_old;f_new = f_old;% 计算材料属性掩膜is_water = f_old >= 1;% 情况1: 冰升温（T < 0）ice_mask = T_old < 0 & ~is_water;delta_T_ice = delta_Q ./ (mass * c_ice);T_new(ice_mask) = T_old(ice_mask) + delta_T_ice(ice_mask);% 处理过零冰节点over_zero = ice_mask & T_new >= 0;Q_used = (0 - T_old(over_zero)) * mass * c_ice;Q_remaining = delta_Q(over_zero) - Q_used;Q_remaining(Q_remaining < 0) = 0;delta_f = Q_remaining / (L * mass);T_new(over_zero) = 0;f_new(over_zero) = delta_f;% 情况2: 相变中（T == 0且f < 1）melting_mask = (T_old == 0) & (f_old < 1);delta_f_melt = delta_Q(melting_mask) / (L * mass);f_new(melting_mask) = f_old(melting_mask) + delta_f_melt;% 情况3: 水升温（T >= 0且f >= 1）water_mask = ~ice_mask & ~melting_mask;delta_T_water = delta_Q(water_mask) ./ (mass * c_water);T_new(water_mask) = T_old(water_mask) + delta_T_water(water_mask);% 处理完全融化full_melt = f_new > 1;excess = f_new(full_melt) - 1;T_new(full_melt) = excess * L / c_water;f_new(full_melt) = 1;% 更新材料属性k = k_water * full_melt + k_ice * ~full_melt;c = c_water * full_melt + c_ice * ~full_melt;% 更新场变量T = T_new;f = f_new;% ========== 可视化更新 ========== %if mod(step, 10) == 0set(h_plot, 'CData', T);title(sprintf('Time: %.2f s', step*dt));drawnow;end
end% 显示最终结果
figure;
pcolor(T);
shading interp;
colormap(jet(1024));
colorbar;
axis equal tight;
title('Final Temperature Distribution (℃)');
xlabel('X');
ylabel('Y');