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

Mid Term Project - Solving the Steady and Unsteady 2D heat conduction problem

AArvind Eashwar Ramesh Babu5 Dec 2019

REPORT:

AIM:

To solve the conduction equation using a Transient solver and a Steady state solver using Iterative techniques (Jacobi,Gauss Seidal,SOR).

EXPLANATION:

  • For a steady state, the limit beyond which there will be no more changes obtained in its value is found.
  • In transient state, \'dt\' value changes, so their will be change in values.
  • In a steady state, there\'s no need for iterative solvers because we only use one single equation to solve.
  • Therefore the fact is that, the answers are going to be the same, but there are three different types involved.

SOLUTION:

I. STEADY STATE FLOW:

The 2 dimentional temperature equation:

`(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`

Discretizing the equation using Central Differencing Method

a. w.r.t \'x\'

`(del^2T)/(delx^2)=(T(i-1,j)-2T(i,j)+T(i+1,j))/(Deltax^2)`

b. w.r.t \'y\'

`(del^2T)/(delx^2)=(T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)`

CODE:

clear all
close all
clc

% Input Values
% Dpmain isa unit square
L = 1;
% nx = ny (given condition)
% Assuming value for nx
nx = 30;
ny = nx;
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; 
alpha = 1.8;

% Boundary Conditions (given)
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 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(nx,2:ny-1) = T_B;

% Calculating Avg. Temperature
T(1,1) = (T_T + T_L)/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;

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

iterative_solver = 1:3;
% Jacobi Method
for iterative_solver = 1
  tic;
   if iterative_solver == 1
      jacobi_iter = 1;
      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;
      end
   end
F = toc;

% PLOTTING
figure(1)
[M,N] = contour(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar    
title_text = sprintf(\'Jacobi Iteration No. = %d, Total Time F = %fs\', jacobi_iter, F);
title(title_text)
xlabel(\'X axis\')
ylabel(\'Y axis\')
fprintf(\'No. of iteration in jacobi = %d\\n\', jacobi_iter);
end

% GS ITERATION
for iterative_solver = 2
  L = 1;
  nx = 30;
  ny = nx;
  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;
  
  % Boundary Conditions (given)
  T_L = 400;
  T_T = 600;
  T_R = 800;
  T_B = 900;
  T = 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(nx,2:ny-1) = T_B;

  % Calculating Avg. Temperature
  T(1,1) = (T_T + T_L)/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;

  [X,Y] = meshgrid(x,y);
  Told = T;
  
  tic;
   if iterative_solver == 2
      gs_iter = 1;
      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;
      end
   end
F = toc;

% Plotting
figure (2)
[M,N] = contour(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar    
title_text = sprintf(\'GS Iteration No. = %d, Total Time F = %fs\', gs_iter, F);
title(title_text)
xlabel(\'X axis\')
ylabel(\'Y axis\')
fprintf(\'No. of iteration in GS = %d\\n\', gs_iter);
end

% SOR ITERATION
for iterative_solver = 3
  L = 1;
  nx = 30;
  ny = nx;
  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;
  
  % Boundary Conditions (given)
  T_L = 400;
  T_T = 600;
  T_R = 800;
  T_B = 900;
  T = 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(nx,2:ny-1) = T_B;

  % Calculating Avg. Temperature
  T(1,1) = (T_T + T_L)/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;

  [X,Y] = meshgrid(x,y);
  Told = T;
  
  tic;
   if iterative_solver == 3
      sor_iter = 1;
      while(error > tol)
         for i=2:nx-1
           for j = 2:ny-1
             T(i,j) = ((1-alpha)*Told(i,j) + (alpha*0.25*((T(i,j-1) + Told(i,j+1) + T(i-1,j) + Told(i+1,j)))));
           end
         end
       error = max(max(abs(Told - T)));
       Told = T;
       sor_iter = sor_iter + 1;
      end
   end
F = toc;

% Plotting
figure (3)
[M,N] = contour(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar    
title_text = sprintf(\'SOR Iteration No. = %d, Total Time F = %fs\', sor_iter, F);
title(title_text)
xlabel(\'X axis\')
ylabel(\'Y axis\')
fprintf(\'No. of iteration in SOR = %d\\n\', sor_iter);
    
end

OUTPUT:

1. Jacobi Iteration:

2. GS Iteration:

3. SOR Iteration:

II. UNSTEADY STATE FLOW:

1. Discretizing the 2D equation by Central Differencing:

a. w.r.t \'x\'

`(del^2T)/(delx^2)=(T(i-1,j)-2T(i,j)+T(i+1,j))/(Deltax^2)`

b. w.r.t \'y\'

`(del^2T)/(delx^2)=(T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)`

2. Discritizing Time Gradient Equation by Forward Differencing Method:

`(delT)/(delt)=(T^(n+1)(i,j) - Tn(i,j))/(Deltat)`

3. Constant k:

`k1 = alpha(Deltat)/(Deltax^2)`

`k2 = alpha(Deltat)/(Deltay^2)`

4. As we assumed nx=ny, we take

`(Deltax^2)/(Deltay^2)`

Therefore, we get,

`T^(n+1)(i,j)=(T^n(i,j)-k1(T(i-1,j)+T(i+1,j))^(n+1) + k2(T(i,j-1))+T(i,j+1))^(n+1)/(1+2K^4+2K^5)`

CODE(EXPLICIT):

% Explicit Solution
clear all
close all
clc

% Input Values
L = 1;
nx = 30;
ny = nx;
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; 
dt = 10;
alpha = 1.88e-5;

% Boundary Conditions (given)
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 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(nx,2:ny-1) = T_B;

% Calculating Avg. Temperature
T(1,1) = (T_T + T_L)/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;

% Constant
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

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

% Explicit Method:
iterative_solver = 1;
tic;
if iterative_solver == 1
  solver_counter = 1;
  while(error > tol)
    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 = k2*(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;
    solver_counter = solver_counter + 1;
  end  
endif
F = toc;

% Plotting
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'Transient iteration number(Explicit)=%d. Total Time(F)=%fs\', solver_counter, F);
title(title_text);
xlabel(\'x axis\')
ylabel(\'y axis\')
fprintf(\'No. of iteration(Explicit)=%d\\n\', solver_counter);

OUTPUT:

III. IMPLICIT METHOD:

% Implicit Solution
close all
clear all
clc

% Input Values
L = 1;
nx = 30;
ny = nx;
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;
dt = 2000;

% Thermal diffucivity
% For steel
alpha = 1.88e-5;
% For SOR
alpha_1 = 1.8;

% Boundary Conditions (given)
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 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(nx,2:ny-1) = T_B;

% Calculating Avg. Temperature
T(1,1) = (T_T + T_L)/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;

% Constant
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

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

% Implicit Method
iterative_solver = 1:3;
% Jacobi Iteration
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

for iterative_solver = 1
tic;
 if iterative_solver == 1
   jacobi_iter = 1;
    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(i,j) = term1*T_Prev_dt(i,j) + term2*H + term3*V;
        end
      end
      error = max(max(abs(Told-T)));
      Told = T;
      jacobi_iter = jacobi_iter+1;
    end
 end
F=toc;

% Plotting
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'Jacobi Unsteady(Implicit) iteration number =%d ,Total Time F =%fs\',jacobi_iter, F);
title(title_text)
xlabel(\'x axis\')
ylabel(\'y axis\')
fprintf(\'No. of iteration in jacobi = %d\\n\', jacobi_iter);
end

% GS iteration
% Input Values
L = 1;
nx = 30;
ny = nx;
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;
dt = 2000;
alpha = 1.88e-5;

% Boundary Conditions (given)
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 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(nx,2:ny-1) = T_B;

% Calculating Avg. Temperature
T(1,1) = (T_T + T_L)/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;

% Constant
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

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

term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

for iterative_solver = 2
tic;
 if iterative_solver == 2
  GS_iter = 1;
    while(error > tol)
      for i=2:nx-1
        for j=2:ny-1
          H = (T(i-1,j) + Told(i+1,j));
          V = (T(i,j-1) + Told(i,j+1));
          T(i,j) = term1*T_Prev_dt(i,j) + term2*H + term3*V;
        end
      end
      error = max(max(abs(Told-T)));
      Told = T;
      GS_iter = GS_iter+1;
    end
 end
F = toc;

%Plotting
figure(2)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'GS Transient(Implicit) iteration number = %d , Total Time F=%fs\', GS_iter, F);
title(title_text)
xlabel(\'xaxis\')
ylabel(\'yaxis\')
fprintf(\'No of Iteration in GS =%d\\n\',GS_iter);
end

%SOR
% Input Values
L = 1;
nx = 30;
ny = nx;
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;
dt = 2000;
alpha1 = 1.88e-5;

% Boundary Conditions (given)
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 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(nx,2:ny-1) = T_B;

% Calculating Avg. Temperature
T(1,1) = (T_T + T_L)/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;

% Constant
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

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

term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

for iterative_solver = 3
tic;
 if iterative_solver == 3  
   SOR_iter = 1;
    while(error > tol)
      for i=2:nx-1
        for j=2:ny-1
          H = (alpha1*T(i-1,j) + (1-alpha1)*Told(i-1,j) + Told(i+1,j));
          V = (alpha1*T(i,j-1) + (1-alpha1)*Told(i,j-1) + Told(i,j+1));
          T(i,j) = T_Prev_dt(i,j)*term1 + (term2*H + term3*V);
        end 
      end
      error = max(max(abs(Told-T)));
      Told = T;
      SOR_iter = SOR_iter+1;
    end
 end 
F=toc;
   
%Plotting
figure(3)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'SOR Transient(Implicit) iteration number =%d ,Total Time F=%fs\', SOR_iter, F);
title(title_text)
xlabel(\'xaxis\')
ylabel(\'yaxis\')
fprintf(\'No of Iteration in SOR =%d\\n\', SOR_iter);
end

OUTPUT:

1. Jacobi Method:

2. GS Method:

3. SOR METHOD:

IV. RESULTS:

1. EXPLICIT:

dx = 0.034483
dy = 0.034483
No. of iteration(Explicit)=2858

2. IMPLICIT:

No. of iteration in jacobi = 806
No of Iteration in GS =433
No of Iteration in SOR =806

From the obtained results, SOR(Successive Over Relaxation) method is time effective compared to other two methods (Jacobi and Gauss-Seidel Method)

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.