APPENDIX: MATLAB Source Code
%% APPENDIX
disp('2D Transient Heat Conduction in Square Ingot');
k=input('Enter thermal conductivity k (W/mK): ');
rho=input('Enter density rho (kg/m3): ');
cp=input('Enter specific heat cp (J/kgK): ');
W=input('Enter total width of ingot (m): ');
Ti=input('Enter initial temperature (K): ');
Tinf=input('Enter ambient temperature (K): ');
h=input('Enter heat transfer coefficient (W/m2K): ');
alpha=k/(rho*cp);
L=W/2;
Bi=h*L/k;
disp(['alpha = ',num2str(alpha),' m2/s Bi = ',num2str(Bi)]);
N=20;
lambda=zeros(1,N);
for m=1:N
if m==1
guess=5;
else
guess=(m-1)*pi/L+1;
end
lambda(m)=fzero(@(x) x*tan(x*L)-h/k,guess);
end
C=zeros(1,N);
A=zeros(N,N);
for m=1:N
denom=lambda(m)*L+sin(lambda(m)*L)*cos(lambda(m)*L);
C(m)=4*sin(lambda(m)*L)/denom;
end
for m=1:N
for n=1:N
A(m,n)=C(m)*C(n)*(Ti-Tinf);
end
end
nx=31;
ny=31;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
[X,Y]=meshgrid(x,y);
Fo_vals=[0.01,0.05,0.1,0.5,2.0];
tvals=Fo_vals*L^2/alpha;
T_ana=zeros(ny,nx,5);
for kt=1:5
t=tvals(kt);
for i=1:nx
for j=1:ny
s=0;
for m=1:N
for n=1:N
s=s+A(m,n)*cos(lambda(m)*x(i))*cos(lambda(n)*y(j))*...
exp(-alpha*(lambda(m)^2+lambda(n)^2)*t);
end
end
T_ana(j,i,kt)=Tinf+s;
end
end
end
MATLAB APPENDIX | Page 1
dt=0.2*dx^2/alpha;
Fo_x=alpha*dt/dx^2;
Fo_y=Fo_x;
disp(['Stability: Fo_x+Fo_y = ',num2str(Fo_x+Fo_y)]);
T=Ti*ones(ny,nx);
T_num=zeros(ny,nx,5);
steps=round(tvals/dt);
for step=1:max(steps)
Tnew=T;
for i=2:nx-1
for j=2:ny-1
Tnew(j,i)=(1-2*Fo_x-2*Fo_y)*T(j,i)+Fo_x*(T(j,i-1)+T(j,i+1))+Fo_y*(T(j-1,i)+T(j+1,i));
end
end
Tnew(:,1)=Tnew(:,2);
Tnew(1,:)=Tnew(2,:);
Bi_dx=h*dx/k;
Tnew(:,nx)=(T(:,nx-1)+Bi_dx*Tinf)/(1+Bi_dx);
Tnew(ny,:)=(T(ny-1,:)+Bi_dx*Tinf)/(1+Bi_dx);
Tnew(ny,nx)=(T(ny-1,nx)+T(ny,nx-1)+2*Bi_dx*Tinf)/(2+2*Bi_dx);
T=Tnew;
for k=1:5
if step==steps(k)
T_num(:,:,k)=T;
end
end
end
disp('Fo Time Center_A Center_N Corner_A Corner_N');
for kt=1:5
fprintf('%.2f %6.2f %.1f %.1f %.1f %.1f\n',Fo_vals(kt),tvals(kt),...
T_ana(1,1,kt),T_num(1,1,kt),T_ana(ny,nx,kt),T_num(ny,nx,kt));
end
disp('Fo MaxError MeanError');
for kt=1:5
err=abs(T_ana(:,:,kt)-T_num(:,:,kt));
fprintf('%.2f %.2f %.2f\n',Fo_vals(kt),max(err(:)),mean(err(:)));
end
figure(1);
for kt=1:5
subplot(2,5,kt);
contourf(X*100,Y*100,T_ana(:,:,kt),20);
colorbar;
axis equal tight;
title(['Ana Fo=',num2str(Fo_vals(kt))]);
subplot(2,5,kt+5);
contourf(X*100,Y*100,T_num(:,:,kt),20);
colorbar;
axis equal tight;
title(['Num Fo=',num2str(Fo_vals(kt))]);
end
figure(2);
tfine=linspace(0,tvals(end),200);
Tc_ana=zeros(size(tfine));
for p=1:length(tfine)
s=0;
for m=1:N
for n=1:N
MATLAB APPENDIX | Page 2
s=s+A(m,n)*exp(-alpha*(lambda(m)^2+lambda(n)^2)*tfine(p));
end
end
Tc_ana(p)=Tinf+s;
end
Tc_num=squeeze(T_num(1,1,:));
plot(tfine,Tc_ana,'b-',tvals,Tc_num,'ro','MarkerFaceColor','r');
xlabel('Time');
ylabel('Center Temperature');
legend('Analytical','Numerical');
grid on;
Tavg_ana=mean(T_ana(:,:,end),'all');
Tavg_num=mean(T_num(:,:,end),'all');
E_removed_ana=rho*cp*W^2*(Ti-Tavg_ana);
E_removed_num=rho*cp*W^2*(Ti-Tavg_num);
E_error=abs(E_removed_ana-E_removed_num)/E_removed_ana*100;
disp(['Energy Error: ',num2str(E_error),' %']);
disp('Grid Independence:');
grids=[21,31,41,51];
Tc_g=zeros(1,4);
for g=1:4
ng=grids(g);
xg=linspace(0,L,ng);
yg=linspace(0,L,ng);
dxg=xg(2)-xg(1);
Fo_xg=alpha*dt/dxg^2;
T_g=Ti*ones(ng,ng);
steps_g=round(tvals(4)/dt);
for step=1:steps_g
T_new_g=T_g;
for i=2:ng-1
for j=2:ng-1
T_new_g(j,i)=(1-2*Fo_xg-2*Fo_xg)*T_g(j,i)+Fo_xg*(T_g(j,i-1)+T_g(j,i+1)+T_g(j-1,i)+T_g(j+1,i)
end
end
T_new_g(:,1)=T_new_g(:,2);
T_new_g(1,:)=T_new_g(2,:);
Bi_dxg=h*dxg/k;
T_new_g(:,ng)=(T_g(:,ng-1)+Bi_dxg*Tinf)/(1+Bi_dxg);
T_new_g(ng,:)=(T_g(ng-1,:)+Bi_dxg*Tinf)/(1+Bi_dxg);
T_g=T_new_g;
end
Tc_g(g)=T_g(1,1);
if g==1
fprintf('%dx%d %.1f -\n',ng,ng,Tc_g(g));
else
fprintf('%dx%d %.1f %.2f\n',ng,ng,Tc_g(g),abs((Tc_g(g-1)-Tc_g(g))/Tc_g(g-1)*100));
end
end
disp('Time Step Independence:');
dt_vals=[0.03,0.025,0.02,0.015];
Tc_dt=zeros(1,4);
MATLAB APPENDIX | Page 3
for d=1:4
dtv=dt_vals(d);
Fo_xt=alpha*dtv/dx^2;
steps_t=round(tvals(4)/dtv);
T_t=Ti*ones(ny,nx);
for step=1:steps_t
T_new_t=T_t;
for i=2:nx-1
for j=2:ny-1
T_new_t(j,i)=(1-2*Fo_xt-2*Fo_xt)*T_t(j,i)+Fo_xt*(T_t(j,i-1)+T_t(j,i+1)+T_t(j-1,i)+T_t(j+1,i)
end
end
T_new_t(:,1)=T_new_t(:,2);
T_new_t(1,:)=T_new_t(2,:);
T_new_t(:,nx)=(T_t(:,nx-1)+Bi_dx*Tinf)/(1+Bi_dx);
T_new_t(ny,:)=(T_t(ny-1,:)+Bi_dx*Tinf)/(1+Bi_dx);
T_t=T_new_t;
end
Tc_dt(d)=T_t(1,1);
if d==1
fprintf('%.3f %.1f -\n',dtv,Tc_dt(d));
else
fprintf('%.3f %.1f %.2f\n',dtv,Tc_dt(d),abs((Tc_dt(d-1)-Tc_dt(d))/Tc_dt(d-1)*100));
end
end
disp('SIMULATION COMPLETED');
MATLAB APPENDIX | Page 4