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;

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