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

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

AAthiyaman R6 Jan 2021

STEADY STATE:

we have to solve the study state equation using different iterative solvers.they are,

  • Jacobi method.
  • Gauss-seidel method.
  • Successive over relaxation method.

Heat Conduction equation is ,

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

by discretize the above equation we get,

`T_P^(n+1)=1/K[(T_L+T_R)/(Deltax^2)+(T_B+T_T)/(Deltay^2)]^n`--------(a)

here , k value is

`K=2**[(Deltax^2+Deltay^2)/(Deltax^2Deltay^2)]`

equation (a) is used for the jacobi method .

For gauss-seidal method ,the equation is

`T_P^(n+1)=1/K[(T_L+T_R)/(Deltax^2)+(T_B+T_T)/(Deltay^2)]^(n+1)`

For successive over relaxation method, the equation is

`T_P^(n+1)=T_P^n(1-lambda)+lamda**[(T_L+T_R)/(Deltax^2)+(T_B+T_T)/(Deltay^2)]^(n+1)`

STEADY STATE HEAT CONDUCTION PROGRAM:

MAIN PROGRAM:

%Solving the steady state equation

clear all
close all
clc

%inputs
L=1;
nx=ny=10;
x=linspace(0,L,nx);
y=linspace(0,L,nx);
dx=dy=x(2)-x(1);
[xx yy] = meshgrid(x,y);

%B.C
T=ones(nx,ny);
T(1,:)=900;%B.s
T(end,:)=600;%T.S
T(:,1)=400;%L.S
T(:,end)=800;%R.S
Told=T;
tol=1e-4;
error=9e9;
k=2*((dx^2+dy^2)/(dx^2*dy^2));
%defining the solvers
for iterative_solver=[1 2 3]
if iterative_solver==1
    tic;
    T_jacobi= jacobi_iter(error,tol,xx,yy,nx,ny,dx,dy,T,k,Told);
    time_jacobi_in_sec=toc
elseif iterative_solver==2
    tic;
    T_gauss= gauss_seidal_iter(error,tol,x,y,nx,ny,dx,dy,k,T,Told);
    time_gauss_in_sec=toc
elseif iterative_solver==3
    tic;
    T_sor= sor_iter(error,tol,x,y,nx,ny,dx,dy,k,T,Told);
    time_sor_in_sec=toc
end
endfor

FUNCTION PROGRAMS FOR SOLVING STEADY STATE EQUATION:

1.JACOBI METHOD:

%function for jacobi method to solve to steady state eqauation
function T= jacobi_iter(error,tol,xx,yy,nx,ny,dx,dy,T,k,Told)
     
     jacobi_iter=1;
     while(error > tol)
     for i=2:nx-1
         for j=2:ny-1
          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;
     jacobi_iter=jacobi_iter+1;
     figure(1)  
     [c,d]=contourf(xx,yy,T);
     colormap(jet)
     colorbar
     title_text=sprintf('Jacobi iteration number=%d @steady state',jacobi_iter);
     title(title_text)
     xlabel('X-axis')
     ylabel('y-axis')
     clabel(c,d)
     pause(0.03)
     end
     fprintf('No of Iteration in jacobi=%d',jacobi_iter);
end

2.GAUSS SEIDAL METHOD:

%function for gauss_seidal method to solve steady state method
function gauss_seidal= gauss_seidal_iter(error,tol,x,y,nx,ny,dx,dy,k,T,Told)
     GS_iter=1;
   
     while(error > tol)
     for i=2:nx-1
         for j=2:ny-1
         T(i,j)=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
         end
     end
     error=max(max(abs(Told-T)));
     Told=T;
     GS_iter=GS_iter+1;
     gauss_seidal=figure(2)   
     [c,d]=contourf(x,y,T);
     colormap(jet)
     colorbar
     title_text=sprintf('GS iteration number=%d @steady state',GS_iter);
     title(title_text)
     xlabel('X-axis')
     ylabel('Y-axis')
     clabel(c,d)
     pause(0.03)
     end
     fprintf('No of Iteration in Gs=%d',GS_iter);
end

3.SUCCESSIVE OVER RELAXATION METHOD:

%function for SOR method to solve steady state method
function T= sor_iter(error,tol,x,y,nx,ny,dx,dy,k,T,Told)
     sor_iter=1;
     w=1.7;
     while(error > tol)
     for i=2:nx-1
         for j=2:ny-1
         T_gs=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
         T(i,j)= w*(T_gs)+(1-w)*(Told(i,j));
         end
     end
     error=max(max(abs(Told-T)));
     Told=T;
     sor_iter=sor_iter+1;
     figure(3)   
     [c,d]=contourf(x,y,T);
     colormap(jet)
     colorbar
     title_text=sprintf('SOR iteration number=%d @steady state',sor_iter);
     title(title_text)
     xlabel('X-axis')
     ylabel('Y-axis')
     clabel(c,d)
     pause(0.03)
     end
     fprintf('No of Iteration in SOR=%d',sor_iter);
end

STEADY STATE HEAT CONDUCTION RESULTS:

1.JACOBI METHOD:

2.GAUSS SEIDEL :

3.SUCCESSIVE OVER RELAXATION :

TRANSIENT/UNSTEADY STATE:

We have to solve the transient heat conduction equation using the given iterative methods for both impicit and explicit forms.

Transient state heat conduction equation is ,

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

1.IMPLICIT FORM:

Implicit form of jacobi method for transient heat conduction equation is ,

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

Implicit form of Gauss seidel method for transient heat conduction equation is ,

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

Implicit form of Successive over relaxation method for transient heat conduction equation is ,

`T_P^(n+1)=T_P^n**(1-lambda)+lambda**[T_P^n+(k1**(T_L^(n+1)+T_R^n)+k2(T_T^n+T_B^(n+1)))/(1+2**k1+2**k2)]`

IMPLICIT METHOD MAIN PROGRAM:

%implicit methods
clear all
close all
clc

% Initial Conditions
L = 1; % Unit square length.
nx = ny = 10; % grid points along X axis and y axis

x = linspace(0, L, nx);
dx = L/(nx-1);
y = linspace(0, L, ny);
dy = L/(ny-1);

error = 9e9;
tol = 1e-5;
t = 5;
dt = 1e-2;
alpha = 0.4;

% Assumption for Unsteady state constants
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);

% Boundary Conditions
T = ones(nx, ny);
T(1,:)=900;%B.s
T(end,:)=600;%T.S
T(:,1)=400;%L.S
T(:,end)=800;%R.S

% Edge Refining
T(1,1) = (400+900)/2;
T(1,ny) = (900+800)/2;
T(nx,1) = (400+600)/2;
T(nx,ny) = (600+800)/2;

for iter_method = [1 2 3]
    if iter_method == 1
       T_jacobi = implicit_1_jacobi( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 );
    else if iter_method == 2
            T_gs = implicit_gauss_1_seidel( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 );
        else if iter_method == 3
                time = implicit_1_sor( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 );
            end
        end
    end
end

FUNCTION PROGRAMS FOR SOLVING TRANSIENT STATE EQUATION IN IMPLICIT FORM:

1.JACOBI METHOD:

function T = implicit_1_jacobi( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 )
% Implicit Jacobi method for Transient state

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
jac_iter = 1;


tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (T_old(i-1, j) + T_old(i+1, j));
                term2 = (T_old(i, j-1) + T_old(i, j+1));
                T(i,j) = (To(i,j) + (k1*term1) + (k2*term2))/(1 + 2*k1 + 2*k2);
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        jac_iter = jac_iter + 1;
    end
    To = T;
end
iter_time = toc;
figure(1)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X axis');
ylabel('Y axis');
colorbar
title_text = sprintf('Jacobi Method Iterations = %d @eimplicit unsteady ', jac_iter);
title(title_text)
end

2.GAUSS SEIDEL METHOD:

function T = implicit_gauss_1_seidel( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 )
% Implicit Gauss Seidel method for Transient state

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
gs_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = (To(i,j) + (k1*term1) + (k2*term2))/(1 + 2*k1 + 2*k2);
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        gs_iter = gs_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(2)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X axis');
ylabel('Y axis');
colorbar
title_text = sprintf('Gauss seidel Method Iterations = %d @implicit unsteady ', gs_iter);
title(title_text)
end

3.SUCCESIVE OVER RELAXATION METHOD:

function T = implicit_1_sor( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 )
% Implicit Successive Over-Relaxation Method for Transient State

lambda = 1.09; % Over Relaxation factor
meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
sor_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = (1-lambda)*T_old(i,j) + (lambda*((To(i,j) + (k1*term1) + (k2*term2))/(1 + 2*k1 + 2*k2)));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        sor_iter = sor_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(3)
[a, b] = contourf(x, y, T);
colormap(jet)
colorbar
clabel(a, b);
xlabel('X axis');
ylabel('Y axis');
title_text=sprintf('Successive Over-Relaxation Method Iterations = %d  @implicit unsteady', sor_iter);
title(title_text)
end

IMPLICIT METHOD RESULTS:

1.JACOBI:

2.GAUSS SEIDEL:

3.SUCCESSIVE OVER RELAXATION

2.EXPLICIT FORM:

Explicit form of jacobi method for transient hetat conduction equation is ,

`T_P^(n+1)=T_P^n**(1-2K_1-2K_2)+(K_1(T_L^n+T_R^n)+K_2(T_B^n+T_T^n))`

Explicit form of Gauss seidel method for transient hetat conduction equation is ,

`T_P^(n+1)=T_P^n**(1-2K_1-2K_2)+(K_1(T_L^(n+1)+T_R^n)+K_2(T_B^(n+1)+T_T^n))`

Explicit form of Successive over relaxation method for transient hetat conduction equation is ,

`T_P^(n+1)=T_P^n**(1-lambda)+lambda**[T_P^n**(1-2k_1-2k_2)+(k1**(T_L^(n+1)+T_R^n)+k2(T_T^n+T_B^(n+1)))]`

EXPLICIT METHOD MAIN PROGRAM:

%explicit method
clear all
close all
clc

% Initial Conditions
L = 1; % Unit square length domain
nx = ny = 10;
% Grid spacing

x = linspace(0, L, nx);
dx = L/(nx-1);
y = linspace(0, L, ny);
dy = L/(ny-1);

error = 9e9;
tol = 1e-5;
t = 5;
dt = 1e-2;
alpha = 0.3;

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

% Boundary Conditions
T = ones(nx, ny);
T(1,:)=900;%B.s
T(end,:)=600;%T.S
T(:,1)=400;%L.S
T(:,end)=800;%R.S

% Edge Refining
T(1,1) = (400+900)/2;
T(1,ny) = (900+800)/2;
T(nx,1) = (400+600)/2;
T(nx,ny) = (600+800)/2;

%Explicit Approach
for iter_method = [1 2 3]
    if iter_method == 1
        T_jacobi = explicit_jacobi_us( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 );
    else if iter_method == 2
            T_gs = explicit_gauss_seidel_us( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 );
        else if iter_method == 3
               T_sor = explicit_sor_us( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 );
            end
        end
    end
end

FUNCTION PROGRAMS TO SOLVE THE EXPLICT FOR TRANSIENT STATE EQUATION:

1.JACOBI METHOD:

function T = explicit_jacobi_us( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 )
% Explicit Jacobi method for Transient state

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
jac_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (T_old(i-1, j) + T_old(i+1, j));
                term2 = (T_old(i, j-1) + T_old(i, j+1));
                T(i,j) = (To(i,j)*(1 - 2*k1 - 2*k2) + (k1*term1) + (k2*term2));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        jac_iter = jac_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(1)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X axis');
ylabel('Y axis');
colorbar
title_text = sprintf('Jacobi Method Iterations = %d @explicit unsteady ', jac_iter);
title(title_text)
end

2.GAUSS SEIDEL METHOD:

function T = explicit_gauss_seidel_us( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 )
% Explicit Gauss Seidel method for Transient State

meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
gs_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = ((To(i,j)*(1 - 2*k1 - 2*k2)) + (k1*term1) + (k2*term2));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        gs_iter = gs_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(2)
[a, b] = contourf(x, y, T);
colormap(jet)
clabel(a, b);
xlabel('X axis');
ylabel('Y axis');
colorbar
title_text = sprintf('Gauss seidel Method Iterations = %d @explicit unsteady ', gs_iter);
title(title_text)
end

3.SUCCESSIVE OVER RELAXATION METHOD :

function T = explicit_sor_us( nx, ny, x, y, dx, dy, T, t, dt, tol, k1, k2 )
% Explicit Successive Over-Relaxation Method for Transient State

lambda = 1.3; % Over Relaxation factor
meshgrid(x,y);

T_old = T;
To = T;
nt = t/dt;
sor_iter = 1;

tic
for m = 1:nt
    error = 9e9;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = T(i-1, j) + T_old(i+1, j);
                term2 = T(i, j-1) + T_old(i, j+1);
                T(i,j) = (1-lambda)*T_old(i,j) + (lambda*((To(i,j)*(1 - 2*k1 - 2*k2)) + (k1*term1) + (k2*term2)));
            end
        end
        error = max(max(abs(T_old - T)));
        T_old = T;
        sor_iter = sor_iter + 1;
    end
    To = T;
end
iter_time = toc;

figure(3)
[a, b] = contourf(x, y, T);
colormap(jet)
colorbar
clabel(a, b);
xlabel('X axis');
ylabel('Y axis');
title_text=sprintf('Successive Over-Relaxation Method Iterations = %d  @explicit unsteady', sor_iter);
title(title_text)
end

EXPLICIT METHOD RESULTS:

1.JACOBI:

2.GAUSS SEIDEL :

3.SUCCESIVE OVER RELAXATION METHOD:

CONCLUSION:

From the results of all methods and comparing thier number of iterations , we can conclude that.

  • jaccobi > Gauss-seidel > SOR according to the no. of iterations computed for all methods and schemes.
  • while comparing steady state and unsteady state , it is proven that the steady state is the most successful one. because it takes less number of iteration and less time to compute the equation
  • In all the iteration method SOR method is the effective method because, it takes less number of iterations for converging solution in both the schemes ( explicit & implicit ) compared to all other iteration methods
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.