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

Solving the steady and unsteady 2D heat conduction problem(MATLAB)

SSaiyan Kumar Behera Behera3 Aug 2020

Objective-To solve the 2D conduction equation using a Transient solver and a Steady state solver using Iterative techniques (Jacobi,Gauss Seidal,SOR).

1. Assume that the domain is a unit square.

2. Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

3. Boundary conditions for steady and transient case

1. Left:400K

2. Top:600K

3. Right:800K

4. Bottom:900K

Steady state analysis using(Jacobi,Gauss Seidal,SOR)

The 2D steady state heat conduction equation is given as:

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

Discertizing the equation using central difference scheme we get

`(T_(i-1,j)^n-2*T_(i,j)^n+T_(i+1,j)^n)/(Deltax^2)+(T_(i,j-1)^n-2*T_(i,j)^n+T_(i,j+1)^n)/(Deltay^2)=0`

Rearraning the above equation we get

`T_(i,j)=C1(T_(i+1,j)+T_(i-1,j))+C2(T_(i,j+1)+T_(i,j-1))`

C1=`(Deltay^2)/(2*(Deltax^2+Deltay^2)` C2=`(Deltax^2)/(2*(Deltax^2+Deltay^2)`

For Jacobi iteration the formula becomes,

`T_(i,j)^(n+1)=(T_(i+1,j)^n+T_(i-1,j)^n+T_(i,j+1)^n+T_(i,j-1)^n)/4`

For Gauss Seidal iteration the formula becomes,

`T_(i,j)^(n+1)=(T_(i-1,j)^(n+1)+T_(i+1,j)^n+T_(i,j-1)^(n+1)+T_(i,j+1)^n)/4`

For SOR iteration the formula becomes,

`T_(i,j)^(n+1)=(1-w)T_(i,j)^n+w*T_(jacobi)`

Transient analysis using(Jacobi,Gauss Seidal,SOR)

The 2D transient heat conduction equation is given as:

`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`

Discertizing the equation using central difference scheme in space and FDS in time we get

`(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)=(T_(i-1,j)^n-2*T_(i,j)^n+T_(i+1,j)^n)/(Deltax^2)+(T_(i,j-1)^n-2*T_(i,j)^n+T_(i,j+1)^n)/(Deltay^2)`

Rearraning the above equation we get (Explicit scheme)

`T_(i,j)^(n+1)=(1-2k1-2k2)T_(i,j)^n+k1(T_(i-1,j)^n+T_(i+1,j)^n)+k2(T_(i,j-1)^n+T_(i,j+1)^n)`

k1=`alpha(Deltat)/(Deltax^2)` k2=`alpha(Deltat)/(Deltay^2)`

Similarly for Implicit Scheme

`T_(i,j)^(n+1)=(T_(i,j)^n+k1(T_(i-1,j)^(n+1)+T_(i+1,j)^(n+1))+k2(T_(i,j-1)^(n+1)+T_(i,j+1)^(n+1)))/(1+2k1+2k2)`

Matlab code for Steady state analysis using Jacobi.

% steady state 2D heat conduction
close all
clear all
clc
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=100*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
k1=0.5*(dy^2/(dy^2+dx^2));
k2=0.5*(dx^2/(dy^2+dx^2));
iterative_solver=1;
error=9e9;
tol=1e-4;
%jacobi
if iterative_solver==1
    jacobi=1;
    while(error>tol)
      for j=2:ny-1  
    for i=2:nx-1
      T(i,j)=k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1));
    end
      end
    error=max(max(abs(Told-T)));
    Told=T;
    jacobi=jacobi+1;
    contourf(x,y,T)
    [C,h]=contourf(x,y,T)
    clabel(C,h)
  
    xlabel('X axis')
    ylabel('Y axis')
    title_text=sprintf("jacobi iteration number=%d",jacobi)
    title(title_text);
    colorbar
      pause(0.03)
    
    end
end

Matlab code for Steady state analysis using Gauss Seidal.

% steady state 2D heat conduction
close all
clear all
clc
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=100*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
k1=0.5*(dy^2/(dy^2+dx^2));
k2=0.5*(dx^2/(dy^2+dx^2));
iterative_solver=2;
error=9e9;
tol=1e-4;
% Gauss seidal
if iterative_solver==2
    gs_iter=1;
    while(error>tol)
      for j=2:ny-1  
    for i=2:nx-1
      T(i,j)=k1*(T(i-1,j)+Told(i+1,j))+k2*(T(i,j-1)+Told(i,j+1));
    end
      end
    error=max(max(abs(Told-T)));
    Told=T;
    gs_iter=gs_iter+1;
    contourf(x,y,T);
   
    [C,h]=contourf(x,y,T);
        clabel(C,h);
          
    xlabel('X axis')
    ylabel('Y axis')
    title_text=sprintf("Gauss iteration number=%d",gs_iter);
    title(title_text);
    colorbar
      pause(0.03)
    
    end
end

Matlab code for Steady state analysis using SOR.

% steady state 2D heat conduction
close all
clear all
clc
nx=10;
ny=10;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=100*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
k1=0.5*(dy^2/(dy^2+dx^2));
k2=0.5*(dx^2/(dy^2+dx^2));
iterative_solver=3;
error=9e9;
tol=1e-4;
w=1.5;
% SOR
if iterative_solver==3
    SOR_iter=1;
    while(error>tol)
      for j=2:ny-1  
    for i=2:nx-1
      Tgauss(i,j)=k1*(T(i-1,j)+Told(i+1,j))+k2*(T(i,j-1)+Told(i,j+1));
      T(i,j)=(1-w)*Told(i,j)+w*Tgauss(i,j);
    end
      end
    error=max(max(abs(Told-T)));
    Told=T;
    SOR_iter=SOR_iter+1;
    contourf(x,y,T);
   
    [C,h]=contourf(x,y,T);
        clabel(C,h);
          
    xlabel('X axis')
    ylabel('Y axis')
    title_text=sprintf("SOR iteration number=%d",SOR_iter);
    title(title_text);
    colorbar
      pause(0.03)
    
    end
end

Results:

Steady state analysis using Jacobi.

Steady state analysis using Gauss seidal.

Steady state analysis using SOR.

Conclusion-The no of grid points taken for steady state analysis was 10. So from all the above graph it can be concluded that,

1. For steady state analysis jacobi method took highest no of iteration(215) to converge followed by gauss-seidal(115) and SOR(30).

2. We can see that SOR method took very less iteration to converge.This is mainly due to introduction of ralaxation factor(`w`). `w<1`is called under-relaxation and `w>1` is called over-relaxation.

3. For SOR method to converge in 30 iteration the relaxation factor used was 1.5(over-relaxation).However w>1.5 increased the no of iteration and below 1.5 also increased the no of iterations.So optimum w should be used for faster convergence.

Matlab code for transient analysis using explicit scheme.

%2D transient conduction explicit
close all
clear all
clc
nx=20;
ny=20;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=100*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
alpha=1.25;
 t=0.1;
dt=5e-4;
nt = t/dt;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
iterative_solver=1;
error=9e9;
tol=1e-4;
time=0;
    for k=1:nt
    
      for j=2:ny-1  
    for i=2:nx-1
      T(i,j)=(1-2*k1-2*k2)*Told(i,j)+k1*(Told(i+1,j)+Told(i-1,j))+k2*(Told(i,j-1)+Told(i,j+1));
    end
      end
 time=time+dt;
   Told=T;
  contourf(x,y,T)
    [C,h]=contourf(x,y,T);
    clabel(C,h)
  
    xlabel('X axis')
    ylabel('Y axis')
     title_text=sprintf("total time=%g",time);
    title(title_text);
    colorbar
      pause(0.03) 
    
    end

Matlab code for transient analysis using Jacobi.

%2D transient heat conduction implicit
close all
clear all
clc
nx=20;
ny=20;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=300*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
Tprev=T;
alpha=1.25;
% t=0.4;
dt=1e-3;
nt = 100;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1+2*k1+2*k2)^-1;
iterative_solver=1;
% error=9e9;
tol=1e-4;
w=1.1;

%jacobi
if iterative_solver==1
     jacobi_iter=1;
time=0;

    for k=1:nt
       error=9e9;
    while(error>tol)
      for j=2:ny-1  
    for i=2:nx-1
      T(i,j)=term1*Tprev(i,j)+k1*term1*(Told(i+1,j)+Told(i-1,j))+k2*term1*(Told(i,j-1)+Told(i,j+1));
    end
      end
    error=max(max(abs(Told-T)));
    Told=T;
  jacobi_iter=jacobi_iter+1;

  
    
    end
   Tprev=T;
    time=time+dt;
   contourf(x,y,T);
     [C,h]=contourf(x,y,T);
    clabel(C,h)
   
    xlabel('X axis')
    ylabel('Y')
    title_text=sprintf('total time=%g and no of iter=%d',time, jacobi_iter);
    title(title_text);
    colorbar
      pause(0.03) 
   
    end
end

Matlab code for transient analysis using Gauss Seidal.

%2D transient heat conduction implicit
close all
clear all
clc
nx=20;
ny=20;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=300*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
Tprev=T;
alpha=1.25;
% t=0.4;
dt=1e-3;
nt = 100;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1+2*k1+2*k2)^-1;
iterative_solver=2;
% error=9e9;
tol=1e-4;
w=1.1;

 %Gauss Seidal
if iterative_solver==2
     gs_iter=1;
time=0;

    for k=1:nt
       error=9e9;
    while(error>tol)
      for j=2:ny-1  
    for i=2:nx-1
      T(i,j)=term1*Tprev(i,j)+k1*term1*(T(i+1,j)+Told(i-1,j))+k2*term1*(T(i,j-1)+Told(i,j+1));
    end
      end
    error=max(max(abs(Told-T)));
    Told=T;
  gs_iter=gs_iter+1;

     
    end
   Tprev=T;
    time=time+dt;
   contourf(x,y,T);
     [C,h]=contourf(x,y,T);
    clabel(C,h)
   
    xlabel('X axis')
    ylabel('Y')
    title_text=sprintf('total time=%g and no of iter=%d',time, gs_iter);
    title(title_text);
    colorbar
      pause(0.03) 
   
    end
  
       
end

Matlab code for transient analysis using SOR.

%2D transient heat conduction implicit
close all
clear all
clc
nx=20;
ny=20;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
T=300*ones(nx,ny);
T(:,1)=400;
T(:,nx)=800;
T(1,:)=900;
T(ny,:)=600;
Told=T;
Tprev=T;
alpha=1.25;
% t=0.4;
dt=1e-3;
nt = 100;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1+2*k1+2*k2)^-1;
iterative_solver=3;
% error=9e9;
tol=1e-4;
w=1.1;

%SOR
if iterative_solver==3
     SOR_iter=1;
time=0;

    for k=1:nt
       error=9e9;
    while(error>tol)
      for j=2:ny-1  
    for i=2:nx-1
         Tgauss(i,j)=term1*(Tprev(i,j)+k1*(T(i+1,j)+Told(i-1,j))+k2*(T(i,j-1)+Told(i,j+1)));
      T(i,j)=(1-w)*Tprev(i,j)+w*Tgauss(i,j);
    end
      end
    error=max(max(abs(Told-T)));
    Told=T;
  SOR_iter=SOR_iter+1;

      end
   Tprev=T;
    time=time+dt;
   contourf(x,y,T);
     [C,h]=contourf(x,y,T);
    clabel(C,h)
   
    xlabel('X axis')
    ylabel('Y')
    title_text=sprintf('total time=%g and no of iter=%d',time, SOR_iter);
    title(title_text);
    colorbar
      pause(0.03) 
   
    end
end

Result:

Transient analysis using explicit scheme.

Transient analysis using Jacobi.

Transient analysis using GS:

.

Transient analysis using SOR.

Conclusion-

1. For transient analysis 20 grid points were taken.The simulation was run for 0.1 seconds.

2. In case of explicit scheme the solution was obtained by time marching.The convergence criteria was not checked.

3. It was seen that the explicit scheme was dominated by Courant-Friedrichs-Lewy(CFL) condition which determines the stability of a solution.

CFL condition for 2D heat equation is:

`alpha(Deltat)/(Deltax^2)+alpha(Deltat)/(Deltay^2)<=0.5`

or `alpha(Deltat)/(Deltax^2)<=0.25` when `Deltax=Deltay`

4. For the above simulation dt=0.0005(CFL no=0.22) which is within the limit so the solution is stable.

CFL no>0.25 showed unstable solution and the solution was blown off.

4. In case of implicit scheme the highest iteration was seen in jacobi(2281) followed by SOR(2258) and Gauss Siedal(1890).In transient analysis the SOR scheme showed faster convergence for under relaxed condition.CFL number had no impact on stability incase of implicit scheme.

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.