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

SOLVING THE STEADY AND UNSTEADY 2D HEAT CONDUCTION PROBLEM

SSunu Cherian2 Mar 2020

Aim

To solve 2D steady and unsteady heat equation.

Program & explanation

2D steady state equation

`(d^2T)/(dx^2) + (d^2t)/(dy^2) = 0`

In this challenge we are trying to solve 2D steady & transient heat conduction over a square, for this we use different iterative methods to find the temperature across the square and we also compare the number of iterations and the time consumed by each method.

clear all
close all;clc

%given values
L = 1;
nx = ny = 10;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;

%assigning values of temperature across the square
T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T_M = ones(nx,ny);
T_M(2:end,ny) = T_T;
T_M(2:end,1)= T_B;
T_M(1,2:end)= T_L;
T_M(nx,2:end)= T_R;



%Assigning corners 
T_M(1,1) = (T_T + T_L)/2;
T_M(1,ny) = (T_T + T_R)/2;
T_M(nx,1) = (T_L + T_B)/2;
T_M(nx,ny) = (T_R + T_B)/2



% required method of solving
% iterative_solver = 1 for Jacobi_iter
% iterative_solver = 2 for Gauss siedel
% iterative_solver = 3 for SOR
iterative_solver = 1 ;


%%Temperature matrix, for iterative solving
Told = T_M;

k_1 = (dx^2)/(2*(dx^2+dy^2));
k_2 = (dy^2)/(2*(dx^2+dy^2));

graphics_toolkit("gnuplot")

   % STEADY STATE, JACOBI

   
if iterative_solver == 1
    jacobi_iter = 1;
 tic   
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   T_M(i,j) = k_1*(Told(i,j-1)+Told(i,j+1))+k_2*(Told(i-1,j)+Told(i+1,j));
   end
  end
     error = max(max(abs(Told-T_M)));
     Told = T_M; 
     jacobi_iter = jacobi_iter + 1; 
 end
 time = toc   
     
     contourf(x, y, T_M);
     [c,h] = contourf(x, y, T_M);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', jacobi_iter);
     title(title_text)
     file_text = sprintf("steady_%d.png",iterative_solver);
     saveas(gca,file_text);
end

    
    
    % STEADY STATE, GS

    
if iterative_solver == 2
  GS_iter = 1;
  
 tic 
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   T_M(i,j) = k_1*(T_M(i,j-1)+T_M(i,j+1))+k_2*(T_M(i-1,j)+T_M(i+1,j));
   end
  end
     error = max(max(abs(Told-T_M)));
     Told = T_M; 
     GS_iter = GS_iter + 1;   
 end
  time = toc  
    
     contourf(x, y, T_M);
     [c,h] = contourf(x, y, T_M);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', GS_iter);
     title(title_text) 
     file_text = sprintf("steady_%d.png",iterative_solver);
     saveas(gca,file_text); 
end
    
    
    
     % STEADY STATE SOR

  
if iterative_solver == 3
   opt_omega = 2/(1+sin(pi*dx));
   SOR_iter = 1;
 tic  
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   T_M(i,j) =opt_omega*(k_1*(T_M(i,j-1)+T_M(i,j+1))+k_2*(T_M(i-1,j)+T_M(i+1,j)))+(1-opt_omega)*Told(i,j);
   end
  end
     error = max(max(abs(Told-T_M)));
     Told = T_M; 
     SOR_iter = SOR_iter + 1; 
 end
  time = toc   
      
     contourf(x, y, T_M);
     [c,h] = contourf(x, y, T_M);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', SOR_iter);
     title(title_text)
     file_text = sprintf("steady_%d.png",iterative_solver);
     saveas(gca,file_text);
end

Output

1.Jacobi

2.Gauss Siedel

3.SOR

Unsteady(Transient) 2D Heat Conduction Equation

`(d^2T)/(dt^2)=1/(alpha)((d^2T)/(dx^2)+(d^2T)/(dy^2))`

clear all
close all;clc

%given values
L = 1;
nx = ny = 10;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
dt = 32
error = 9e9;
tol = 1e-4;

%assigning values of temperature across the square
T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T_M = ones(nx,ny);
T_M(2:end,ny) = T_T;
T_M(2:end,1)= T_B;
T_M(1,2:end)= T_L;
T_M(nx,2:end)= T_R;



%Assigning corners 
T_M(1,1) = (T_T + T_L)/2;
T_M(1,ny) = (T_T + T_R)/2;
T_M(nx,1) = (T_L + T_B)/2;
T_M(nx,ny) = (T_R + T_B)/2

Told = T_M


% required method of solving
% iterative_solver = 1 for Jacobi_iter
% iterative_solver = 2 for Gauss siedel
% iterative_solver = 3 for SOR
iterative_solver = 3 ;


    
     % UNSTEADY IMPLICIT JACOBI
     
      alpha = 1e-4;
 k_3 = alpha*dt/(dx^2);
 k_4 = alpha*dt/(dy^2);  
     
     term1 = (1+2*k_3+2*k_4)^-1;
     term2 = k_3*term1;
     term3 = k_4*term1;

  
 
    
if iterative_solver == 1
tic

% Copy of the Temperature matrix 
  T_prev_dt = Told;
  jacobi_iter = 1;

% Time loop  
for k = 1:200  
 % Convergence loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (Told(i-1,j)+Told(i+1,j));
   V = (Told(i,j-1)+Told(i,j+1));
  T_M(i,j) = (T_prev_dt(i,j)*term1)+(term2*H)+(term3*V);
    end
   end
     error = max(max(abs(Told-T_M)));
     Told = T_M; 
     jacobi_iter = jacobi_iter + 1; 
 end
    
    T_prev_dt = T_M;
    error = 9e9;
end
time = toc
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_M);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', jacobi_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
    
    
    
   
    
     % UNSTEADY IMPLICIT GAUSS SIEDEL
     
     
     term1 = (1+2*k_3+2*k_4)^-1;
     term2 = k_3*term1;
     term3 = k_4*term1;
if iterative_solver == 2
tic
% Copy of the Temperature matrix 
  T_prev_dt = Told;
  GS_iter = 1;
  
% Time loop
 for k = 1:200
 % Convergence loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (T_M(i-1,j)+Told(i+1,j));
   V = (T_M(i,j-1)+Told(i,j+1));
  T_M(i,j) = (T_prev_dt(i,j)*term1)+(term2*H)+(term3*V);
    end
   end
     error = max(max(abs(Told-T_M)));
     Told = T_M; 
     GS_iter = GS_iter + 1; 
  end
   
    T_prev_dt = T_M;
    error = 9e9;
    
 end
time = toc 
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_M);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', GS_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
    
    
    
    
    
     % UNSTEADY IMPLICIT GAUSS SIEDEL WITH OVER RELAXATION
     
     
     term1 = (1+2*k_3+2*k_4)^-1;
     term2 = k_3*term1;
     term3 = k_4*term1;
     
if iterative_solver == 3
tic
  opt_omega = 2/(1+sin(pi*dx));
  
% Copy of the Temperature matrix  T_prev_dt = Told;
  GS_iter = 1;
  
% Time loop
 for k = 1:200
 % Convergence loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (T_M(i-1,j)+Told(i+1,j));
   V = (T_M(i,j-1)+Told(i,j+1));
  T_M(i,j) = opt_omega*((T_prev_dt(i,j)*term1)+(term2*H)+(term3*V))+(1-opt_omega)*T_prev_dt(i,j);
    end
   end
     error = max(max(abs(Told-T_M)));
     Told = T_M; 
     GS_iter = GS_iter + 1; 
  end
   
    T_prev_dt = T_M;
    error = 9e9;
    
 end 
 time = toc
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_M);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', GS_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end

Outputs

1.Jacobi Method

2.GS Method

3.SOR

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.