## MATLAB Assignment Solution on Finite Difference Project

• 5th Oct, 2021
• 19:07 PM
```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;

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)))
```