Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

%code for steady state
n=11
x=linspace(0,1,n)
y=linspace(0,1,n)
x=y
Told=zeros(n,n)
Told(:,1)=400
Told(:,n)=800
Told(n,:)=900
Told(1,:)=600
T2old=zeros(n,n)
T2old(:,1)=400
T2old(:,n)=800
T2old(n,:)=900
T2old(1,:)=600
T3old=zeros(n,n)
T3old(:,1)=400
T3old(:,n)=800
T3old(n,:)=900
T3old(1,:)=600
T=Told,T2=T2old,T3=T3old
error1=9e9,error2=9e9,error3=9e9
gauss_jacobi_iteration=0;gauss_siedel_iteration=0;sor_iteration=0;
m=3
b=2*cos(pi/n)
a=4/(2+sqrt(4-b^2))
for method=1:m
    if method==1
while error1>1e-4
    gauss_jacobi_iteration=gauss_jacobi_iteration+1

    for j=2:n-1
        for i=2:n-1
            T(i,j)=(Told(i+1,j)+Told(i-1,j)+Told(i,j+1)+Told(i,j-1))/4
            
       
      end
         
    end
     
         error1=max(max(abs(T-Told)))
        Told=T
 
end

subplot(3,3,1)
surf(x,-y,Told)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')
title(['gauss jacobi iteration=' num2str(gauss_jacobi_iteration)])
subplot(3,3,2)
contour(x,-y,Told)
xlabel("position(m)")
ylabel("position(m)")
zlabel("temperature(k)")
subplot(3,3,3)
pcolor(x,-y,Told)
xlabel("position(m)")
ylabel("position(m)")
zlabel("temperature(k)")




    
    elseif method==2
        while error2>1e-4
    gauss_siedel_iteration=gauss_siedel_iteration+1

    for j=2:n-1
        for i=2:n-1
            T2(i,j)=(T2old(i+1,j)+T2(i-1,j)+T2old(i,j+1)+T2(i,j-1))/4
       
        end
    end
     error2=max(max(abs(T2-T2old)))
T2old=T2

end

subplot(3,3,4)
surf(x,-y,T2old)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')
title(['gauss siedel iteration=' num2str( gauss_siedel_iteration)])
subplot(3,3,5)
contour(x,-y,T2old)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')
subplot(3,3,6)
pcolor(x,-y,T2old)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')



    

    else
          while error3>1e-4
    sor_iteration=sor_iteration+1

    
     for j=2:n-1
        for i=2:n-1
        T3(i,j)=a*(T3old(i+1,j)+T3(i-1,j)+T3old(i,j+1)+T3(i,j-1))/4+(1-a)*(T3(i,j))
            
            
       
        end
    end
     error3=max(max(abs(T3-T3old)))
T3old=T3
    end
end
end

subplot(3,3,7)
surf(x,-y,T3old)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')
title(['Sor iteration=' num2str(sor_iteration)])
subplot(3,3,8)
contour(x,-y,T3old)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')
subplot(3,3,9)
pcolor(x,-y,T3old)
xlabel('position(m)')
ylabel('position(m)')
zlabel('temperature(k)')
% code for explicit unsteady
n=11
x=linspace(0,1,n)
y=linspace(0,1,n)
dx=x(2)-x(1)
dy=dx
x=y
Told=zeros(n,n)
Told(:,1)=400
Told(:,n)=800
Told(n,:)=900
Told(1,:)=600


T=Told


alpha=0.001
dt=1
t=820
for k=1:(t/dt)

    for j=2:n-1
        for i=2:n-1
            T(i,j)=Told(i,j)+((alpha*dt)/((dx)^2))*(Told(i,j+1)+Told(i,j-1)+Told(i-1,j)+Told(i+1,j)-4*Told(i,j))
            
       
      end
         
    end
     
  
        Told=T
 
end
figure(1)
surf(x,-y,Told)
xlabel("position(m)")
ylabel("position(m)")
title("temperature(K)")
figure(2)
contour(x,-y,Told)
xlabel("position(m)")
ylabel("position(m)")
figure(3)
pcolor(x,-y,Told)
xlabel("position(m)")
ylabel("position(m)")
% code for implicit unsteady
n=11
x=linspace(0,1,n)
y=linspace(0,1,n)
x=y
dx=x(2)-x(1)
dy=dx
Told=zeros(n,n)
Told(:,1)=400
Told(:,n)=800
Told(n,:)=900
Told(1,:)=600
T2old=zeros(n,n)
T2old(:,1)=400
T2old(:,n)=800
T2old(n,:)=900
T2old(1,:)=600
T3old=zeros(n,n)
T3old(:,1)=400
T3old(:,n)=800
T3old(n,:)=900
T3old(1,:)=600
T_1=zeros(n,n)
T_1(:,1)=400
T_1(:,n)=800
T_1(n,:)=900
T_1(1,:)=600
T_2=zeros(n,n)
T_2(:,1)=400
T_2(:,n)=800
T_2(n,:)=900
T_2(1,:)=600
T_3=zeros(n,n)
T_3(:,1)=400
T_3(:,n)=800
T_3(n,:)=900
T_3(1,:)=600
T=Told,T2=T2old,T3=T3old
error1=9e9,error2=9e9,error3=9e9
gauss_jacobi_iteration=0;gauss_siedel_iteration=0;sor_iteration=0;
m=3
alpha=0.001
dt=10
t=200
k=((alpha*dt)/(dx)^2)
a=1.2
for method=1:m
    if method==1
        
        for r=1:(t/dt)
            error1=9e9
            while error1>1e-4
                gauss_jacobi_iteration=gauss_jacobi_iteration+1
            
                
            for j=2:n-1
                for i=2:n-1
        T(i,j)=(k/(1+4*k))*(Told(i+1,j)+Told(i-1,j)+Told(i,j+1)+Told(i,j-1))+(T_1(i,j)/(1+4*k))
                end
            end
            error1=max(max(abs(T-Told)))
            Told=T
        
            
            end
            T_1=T
        end
        subplot(3,3,1)
        surf(x,-y,Told)
        xlabel("position(m)")
        ylabel("position(m)")
        zlabel("temperature(k)")
        title(['gauss jacobi iteration= ' num2str(gauss_jacobi_iteration)])
        subplot(3,3,2)
        contour(x,-y,Told)
         xlabel("position(m)")
        ylabel("position(m)")
        subplot(3,3,3)
        pcolor(x,-y,Told)
         xlabel("position(m)")
        ylabel("position(m)")
        
        
      
    elseif method==2
          
        for r=1:(t/dt)
            error2=9e9
            while error2>1e-4
                gauss_siedel_iteration=gauss_siedel_iteration+1  
           
                
            for j=2:n-1
                for i=2:n-1
        T2(i,j)=(k/(1+4*k))*(T2old(i+1,j)+T2(i-1,j)+T2old(i,j+1)+T2(i,j-1))+(T_2(i,j)/(1+4*k))
                end
            end
            error2=max(max(abs(T2-T2old)))
            T2old=T2
            
           
            end
         T_2=T2
         end
         subplot(3,3,4)
         surf(x,-y,T2old)
        xlabel("position(m)")
        ylabel("position(m)")
        zlabel("temperature(k)")
        title(['gauss siedel iteration=' num2str(gauss_siedel_iteration)])
        subplot(3,3,5)
        contour(x,-y,T2old)
        xlabel("position(m)")
        ylabel("position(m)")
        subplot(3,3,6)
        pcolor(x,-y,T2old)
        xlabel("position(m)")
        ylabel("position(m)")
   
    else
        
         for r=1:(t/dt)
             error3=9e9
             while error3>1e-4
                sor_iteration=sor_iteration+1
            
                
            for j=2:n-1
                for i=2:n-1
        T3(i,j)=a*((k/(1+4*k))*(T3old(i+1,j)+T3(i-1,j)+T3old(i,j+1)+T3(i,j-1))+(T_3(i,j)/(1+4*k)))+(1-a)*T3(i,j)
                end
            end
            error3=max(max(abs(T3-T3old)))
            T3old=T3
            
          
             end
           T_3=T3
        end
        subplot(3,3,7)
        surf(x,-y,T3old)
         xlabel("position(m)")
        ylabel("position(m)")
        zlabel("temperature(k)")
        title(['sor iteration=' num2str(sor_iteration)])
        subplot(3,3,8)
        contour(x,-y,T3old)
         xlabel("position(m)")
        ylabel("position(m)")
        subplot(3,3,9)
        pcolor(x,-y,T3old)
        xlabel("position(m)")
        ylabel("position(m)")
        
         
    end
end 

number of grids taken for simultaion is 11.n=11 so dx=(l/n-1)=0.1

according to the boundary conditions  T at top surface=600K

bottom surface=900K

left surface=400K

right surface=800K

These conditions have been considered in the code posted and taken care of.

The intitial error has be taken as 9e9 which is an arbitrary value and a while loop was used for a converged matrice and the minimum error was 1e-4.There were 3 schemes that were run for the steady state problem.

for steady state after solving the 2d diffusion equaton we arrived at

T(i,j)=(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))/4

scheme 1 gauss jacobi method

scheme 2 gauss siedel method

scheme 3 successive over relaxtive method

firstly the matrices of T were given zero(n,n) so as to consider the values for first iteration and then the scheme was implemented accordingly in the code. The gauss jacobi method took 264 iterations to convege.The gauss siedel method convereged for less number of iterations which is 140 as in gauss siedel the value is taken form the same iteration but the previous equation which makes the converegene easy.SOR took 35 iterations to converge as the value of a(over relaxive factor ) was 1.2 which gave a quicker convergence.NOTE: as the value of a was increased close to 2 say 1.99 the solution started to diverge thats why the value close to 1 was chosen. 

NOTE: the plots that are made w.r.t x and (-y) because the origin is supposed to be the top left corner of the 1 by 1 unit square. so the values of y ranges form (0 to -1).

so the surface temperature at the bottom is 900K so heat transfer is from bottom surfce to top surface and from right to left which is clearly visible by the plots.

Explicit unsteady:

T(p+1)=T(p)+(dT/dt)(p)

In case of explicit unsteady the value of a flow field is know at a previous time step and the value at the next time step is to be computed.

alpha(therma diffusivity is taken as 0.001) time step=1sec total time =820 second

Implicit unsteady:

in implicit unsteady the value at requiered tim step is calculated and the derivative at the required time step is know.

T(p+1)=T(p)+(dT/dt)(p+1)

so this was incorporated in the code and was solved with 3 different schemes

the error tolerance was given as 1e-4 and the intital error was 9e-9.

there were  new matrices  intorduced which were T_1,2,3 which was done to update the Tempertaure matrix w.r.t to time. 

Gauss jacobi iteration=853

gausssiedel iteration=500

sor iteration=320

over relaxation factor=1.2