clc
clear all
close all
%% 1-Inputs section
name=('Titanium');
denisty = 4500;
conductivity = 22;
spacific_heat = 0.54;
Lx = 8;
Ly = 10;
Nx = 10;
Ny = 10;
T_initial= 0 ;
T_east = 0 ;
T_west = 0 ;
T_north = 550 ;
T_south = 550 ;
t_end=550;
dt = 10;
tolerence = 100;
tolerence_ss=0.1;
%% 2- Constants Section
k=1;
err_SS_max(k)=1;
err_SS_min(k)=1;
dx=Lx/Nx;
dy=Ly/Ny;
n_time=round(t_end/dt);
alpha = conductivity/(spacific_heat*denisty);
T_max=max([T_east T_west T_north T_south T_initial]);
T_min=min([T_east T_west T_north T_south T_initial]);
if dt<= 1/(2*alpha*((1/dx^2)+(1/dy^2)))
else
fprintf('Error, the stability condition is not met\nPlease return to "Inputs Section" and choose a "dt" smaller than %f \n',1/(2*alpha*((1/dx^2)+(1/dy^2))))
return
end
% ----------------- Initial Conditions for finite difference section ---------------
T=zeros(Nx+2,Ny+2,75000);
T(:,1,:)=T_south;
T(:,Ny+1,:)=T_north;
T(:,Ny+2,:)=T_north;
T(Nx+1,:,:)=T_east;
T(Nx+2,:,:)=T_east;
T(1,:,:)=T_west;
T(:,:,1)=T_initial;
% ------------------- Initial Conditions for steady state section -------------------
Tss=zeros(Nx+2,Ny+2); Tss2=zeros(Nx+2,Ny+2);
Tss(:,1)=T_south; Tss2(:,1)=T_south;
Tss(:,Ny+1)=T_north; Tss2(:,Ny+1)=T_north;
Tss(:,Ny+2)=T_north; Tss2(:,Ny+2)=T_north;
Tss(Nx+1,:)=T_east; Tss2(Nx+1,:)=T_east;
Tss(Nx+2,:)=T_east; Tss2(Nx+2,:)=T_east;
Tss(1,:)=T_west; Tss2(1,:)=T_west;
%% 3- Steady-State section
while err_SS_max(k)>=tolerence_ss || err_SS_min(k)>=tolerence_ss
for i=2:Nx
for j=2:Ny
Tss2(i,j)=0.25*( Tss(i+1,j) + Tss(i,j+1) + Tss(i-1,j) + Tss(i,j-1) );
end
end
k=k+1;
err_SS_max(k)=abs(max(max(Tss2-Tss)));
err_SS_min(k)=abs(min(min(Tss2-Tss)));
Tss=Tss2;
end
%% 4- Finite difference section
k=1;
err_E_max(k)=100;
err_E_min(k)=100;
while err_E_max(k)>=tolerence || err_E_min(k)>=tolerence
for i=2:Nx
for j=2:Ny
T(i,j,k+1) =T(i,j,k)+dt*alpha*(((T(i-1,j,k)-2*T(i,j,k)+T(i+1,j,k))/dx^2)+((T(i,j-1,k)-2*T(i,j,k)+T(i,j+1,k))/dy^2));
end
end
k=k+1;
err_E_max(k)=abs(max(max(T(:,:,k)-Tss)));
err_E_min(k)=abs(min(min(T(:,:,k)-Tss)));
if round(err_E_max(k),5)==round(err_E_max(k-1),5) && err_E_max(k)~= 0
errordlg('The solution is not converging, Please choose a larger tolerence','Tolerence Error');
close(message)
return
end
if round(err_E_min(k),5)==round(err_E_min(k-1),5) && err_E_min(k)~= 0
errordlg('The solution is not converging, Please choose a larger tolerence','Tolerence Error');
close(message)
return
end
end
T=T(:,:,1:k);
SStime=k*dt;
%% 6- Plotting section
x=zeros(1,Nx+2);y=zeros(1,Ny+2);
for i = 1:Nx+2
x(i) =(i-1)*dx;
end
for i = 1:Ny+2
y(i) =(i-1)*dy;
end
% -------------- Constant plot ----------------
figure
% suptitle('TEMPERATURE DISTRIBUTION USING FINITE DIFF METHOD');
subplot(2,1,1)
surf(x,y,Tss)
title(sprintf('TEMPERATURE DISTRIBUTION USING FINITE DIFFERENCE METHOD\n\n Temperature at steady state time : %i seconds ',round(SStime)))
caxis([T_min T_max]);
view(90,-90);
xlim([0 Lx+dx]); xlabel('Length');
ylim([0 Ly+dy]); ylabel('Width');
zlim([T_min T_max]); zlabel('Temprature');
drawnow
% ------------ Animated plot ----------
for j=1:n_time
subplot(2,1,2)
surf(x,y,T(:,:,j))
title('Animation')
caxis([T_min T_max]);
view(90,-90);
xlim([0 Lx+dx]); xlabel('Length');
ylim([0 Ly+dy]); ylabel('Width');
zlim([T_min T_max]); zlabel('Temprature');
drawnow
end
surf(x,y,Tss)
view(90,-90);
xlim([0 Lx+dx]); xlabel('Length');
ylim([0 Ly+dy]); ylabel('Width');
zlim([T_min T_max]); zlabel('Temprature');
title(sprintf('Temperature at time : %i seconds ',round(SStime)))
