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

Solving 2D Heat Conduction Equation in both steady and transient case using iterative solvers

RRathish Gupta10 Jan 2020

Thie task is to solve 2D Heat Conduction equation

`(del T)/(del t) = ((del^2T)/(del x^2) +(del^2T)/(del x^2))*K`

where K is the thermal diffusivity

Assumptions in the equation: Internal heat generation is zero which implies a change of heat flow from one to another side is converted into internal energy.

Abstract: This challenge will focus on solving the equation with three solvers jacobian, gauss-seidel, and SOR iteration for two cases steady and transient. While the explicit method of solving is used to solve steady PDE and both explicit and implicit methods are used to solve transient PDE.

Equation of steady-state- since, the temperature doesn\'t change with time the equation reduces to

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

Steps

Discretization: The discretized equation of the steady-state is

`(T_(i+1,j) - 2* T_(i.j)+T_(i-1,,j))/dx^2 + (T_(i,j+1) - 2* T_(i.j)+T_(i,,j-1))/dy^2`

Input Variable:

  1. i) Nodes: we are going to give as a user input but the answers are shown for n=20
  2. ii) Length of the domain: This is taken as 1 along x and y-axis
  3. ii) Tolerance assumed 1e-4

Generating a 2D Plate profile

I have generated a 2D profile of temperature which gets analyzed. The temperature at 4 boundaries are given as below

Top=400K

Bottom=900k

Left=600K

Right=800K

Solving

  1. i) Jacobian Solver

Here the algebraic equation is solved such that the initial conditions of T is used to solve every step

  1. ii) Gauss Siedel Solver

The newly discovered value of T is used to find out the temperature at the next node

iii) SOR: Here an extra term called relaxation parameter is multiplied based on the slope and gauss seidel generated term to solve the equation faster

Code for Steady State:

  clear all
close all
clc

%intial Conditions
nx = input(\'please give your node number\');
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;

error = 1e8;
tol = 1e-4;

%temp matrix
T= 300*ones(nx,ny)

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

T(:,1)= T_L;
T(1,:)= T_T;
T(:,ny)= T_R;
T(nx,:)= T_B;

%edge refining 

T(1,1) = (600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (400+900)/2;
T(end,end) = (800+900)/2;

Told=T;

solver = input(\'enter the desired number\');

tic;
%% Gauss Jacobi Iteration Scheme
if solver == 1
 iter_count_jac = 1;
 while(error>tol)
  for i=2:nx-1
      for j=2:ny-1
      term1 = (Told(i-1,j)+Told(i+1,j))/4;
      term2 = (Told(i,j-1)+Told(i,j+1))/4;
      T(i,j)= term1+term2;
      end
  end
  error = max(max(abs(Told-T)));
  Told=T;
  iter_count_jac = iter_count_jac + 1;
end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d n computation time: %0.3f s\',iter_count_jac,time_required)
title(title_text);
end

tic;
%% Gauss Siedel Iteration Scheme
if solver == 2
 iter_count_GS = 1;
 while(error>tol)
  for i=2:nx-1
      for j=2:ny-1
      term1 = (T(i-1,j)+Told(i+1,j))/4;
      term2 = (T(i,j-1)+Told(i,j+1))/4;
      T(i,j)= term1+term2;
      end
      end  
  error = max(max(abs(Told-T)));
  Told=T;
  iter_count_GS = iter_count_GS + 1;
end
time_required=toc;
figure(2)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d n computation time: %0.3f s\',iter_count_GS,time_required)
title(title_text);
end


%% omega = 0.5;
%% SOR Iteration Scheme
omega = input(\'enter any of your fav number in between 0 to 2 your wish\');
if solver == 3
 iter_count_SOR = 1;
 while(error>tol)
  for i=2:nx-1
      for j=2:ny-1
      term1 = (T(i-1,j)+Told(i+1,j))/4;
      term2 = (T(i,j-1)+Told(i,j+1))/4;
      T(i,j)=(1-omega)* Told(i,j)+ omega*(term1+term2);  
      end
      end  
  error = max(max(abs(Told-T)));
  Told=T;
  iter_count_SOR = iter_count_SOR + 1;
end
time_required=toc;
figure(2)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d n computation time: %0.3f s\',iter_count_SOR,time_required)
title(title_text);
end

Results:

Jacobian:

Gauss-Seidel:

SOR:

As we move ahead with the method the iteration number gets reduced but the time of computing gets increased. The values remained the same in all cases. In SOR the computation has been carried out for 10 values between 0.5 to 2.0. The optimized value comes out to be 1.6

Transient - Implicit:

`((T_(i+1,j) - 2* T_(i.j)+T_(i-1,,j))/dx^2 + (T_(i,j+1) - 2* T_(i.j)+T_(i,,j-1))/dy^2)^n*K = (T_(i,j)^(n+1)-T_(i,j)^(n))/dt^2`

Here the unknowns of this equation is 1

rearranging it we can get the unknown T

The only difference with that of steady state is that here we are adding an extra loop of time.

clear all
close all
clc

%intial Conditions
nx = 10;
ny = nx;
nt= 50;
a = 0.001;
dt = 1;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;

error = 1e8;
tol = 1e-4;

%temp matrix
T= 400*ones(nx,ny);

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

T(:,1)= T_L;
T(1,:)= T_T;
T(:,ny)= T_R;
T(nx,:)= T_B;

%edge refining 

T(1,1) = (600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (400+900)/2;
T(end,end) = (800+900)/2;

Told=T;
solver = input(\'give a number\')
if solver == 1
tic;
%% Gauss Jacobi Iteration Scheme
iter_count_jac = 1;
k1 =  a*dt/dx^2;
for k = 1:nt
  for i=2:nx-1
      for j=2:ny-1
      term1 = Told(i,j);
      term2 = k1*(Told(i-1,j)- 2*Told(i,j)+Told(i+1,j)); 
      term3 = k1*(Told(i,j-1)- 2*Told(i,j)+Told(i,j+1)); 
      T(i,j)= term1+term2+term3;
      end
  end
  error = max(max(abs(Told-T)));
  Told=T;
  iter_count_jac = iter_count_jac + 1;
  end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d, computation time: %0.3f s,k = %d\',iter_count_jac,time_required,k)
title(title_text);
end



if solver == 2
tic;
%% Gauss Seidel Iteration Scheme
iter_count_GS = 1;
k1 =  a*dt/dx^2;
for k = 1:nt
  for i=2:nx-1
      for j=2:ny-1
      term1 = Told(i,j)- 2*k1*Told(i,j)- 2*k1*Told(i,j);
      term2 = k1*(T(i-1,j)+Told(i+1,j)); 
      term3 = k1*(T(i,j-1)+Told(i,j+1)); 
      T(i,j)= term1+term2+term3;
      end
  end
  error = max(max(abs(Told-T)));
  Told=T;
  iter_count_GS = iter_count_GS + 1;
  end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d, computation time: %0.3f s,k = %d\',iter_count_GS,time_required,k)
title(title_text);
end

if solver == 3
tic;
%% SOR Iteration Scheme
omega=1.5
iter_count_SOR = 1;
k1 =  a*dt/dx^2;
for k = 1:nt
  for i=2:nx-1
      for j=2:ny-1
      term1 = Told(i,j)- 2*k1*Told(i,j)- 2*k1*Told(i,j);
      term2 = k1*(T(i-1,j)+Told(i+1,j)); 
      term3 = k1*(T(i,j-1)+Told(i,j+1)); 
      T(i,j)=(1-omega)* Told(i,j)+ omega*(term1+term2+term3);
      end
  end
  error = max(max(abs(Told-T)));
  Told=T;
  iter_count_SOR = iter_count_SOR + 1;
  end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d, computation time: %0.3f s,k = %d\',iter_count_SOR,time_required,k)
title(title_text);
end

Results:

The values come out to be the same in all cases. Because the number of iterations is the same as the number of steps so there is no significance of using iterative solvers here

Transient Implicit:

clear all
close all
clc

%intial Conditions
nx = 10;
ny = nx;
nt= 50;
a = 1.2;
dt = 1e-2;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;

error = 1e8;
tol = 1e-4;

%temp matrix
T= 400*ones(nx,ny);

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

T(:,1)= T_L;
T(1,:)= T_T;
T(:,ny)= T_R;
T(nx,:)= T_B;

%edge refining 

T(1,1) = (600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (400+900)/2;
T(end,end) = (800+900)/2;

Told=T;
T_prev = T;
solver=input(\'give a number\')

if solver == 1
tic;
%% Gauss Jacobi Iteration Scheme
iter_count_jac = 1;
k1 =  a*dt/dx^2;
cons1 = 1+2*k1+2*k1;
for k = 1:nt
  while(error>tol)
   for i=2:nx-1
      for j=2:ny-1
      term1 = T_prev(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)/cons1);
      end
   end
   error = max(max(abs(Told-T)));
   Told=T;
   iter_count_jac = iter_count_jac + 1;
  end
  T_prev = T;
end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d, computation time: %0.3f s,k = %d\',iter_count_jac,time_required,k)
title(title_text);
end

if solver == 2
tic;
%% Gauss Seidel Iteration Scheme
iter_count_GS = 1;
k1 =  a*dt/dx^2;
cons1 = 1+2*k1+2*k1;
for k = 1:nt
  while(error>tol)
   for i=2:nx-1
      for j=2:ny-1
      term1 = T_prev(i,j);
      term2 = k1*(T(i-1,j)+Told(i+1,j)); 
      term3 = k1*(T(i,j-1)+Told(i,j+1)); 
      T(i,j)= ((term1+term2+term3)/cons1);
      end
   end
   error = max(max(abs(Told-T)));
   Told=T;
   iter_count_GS = iter_count_GS + 1;
  end
  T_prev = T;
end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d, computation time: %0.3f s,k = %d\',iter_count_GS,time_required,k)
title(title_text);
end

if solver == 3
tic;
%% Gauss Seidel Iteration Scheme
iter_count_SOR = 1;
omega = 0.5
k1 =  a*dt/dx^2;
cons1 = 1+2*k1+2*k1;
for k = 1:nt
  while(error>tol)
   for i=2:nx-1
      for j=2:ny-1
      term1 = T_prev(i,j);
      term2 = k1*(T(i-1,j)+Told(i+1,j)); 
      term3 = k1*(T(i,j-1)+Told(i,j+1)); 
      T(i,j)= (1-omega)*T_prev(i,j) + omega*((term1+term2+term3)/cons1);
      end
   end
   error = max(max(abs(Told-T)));
   Told=T;
   iter_count_SOR = iter_count_SOR + 1;
  end
  T_prev = T;
end
time_required=toc;
figure(1)
contourf(x,fliplr(y),T);
colorbar;
title_text = sprintf(\'iteration number=%d, computation time: %0.3f s,k = %d\',iter_count_SOR,time_required,k)
title(title_text);
end

Results:

The same results are found in steady-state although the nature of using iterative solvers is different here. The no.of iterations get down as we put iterative solvers. The optimum value of omega found out to be around 0.9 after trying values from 0.5 to 2.0.

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.