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

Solving the steady and unsteady 2D heat conduction problem

PPradyumna Kodancha2 Oct 2020

Aim: To solve the heat conduction Equations and plot graphs.

Objective:

Solving 2D heat conduction equation using different iterative techniques.

1) Explicit Solver (unsteady state)

2) Implicit Solver (steady and unsteady state)

(a) Jacobi method

(b) Gauss Siedel method

(c) Successive Over Relaxation method.

Problem description:

1. Assume domain is unit square.

2. Gridpoint is the same in both X & Y direction.

3. Boundary condition

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

The absolute error criteria is 1e-4, and the initial domain condition at t=0 temperature=300K.

Unsteady / Transient state solution:

Equation is

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

After applying forward difference scheme to time derivative, central difference scheme to spaciatial derivative, we have

`(T^(n+1))underset(i,j)=((T^n)underset(i,j)+(K1((T^n)underset(i-1,j)+(T^n)underset(i+1,j))+K2((T^n)underset(i,j-1)+(T^n)underset(i,j+1)))^(n+1))/(1+2K1+2K2)`

Where

K1=`alpha(Deltat)/(Deltax)^2`

K2=`alpha(Deltat)/(Deltay)^2`

`Deltax=Deltay`

Initial Temperature profile is attached below.

1) Explicit Scheme program

clear all
close all
clc

nx=50;%grid
ny=nx;
L=1;%length of domain
x=linspace(0,1,nx);%grid spacing
y=linspace(0,1,ny);%grid spacing
dx=L/(nx-1);
dy=L/(ny-1);
alpha=1e-3;%conduction coefficient
t=500;%total time
dt=0.001;%time step
nt=t/dt;%time for loop
conv=1e-4;%convergence criteria
T=300*ones(nx,ny);%initialization of domain
%boundary condition
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
%boundary condition for corner nodes
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
%cfl number or stability criteria
v1=alpha*(dt/(dx^2));
v2=alpha*(dt/(dy^2));
cfl=v1+v2;
%mesh grid 
[xx,yy]=meshgrid(x,y);
%initilization
Told=T;
tic; %time calculator
c=0;
t1=0;
A1=0;

for k=1:nt  %time loop
    for i=2:nx-1 %space loop
        for j=2:ny-1 %space loop
            term1=Told(i,j);
            term2=v1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
            term3=v2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
            T(i,j)=term1+term2+term3;
        end
    end
    A1=T;
    %convergence checker
    t1 = t1+1;
    if (max(max(abs(A1-Told)))<conv)
       break;
    end
    Told=T;
    c=c+1;
end
%plotting graph
 figure(1)
    [G,A]=contourf(xx,yy,T);
    time=toc;
    colorbar;
    colormap(jet);
    clabel(G,A);
    title(sprintf('2D Unsteady Explicit heat conduction n Number of iterations:%d(time=%0.4fs)n Computation Time = %.4fs n CFL number = %.4f',t1,k*dt,time,cfl));
    xlabel('xx');
    ylabel('yy');

The results are attached.

2) Implicit Scheme program

clear all
close all
clc

nx=50;%grid
ny=nx;
L=1;%length of domain
x=linspace(0,1,nx);%grid spacing
y=linspace(0,1,ny);%grid spacing
dx=L/(nx-1);
dy=L/(ny-1);
alpha=1e-2;%conduction coefficient
dt=0.1;%time step
tol=1e-4;

T=300*ones(nx,ny);%initialization of domain
%boundary condition
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
%boundary condition for corner nodes
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
omega= (2/(1+sin(pi*dx))); %relaxation value
%cfl number or stability criteria
v1=alpha*(dt/(dx^2));
v2=alpha*(dt/(dy^2));
%mesh grid 
[xx,yy]=meshgrid(x,y);
%initilization
Told=T;
T_prev=Told;
c=0;
t1=0;
solver= input('Enter the number for the solver method :n 1) Jacobi n 2) Gauss siedel n 3) SORn');
tic;%time calculator
for k=1:400000 %time loop
   if solver==1
        term1=(1+2*v1+2*v2)^(-1);
        term2=v1*term1;
        term3=v2*term1;
        error=9e9;
        while(error>tol) %iteration conversion checker
            for i=2:nx-1 %space loop
                for j=2:ny-1 %space loop
                    h=(Told(i-1,j)+Told(i+1,j));
                    v=(Told(i,j-1)+Told(i,j+1));
                    T(i,j)= term1*T_prev(i,j)+term2*h+term3*v;
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            c=c+1;
        end
    end
    if solver==2
        term1=(1+2*v1+2*v2)^(-1);
        term2=v1*term1;
        term3=v2*term1;
        error=9e9;
        while(error>tol) %iteration convergence checker
            for i=2:nx-1 %space loop
                for j=2:ny-1 %space loop
                    h=(T(i-1,j)+Told(i+1,j));
                    v=(T(i,j-1)+Told(i,j+1));
                    T(i,j)= term1*T_prev(i,j)+term2*h+term3*v;
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            c=c+1;
        end
    end
    if solver==3
        term1=(1+2*v1+2*v2)^-1;
        term2=v1*term1;
        term3=v2*term1;
        error=9e9;
        while(error>tol) %iteration conversion checker
            for i=2:nx-1 %space loop
                for j=2:ny-1 %space loop
                    h=(omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
                    v=(omega*T(i,j-1)+(1-omega)*Told(i,j-1)+Told(i,j+1));
                    T(i,j)= term1*T_prev(i,j)+(term2*h+term3*v);
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            c=c+1;
        end
   end
  
   conv = max(max(abs(T-T_prev)));
   T_prev=T;
   t1 = t1+1;
   if(conv<tol) %time convergence checker
       break;
   end
  
end
%plotting graph
 figure(2);
   str_array=["Jacobi","Gauss-siedel","SOR"];
   [G,A]=contourf(xx,yy,T);
   Simulation_time=toc;
   colorbar;
   colormap(jet);
   clabel(G,A);
   title(sprintf('2D Unsteady Implicit heat conduction solver method = %s n  Iteratons = %d n Computation Time = %.4fs n convergence step = %d(time=%.4fs)',str_array(solver),c,Simulation_time,t1,t1*dt));
   xlabel('xx');
   ylabel('yy');

The results are attached.

Steady-state solution:

important modification is `(delT)/(delt)=0`

Equation reduced to

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

By `Deltax=Deltay`

`Tunderset(i,j)=1/4(Tunderset(i-1,j)+Tunderset(i+1,j)+Tunderset(i,j-1)+Tunderset(i,j+1))`

Above equation is also called 5-point formula.

steady-state implicit program

clear all
close all
clc

nx=80;%grid
ny=nx;
L=1;%length of domain
x=linspace(0,1,nx);%grid spacing
y=linspace(0,1,ny);%grid spacing
dx=L/(nx-1);
dy=L/(ny-1);

error=9e9;
tol=1e-4;
k=2*(1/(dx^2)+1/(dy^2)); %constant 

T=300*ones(nx,ny);%initialization of domain
%boundary condition
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
%boundary condition for corner nodes
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
omega= (2/(1+sin(pi*dx))); %relaxation value
%mesh grid 
[xx,yy]=meshgrid(x,y);
%initilization
Told=T;
c=0;
solver= input('Enter the number for the solver method :n 1) Jacobi n 2) Gauss siedel n 3) SORn');
tic;%time calculator

   if solver==1
 
        while(error>tol)%iteration conversion checker
            for i=2:nx-1 %space loop
                for j=2:ny-1 %space loop
                    term1=(Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
                    term2=(Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
                    T(i,j)= term1+term2;
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            c=c+1;
        end
    end
    if solver==2
       
        while(error>tol)%iteration conversion checker
            for i=2:nx-1 %space loop
                for j=2:ny-1 %space loop
                    term1=(T(i-1,j)+Told(i+1,j))/(k*(dx^2));
                    term2=(T(i,j-1)+Told(i,j+1))/(k*(dy^2));
                    T(i,j)= term1+term2;
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            c=c+1;
        end
    end
    if solver==3
        
        while(error>tol)%iteration conversion checker
            for i=2:nx-1 %space loop
                for j=2:ny-1 %space loop
                    term1=(T(i-1,j)+Told(i+1,j))/(k*(dx^2));
                    term2=(T(i,j-1)+Told(i,j+1))/(k*(dy^2));
                    T(i,j)= (1-omega)*Told(i,j)+omega*(term1+term2);
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            c=c+1;
        end
   end
%plotting

figure(1);
str_array=["Jacobi","Gauss-siedel","SOR"];
[G,A]=contourf(xx,yy,T);
Simulation_time=toc;
colorbar;
colormap(jet);
clabel(G,A);
title(sprintf('2D Unsteady Implicit heat conduction n solver method = %s n  Iteratons = %d n Computation Time = %.4fs',str_array(solver),c,Simulation_time));
xlabel('xx');
ylabel('yy');

The results are attached.

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.