2D Heat Conduction Equation

1.) STEADY STATE PROGRAM

For Solving 2D Heat conduction problem in steady state. First we will take a steady state heat conduction equation:

`frac{delta^2T}{deltax^2}+frac{delta^2T}{deltay^2}=0` .................(1)

Approximating the first two terms with the central difference scheme:

`frac{delta^2T}{deltax^2}=frac{Tunderset(i-1,j)( )-2Tunderset(ij)( )+Tunderset(i+1,j)()}{Deltax^2}` .................(2)

`frac{delta^2T}{deltay^2}=frac{Tunderset(i,j-1)( )-2Tunderset(ij)( )+Tunderset(i,j+1)()}{Deltay^2}` ...................(3)

As, given in the challenge that nx=ny and thus `Deltax=Deltay`

Substituting the equation number 2 and 3 in equation 1 and after rearranging the terms the final equation, thus the final becomes:

`Tunderset(ij)( )=frac{1}{4}(Tunderset(i-1,j)( )+Tunderset(i+1,j)( )+Tunderset(i,j-1)( )+Tunderset(i,j+1)( ))`

As this equation has no time derivative the time marching process will not be needed

we will be solving this equation by three iterative methods as

1.) Jacobi

2.) Guass-Seidel

3.) Successive Over Relaxation (SOR)

For solving the equation with these methods first we will make a matrix of T with the guess terms and putting the boundary condition and in time marching we will call the previous iteration value and this process will be repeated until the error value comes under a prescribed value.

In Guass-Seidel method of approximation while calculating t(i,j), the t(i-1,j) and t(i,j-1) will already be calculated thus we will replace the terms from previous time step iterations to the newer one thus keeping the number of iterations to lower value and the thus taking less time to solve the equation.

In SOR method a factor `alpha`is multiplied to either the equation of jacobi and gauss-seidel which helps us to speed us the method of iteration

For steady state equation the program can be seen below:

close all
clear all
clc
%defining constants
nx=20;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=dx;
tol=1e-4;
%defining boundary conditions
t=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;
%defining corner points by taking average of bounday conditions
t(1,1)=500;
t(1,ny)=700;
t(nx,1)=650;
t(nx,ny)=850;
t_initial=t;
%loop for different methods
for method=1:3
  error=2000;  
  if method==1
    %jacobi iteration
    told=t;
    jacobi_iter=0;
    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(method)
    contourf(x,y,t)
    clabel(contourf(x,y,t))
    xlabel(\'x\')
    ylabel(\'y\')
    colorbar
    title(sprintf(\'Jacobi No. of Iter %d\',jacobi_iter))
    pause(0.000001)
    end
  

   elseif method==2
      %gs iteration
      told=t_initial;
      gs_iter=0
      error=2000;
      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(method)
      contourf(x,y,t)
      clabel(contourf(x,y,t))
      xlabel(\'x\')
      ylabel(\'y\')
      colorbar
      title(sprintf(\'Gauss Seidel No. of Iter %d\',gs_iter))
      pause(0.0000001)
      end
      
    
    elseif method==3
    %sor iteration
      told=t_initial;
      sor_iter=0;
      error=2000;
      %sor factor calculation
      z=cos(pi/nx)+cos(pi/ny);
      omega=4/(2+sqrt(4-z^2));
      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(method)
    contourf(x,y,t)
    clabel(contourf(x,y,t))
    xlabel(\'X\')
    ylabel(\'Y\')
    colorbar
    title(sprintf(\'Successive Over Relaxation No. of Iter %d\',sor_iter))
    pause(0.0000001)
      end
    end
end

The results for the steady state program for the 2D heat conduction equation can be seen below:

1.) Steady State

Method used= Jacobi Approximation

Number of Iterations=868

2.) Steady State

Method Used = Gauss Seidel

No. of Iteratons = 462

3.) Steady State

Method Used = Successive Over Relaxation (SOR)

No. of Iteratons = 60

Here we can see that the all the four graphs have converged and the steady state achieved can be seen in the figure. Temperature at different boundary points can be seen written in the graph. In jacobi method it can be seen that the number of iteration taken during the course are 868 whereas in the gauss seidel program it can be seen that the number of iterations taken are 462. This can be explained due to choosing the updated value rather than the older values as explained above. For successive over relaxation the alpha is taken as 1.2 and this factor can be seen as it boosted our results and thus solved our equation in just 60 iterations.

2.) TRANSIENT EXPLICIT

In unsteady state the the heat convection equation becomes:

`frac{dT}{dt}=alpha(frac{delta^2T}{deltax^2}+frac{delta^2t}{deltay^2})` ...................1

Approximating the time differencial with backward differencing scheme:

`frac{dT}{dt}=frac{T^(n+1)underset(ij)()-T^(n)underset(ij)( )}{Deltat}` ......................2

Approximating the spacial terms by central differencing:

`frac{delta^2T}{deltax^2}=frac{Tunderset(i+1,j)()-2Tunderset(i,j)()+Tunderset(i-1,j)()}{Deltax^2}`..........................3

`frac{delta^2T}{deltay^2}=frac{Tunderset(i,j+1)()-2Tunderset(i,j)()+Tunderset(i,j-1)()}{Deltay^2}`..........................4

For the case of explicit condition the spacial terms are solved at the last time step solved which in this case will be n.

substituting the equation 2,3 and 4 in equation 1 and further rearranging:

`Tunderset(i,j)()^(n+1)=Tunderset(i,j)()^(n)+alphaDeltat(frac{Tunderset(i+1,j)()^(n)-2Tunderset(i,j)()^(n)+Tunderset(i-1,j)()^(n)}{Deltax^2}+frac{Tunderset(i,j+1)()^(n)-2Tunderset(i,j)()^(n)+Tunderset(i,j-1)()^(n)}{Deltay^2})`

As for this case the nx=ny thus dx=dy putting this into equation and further simplifying we get:

`Tunderset(i,j)()^(n+1)=Tunderset(i,j)()^(n)+frac{alphaDeltat}{Deltax^2}(Tunderset(i+1,j)()^(n)-2Tunderset(i,j)()^(n)+Tunderset(i-1,j)()^(n)+Tunderset(i,j+1)()^(n)-2Tunderset(i,j)()^(n)+Tunderset(i,j-1)()^(n))`

Replacing `frac{alphaDeltat}{Deltax^2}`with k and rearranging we get:

`Tunderset(i,j)()^(n+1)=(1-4k)Tunderset(i,j)()^(n)+k(Tunderset(i+1,j)()^(n)+Tunderset(i-1,j)()^(n)+Tunderset(i,j+1)()^(n)+Tunderset(i,j-1)()^(n))`

Now for this equation we will write a program to solve the 2D heat conduction equation explicitly. For the solution to be stable in CFD we should follow CFL criteria. For the case of diffusion the CFL Criteria is written below:

`frac{alphaDeltat}{MIN(Deltax^2,Deltay^2)}<=0.5`

As we are solving the current equation for the case where dx=dy. So, for thi special case, the criteria becomes:

`frac{alphaDeltat}{Deltax^2}<=0.25`

We can see the program attached below for explicit 2D heat Conduction equation:

clear all
close all
clc

nx=10;
ny=nx;
x=linspace(0,1,nx);
y=x;
nt=300;
dx=x(2)-x(1);
dy=dx;
dt=1e-3;
alpha=2;
%boundary conditions
T=ones(nx,ny);
T(1,2:ny-1)=600;
T(nx,2:ny-1)=900;
T(2:nx-1,1)=400;
T(2:nx-1,ny)=800;
%average for corner points
T(1,1)=500;
T(1,ny)=700;
T(nx,1)=650;
T(nx,ny)=850;

%CFL Number
m=alpha*dt/(dx^2);
k1=m;
Told=T;
iter=0;
%time loop
for k=1:nt
  %spacial loop
  for i=2:nx-1
    for j=2:ny-1
      term1=(1-4*k1)*Told(i,j);
      term2=k1*(Told(i-1,j)+Told(i+1,j));
      term3=k1*(Told(i,j-1)+Told(i,j+1));
      T(i,j)=term1+term2+term3;
    endfor
  endfor
   %error=max(max(abs(T-Told)));
   Told=T;
   iter=iter+1;
   figure(1)
   [M,N]=contourf(x,y,T);
   clabel(M,N)
   colormap(jet)
   colorbar
   t_text=sprintf(\'Number of Iterations = %d\',iter);
   title(t_text)
   pause(0.00001)
   xlabel(\'X\')
   ylabel(\'Y\')
   fprintf(\'CFL Number=%d\',m);
endfor

Here in this program the alpha is 2. the value of dt is 0.001 and that of dx or dy is 0.1111.

For this case we can solve the

`frac{2xx0.001}{0.11111^2}=0.16`

Here we can see that our CFL Number is less than that of desired criteria thus our solution should be stable.

The graph for the results can be seen below:

In this case we dont need to solve this equation with jacobi, gauss seidel or SOR methods because in explict method we solve the spacial frame only at previous time step, thus we don\'t need any other iterative schemes. Here in this graph for alpha =2, dt as 0.001for about 300 time steps. The number of nodes along x or y axis were taken as 10. We reached steady state at about 300 iterations.

3.) TRANSIENT IMPLICIT SCHEME

Similarly in the case of implicit scheme the equation becomes the same as that in case of the explicit but with a little change as shown below:

`Tunderset(i,j)()^(n+1)=(1-4k)Tunderset(i,j)()^(n)+k(Tunderset(i+1,j)()^(n+1)+Tunderset(i-1,j)()^(n+1)+Tunderset(i,j+1)()^(n+1)+Tunderset(i,j-1)()^(n+1))`

As this equation is solved at the same time step at which it is being calculated, so the result obained from this case is typically stable so we dont require to fulfill CFL criteria:

This case can be solved under three iterative schemes as explained above that are:

1. Jacobi
2. Gauss Seidel
3. Successive Over Relaxation

In Jacobi iteration the solving equation becomes :

`Tunderset(i,j)()^(n+1)=(1-4k)Tunderset(i,j)()^(n)+k(Tunderset(i+1,j)()^(n)+Tunderset(i-1,j)()^(n)+Tunderset(i,j+1)()^(n)+Tunderset(i,j-1)()^(n))`

as the terms taken in this iterations are less updated thus in gauss seidel method of iterations we take more updated terms as shown below:

`Tunderset(i,j)()^(n+1)=(1-4k)Tunderset(i,j)()^(n)+k(Tunderset(i+1,j)()^(n)+Tunderset(i-1,j)()^(n+1)+Tunderset(i,j+1)()^(n)+Tunderset(i,j-1)()^(n+1))`

In Case of SOR we take a multiplying factor omega and then the factor is multiplied by the equation thus our final equation becomes:

`Tunderset(i,j)()^(n+1)=(1-Omega)Tunderset(i,j)()^(n)+Omegaxx(1-4k)Tunderset(i,j)()^(n)+k(Tunderset(i+1,j)()^(n)+Tunderset(i-1,j)()^(n)+Tunderset(i,j+1)()^(n)+Tunderset(i,j-1)()^(n))`

in this case the omega is generally taken in between 1 to 1.5 and it is generally taken by hit and trail method. If the value is taken above a permissible limit the result can explode. This multiplication factor is generally taken to boost our results such that at less number of iterations we could converge the results.

The program for the implicit scheme can be seen below:

close all
clear all
clc
%defining constants
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=dx;
tol=1e-4;
t=2;
dt=1e-2;
nt=t/dt;
%defining boundary conditions
T=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;
%defining corner points by taking average of bounday conditions
T(1,1)=500;
T(1,ny)=700;
T(nx,1)=650;
T(nx,ny)=850;
alpha=10;
k1=alpha*dt/dx^2;
t_initial=T;
t_prev = T;
told = T;
method=1;
%loop for different methods
for method=1:3 
  if method==1
    %jacobi iteration
    told=t_initial;
    t_prev=t_initial;
    jacobi_iter=0;
     for k=1:nt  
     %disp(k);
     error=2000;  
      while (error>tol)
        for i=2:nx-1
          for j=2:ny-1
            T(i,j)=(t_prev(i,j)+k1*(told(i-1,j)+told(i+1,j)+told(i,j-1)+told(i,j+1)))/(1+4*k1);
          end
        end
        error = max(max(abs(told-T)));
        told=T;
        jacobi_iter=jacobi_iter+1;
      end
      t_prev=T;
     end
      figure(method)
      [C,H]=contourf(x,y,T);
        clabel(C,H)
      xlabel(\'x\')
      ylabel(\'y\')
      colorbar
      title(sprintf(\'Jacobi No. of Iter %d\',jacobi_iter))
      pause(0.001)
  elseif method==2
      %gs iteration
      told=t_initial;
      t_prev = t_initial;
      gs_iter=0;
       error=2000; 
       for k=1:nt 
        error=2000;
        while (error>tol)
          for i=2:nx-1
            for j=2:ny-1
              T(i,j)=(t_prev(i,j)+k1*(T(i-1,j)+told(i+1,j)+T(i,j-1)+told(i,j+1)))/(1+4*k1);
            end
          end
        error = max(max(abs(told-T)));
        told=T;
        gs_iter=gs_iter+1;
        figure(method)
        [C,H]=contourf(x,y,T);
        clabel(C,H)
        xlabel(\'x\')
        ylabel(\'y\')
        colorbar
        title(sprintf(\'Gauss Seidel No. of Iter %d\',gs_iter))
        pause(0.001)
        end
        t_prev=T;
      end
    
    elseif method==3
    %sor iteration
      told=t_initial;
      t_prev = t_initial;
      sor_iter=0;
       error=2000; 
      %sor factor calculation
      %z=cos(pi/nx)+cos(pi/ny);
      omega=1.3;
       for k=1:nt 
          error=2000; 
        while (error>tol)
          for i=2:nx-1
            for j=2:ny-1
              T(i,j)=(1-omega)*told(i,j) + omega*((t_prev(i,j)+k1*(T(i-1,j)+told(i+1,j)+T(i,j-1)+told(i,j+1)))/(1+4*k1));
            end
          end
          error = max(max(abs(told-T)));
          told=T;
          sor_iter=sor_iter+1;
          figure(method)
      [C,H]=contourf(x,y,T);
        clabel(C,H)
      xlabel(\'X\')
      ylabel(\'Y\')
      colorbar
      title(sprintf(\'Successive Over Relaxation No. of Iter %d\',sor_iter))
      pause(0.00001)
        end
        t_prev=T;
       end   
end
end

The result for the jacobi scheme for the plot can be seen below:

No. of Iterations: 1177

Length of Time Step:0.001

alpha=10

The result for the Gauss Seidel scheme for the plot can be seen below:

No. of Iterations: 747

Length of Time Step:0.001

alpha=10

The result for the Successive Over Relaxation scheme for the plot can be seen below

No. of Iterations: 497

Length of Time Step:0.001

alpha=10

omega = 1.3

Finally here it can be seen that the in the jacobi method the number of iteration taken to solve are much more larger than that of the gauss seidel and no. of iteration taken to solve the gauss seidel are much more larger than that of SOR. Thus Sor can be used for better covergence of results and a better solution thus helping to solve the equation quicker without using less computing resources.