模拟圆柱体的热传导.txt

问题描述: 一杯沸水冷却

% 一杯沸水冷却，圆柱体模型，底面半径0.05m，高0.1m，周围温度20度，初始水温100度
% 方程是四维输运方程（常数a^2=k/(c*p)，k是热传导系数0.6006焦/（米*秒*度））
% 初始条件：t=0时水等于100度
% 边界条件：1.上下壁都是自由冷却，第三类边界条件，周围温度保持在20度（H=k/h,h取1）
%           2.杯壁绝热，第二类边界条件
% du/dt=a^2*{d[ρ*du/dρ]/ρ+d^2u/dz^2}
% 

clc;clear;close all;
N=200;
H=0.6006;
a2=H;
figure;
% \copyright: zjliu
% Author's email: [e-mail] 
Nh=round(N+N/5);
u0=100*ones(N);
rs=linspace(-0.05,0.05,N);
zs=linspace(0,0.1,N);
[rho,z]=meshgrid(rs,zs);
hp=image(rs,zs,u0);
set(gca,'ydir','reverse');
set(gcf,'DoubleBuffer','on');
shading('interp');
xlim([-0.1,0.1]);
ylim([-0.01,0.11]);
colorbar;
t=0;
dt=0.04;
drho=rho(2,2)-rho(2,1);
dz=z(1,2)-z(1,1);
tx=title('圆柱的热传导问题: t=0');
xlabel('\rho (/m)','fontsize',16);
ylabel('z (/m)','fontsize',16);
while t<300;
    du2r=diff(diff(u0,1,2).*rho(:,2:end),1,2)./rho(:,2:end-1);
    du2r=[zeros(N,1),du2r,zeros(N,1)];
    du2z=diff(u0,2,1);
    du2zu=-[u0(1,:)-20]/H^2;
    du2zd=-[u0(end,:)-20]/H^2;
    du2z=[du2zu;du2z;du2zd];
    u0=u0+dt*a2*[du2r+du2z];
    set(hp,'CData',u0);
    t=t+dt;
    set(tx,'string',['圆柱的热传导问题: t=',num2str(t)]);
    pause(0.02);
end