Steady and Transient analysis of 2D Heat Equation using Jacobi Gauss Siedel and SOR Iterative Solver

Objective:-

To perform steady and Transient state analysis of a 2D heat conduction equation using iterative solver like Jacobi, Gauss Siedel and SOR.

Input Parameter

Length=1 m

`alpha`=1

Boundary condition

Stability Criteria

Time step can be calculated using CFL condition

`alpha{Deltat}/{Deltax^2}<=1/4`

spatial step dx=0.1111, dt=0.003085

So, we have taken dt=0.001 during calculation

Steady State Analysis

2D steady state heat equation is given by

`{del^2T}/{delx^2}+{del^2T}/{dely^2}=0`

Since the above equation is elliptic in nature, it requires data from all neighbouring points. So, we should use central difference scheme.

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

solving, we get

`T_(i,j)={Deltax^2Deltay^2}/{2(Deltax^2+Deltay^2)}({T_(i+1,j)+T_(i-1,j)}/{Deltax^2}+{T_(i,j+1)+T_(i,j-1)}/{Deltay^2})`

Jacobi Iteration Formula:-

`T_(i,j)^(n+1)={dy^2}/{2(dx^2+dy^2)}(T_(i+1,j)^(n)+T_(i-1,j)^(n))+{dx^2}/{2(dx^2+dy^2)}(T_(i,j+1)^(n)+T_(i,j-1)^(n))`

In this method, T at n+1 iteration are calculated with the values of T at n interation.

Gauss Siedel Iteration Formula:-

`T_(i,j)^(n+1)={dy^2}/{2(dx^2+dy^2)}(T_(i+1,j)^(n)+T_(i-1,j)^(n+1))+{dx^2}/{2(dx^2+dy^2)}(T_(i,j+1)^(n)+T_(i,j-1)^(n+1))`

In this method, T at n+1 iteration are calculated with the already calculated values of T at n+1 interation. The values of T that have not been calculated at n+1 are taken from nth iteration.

Successive Over Relaxation:-

`T_(i,j)^(K+1)=(1-w)*T_(i,j)^K+w*T_(gauss)(i,j)`

Where T_gasuss is the result of Gauss Siedel method. Omega values varies from 0 to 1.

MATLAB main code

%Steady State Analysis
clear all;
close all;
clc;

%no of nodes
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(1,0,ny);

%grid length
dx=x(2)-x(1);
dy=y(1)-y(2);
T=100*ones(10,10);

%boundary condition
T(1,:)=600;
T(:,1)=400;
T(:,10)=800;
T(10,:)=900;

tol=1e-4;
T_old=T;
prompt=[\'press\\n\'...
    \'1 for jacobi iterative\\n\'...
    \'2 for gauss siedel iterative\\n\' ...
    \'3 for SOR iterative\\n\'];
iter=input(prompt);
error=9e9;
count=0;
w=1.5;

%Steady State Analysis
%Jacobi method
if iter==1
    while error>=tol
        T=jacobi_iterative(T,dx,dy);
        error=max(max(abs(T_old-T)));
        T_old=T;
        count=count+1
    end
    [c,h]=contourf(x,y,T);
    clabel(c,h);
    colorbar
    text=sprintf([\'Solving 2D steady heat equation using jacobi\'... 
        \' iterative solver\\n\'...
        \'No of grid points=%d;\'... 
        \' Total no of iteration=%d\'],nx,count);
    title(text);
end

%Gauss Siedel method
if iter==2
    while error>tol
        T=gauss_iterative(dx,dy,T);
        error=max(max(abs(T_old-T)));
        T_old=T;
        count=count+1
    end
[c,h]=contourf(x,y,T);
clabel(c,h);
colorbar;
text=sprintf([\'Solving 2D steady heat equation using Gauss Siedel\'...
    \' iterative solver\\n\'...
    \'No of grid points=%d;\'... 
    \' Total no of iteration=%d\'],nx,count);
title(text);
end

%Successive over relaxation
if iter==3
    K1=dy^2/(2*(dx^2+dy^2));
    K2=dx^2/(2*(dx^2+dy^2));
    while error>tol
        T=SOR_iterative(T,K1,K2);
        error=max(max(abs(T_old-T)));
        T_old=T;
        count=count+1
    end
[c,h]=contourf(x,y,T);
clabel(c,h);
colorbar;
text=sprintf([\'Solving 2D steady heat equation using SOR\'...
    \' iterative solver\\n\'...
    \'No of grid points=%d;\'... 
    \' Total no of iteration=%d\'],nx,count);
title(text);
end

%Jacobi Iterative Solver
function[P]=jacobi_iterative(t,x,y)
dx=x;
dy=y;
%K1=dy^2/2(dx^2+dy^2)
%K2=dx^2/2(dx^2+dy^2)
K1=0.025;
k2=0.025;
T_old=t;
for j=2:9
    for i=2:9
        t(i,j)=(T_old(i+1,j)+T_old(i-1,j)+T_old(i,j+1)+T_old(i,j-1))/4;
    end
end
P=t;
end

%Gauss Siedel Iterative solver
function[P] =gauss_iterative(x,y,t)
dx=x;
dy=y;
%K1=dy^2/2(dx^2+dy^2)
%K2=dx^2/2(dx^2+dy^2)
K1=dy^2/(2*(dx^2+dy^2));
k2=dx^2/(2*(dx^2+dy^2));
T_old=t;
for j=2:9
    for i=2:9
        t(i,j)=(T_old(i+1,j)+t(i-1,j)+T_old(i,j+1)+t(i,j-1))*0.25;
    end
end
P=t;
end

%SOR Iterative solver
function[P]=SOR_iterative(T,K1,K2)
w=1.5;
T_old=T;
for i=2:9
    for j=2:9
        t(i,j)=K1*(T_old(i+1,j)+T(i-1,j))+...
            K2*(T_old(i,j+1)+T(i,j-1));
        T(i,j)=(1-w)*T_old(i,j)+w*t(i,j);
    end
end
P=T;
end

Output

Conclusion

From above results, we can observe that Jacobi iterative solver takes more number of iteration. Jacobi iterative solver takes 214 number of iteration to reach convergence while Gauss siedel and successive over relaxation iterative solver takes 114 and 29 number of iteration respectively. The number of iteration in SOR depends on the selection of omega. In the above graph, convergence is achieved in 29 iteration when omega is 1.5. This value of omega was found to be optimum. The value greater or smaller than 1.5 slows down the rate of convergence. When omega is 2 or more than 2, solution will not converge at all.

Transient State Ananlysis

The 2D heat equation for transient state is given by

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

Since the above equation hyperbolic part and elliptic part, we should use forward difference scheme for time derivative and central difference scheme for spatial derivative.

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

On solving, we get

`T_(i,j)^(n+1)=(1-2k_1-2k_2)T_(i,j)^n+k_1(T_(i-1,j)^n+T_(i+1,j)^n)+k_2(T_(i,j-1)^n+T_(i,j+1)^n)`

The above equation is formula for explicit scheme because we don\'t need system of linear equation to solve.

where

`k_1=alpha{Deltat}/{Deltax^2}`

`k_2=alpha{Deltat}/{Deltay^2}`

Implicit Scheme for 2d heat equation for transient state is given by

`T_(i,j)^(n+1)={T_(i,j)^n+k_1(T_(i-1,j)^(n+1)+T_(i+1,j)^(n+1))+k_2(T_(i,j-1)^(n+1)+T_(i,j+1)^(n+1))}/{1+2k_1+2k_2}`

Implicit scheme is solved by solving system of coupled linear equation but iterative solver helps to solve these equation by direct method. Jacobi, Gauss Siedel and SOR iterative solver are used to solve implicit scheme.

MATLAB main code

%Transient State Analysis
clear all;
close all;
clc;
%no of grid
nx=10;
ny=10;
T=ones(10,10);
%boundary condition
T(1,:)=600;
T(:,1)=400;
T(:,10)=800;
T(10,:)=900;
error=1;
prompt=[\'press\\n\'...
    \'1 for explicit\\n\'...
    \'2 for jacobi implicit\\n\'...
    \'3 for gauss siedel implicit\\n\' ...
    \'4 for SOR implicit\\n\'];
iter=input(prompt);
%unsteady state 
%explicit method
if iter==1
    [T_final,error,c]=unsteady_explicit(nx,ny,T);
end
%implicit jacobi method
if iter==2
    [T_final,error,c]=implicit_jacobi(nx,ny,T);
end
%implicit gauss siedel method
if iter==3
    [T_final,error,c]=implicit_gauss(nx,ny,T);
end
%implicit SOR method
if iter==4
    [T_final,error,c]=implicit_sor(nx,ny,T);
end

%Explicit Scheme
function[T,error,count]=unsteady_explicit(nx,ny,t)
x=linspace(0,1,nx);
y=linspace(1,0,ny);
dx=x(2)-x(1);
dy=y(1)-y(2);
error=9e9;
count=0;
T_old=t;
dt=1/1000;
alpha=1;
time=0.4;
nt=time/dt;
K1=alpha*dt/(dx^2);
K2=alpha*dt/(dy^2);
T=t;
for k=1:nt
    for j=2:9
        for i=2:9
            term_1=T_old(i+1,j)+T_old(i-1,j);
            term_2=T_old(i,j+1)+T_old(i,j-1);
            T(i,j)=(1-2*K1-2*K2)*T_old(i,j)+K1*term_1+K2*term_2;
        end
    end
    error=max(max(abs(T_old-T)));
    T_old=T;
    count=count+1;
end
figure(1)
[c,h]=contourf(x,y,T);
clabel(c,h);
colorbar;
text=sprintf([\'Solving 2D unsteady heat equation using\\n\'... 
    \'explicit\'... 
    \' iterative solver\\n\'...
    \'No of grid points=%d;\'... 
    \' Total no of iteration=%d\'],nx,count);
title(text);
end

%Implicit Scheme using Jacobi
function[T,error,count]=implicit_jacobi(nx,ny,t)
x=linspace(0,1,nx);
y=linspace(1,0,ny);
dx=x(2)-x(1);
dy=y(1)-y(2);
tol=1e-4;
error=9e9;
count=0;
T_old=t;
dt=1/1000;
alpha=1;
time=0.4;
nt=time/dt;
K1=alpha*dt/(dx^2);
K2=alpha*dt/(dy^2);
T=t;
T_prev=t;
for k=1:nt
    error=1;
    while error>tol
        for j=2:9
            for i=2:9
                term_1=T_old(i+1,j)+T_old(i-1,j);
                term_2=T_old(i,j+1)+T_old(i,j-1);
                T(i,j)=(T_prev(i,j)+K1*term_1+K2*term_2)/(1+2*K1+2*K2);
            end
        end
    error=max(max(abs(T_old-T)));
    T_old=T;
    count=count+1;
    end
    T_prev=T;
end
figure(1)
[c,h]=contourf(x,y,T);
clabel(c,h);
colorbar;
text=sprintf([\'Solving 2D unsteady heat equation using\\n\'... 
    \'implicit jacobi\'... 
    \' iterative solver\\n\'...
    \'No of grid points=%d;\'... 
    \' Total no of iteration=%d\'],nx,count);
title(text);
end
            

%Implicit Scheme using Gauss Siedel
function[T,error,count]=implicit_gauss(nx,ny,t)
x=linspace(0,1,nx);
y=linspace(1,0,ny);
dx=x(2)-x(1);
dy=y(1)-y(2);
tol=1e-4;
error=9e9;
count=0;
T_old=t;
dt=1/1000;
alpha=1;
time=0.4;
nt=time/dt;
K1=alpha*dt/(dx^2);
K2=alpha*dt/(dy^2);
T=t;
T_prev=t;
for k=1:nt
    error=1;
    while error>tol
        for j=2:9
            for i=2:9
                term_1=T_old(i+1,j)+T(i-1,j);
                term_2=T_old(i,j+1)+T(i,j-1);
                T(i,j)=(T_prev(i,j)+K1*term_1+K2*term_2)/(1+2*K1+2*K2);
            end
        end
    error=max(max(abs(T_old-T)));
    T_old=T;
    count=count+1;
    end
    T_prev=T;
end
figure(1)
[c,h]=contourf(x,y,T);
clabel(c,h);
colorbar;
text=sprintf([\'Solving 2D unsteady heat equation using\\n\'... 
    \'implicit gauss siedel\'... 
    \' iterative solver\\n\'...
    \'No of grid points=%d;\'... 
    \' Total no of iteration=%d\'],nx,count);
title(text);
end

%Implicit Scheme using SOR
function[T,error,count]=implicit_sor(nx,ny,t)
x=linspace(0,1,nx);
y=linspace(1,0,ny);
dx=x(2)-x(1);
dy=y(1)-y(2);
tol=1e-4;
error=9e9;
count=1;
T_old=t;
dt=1/1000;
alpha=1;
time=0.4;
nt=time/dt;
K1=alpha*dt/(dx^2);
K2=alpha*dt/(dy^2);
T=t;
T_prev=t;
w=1.5;
T_gauss=ones(10,10);
for k=1:nt
    error=1;
    while error>tol
        for j=2:9
            for i=2:9
                term_1=T_old(i+1,j)+T(i-1,j);
                term_2=T_old(i,j+1)+T(i,j-1);
                T_gauss(i,j)=(T_prev(i,j)+K1*term_1+K2*term_2)/(1+2*K1+2*K2);
                T(i,j)=(1-w)*T_prev(i,j)+w*T_gauss(i,j);
            end
        end
    error=max(max(abs(T_old-T)));
    T_old=T;
    count=count+1;
    end
    T_prev=T;
end
figure(1)
[c,h]=contourf(x,y,T);
clabel(c,h);
colorbar;
text=sprintf([\'Solving 2D unsteady heat equation using\\n\'... 
    \'implicit SOR\'... 
    \' iterative solver\\n\'...
    \'No of grid points=%d;\'... 
    \' Total no of iteration=%d\'],nx,count);
title(text);
end

Output

Conclusion

Explicit scheme is the fastest to solve 2D transient heat equation because next iteration is based upon previous state which has been already calculated.

Implicit scheme requires guess values and iteration for convergence as it is based upon future values which are unknown at time of calculation. After each time step, it checks convergence.

Among all the implicit methods, SOR method has fastest convergence rate (in 2083 iterations) followed by Gauss Siedel method (in 2238 iterations) and then Jacobi method (in 2796 iterations).