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

2D heat conduction equation Using Point Iterative techniques

AAditya Soni8 Oct 2018

Given Problem Statement- Solve the 2D heat conduction equation by using the point iterative techniques:-

Implement the following methods

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Absolute error criteria is 1e-4

Problem :
1. The domain is a unit square.
2. The number of points along the x-direction is equal to the number of points along the y-direction.
3. Following are the boundary conditions used:
I. Left:400K II. Top:600K III. Right:800K IV. Bottom:900K

1. Jacobi Method=

Jacobi Method:

then successive approximations are:

The process can be continued with the scheme

If the vectors u(m) approach a limit and the difference of successive vectors approaches the zero vector then it seems clear that a solution of the system will be found.

Snippet for Steady state case

T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));​

Snippet for Unstaedy Implicit case

term1 = (1+(2*k1)+(2*k2))^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                H = (Told(i-1,j)+Told(i+1,j));
                V = (Told(i,j-1)+Told(i,j+1));
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);

2. Gauss-siedel=n this method, rather than computing the new vectors u(m) entirely from the previous vector u(m-1) , we could use a new element as soon as we have it because if convergence takes place the new elements should be more accurate than the old ones .

The Gauss-Seidel algorithm usually converges much faster than the Jacobi method. Even though we can define both methods for any matrix, the method may or may not converge. Even if it does converge, it may do so very slowly and prove to be of no practical use. Fortunately, for many problems similar to heat conduction, these methods and their variation do effectively

converge.

Snippet for Steady state case

T(i,j) = 0.25*(T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1));

Snippet for unsteady implicit case

 term1 = (1+(2*k1)+(2*k2))^-1;    %Gauss siedel
                term2 = k1*term1;
                term3 = k2*term1;
                H = (T(i-1,j)+T(i+1,j));
                V = (T(i,j-1)+T(i,j+1));
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);

3. Successive Over realaxation= Its a variation of gauss siedel method which has an acceleration parameter (Omega)rf in our case. Here the choice of the parameter should be between 1 and 2 so that convergence is as rapid as possible. For very special matrices there are formulae for such optimal, but generally the optimal values are approximated by virtue of experience. Also the initial guess should be as close as possible to the solution, and in this case one may rely on the nature of the solution as dictated by the physical problem that is being modelled.

Snippet for Steady state case

T(i,j) = rf*(0.25*(T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1))) + (1-rf)*Told(i,j);

Snippet for Unsteady Implicit case

 term1 = (1+(2*k1)+(2*k2))^-1;     %SOR
                term2 = k1*term1;
                term3 = k2*term1;
                H = (T(i-1,j)+T(i+1,j));
                V = (T(i,j-1)+T(i,j+1));
                T(i,j) = rf+(1-rf)*(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);

4. Unsteady explicit equation:-

Solving explicit equation requires less time, than solving the implicit equation.

snippet

  T(i,j) = Told(i,j)+k3*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j))+k4*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));

Hereby attaching The plots with Information:-

Series of attachment:-

1. Jacobi Method Steady state

2. Gauss siedel steady state

3. SOR explicit steady state

4. Jacobi Method Implicit Form

5. Gauss Siedel Implicit form

6. SOR implicit form

7. Unsteady explicit form

7. steady state form code

8. Unsteady form code

Elapsed time is 343.062 seconds.

No of iteration taken is 869

Time taken for solving is 343.063193

Elapsed time is 201.556 seconds.

No of iteration taken is 462

Time taken for solving is 201.557353

Elapsed time is 190.694 seconds.

No of iteration taken is 461

Time taken for solving is 190.695314

Result= 1. SOR method is the fastest as it is taking very less time comapared to other two methos.

2. SOR is just a slightly modification of Gauss siedel method so the no. of iterations in this example are nearly same but SOR is taking less time comapre to gauss siedel.

Implicit Form:-

IMPLICIT - Influence of both time loop and convergence loop is experienced. The rate of convergence is modified by using the above mentioned three iterative methods(jacobi, GS, SOR).

Elapsed time= 6.9889 seconds

Number of iterations taken= 2209

Time taken for solving it= 7.4714 seconds

Elapsed time= 5.24119 seconds

Number of iterations taken= 1507

Time taken for solving it= 5.802792 seconds

Elapsed time=3.49795 seconds

Number of iterations taken=788

Time taken for solving it= 4.32323 seconds

Results- 1. SOR is again fastest method to solve the heat conduction equation over a given period of time step.

UNSTEADY STATE EQUATION:

1. EXPLICIT - It is solved using the concept of time marching, and doesn\'t require formation of any convergence loop.

Elapsed time is 1.60286 seconds.

No of iteration taken is 1

Time taken for solving is 4.386900

Script for Steady state condition:-

%--------------------------------Solving 2D Heat Conduction equation-------------------------------------%

%---------------------------------Using explicit Method--------------------------------------------------%

%---------------------------------Steady State analysis--------------------------------------------------%

%Given equation= dT/dt-(alpha)((d^2T/dx^2)+(d^2T/dy^2))=0

clc;
close all;
clear all;

%Inputs from question
%nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
% Assume that the domain is a unit square. 
% Boundary conditions for steady and transient case:-

    %1. Left:400K

    %2. Top:600K

    %3. Right:800K

    %4. Bottom:900K
    
Lx=10;    %TAking length of grid in X direction

Ly=Lx;    %Taking Length of grid in  y direction

nx=20;    %taking no. of nodes in x direction

ny=nx;    %taking no. of nodes in y direction

nt=40;




tic

%Taking Step sizes now:-

dx=Lx/(nx-1); %for x coordinate

dy=Ly/(ny-1); %for x coordinate

%setting loop for iteration solution:-

I=3; %or 1 and 2

%Entering other necessary values:-

tol=1e-4;   %Provided in question

error=9e9;  %Initialising the error for running iterations

%Grids formation for 1D and 2D

x=linspace(0,Lx,nx);

y=linspace(0,Ly,ny);


dx = x(2)-x(1);
dy = y(2)-y(1);

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


%combining grids to create 2d grid

[X Y]=meshgrid(x,y);
 
%assisgning boundary conditions

T=ones(nx,ny);  %Temp array

%Boundary Condition, 900k=900,000

T(1,:)=900;       %T at bottom phase
T(end,:)=600;       %T at top phase
T(:,1)=400;       %T at left phase
T(:,end)=800;       %T at right phase

%copy of temperature array
Told=T;

K=1;  %Iteration starting

rf=2/(1+sin(pi*dx));    %rf=optimum relaxation factor


%starting of loops-----------------------------------------------------------------------------------%




for f=1:nt 

 while error>tol       %condition for loop run until error is greater than tolerance.
  
  for i=2:nx-1;   %X value loop
    
    for j=2:ny-1; %Y value loop
      
      
      
      %Iterative Solving start using all 3 methods. USing if loop to run all three method using all 3 values of I array 
      
      if I==1     %jacobi
       
        
      
         T(i,j)=  k1*(Told(i,j-1)+Told(i,j+1))+k2*(Told(i-1,j)+Told(i+1,j));
        
      elseif I==2     %gauss-siedel
       
        
        
         T(i,j) = k1*(T(i,j-1)+T(i,j+1))+k2*(T(i-1,j)+T(i+1,j));
          
      elseif I==3     %Successive over relaxation
       
        
        
          
       T(i,j) =rf*(k1*(T(i,j-1)+T(i,j+1))+k2*(T(i-1,j)+T(i+1,j)))+(1-rf)*Told(i,j);     
          
        end
        
      end
      
    end

error1=max(abs(Told-T));  %error calculation for loop condition satisfaction


error=max(error1);   %error between old array and new array of temperature T.

Told=T;  %(Updating temperature after loop)

K=K+1;  %incrementation of k for all 3 methods

#figure(1)
#contourf(T)
#pause(0.003)
%Creating Contour plot with temp colorbar
[C,h]=contourf(X,Y,T);            %Contur function 

set(gca,\'Xtick\',(0:dx:1))

set(gca,\'Ytick\',(0:dy:1))

clabel(C,h,\'FontSize\',8,\'FontWeight\',\'bold\',\'color\',\'k\')

xlabel(\'X co-ordinate\',\'FontSize\',12)

ylabel(\'Y co-ordinate\',\'FontSize\',12)

title(sprintf(\'Temperature distribution on %d x %d grid using SOR method at Iteration no. %d\',Lx,Ly,K),\'FontSize\',15);

pause(0.003)



  end

end

toc





      





  
% % Iteration and Time display
display(\'2D heat conduction explicit method\');
fprintf(\'No of iteration taken is %d\\n\',K);
fprintf(\'Time taken for solving is %f\\n\',toc);


if I==1
  ct_j=toc;
  fprintf(\'%d Jacobi method iteration for computational time of %0.3 second\',K,ct_j)
  
  elseif I==2
  ct_gs=toc;
  fprintf(\'%d Gauss Siedel method iteration for computational time of %0.3 second\',K,ct_gs)
  
  elseif I==3
  ct_sor=toc;
  fprintf(\'%d SOR method iteration for computational time of %0.3 second\',K,ct_sor)
  
end

For unsteady state

%--------------------------------Solving 2D Heat Conduction equation-------------------------------------%

%---------------------------------Using explicit and Implicit Method--------------------------------------------------%

%---------------------------------Unsteady State analysis--------------------------------------------------%

%Given equation= dT/dt-(alpha)((d^2T/dx^2)+(d^2T/dy^2))=0


clc
clear all
close all
% Space Initialization 
I=3; %or 1 or 2
Lx = 1; Ly = 1;
nx = 11;
ny = 11;
dx = Lx/(nx-1); dy = Ly/(ny-1);
 % Time Initialization
time = 0.2;
k = 1.4;             % Thermal Diffusivity
dt = 1e-2;
nt = 500;
k1 = k*dt/(dx^2);
k2 = k1;
 % Mesh
x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
[X, Y] = meshgrid(x,y);
 % Limits
tol = 1e-4;
error = 9e9;
 % Temp Initialization
T = ones(nx,ny);
 % Boundary conditions
 


T(:,1) = 400; % Left
T(:,end) = 800; % Right
T(1,2:end-1) = 600; % Top
T(end,2:end-1) = 900; % Bottom
 % Copying T
Told = T;
 % Implicit method
iter = 0;
T_prev_dt = Told;
rf=2/(1+sin(pi*dx)); %SOR const.
tic
 % Loop
for k = 1:nt
  error=9e9;
 while(error > tol)
 
 
        for i = 2:nx-1
            for j = 2:ny-1
              
                 if I==1     %jacobi
                term1 = (1+(2*k1)+(2*k2))^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                H = (Told(i-1,j)+Told(i+1,j));
                V = (Told(i,j-1)+Told(i,j+1));
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
                
                    
                elseif I==2     %gauss-siedel
        
          
                term1 = (1+(2*k1)+(2*k2))^-1;    %Gauss siedel
                term2 = k1*term1;
                term3 = k2*term1;
                H = (T(i-1,j)+T(i+1,j));
                V = (T(i,j-1)+T(i,j+1));
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
                
                  elseif I==3     %Successive over relaxation
                term1 = (1+(2*k1)+(2*k2))^-1;     %SOR
                term2 = k1*term1;
                term3 = k2*term1;
                H = (T(i-1,j)+Told(i+1,j));
                V = (T(i,j-1)+Told(i,j+1));
                T(i,j) = rf+(1-rf)*(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
                  
                

                end
            end
        end          
            error = max(max(abs(Told - T)));
            Told = T;
            iter = iter+1;
            
           
  
            
  end
  

 % Plotting
 
 
 
 
 %Looping the previous value
##    Told = T;
    T_prev_dt = T;

end


toc
 #figure(1)
   %contourf(T);
   %pause(0.003)
 %Plotting
  [C,h]=contourf(X,Y,T);            %Contur function 

set(gca,\'Xtick\',(0:dx:1))

set(gca,\'Ytick\',(0:dy:1))

clabel(C,h,\'FontSize\',8,\'FontWeight\',\'bold\',\'color\',\'k\')

xlabel(\'X co-ordinate\',\'FontSize\',12)

ylabel(\'Y co-ordinate\',\'FontSize\',12)

title(sprintf(\'Temperature distribution on %d x %d grid using SOR method at Iteration no. %d\',Lx,Ly,iter),\'FontSize\',15);

pause(0.003)


 

  




  
  
% % Iteration and Time display
display(\'2D heat conduction Implicit method\');
fprintf(\'No of iteration taken is %d\\n\',iter);
fprintf(\'Time taken for solving is %f\\n\',toc);





 
 %--------------------------------Explicit form---------------------------------------


clc
clear all
close all
% Space Initialization 
I=1;
Lx = 1; Ly = 1;
nx = 11;
ny = 11;
dx = Lx/(nx-1); dy = Ly/(ny-1);
 % Time Initialization
time = 0.2;
            
dt=1e-2;
nt = 500;

 % Mesh
x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
[X, Y] = meshgrid(x,y);


 alpha = 1e-2;
 k3 = alpha*dt/(dx^2);
 k4 = alpha*dt/(dy^2);
 % Limits
tol = 1e-4;
error = 9e9;
 % Temp Initialization
T = ones(nx,ny);
 % Boundary conditions
 


T(:,1) = 400; % Left
T(:,end) = 800; % Right
T(1,2:end-1) = 600; % Top
T(end,2:end-1) = 900; % Bottom
 % Copying T
Told = T;
 % Implicit method
iter = 1;
T_prev_dt = Told;

tic
% UNSTEADY EXPLICIT 
   
if I == 1
tic
% Time marching solution 
% With dt = 32, the steady state is obtained at less no.of iterations @ k = 200, with the negligible error of 1.7071e-004.
% Therefore the solution resembles the solution of the steady state solution
 for k = 1:nt
  for i = 2:nx-1
   for j = 2:ny-1
  T(i,j) = Told(i,j)+k3*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j))+k4*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
   end
  end
     Told = T;  
 end
 


end
    toc
    
    [C,h]=contourf(X,Y,T);            %Contur function 

set(gca,\'Xtick\',(0:dx:1))

set(gca,\'Ytick\',(0:dy:1))

clabel(C,h,\'FontSize\',8,\'FontWeight\',\'bold\',\'color\',\'k\')

xlabel(\'X co-ordinate\',\'FontSize\',12)

ylabel(\'Y co-ordinate\',\'FontSize\',12)

title(sprintf(\'Temperature distribution on %d x %d grid using unsteady explicit method at Iteration no. %d\',Lx,Ly,iter),\'FontSize\',15);

pause(0.003)


 

  




  
  
% % Iteration and Time display
display(\'2D heat conduction Implicit method\');
fprintf(\'No of iteration taken is %d\\n\',iter);
fprintf(\'Time taken for solving is %f\\n\',toc);
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.