To solve 2D heat conduction equation by using the various point iterative techniques in MATLAB

AIM:

To solve the 2D heat conduction equation by using the following point iterative techniques in MATLAB

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Explicit Scheme vs Implicit scheme

Numerical solution schemes are often referred to as being explicit or implicit. When a direct computation of the dependent variables can be made in terms of known quantities, the computation is said to be explicit. When the dependent variables are defined by coupled sets of equations, and either a matrix or iterative technique is needed to obtain the solution, the numerical method is said to be implicit.

Explicit Scheme: Is one in which the differential equation is discretized in such a way that there is only one unknown (at new time level n+1) on the left hand side (LHS) of the difference equation and it is computed in terms of all other terms on the RHS which are known (at previous time level n+1). Let us write an Explicit Scheme for eqn

Tt = Txx ............................................................................…(1)

Discretizing time derivative using forward difference and the second order derivative using central difference we get:

(Tjn+1 - Tjn )/ ∆t = (Tj+1n – 2 Tjn + Tj-1n) / ∆x2

This can be simplified as:

Tjn+1 = (1-2r) Tjn + r(Tj-1n + Tj+1n) ………………………………………..(2)

where r = ∆t/ ∆x2

This is an explicit scheme called FTCS (Forward differencing in Time and Central differencing in Space at time level n) for solving a 1-D heat equation. Explicit as we see one unknown on LHS (Tjn+1) being calculated in terms on all the term on RHS which are known as they are at previous time level n.

Explicit Schemes have to satisfy stability conditions in order for solution to converge. For example above Explicit Scheme in eqn (2) will be stable if r ≤ ½.

Implicit Scheme: Is one in which the differential equation is discretized in such a way that there are multiple unknowns at n+1 time level on the LHS of the equation and the terms on RHS are known ones at n time level. Let us write Implicit Scheme for eqn (1):

(Tjn+1 - Tjn )/ ∆t = (Tj+1n+1 – 2 Tjn+1 + Tj-1n+1) / ∆x2

(discretizing 2nd order derivative at time level n+1 unlike Explicit Scheme where we discretized it at n level). We simplify it to:

(1+2r) Tjn+1 - rTj-1n+1 - rTj+1n+1 = Tjn …………………………………..(3)

Here we discretized time derivative term in backward differencing (Backward in Time) and discretized 2nd order term in Central differencing in Space both at time level n+1. Therefore this implicit scheme is called BTCS scheme.

Implicit Schemes result in multiple simultaneous algebraic equations (here we’ll have 3 equations corresponding to 3 unknowns at time n+1 level) that need to be solved as a system of equations of the form[K]{T}n+1 = {F}n . In Our simple case here we have [K]3x3 {T}3n+1 = {F}3n

Good thing about Implicit Schemes is that they don’t have to satisfy any stability conditions. They are unconditionally stable. But they are computationally expensive as you need more storage space for [K] matrix and need sophisticated algorithms to solve.

The 2D Heat Conduction Equation is given by the following equation

``(delT)/dt+alpha((del^2T)/(del^2x)+(del^2T)/(del^2y)+(del^2T)/(del^2z))=0`

Assumptions made to solve the problem are as below-

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

4.Tolerance=1e-4 i.e error criteria considered

For Steady State Analysis we have `(delT)/dt=0``

So our Equation can be rewritten as-

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

Applying Second order CDS scheme on the above equation we get-

`alpha((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`

as `(Deltax^2)=(Deltay^2)` hence final equation becomes

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

Solving the above equation by iterative methods

1.In Jacobian method we have

`T_(i,j)=((T_(i-1,j)+T_(i+1,j)+T_(i,j-1)+T_(i,j+1))_(old))/4` , i.e we use previous or known values to compute the dependent variable.

2.In Gauss Seidal Method we have

`T_(i,j)=((T_(i-1,j)+(T_(i+1,j))_(old)+T_(i,j-1)+(T_(i,j+1))_(old))/4` , i.e we use previous or known values to compute the n+1 terms as well as use the current timestep values to compute n-1 terms simulateously

3.In Successive Over Relaxation (SOR) method

Over Relaxation factor f is given by

f=2/(1+sin*pie*dx)

hence `T_(i,j)=f*((T_(i-1,j)+(T_(i+1,j))_(old)+T_(i,j-1)+(T_(i,j+1))_(old))/4+(1-f)*(T_(old))``

MATLAB CODE

%Steady State Conduction By Different Iterative Methods
clear all
close all
clc
nx=100;%No of grid points
ny=nx;
u=298*ones(nx,ny);%Initializing Temperature Matrix
%Temperature at Boundary
u(:,ny)=800;
u(1,:)=600;
u(:,1)=400;
u(nx,:)=900;
uold=u;
%Grid Preparation
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=x;
iter=input(\'Enter input 1 -For JACOBI, 2 -For GAUSS-SEIDAL and 3-For SOR GAUSS-SEIDAL\\n\')
tol=1e-4;  %Limiter for iterations
err=1;     %Error for starting iterations
[xx,yy]=meshgrid(x,y); %2D space for plotting
tic;%Timer Starts
%JACOBI METHOD
if(iter==1)
    jacobi_iter=1;
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=0.25*(uold(i,j-1)+uold(i,j+1)+uold(i-1,j)+uold(i+1,j));
              end
        end
        err=max(abs(uold(:)-u(:)));
        uold=u;
       jacobi_iter=jacobi_iter+1;
    end
end

%GAUSS SIEDEL METHOD
if(iter==2)
    gs_iter=1;
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=0.25*(u(i,j-1)+uold(i,j+1)+u(i-1,j)+uold(i+1,j));
              end
        end
        err=max(abs(uold(:)-u(:)));
        uold=u;
       gs_iter=gs_iter+1;
    end
end


    %SUCCESSIVE OVER RELAXATION METHOD
    f=2/(1+sin(pi*dx));
if(iter==3)
    sor_iter=1;
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=(1-f)*uold(i,j)+(f*0.25*(u(i,j-1)+uold(i,j+1)+u(i-1,j)+uold(i+1,j)));
              end
        end
        err=max(abs(uold(:)-u(:)));
        uold=u;
       sor_iter=sor_iter+1;
    end
end
time=toc;%Timer Stops
figure(1)
[C,h]=contourf(xx,yy,u);
colorbar;
colormap(jet);
clabel(C,h);
if(iter==1)
    title(sprintf(\'Temperature profile for steady state heat conduction using Gauss Jacobi method \\n No. of Iterations:%d \\n Computation Time:%0.3f s\',jacobi_iter,time));
end
if(iter==2)
    title(sprintf(\'Temperature profile for steady state heat conduction using Gauss Seidal Method \\n No. of Iterations:%d \\n Computation Time:%0.3f s\',gs_iter,time));
end
if(iter==3)
    title(sprintf(\'Temperature profile for steady state heat conduction using SOR method \\n No. of Iterations:%d \\n Computation Time:%0.3f s\',sor_iter,time)); 
end
xlabel(\'xx\');
ylabel(\'yy\');

Thermal Contour Plots For Different Iterative methods

OBSERVATIONS:

The Jacobian method takes total 15951 iterations and computational time of 3.206 Secs to converge.

The Gauss Seidal method takes total 8677 iterations and computational time of 1.785 Secs to converge.

In SOR Gauss Seidal Method by providing a suitable relaxation factor takes total 289 iterations and computaional time of .072 s which is least of the other two methods at an error allowance of 1e-4 which is reasonably accurate as compared to other two methods.

The reason for the Jacobian method to take higher number of iterations is that it takes older value as the solution proceeds and hence it takes higher computational time to converge to the result.

While the other two method such as Gauss Seidal and SOR Gauss Seidal takes the updated value as the solution proceeds and hence take even lesser computational time which gives an upper hand to this type of iterative solver to get numerical solution.

For Transient State Analysis we have

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

Explicit Scheme:-

Applying Forward Differencing Scheme on the time derivative and Second order CDS scheme on the spatial derivative in the above equation we get-

`(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))=0`

For our convenience lets say

`T_(i,j)^(n+1)=T_(p)^(n+1),T_(i,j)^(n)=T_(p)^(n),T_(i+1,j)^(n)=T_(R)^(n),T_(i-1,j)^(n)=T_(L)^(n),T_(i,j+1)^(n)=T_(T)^(n) and T_(i,j-1)^(n)=T_(B)^(n)`

Now our final equation is-

`T_(p)^(n+1)=T_(p)^(n)(1-2k_1-2k_2)+k_1T_(L)^(n)+k_1T_(R)^(n)+k_2T_(T)^(n) +k_2T_(B)^(n)`

where

`k_1=alpha*(Deltat)/(Deltax^2) and k_2=alpha*(Deltat)/(Deltay^2)`

Implicit Scheme:-

Applying Backward Differencing Scheme on the time derivative and Second order CDS scheme on the spatial derivative in the above equation we get-

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

For our convenience lets say

`T_(i,j)^(n+1)=T_(p)^(n+1),T_(i,j)^(n)=T_(p)^(n),T_(i+1,j)^(n+1)=T_(R)^(n+1),T_(i-1,j)^(n+1)=T_(L)^(n+1),T_(i,j+1)^(n+1)=T_(T)^(n+1) and T_(i,j-1)^(n+1)=T_(B)^(n+1)`

Now our final equation is-

`T_(p)^(n+1)=(T_(p)^(n)+k_1T_(L)^(n+1)+k_1T_(R)^(n+1)+k_2T_(T)^(n+1) +k_2T_(B)^(n+1))/(1+2k_1+2k_2)`

where

`k_1=alpha*(Deltat)/(Deltax^2) and k_2=alpha*(Deltat)/(Deltay^2)`

Explicit Method Code

%Explicit Method for Transient State Heat Conduction
clear all
close all
clc

nx=100;%No of grid points
ny=nx;%Considering square profile
nt=600;%No of time steps
c=1.6563e-4;%Silver
x=linspace(0,1,nx);
y=x;
dx=1/nx;
dy=1/ny;
dt=0.1;
k=(c*dt)/(dx^2);
u=298*ones(nx,ny);%Initializing Temperature at all nodes 
%Temperature at Boundary
u(:,ny)=800;
u(1,:)=600;
u(:,1)=400;
u(nx,:)=900;
uold=u;
[xx,yy]=meshgrid(x,y);%2D Space for plotting
tic;
tol=1e-4;%limiter for iterations
err=1;%error allowance for starting iterations
tic;%Timer Starts
iter=input(\'Enter input 1 -For JACOBI, 2 -For GAUSS-SEIDAL and 3-For SOR GAUSS-SEIDAL\\n\')

%JACOBI METHOD
if(iter==1)
    jacobi_iter=1;
    for p=1:nt
        err=1;
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=k*(uold(i-1,j)+uold(i,j-1)+uold(i+1,j)+uold(i,j+1))+(1-4*k)*uold(i,j);
              end
        
 end
        err=max(abs(uold(:)-u(:)));
        uold=u;
       jacobi_iter=jacobi_iter+1;
    end
   
   end
end

%GAUSS SIEDEL METHOD
if(iter==2)
    gs_iter=1;
    for p=1:nt
        err=1;
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                   u(i,j)=k*(u(i-1,j)+u(i,j-1)+uold(i+1,j)+uold(i,j+1))+(1-4*k)*uold(i,j);
              end
        end
        err=max(abs(uold(:)-u(:)));
        uold=u;
       gs_iter=gs_iter+1;
    end
    
end
end

    %SUCCESSIVE OVER RELAXATION METHOD
    f=2/(1+sin(pi*dx));
if(iter==3)
    sor_iter=1;
     for p=1:nt
        err=1;
    
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=(1-f)*uold(i,j)+(f*(k*(u(i-1,j)+u(i,j-1)+uold(i+1,j)+uold(i,j+1))+(1-4*k)*uold(i,j)));
              end
        end
          err=max(abs(uold(:)-u(:)));
        uold=u;
       sor_iter=sor_iter+1;
    
    end
    
     end
end
time=toc;%Timer Stops
figure(1)
[C,h]=contourf(xx,yy,u);%Plotting of Transient State
colorbar;
colormap(jet);
clabel(C,h);
if(iter==1)
    title(sprintf(\'Temperature profile for Transient state heat conduction using Gauss Jacobi method \\n No. of Iterations:%d \\n Computation Time:%0.3f s\',jacobi_iter,time));
end
if(iter==2)
    title(sprintf(\'Temperature profile for  Transient state heat conduction using Gauss Seidal Method \\n No. of Iterations:%d \\n Computation Time:%0.3f s\',gs_iter,time));
end
if(iter==3)
    title(sprintf(\'Temperature profile for Transient state heat conduction using SOR method \\n No. of Iterations:%d \\n Computation Time:%0.3f s\',sor_iter,time)); 
end
xlabel(\'xx\');
ylabel(\'yy\');

Implicit Method Code

%Implicit Method For Transient State Heat Conduction

clear all
close all
clc

nx=100;%No of grid points in X direction
ny=nx; %No of grid points in Y direction
nt=600;%No of time steps
c=1.6563e-4;%Silver
x=linspace(0,1,nx);
y=x;%Considering a square grid
dx=1/nx;
dy=1/ny;
dt=0.1;
k1=(c*dt)/(dx^2);
k2=(c*dt)/(dy^2);
K=((1+(2*k1)+(2*k2))^-1);

u=298*ones(nx,ny);%Initializing Temperature at all nodes
%Temperature at Boundary
u(:,ny)=800;
u(1,:)=600;
u(:,1)=400;
u(nx,:)=900;
uold=u;
uinitial=u;
[xx,yy]=meshgrid(x,y);
tol=1e-4;%limiter for iterations
 err=1;%error allowance
tic;
iter=input(\'Enter input 1 for JACOBI 2 for GAUSS-SEIDAL and 3 for SOR GAUSS-SEIDAL\\n\')

%JACOBI METHOD
if(iter==1)
    jacobi_iter=1;
    for p=1:nt
        err=1
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=(uold(i,j)+k1*(uinitial(i-1,j)+uinitial(i+1,j))+k2*(uinitial(i,j-1)+uinitial(i,j+1)))*K;
              end
        
 end
        err=max(abs(uinitial(:)-u(:)));
        uinitial=u;
       jacobi_iter=jacobi_iter+1;
    end
    uold=u
   end
end

%GAUSS SIEDEL METHOD
if(iter==2)
    gs_iter=1;
    for p=1:nt
        err=1
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=(uold(i,j)+k1*(u(i-1,j)+uinitial(i+1,j))+k2*(u(i,j-1)+uinitial(i,j+1)))*K;
              end
        end
        err=max(abs(uinitial(:)-u(:)));
        uinitial=u;
       gs_iter=gs_iter+1;
    end
     uold=u;
end
end

    %SUCCESSIVE OVER RELAXATION METHOD
    f=2/(1+sin(pi*dx))
if(iter==3)
    sor_iter=1;
     for p=1:nt
        err=1;
    
    while(err>tol)
        for i=2:nx-1
              for j=2:ny-1
                  u(i,j)=uold(i,j)*K+((k1*K)*(u(i-1,j)*f+uinitial(i-1,j)*(1-f)+uinitial(i+1,j)))+((k2*K)*(u(i,j-1)*f+uinitial(i,j-1)*(1-f)+uinitial(i,j+1)));
              end