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

Solving the steady and unsteady 2D heat conduction problem

EEpuri Yoga Narasimha18 Feb 2021

Objective : To solve steady and unsteady in Explicit and implicit form using iterative methods.

Iterative methods:

1. Jacobi method.

2. Gauss siedel method.

3. Point successive over-relaxation using any iterative methods.

4. Line successive over-relaxation using any iterative methods.

5. Alternating Direction Implicit(ADI) Method.

6. Conjugate gradient method.

7. Biconjugate method.

8. Multigrid method.

Requirement of Iterative solvers :

As in general implicit method can be solved by writing discretized equation at all discrete points in the domain , and forming a matrix 'AX = B' ,

It would be difficult to form equations for matrix as the discrete points in the physical domain can't be predictable.

So , we use Iterative solvers.

Nodes reference for Discretization :

Note:

At present time step , the blue part will be get executed before T(i,j) , orange grid points won't executed before T(i,j).

Description about Iterative Solvers :

1. Jacobi :

In this method we use , values of previous iteration only.

2. Gauss siedel :

In this method we use , values of previous iteration values of which is not executed at the present time step iteration ,

Present iteration values which is executed at the present time step.

Blue part grid points will be considered the present time step values.

Orange part grid points will be considered the previous time step values.

3. Successive over relaxation method : (`omega` > 1)

`T_(i,j)^(k+1) = (1 - omega)*T_(i,j)^k + omega*(T_(i,j)^ast)`

PDE considered 2D Heat conduction equation with No heat geneation:

Transient :

`frac{delT}{delt} = alpha*frac{del^2T}{delx^2} + alpha*frac{del^2T}{dely^2}`

Steady :

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

Discretized UnSteady Explicit form :

`T_(i,j)^(k+1) = T_(i,j) ^k + alpha*dt*(frac{T_(i+1,j)^k - 2*T_(i,j)^k + T_(i-1,j)^k}{Deltax^2} + frac{T_(i,j+1)^k - 2*T_(i,j)^k + T_(i,j-1)^k}{Deltay^2})`

Discretized Unsteady Implicit form:

`T_(i,j)^(k+1) = frac{1}{1+2*k_1 + 2*k_2}*T_(i,j)^k + frac{alpha*Deltat(T_(i+1,j)^(k+1) + T_(i-1,j)^(k+1))}{(1+2*k_1+2*k_2)*Deltax^2} + frac{alpha*Deltat(T_(i,j+1)^(k+1) + T_(i,j-1)^(k+1))}{(1+2*k_1+2*k_2)*Deltay^2}`

Discretized Steady form :

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

For square domain and dirichilet condition :

Relaxation factor (`omega`):

b = `(frac{Deltax}{Deltay})^2`

a = `frac{cos(frac{pi}{n-1}) + b*cos(frac{pi}{n-1}) }{(1 + b^2)^2}`

`omega = frac{2 - sqrt(1 - a)}{a}`

Matlab code for Steady 2D Heat conduction Equation :

%% Clearing commands
close all
clear all
clc

% Initial data
% Thermal diffusivity
alpha = 1.6;    

% Considered 1x1 square grid
% Length of domain
l = 1;

% Number of girds along x and y swwep 100
n = 20;

% grid size along x 
d_x = l/(n-1);

% grid size along y
d_y = l/(n-1);

% Tolerance 
Tol = 1e-4;

% error initial
error = 9e9;

% Creating array for iuput data
T = 300*ones(n,n);

% Boundary conditions

% Bottom
T(1,:) = 900;
 
% Top
T(n,:) = 600;

% Left 
T(:,1) = 400;

% Right 
T(:,n) = 800;

% Corner points
T(1,1) = (T(1,2) + T(2,1))/2;
T(n,1) = (T(n-1,1) + T(n,2))/2;
T(1,n) = (T(1,n-1) + T(2,n))/2;
T(n,n) = (T(n-1,n) + T(n,n-1))/2;

k = 2*((d_x^2 + d_y^2)/(d_x^2*d_y^2));

y1 = 1;

x = linspace(0,l,n);

y = x;

[X,Y] = meshgrid(x,y);

% directing user to choose the method
fprintf('solver 1 : steady n solver 2 : Transient n :')
solver = input('Enter the number to select the type of solver :');
fprintf(' steady methods n method 1 : jacobi n method 2 : point gauss sidel n method 3 : SOR.P n ')
fprintf('method 4 : Transient Explicit  method )
method = input('Enter the number to execute the method :');

if solver == 1
 if method == 1
   disp(" hey jacobi here thanks for choosing me, haha , let's get numbers...Hurrayyyy!!!")
   T_new = T;
   tic
    while error > 1e-4
    %for i=1:2
      % y sweep loop
      for i=2:n-1
          % x sweep loop
          for j = 2:n-1
             T_new(i,j) = (1/k)*(((T(i-1,j) + T(i+1,j))/d_x^2)  + ((T(i,j-1) + T(i,j+1))/d_x^2));
          end         
      end %Iteration loop  completed
      % Checking bound of error
      error = max(max(abs(T_new - T)));
      T = T_new;
      contourf(X,Y,T,'ShowText','on')
      colorbar
      colormap(jet)
      title('Steady Jacobi time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01) 
    end %while loop end   
      ax = gca;
      jacobi = toc;
     video_file = VideoWriter('method_1');
     open(video_file)
     writeVideo(video_file,ac)
     close(video_file)
elseif method == 2
   disp(" hey I am faster than Jacobi, jai gauss-siedel,let's get numbers...Hurrayyyy!!!")
   T_new = T;
   tic
    while error > 1e-4
    %for i=1:2
      % y sweep loop
      for i=2:n-1
          % x sweep loop
          for j = 2:n-1
             T_new(i,j) = (1/k)*((T_new(i-1,j)+T(i+1,j))/d_x^2) + (1/k)*((T_new(i,j-1) + T(i,j+1))/d_y^2); 
          end          
      end
      error = max(max(abs(T_new - T)));
      T = T_new;   
      contourf(X,Y,T,'ShowText','on')
      colorbar
      colormap(jet)
      title('steady Gauss siedel time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01) 
      %{
      figure(2)
      c1 = x;
      c2 = c1;
      imagesc(c1,c2,T_new)
      colorbar
      colormap(jet)
      %}
    end %end while
    g_s = toc;
     video_file = VideoWriter('method_8');
     open(video_file)
     writeVideo(video_file,ac)
     close(video_file)
     ax = gca;
     set(ax,'YDir','reverse')
elseif method == 3
   disp(" If fastest method is your dad , then i am your dad..@KGF-chapter-SOR")
   T_new = T;
   b = (d_x/d_y)^2;
   a = ((cos(pi/(n-1)) + b*cos(pi/(n-1)) )/( 1 + b^2 ))^2;
   w_opt = (2-sqrt(1-a))/a;
   tic
    while error > 1e-4
    %for i=1:2
      % y sweep loop
      for i=2:n-1
          % x sweep loop
          for j = 2:n-1
              % For point SOR-g.s
              % T_new = (1-w)*T + w*(T_g-s)
             T_new(i,j) =  (1-w_opt)*T(i,j) + w_opt*((1/k)*((T_new(i-1,j)+T(i+1,j))/d_x^2) + (1/k)*((T_new(i,j-1) + T(i,j+1))/d_y^2));
          end   
      end
      error = max(max(abs(T_new - T)));
      T = T_new;
      contourf(X,Y,T,'ShowText','on')
      colorbar
      colormap(jet)
      title('Steady SOR_g-s time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01)
      ax = gca;
      set(ax,'YDir','reverse')
      SOR_g_s = toc;
    end % While end loop
     video_file = VideoWriter('method_3');
     open(video_file)
     writeVideo(video_file,ac)
     close(video_file)
     ax = gca;
     set(ax,'YDir','reverse')
 end
end

Results :

Steady Heat conduction equation using 3 iterative solvers :

Jacobi method:

Gauss - siedel method :

Successive-over-relaxation method:

Description About Methods :

As can observe from the above results for to approach steady state , Time taken : Jacobi > Gauss - Siedel > Successive over relaxation using gauss siedel.

So , Succesive over relaxation method is faster than both methods.

So , SOR-g-s is more efficient.

Explicit transient equation :

               *******Explicit transient equation *********

% condition for convergence
% CFL1 + CFL2 < 0.5;
% CFL1 = alpha*dt/d_x^2 ;
% CFL2 = alphs*dt/d_y^2;
dt1 = (1/2)*((d_x^2*d_y^2)/((d_x^2+d_y^2)*alpha));
dt = dt1 - 0.01*dt1;
k1 = (alpha*dt)/(d_x^2);
k2 = (alpha*dt)/(d_y^2);
 if method == 4
        disp(" hey Explicit jacobi here thanks for choosing me, haha , let's get numbers...Hurrayyyy!!!")
   T_new = T;
   tic
    while error > 1e-4
    %for i=1:2
      % y sweep loop
      for i=2:n-1
          % x sweep loop
          for j = 2:n-1
             T_new(i,j) = (1-2*k1-2*k2)*T(i,j) + k1*(T(i+1,j) + T(i-1,j)) + k2*(T(i,j+1) + T(i,j-1));
          end         
      end %Iteration loop  completed
      % Checking bound of error
      error = max(max(abs(T_new - T)));
      T = T_new;
      contourf(X,Y,T,'ShowText','on')
      colorbar
      colormap(jet)
      title('Transient Explicit Method time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01)
      ax = gca;
      jacobi = toc;
    end %while loop end 
    video_file = VideoWriter('method_4');
    open(video_file)
    writeVideo(video_file,ac)
    close(video_file)
 end

Explicit Transient Equation :

Unsteady using Jacobi Method :

                   ****** Unsteady using Jacobi method ******

close all
clear all
clc

nx = 20;
ny = 20;
nt = 1000;

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

d_x = x(2) - x(1);
d_y = y(2) - y(1);

y1 = 1;

error = 9e9;
tolerance = 1e-4;
dt = 1e-4;

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T = 300*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:ny-1) = T_B;

T(1,1) = (T_L + T_T)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

 k1 = (1.6*dt)/d_x^2; 
 k2 = (1.6*dt)/d_y^2; 
 term1 = (1/(1+2*k1+2*k2)); 
 term2 = (k1*term1);
 term3 = (k2*term1); 
 
   disp(" hey I am faster than Jacobi, jai gauss-siedel,let's get numbers...Hurrayyyy!!!")
   T_new = T;
   T_prev = T;
   tic
   for nt = 1:1000
       error = 9e9;
    while error > 1e-4     
      for i=2:nx-1
          for j = 2:ny-1
             T_new(i,j) = ((term1)*T_prev(i,j)) + ((term2)*(T(i+1,j) + T(i-1,j))) + ((term3)*(T(i,j+1) + T(i,j-1)));
          end         
      end %end for loop
      error = max(max(abs(T_new - T)));
      T = T_new;
      contourf(x,y,T,'ShowText','on')
      colorbar
      colormap(jet)
      ax = gca;  
      set(ax,'YDir','reverse')
      title('Transient Implicit Gauss siedel time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01)          
      Transient_g_s = toc;
    end %while loop end
     T_prev = T_new;
   end % Time loop end
    video_file = VideoWriter('method_6');
    open(video_file)
    writeVideo(video_file,ac)
    close(video_file)

Result :

Time : 2782.6s

Unsteady using Gauss-siedel Method :

                  *****  Unsteady Using Gauss siedel *****

close all
clear all
clc

nx = 20;
ny = 20;
nt = 1000;

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

d_x = x(2) - x(1);
d_y = y(2) - y(1);

y1 = 1;

error = 9e9;
tolerance = 1e-4;
dt = 1e-4;

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T = 300*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:ny-1) = T_B;

T(1,1) = (T_L + T_T)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

 k1 = (1.6*dt)/d_x^2; 
 k2 = (1.6*dt)/d_y^2; 
 term1 = (1/(1+2*k1+2*k2)); 
 term2 = (k1*term1);
 term3 = (k2*term1); 
 
   disp(" hey I am faster than Jacobi, jai gauss-siedel,let's get numbers...Hurrayyyy!!!")
   T_new = T;
   T_prev = T;
   tic
   for nt = 1:1000
       error = 9e9;
    while error > 1e-4     
      for i=2:nx-1
          for j = 2:ny-1
             T_new(i,j) = ((term1)*T_prev(i,j)) + ((term2)*(T(i+1,j) + T_new(i-1,j))) + ((term3)*(T(i,j+1) + T_new(i,j-1)));
          end         
      end %end for loop
      error = max(max(abs(T_new - T)));
      T = T_new;
      contourf(x,y,T,'ShowText','on')
      colorbar
      colormap(jet)
      ax = gca;  
      set(ax,'YDir','reverse')
      title('Transient Implicit Gauss siedel time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01)          
      Transient_g_s = toc;
    end %while loop end
     T_prev = T_new;
   end % Time loop end
    video_file = VideoWriter('method_6');
    open(video_file)
    writeVideo(video_file,ac)
    close(video_file)

Result :

Here time : 3100s , it took time for video.

But time of gauss siedel method is less than jacobi method.

Unsteady using SOR - g-s :

        ****** Unstaedy Heat Conduction using SOR-g-s *******

close all
clear all
clc

nx = 20;
ny = 20;
nt = 1000;

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

d_x = x(2) - x(1);
d_y = y(2) - y(1);

y1 = 1;

error = 9e9;
tolerance = 1e-4;
dt = 1e-4;

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T = 300*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:ny-1) = T_B;

T(1,1) = (T_L + T_T)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

 k1 = (1.6*dt)/d_x^2; 
 k2 = (1.6*dt)/d_y^2; 
 term1 = (1/(1+2*k1+2*k2)); 
 term2 = (k1*term1);
 term3 = (k2*term1); 
 
 b = (d_x/d_y)^2;
 a = ((cos(pi/(n-1)) + b*cos(pi/(n-1)) )/( 1 + b^2 ))^2;
 w_opt = (2-sqrt(1-a))/a;
 
   disp(" hey I am faster than Jacobi, jai gauss-siedel,let's get numbers...Hurrayyyy!!!")
   T_new = T;
   T_prev = T;
   tic
   for nt = 1:1000
       error = 9e9;
    while error > 1e-4     
      for i=2:nx-1
          for j = 2:ny-1
             T_new(i,j) = (1 - w_opt)*T(i,j) +  w_opt*(((term1)*T_prev(i,j)) + ((term2)*(T(i+1,j) + T_new(i-1,j))) + ((term3)*(T(i,j+1) + T_new(i,j-1))));
          end         
      end %end for loop
      error = max(max(abs(T_new - T)));
      T = T_new;
      contourf(x,y,T,'ShowText','on')
      colorbar
      colormap(jet)
      ax = gca;  
      set(ax,'YDir','reverse')
      title('Transient Implicit Gauss siedel time @',toc)
      ac(y1) = getframe(gcf);
      y1 = y1 + 1; 
      pause(0.01)          
      Transient_g_s = toc;
    end %while loop end
     T_prev = T_new;
   end % Time loop end
    video_file = VideoWriter('method_6');
    open(video_file)
    writeVideo(video_file,ac)
    close(video_file)

Result :

Time : 1753.0

Method Time (s)

Steady State Equation Jacobi 280.0793

Gauss siedel 139.95

SOR-g-s 43.116

Unsteady state Equation Jacobi 2782.6

Implicit Method Gauss siedel 2782.6< t < 1753.

SOR-g-s 1753.0

Explicit Method 268.9227

Conclusions :

From the above results can say SOR would be effcient method than jacobi and gauss-siedel iterative 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.