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