案例描述

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

02

温度场计算

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

1min	5min
10min	20min

03

matlab源代码附上matlab源代码如下：

% 参数设置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 (℃)');% 主循环for step = 1:n_steps    T_old = T;    f_old = f;    k_old = k;    c_old = c;    for i = 1:nx        for j = 1:ny            % 当前节点属性            current_k = k_old(i,j);            current_f = f_old(i,j);            current_T = T_old(i,j);              % 确定材料属性            if current_f >= 1                current_c = c_water;                current_k_val = k_water;            else                current_c = c_ice;                current_k_val = k_ice;            end            % 计算热流            Q_total = 0;            % 东向            if i < nx                neighbor_k = k_old(i+1,j);                k_east = 2*current_k_val*neighbor_k/(current_k_val + neighbor_k);                Q_total = Q_total + k_east*(T_old(i+1,j) - current_T);            elseif i == nx                Q_total = Q_total + h*(T_env - current_T)*dx;            end            % 西向            if i > 1                neighbor_k = k_old(i-1,j);                k_west = 2*current_k_val*neighbor_k/(current_k_val + neighbor_k);                Q_total = Q_total + k_west*(T_old(i-1,j) - current_T);            elseif i == 1                Q_total = Q_total + h*(T_env - current_T)*dx;            end            % 北向            if j < ny                neighbor_k = k_old(i,j+1);                k_north = 2*current_k_val*neighbor_k/(current_k_val + neighbor_k);                Q_total = Q_total + k_north*(T_old(i,j+1) - current_T);            elseif j == ny                Q_total = Q_total + h*(T_env - current_T)*dx;            end            % 南向            if j > 1                neighbor_k = k_old(i,j-1);                k_south = 2*current_k_val*neighbor_k/(current_k_val + neighbor_k);                Q_total = Q_total + k_south*(T_old(i,j-1) - current_T);            elseif j == 1                Q_total = Q_total + h*(T_env - current_T)*dx;            end            % 计算热量            delta_Q = Q_total * dt;            mass = rho_ice * dx^2;            % 相变处理            if current_T < 0                delta_T = delta_Q / (mass * current_c);                T_new_val = current_T + delta_T;                f_new_val = 0;                if T_new_val >= 0                    Q_used = (0 - current_T) * mass * current_c;                    Q_remaining = delta_Q - Q_used;                    if Q_remaining > 0                        delta_f = Q_remaining / (L * mass);                        f_new_val = delta_f;                    end                end            elseif current_T == 0 && current_f < 1                delta_f = delta_Q / (L * mass);                f_new_val = current_f + delta_f;            else                delta_T = delta_Q / (mass * current_c);                T_new_val = current_T + delta_T;                f_new_val = current_f;            end            % 处理完全融化            if f_new_val >= 1                excess = f_new_val - 1;                T_new_val = excess * L / c_water;                f_new_val = 1;                k(i,j) = k_water;                c(i,j) = c_water;            end            % 更新值            T(i,j) = T_new_val;            f(i,j) = f_new_val;        end    end    % 更新图形    if mod(step, 10) == 0        set(h_plot, 'CData', T);        title(sprintf('Time: %.2f s', step*dt));        drawnow;    endend% 显示最终结果figure;pcolor(T);shading interp;          % 双线性插值colormap(jet(1024));     % 1024级颜色渐变colorbar;axis equal tight;title('Final Temperature Distribution (℃)');xlabel('X');ylabel('Y');

04

matlab源文件下载

https://pan.baidu.com/s/1ydrtAgMV7d9_1DIUdgZNNQ?pwd=xxn5 提取码: xxn5