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

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

PPradeep R27 May 2020

OBJECTIVE:

To solve the traditional 2D heat conduction equation for steady-state and transient conditions using the iterative solvers such as the Jacobi, Gauss seidel and Successive over relaxation.

CONDITIONS FOR A 2D DOMAIN:

Assume that the domain is a unit square.

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

Here we assume the number of nodes as 10.

Boundary conditions for steady and transient case

  • Left:400K

  • Top:600K

  • Right:800K

  • Bottom:900K

At Corners,

  • Left Top =500K

  • Left Bottom = 650K

  • Right Top = 700K

  • Right Bottom = 850K

All other nodes are given an initial temperature of 300K

The logic behind iterative solvers used here is:

No need of time loop for Steady-state.

LOGIC:

STEADY STATE:

Now in heat transfer steady-state means, the temperature of the body does not vary with time. Amount of heat enters from one side of the body the same amount of heat leaves the body from the other side to maintain the temperature constant or steady.

`(del^2 T) / (del x^2) + (del^2 T)/ (del y^2) = 0`

By discretizing the equation above it becomes,

`T_p = 1/k [(T_L + T_R)/( Delta x^2) + (T_T + T_B)/( Delta y^2) ]`

Here the value of k is,

`K = 2((Delta x^2 + Deltay^2) /(Delta x^2 Deltay^2) )`

If we further simply this the equation for code becomes

`0.25 * (T_L + T_R + T_T +T_B)`

which can also be written as,

`T_(i,j) = 0.25* (T_(i-1,j) + T_(i+1,j) + T_(i,j-1) + T_(i,j+1))`

MAIN CODE FOR STEADY STATE JACOBI:

clear all
close all
clc

nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

tol = 1e-4;


t = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;

t_initial = t;
told = t;
error = 9e9;
iterative_solver = 1;
if iterative_solver == 1
  
    jacobi_iter = 1;
    while(error > tol)
        for i = 2:nx-1 
        for    j = 2:ny-1
            
            t(i,j) = 0.25*(told(i-1,j)+ told(i+1,j) + told(i,j-1) + told(i,j+1));
        end
        end
         error = max(max(abs(told-t)));
    told = t;
    jacobi_iter = jacobi_iter + 1;
    figure(1)
    contourf(x,y,t)
    clabel(contourf(x,y,t))
    colorbar
    set(gca, 'ydir', 'reverse')
    title_text = sprintf('No of Jacobi iterations = %d', jacobi_iter);
    title(title_text)
    pause(0.03)  
    

    end
end

RESULT FOR STEADY STATE JACOBI:

MAIN CODE FOR STEADY STATE GAUSS SEIDEL:

clear all
close all
clc

nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

tol = 1e-4;


t = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;

t_initial = t;
told = t;
error = 9e9;
iterative_solver = 2;
if iterative_solver == 2
  
    GS_iter = 1;
    while(error > tol)
        for i = 2:nx-1
           for j = 2:ny-1
            
            t(i,j) = 0.25*(t(i-1,j)+ told(i+1,j) + t(i,j-1) + told(i,j+1));
           end
        end
         error = max(max(abs(told-t)));
    told = t;
    GS_iter = GS_iter + 1;
   
    figure(1)
    contourf(x,y,t)
    clabel(contourf(x,y,t))
    colorbar
    set(gca, 'ydir', 'reverse')
    title_text = sprintf( 'Gauss seidal iterations = %d' ,GS_iter);
    title(title_text)
    pause(0.03)  
    

    end
end

RESULT FOR STEADY STATE GAUSS SEIDEL:

MAIN CODE FOR STEADY STATE SOR:

A factor called Omega is introduced for SOR iterations. When Omega is less than one it is under relaxed when it is greater than 1 it is over relaxed. It can't be 1 because we will again get Jacobi or gauss seidel again depending on which solver we use.

clear all
close all
clc

nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

tol = 1e-4;

omega = 1.4;


t = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;

t_initial = t;
told = t;
error = 9e9;
iterative_solver = 3;
if iterative_solver == 3
  
    sor_iter = 1;
    while(error > tol)
        for i = 2:nx-1
        for j = 2:ny-1
            
            t(i,j) = omega*(0.25*(t(i-1,j)+ told(i+1,j) + t(i,j-1) + told(i,j+1))) + (1 - omega)*told(i,j);
        end
        end
    error = max(max(abs(told-t)));
    told = t;
    sor_iter = sor_iter + 1;
   
    figure(1)
    contourf(x,y,t)
    clabel(contourf(x,y,t))
    colorbar
    set(gca, 'ydir' , 'reverse')
    title_text = sprintf( 'SOR iterations = %d' ,sor_iter);
    title(title_text)
    pause(0.03)  
    
 
    end
end

RESULT FOR STEADY STATE SOR:

EQUATION OF UNSTEADY STATE EXPLICIT FUNCTION:

Transient state: In this, the temperature in a body keeps on changing with time or body is not in equilibrium. This kind of situations arises due to many factors like a sudden change in environmental temperature, changing properties of the body and many others.

There is no need for iterative solvers in explicit form sin all the variables on R.H.S is known.

`(delT )/( delt) = alpha (del^2T) /( del x^2) + (del^2T) / (del y^2 )`

If we discretize the equations it becomes,

`(T_P^(n+1) - T_P^n )/ (Deltat) = alpha (T_L - 2T_P + T_R) /( del x^2) + (T_T - 2T_P + T_B) / (del y^2 )`

If we simplify it further,

`(T_P^(n+1) )= T_P^n + k1 (T_L - 2T_P + T_R)^n + k2(T_T - 2T_P + T_B) ^n`

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

For the better understanding for code,

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

CFL CRITERIA:

For Explicit functions we have criteria such that it reaches a steady solution.

CFL <= 0.5 which means k1 and k2 should be less than 0.5 in order to reach steady state.

MAIN CODE FOR EXPLICIT FUNCTION:

clear all
close all
clc

nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

tol = 1e-4;

dt = 1e-3;

t = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;


told = t;



k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));


tol = 1e-4;
iter = 1;


if  iter == 1

for k = 1:1400
  error = 1;
    while(error>tol)
    for i = 2:nx-1
    for j = 2:ny-1
        
       A = told(i,j);
       B = k1*(told(i+1,j) -2*(told(i,j)) + told(i-1,j));
       C = k2*(told(i,j+1) -2*(told(i,j)) + told(i,j-1));
       
       t(i,j) = A + B + C;
        
        
    end
    end
        error = max(max(abs(told-t)));
        told = t;
        iter = iter + 1;
        
        figure(1)
        contourf(x,y,t);
        clabel(contourf(x,y,t));
        set(gca, 'ydir', 'reverse')
        colormap(jet)
        colorbar
        title_text = sprintf('No of iter = %d' , iter);
        title(title_text)
        pause(0.03)      
         
        
        
    end     
end
end

RESULTS:

FORMULA FOR IMPLICIT FUNCTIONS:

`(delT )/( delt) = alpha (del^2T) /( del x^2) + (del^2T) / (del y^2 )`

If we discretize the equations it becomes,

`(T_P^(n+1) - T_P^n )/ (Deltat) = alpha (T_L - 2T_P + T_R) /( del x^2) + (T_T - 2T_P + T_B) / (del y^2 )`

If we simplify it further,

`(T_P^(n+1) )=(1)/(1+k1+k2) (T_P^n + k1 (T_L + T_R)^(n+1) + k2(T_T + T_B) ^(n+1))`

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

Alpha is the diffusivity term.

For the better understanding for code,

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

MAIN CODE FOR IMPLICIT JACOBI:

clear all
close all
clc

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 = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;

told = t;
t_prev = t;

 dt = 1e-3;

 tol = 1e-4;
 k1 = 1.1*(dt/(dx^2));
 k2 = 1.1*(dt/(dy^2));


    jacobi_iter = 1;
for k = 1:1400
    error = 9e9;
 while(error > tol)
       
    for i = 2:nx-1
    for j = 2:ny-1
        
       term1 = 1/(1+ (2*k1) + (2*k2));
       term2 = k1*term1;
       term3 = k2*term1;
       H = (told(i-1,j)) + (told(i+1,j));
       V = (told(i,j+1)) + (told(i,j-1));
       
       t(i,j) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
    end 
    end 
        error = max(max(abs(told-t)));
        told = t; 
        jacobi_iter = jacobi_iter + 1; 
 end
  t_prev = t; 
 
    
 figure(1)
 contourf(x,y,t)
 clabel(contourf(x,y,t))
 colorbar
 set(gca, 'ydir', 'reverse')
 title_text = sprintf('Jacobi iterations = %d' , jacobi_iter);
 title(title_text)
  
   

end

RESULT FOR IMPLICIT JACOBI:

MAIN CODE FOR IMPLICIT GAUSS SEIDEL:

clear all
close all
clc

nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

t = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;

told = t;
t_prev = t;



dt = 1e-3;
  
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

error = 9e9;
tol = 1e-4;


GS_iter = 1;
for k = 1:1400
error = 1;
    while(error > tol)
       
    for i = 2:nx-1
    for j = 2:ny-1
     
        
       term1 = 1/(1+ (2*k1) + (2*k2));
       term2 = k1*term1;
       term3 = k2*term1;
       H = t(i-1,j) + told(i+1,j);
       V = t(i,j+1) + told(i,j-1);
       
       t(i,j) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
    end 
    end 
        error = max(max(abs(told-t)));
        told = t;        
        GS_iter = GS_iter + 1; 
 end
  
 t_prev = t;
 
 figure(1)
 contourf(x,y,t);
 set(gca, 'ydir', 'reverse')
 clabel(contourf(x,y,t));
 colorbar
 title_text = sprintf('Gauss seidal iterations = %d' , GS_iter);
 title(title_text)
end

RESULT FOR IMPLICIT GAUSS SEIDEL:

MAIN CODE FOR IMPLICIT SOR:

clear all
close all
clc

nx = 10;
ny = nx;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = dx;

t = 300*ones(nx,ny);

t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;


t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;

told = t;
t_prev = t;
t_sor = t;


dt = 1e-3;


k1 = 1.1*(dt)/(dx^2);
k2 = 1.1*(dt)/(dy^2);

omega = 1.1;


tol = 1e-4;


SOR_iter = 0;
for k = 1:1400
    error = 1;
 while(error > tol)
       
    for i = 2:nx-1
    for j = 2:ny-1
       
        
       term1 = 1/(1+ (2*k1) + (2*k2));
       term2 = k1*term1;
       term3 = k2*term1;
       H = t(i-1,j) + told(i+1,j);
       V = t(i,j+1) + told(i,j-1);
       
       t(i,j) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
       t_sor(i,j) = ((1-omega)*told(i,j)) + (t(i,j)*omega);
    end 
    end 
        error = max(max(abs(told-t)));
        told = t_sor;        
        SOR_iter = SOR_iter + 1; 
 end
  
 t_prev = t_sor;
  
 figure(1)
 contourf(x,y,t);
 set(gca, 'ydir', 'reverse')
 clabel(contourf(x,y,t));
 colorbar
 title_text = sprintf('SOR iterations = %d' , SOR_iter);
 title(title_text)

  
end

RESULT FOR IMPLICIT SOR:

CONCLUSIONS:

  • It is clearly evident that the steady-state takes less no of iterations because of the lack of time derivative and the flow properties are constant.
  • The disadvantage of Explicit is, if CFL goes beyond 0.5, the solution becomes unstable. Meanwhile Implicit doesn't have to follow the CFL criteria.
  • Here for SOR, gauss seidel is used for both the cases and it clear that SOR's iteration is less in both steady-state as well as in implicit with a relaxation factor greater than 1 conveying an over-relaxed solution.
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.