mid term project unsteady state analysis

clear all
close all
clc
% explicit scheme 
% plate lenght in x,y direction
lx=1;
ly=1;
% number of grid in x,y direction
nx=10;
ny=10;
dx=lx/(nx-1);
dy=ly/(ny-1);
% temperature

T1=400;
T2=600;
T3=800;
T4=900;
% themal conductivity,density,specific heat 
k=0.72;
rho=1920;
cp=835;
dt=0.01;
alpha=k/(rho*cp);
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;
%n tim step
deltaT=1300;
tol=0.1;

% maximum number of of time step before stopping
lmax=10^4;
% initializing the arrayfun
T=T1*ones(nx+1,ny+1);
% initializing boundry rows and column it does not have any effect in the program but plaus crucial role in 
% in ploting the contour
T(2:nx,ny+1)=T4
T(1,ny+1)=T3
T(nx+1,ny+1)=T4
T(1,1)=T1
% creating node location of x y
x=0:dx:lx;
y=0:dy:ly;

 for L=1:lmax
       Told=T
        for j=2:ny-1
    for i=2:nx-1
     % calculating temprature for new time step
     T(i,j)=Told(i,j)+k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1));
            end
        end
        Tplot(L)=T(nx,ny);
        % diff
        diff=abs((T(nx,ny) - tol)/tol*100);
        fprintf(\'time step = %8.0f - diff. = %10.6f percent\\n\',L, diff);
        % exit iteration loop to reach stesdy stste
        if (diff<tol)
         break
     end
     % number of contour for plot
        
        end
        % number of contour for plot
        nc=50;
        dT=(T2-T1)/nc;
        v=T1:dT:T2;
        colormap(jet)
        contourf(x,y,T\',v, \'linestyle\',\'none\')
        colorbar
        axis equal tight
        title([\'contour plot of temperature in deg. C at time=\', num2str(deltaT*L/3600), \'h\'])
        xlabel(\'x(m)\')
        ylabel(\'y(m)\')
        pause(0.01)

% implicit schemes 2d heat condunction eqn
clear all;
close all
clc;


%initial variables
nx = 20; %grid number for x axis
ny = 20; %grid number for y axis

x = 1; % width in x dimension\',unit is m
y = 1; % width in Y dimension\',unit is m

% initial condition set up

T1 = 400; %constant temperature at left side,unit is K
T2 = 600; %constant temperature at top side,unit is K
T3 =800 ; %constant temperature at right side,unit is K
T4 = 900; %constant temperature at bttom side,unit is K

Ti = 298; %initial temperature of the plate
k = 237; %conduction coefficience,unit is W/m*K
alpha = 9.71e-5; %thermal diffusivity of the material, unit is m^2/s


%width,lenght of the grid
dx = x/(nx-1); 
dy = y/(ny-1); 

%setup dt and the experime
 

delta_t = 0.001; %set the basic dt. make sure to satisfy convergence condition.
%initialize the minimum time equal to experiment time

%initiate the temperature array of the plate
T = ones(nx,ny); %T array holds the temperature caculated from the last iteration
Tnew = Ti*ones(nx,ny);

%    apply boundary conditio
% since we have constant temperature boundary condition, we can apply it
% out of the main loop (boundary temperatures will not change inside the main loop).

%right boundary condition(constant temperature T1)
for j=1:ny
i=nx;
Tnew(i,j)=T1;
end
%top boundary condition (constant temperature T2)
for i=2:nx-1
j=ny;
Tnew(i,j)=T2;
end
%left boundary condition (constant temperature T3)
for j=1:ny
i=1;
Tnew(i,j)=T3;
end
%bottom boundary condition(constant temperature T4)
for i=2:nx-1
j=1;
Tnew(i,j)=T4;
end
%treat the corner points seperately in order to avoid singular gridpoints.
Tnew(1,ny) = 0.5*(T3+T2);
Tnew(nx,ny) = 0.5*(T1+T2);
Tnew(nx,1) = 0.5*(T4+T1);
Tnew(1,1) = 0.5*(T4+T3);
T = Tnew;

% main loop

%state, and calculate the time used to reach the steady state.

while counter_t > 0

%interior node
        for i=2:nx-1
            for j=2:ny-1
             Tnew(i,j) = T(i,j)+alpha*delta_t*((T(i+1,j)-2*T(i,j)+T(i-1,j)/dx^2)+(T(i,j+1-2*T(i,j)+T(i,j-1)/dy^2))
     end
  end
end