BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

SShanaj Ali5 Mar 2021

AIM:- To write code for steady and unsteady equations .

GIVEN:-

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

EXPLANATION:-

1. Steady state:

steady state equation for 2D is `((d^2T)/dx^2)+ ((d^2T)/dy^2)`=0

There are three iterative solvers which can be used to solve equations. Those three are jacobi,gauss-seidel and succesive over relation(SOR). Lets see how they works.

a. code for jacobi method:-

clear all
close all
clc

%(dT^2/dx^2)+(dT^2/dy^2)=0

%initial conditions

%unit square
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);;


T=400*ones(nx,ny);
T(1,:)= 600;     %Top
T(nx,:)= 900;    %Bottom
T(:,1)= 400;     %Left
T(:,ny)= 800;    %Right
Told=T;
tol=1e-4;
error=9e9;
jacobi_iter=1


%convergence loop
while error > tol
  %space loop
  for i=2:nx-1
    for j=2:ny-1
        T(i,j)= 0.25*(Told(i,j-1)+Told(i,j+1)+Told(i-1,j)+Told(i+1,j));
    endfor
  endfor
  error=max(max(abs(Told-T)));
  jacobi_iter=jacobi_iter+1;
  Told=T;
end
 contour(x,y,T,'showText','on')
 xlabel('x')
 ylabel('y')
 title_text1=sprintf('jacobi iteration= %d',jacobi_iter)
 title_text2=sprintf('Temperature field for jacobi steady method')
 title({title_text1,title_text2})

b. Code for Gauss-seidel method:-

clear all
close all
clc

%(dT^2/dx^2)+(dT^2/dy^2)=0

%initial conditions

%unit square
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);;


T=400*ones(nx,ny);
T(1,:)= 600;     %Top
T(nx,:)= 900;    %Bottom
T(:,1)= 400;     %Left
T(:,ny)= 800;    %Right
Told=T;
tol=1e-4;
error=9e9;
gs_iter=1


%convergence loop
while error > tol
  gs_iter=gs_iter+1
  Told=T
  %space loop
  for i=2:nx-1
    for j=2:ny-1
        T(i,j)= 0.25*(T(i,j-1)+T(i,j+1)+T(i-1,j)+T(i+1,j));
    endfor
  endfor
  error=max(max(abs(Told-T)));
end
 contour(x,y,T,'showText','on')
 xlabel('x')
 ylabel('y')
  title_text1=sprintf('gs iteration= %d',gs_iter)
 title_text2=sprintf('Temperature field for gs steady method')
 title({title_text1,title_text2})

c. Code for succesive over relaxation(SOR)method:-

clear all
close all
clc

%(dT^2/dx^2)+(dT^2/dy^2)=0

%initial conditions

%unit square
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);;


T=400*ones(nx,ny);
T(1,:)= 600;     %Top
T(nx,:)= 900;    %Bottom
T(:,1)= 400;     %Left
T(:,nx)= 800;    %Right
Told=T;
tol=1e-4;
error=9e9;
w=1.1;
 sor_iter=1;

%convergence loop
while error > tol
 
  %space loop
  for i=2:nx-1
    for j=2:ny-1
        T(i,j)=(1-w)* Told(i,j) +w*(0.25*(T(i,j-1)+T(i,j+1)+T(i-1,j)+T(i+1,j)));
    endfor
  endfor
  error=max(max(abs(Told-T)));
  sor_iter=sor_iter+1;
  Told=T
end
 contour(x,y,T,'showText','on')
 xlabel('x')
 ylabel('y')
  title_text1=sprintf('SOR iteration= %d',sor_iter)
 title_text2=sprintf('Temperature field for SOR steady method')
 title({title_text1,title_text2})

As you will observe in graphs that all look same but the only difference is number of iterations used. You will observe that SOR take less number of iterations as compare to other two methods.

Now lets look for unsteady /transient equation

2. Unsteady equation=`(dT)/dt+(d^2T)/dx^2+(d^2T)/dy^2`=0

These are two schemes for unsteady equation. One is Explicit and other is implicit.

a. Explicit scheme code:-

clear all
close all
clc

%(dT^2/dx^2)+(dT^2/dy^2)=0

%initial conditions

%unit square
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);


T=400*ones(nx,ny);
T(1,:)= 600;     %Top
T(nx,:)= 900;    %Bottom
T(:,1)= 400;     %Left
T(:,nx)= 800;    %Right

alpha=1.4;   
Told=T;
t=0.4;    %time 
c=0.5;    %courant number
dt=(c*(dx^2))/(2*alpha);
n=t/dt;


%time loop

for k=1:n
  for i= 2:nx-1
    for j=2:ny-1
      H=((T(i,j+1)-2*Told(i,j)+T(i,j-1))/(dx^2));
      V=((T(i+1,j)-2*Told(i,j)+T(i-1,j))/(dy^2));
      T(i,j)=Told(i,j)+alpha*dt*(H+V)
    endfor
  endfor
  error=max(max(abs(Told-T)));
  Told=T;
end


contour(x,y,T,'showtext','on')
colorbar;
xlabel('Length')
ylabel('Width')
title('Temperature field for explicit method')

(NOTE:- For explicit ,make sure to use courant number(c) equals to 0.5 for 2D heat conduction equation)

b. Code for implicit scheme(Use three iterative solvers- Jacobi,gs,SOR):-

clear all
close all
clc

%(dT^2/dx^2)+(dT^2/dy^2)=0

%initial conditions

%unit square
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);


T=400*ones(nx,ny);
T(1,:)= 600;     %Top
T(nx,:)= 900;    %Bottom
T(:,1)= 400;     %Left
T(:,ny)= 800;    %Right
Told=T;
error=1;

%define time loop
t=40;
nt=0.1;
n=t/nt;
alpha=1.4;      %thermal diffusivity

%define constant terms
k1=(alpha*nt)/dx^2;
k2=(alpha*nt)/dy^2;
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
Tprev=T;
tol=1e-4;
%iterative solver selection
iterative_solver=input('iterative_solver(1. jacobi 2. Gauss siedel 3. SOR)=');

%jacobi method
if iterative_solver==1
  jacobi_iter=1
  jacobi_iter_timestep=1;
  
  for k=1:n
    error=9e9;
    while(error>tol)
    for i=2:nx-1
      for j=2:ny-1
        H=(Told(i,j+1)+Told(i,j-1));
        V=(Told(i+1,j)+Told(i-1,j));
        T(i,j)=Tprev(i,j)*term1+H*term2+V*term3;
        
      endfor
      
    endfor
    error=max(max(abs(Told-T)));
    Told=T;
    jacobi_iter=jacobi_iter+1;
   
  end
contour(x,y,T,'showtext','on')
colorbar;
xlabel('length')
ylabel('width')
title_text1=sprintf('Total jacobi iteration= %d',jacobi_iter);
title_text2=sprintf('Total jacobi Timesteps=%d',jacobi_iter_timestep);
title({title_text1,title_text2})
jacobi_iter_timestep=jacobi_iter_timestep+1;
pause(0.05)
   if(T-Tprev< tol)
     break
  end
  Tprev=T;
 end
end

%gauss seidel
if iterative_solver==2
  gs_iter=1;
  gs_iter_timesteps=1;
  for k=1:n
    error=9e9;
    while(error>tol)
    for i=2:nx-1
      for j=2:ny-1
        H=(T(i,j+1)+T(i,j-1));
        V=(T(i+1,j)+T(i-1,j));
        T(i,j)=Tprev(i,j)*term1+H*term2+V*term3;
        
      endfor
    end
   error=max(max(abs(Told-T)));
   Told=T;
   gs_iter=gs_iter+1; 
 end
contour(x,y,T,'showtext','on')
colorbar
xlabel('length')
ylabel('width')
title_text1=sprintf('Total gs iteration= %d',gs_iter);
title_text2=sprintf('Total gs Timesteps= %d',gs_iter_timesteps);
title({title_text1,title_text2})
gs_iter_timesteps=gs_iter_timesteps+1;
pause(0.05)
 if(T-Tprev < tol)
 break
end
Tprev=T;
 end
end

%successive over relaxation
if iterative_solver==3
  SOR_iter=1;
  SOR_iter_timesteps=1;
  for k=1:n
    error=9e9;
    while(error>tol)
    for i=2:nx-1
      for j=2:ny-1
        H=(T(i,j+1)+T(i,j-1));
        V=(T(i+1,j)+T(i-1,j));
        T(i,j)=((1-alpha)*Told(i,j))+alpha*(Tprev(i,j)*term1+H*term2+V*term3);
        
      endfor
    end
   error=max(max(abs(Told-T)));
   Told=T;
   SOR_iter=SOR_iter+1; 
 end
contour(x,y,T,'showtext','on')
colorbar;
xlabel('length')
ylabel('width')
title_text1=sprintf('Total SOR iteration= %d',SOR_iter);
title_text2=sprintf('Total SOR Timesteps=%d',SOR_iter_timesteps);
title({title_text1,title_text2})
SOR_iter_timesteps=SOR_iter_timesteps+1;
pause(0.05)
 if(T-Tprev < 1e-4)
    break
  end
 Tprev=T;
 end
end

Observation:- As one can observe ,explicit is direct method i.e it doesnt use any iterative solvers to solve but on the same time its unstable whereas implicit scheme use three iterative solvers to solve but it is stable.Meanwhile one more point to note here is, SOR is far better than other two methods(jacobi,gauss-seidel).

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.