% 一杯沸水冷卻，圓柱體模型，底面半徑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}
% 
function heat
clc;clear;close all;
N=200;
H=0.6006;
a2=H;
figure;
% \copyright: zjliu
% Author's email: zjliu2001@163.com 
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