Solving the steady and unsteady 2D heat conduction problem

AIM

The aim of this project is to solve the 2 Dimensional equation using iterative solvers.

OBJECTIVE 

The main objective is to write an program to solve 2D steady and unsteady state equations by using following iterative solvers

• Jacobi
• Gauss Seidel
• Successive Over Relaxation

TWO DIMENSIONAL HEAT CONDUCTION (STEADY STATE)

Conduction is the mode of energy transfer as heat due to temperature difference in solids where the mass is contiguous. Steady state heat conduction occurs when the temperature difference driving the conduction of heat is constant. Amount of heat entering and the amount of heat leaving the system should be same for steady state heat conduction. In this problem we are taking a square plate in which there is a steady state heat conduction whose boundary conditions are 

• Left: 400K
• Top: 600K
• Right: 800K
• Bottom: 900K

The two dimensional heat steady state equation is given by the Laplace equation

`frac[del^2T][delx^2]+frac[del^2T][dely^2]=0`

To solve the above equation in a computer program we have to first discretise the two dimensional steady state heat conduction which is given by

`T_p=1/k[frac[T_l+t_r][trianglex^2]]+1/k[frac[T_t+T_b][triangley^2]]`

where,

T_p = temperature at point p which is to be calculated

T_l= temperature at the left node

T_r= temperature at the right node

T_t= temperature at the top node

T_b= temperature at the bottom node

The above equation can be written in the generalised node format 

`T_p=1/k[frac[T_(i-1,j)+T_(i+1,j)][trianglex^2]]+1/k[frac[T_(i,j-1)+T_(i,j+1)][triangley^2]]`

MATLAB CODE

clear all;
close all;
clc;
% Input for the steady state conduction  
L =1; 
nx = 10;
ny = 10;
x =linspace(0,1,nx);
y =linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
  
% Boundaty Conditions 
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T = ones(nx,ny);

%Temperature at boundary
 T(:,1)=T_left;
 T(:,end)=T_right;
 T(end,:)=T_top;
 T(1,:)=T_bottom;
 
%creating a copy of T
[X,Y] = meshgrid(x,y);
Told =T;
  
%Jacobi Iterative Solver
 for iterative_solver =1;
tic;
 if iterative_solver ==1
    jacobi_iterations =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_iterations = jacobi_iterations+1;
end
F = toc;
%Jacobi plot for steady state conduction
figure(1);
[P,Q] = contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('Jacobi iterative solver,Iteration number=%d,computation time F=%fs',jacobi_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf( 'jacobi=%d/n',jacobi_iterations);
pause(0.03);
end
    
%Gauss Seidel Solver
 for iterative_solver =2;
L =1; 
nx = 10;
ny = 10;
x =linspace(0,1,nx);
y =linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
  
%Boundary Conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T = 298*ones(nx,ny);
  
%Temperature at boundary
 T(:,1)=T_left;
 T(:,end)=T_right;
 T(end,:)=T_top;
 T(1,:)=T_bottom;
 
%creating a copy of T
[X,Y] = meshgrid(x,y);
Told =T;
tic;
    if iterative_solver ==2
    gs_iterations =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_iterations = gs_iterations+1;
end
end
F = toc;
% Gauss Seidel plot for steady state conduction
figure(2);
[P,Q] = contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf( 'Gauss seidel iterative solver,Iteration number=%d,computation time F=%fs',gs_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf( 'gs=%d/n',gs_iterations);
pause(0.03);
end
%Successive Over Relaxation
for iterative_solver =3;
L =1; 
nx = 10;
ny = 10;
x =linspace(0,1,nx);
y =linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
  
% Boundary Conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T = ones(nx,ny);
% calculation of relaxation factor
omega = 2/(1+sin(pi*dx));
  
%Temperature at boundary
T(:,1)=T_left;
T(:,end)=T_right;
T(end,:)=T_top;
T(1,:)=T_bottom;
 
%creating a copy of T
[X,Y] = meshgrid(x,y);
Told =T;
tic;

if iterative_solver ==3
    SOR_iterations =1;
    while(error >tol)
    for i=2:nx-1
     for j=2:ny-1
      T(i,j) =((1-omega)*Told(i,j)+(omega*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;
    SOR_iterations = SOR_iterations+1;
end
end
F=toc;
% SOR plot for steady state conduction 
figure(3);
[P,Q] = contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf( 'SOR Iterative solver, Iteration number=%d,computation time F=%fs',SOR_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf( 'SOR=%d/n',SOR_iterations);
pause(0.03);
end
end

  We have to define the square plate of length 1m with 10 nodes so we define some variables in the start of the program.

• Length of the plate (l) = 1m
• number of nodes on x-axis (nx) = 10
• number of nodes on y-axis (ny) = 10
• X is the array of nodes in the x axis
• Y is the array of nodes in y axis.
• error term for temperature= 9e9
• tolerance for the temperature criteria = 1e-4

After defining the square plate we must assign the boundary conditions at which the temperatures at the boundaries is assigned as

• Left: 400K
• Top: 600K
• Right: 800K
• Bottom: 900K

The above condition should be applied to all the nodes in the outer boundary so the next 4 lines indicates that we are assigning the same boundary conditions for all the outer nodes. When we come to the Jacobi method portion of the code, here we ae defining the iterative solver if iterative solver is equal to 1 then, we are basically saying the code to execute Jacobi method. So we have two if statements for the space loop and time loop and the "while" statement tells the condition to repeat the iterations until the tolerance value is reached. Once we are running Jacobi method, we have to update the temperature values and again run the iterations till we reach the convergence in the solution. We then use "countourf" command to execute the plot for the temperature distribution within the square plate.

Similarly, the same procedure is implemented for gauss seidel and successive over relaxation methods.

RESULTS AND DISCUSSION

Jacobi method for solving the 2D heat conduction is as shown in the above figure.Here we can observe that the jacobi method has taken 217 iterations to reach the convergence in the solution and the duration for which the simulation ran is for 0.367446 seconds

Gauss Seidel method for solving the 2D heat conduction is as shown in the above figure. Here we can observe that the Gauss Seidel method has taken 112 iterations to reach the convergence in the solution and the simulation just ran for 0.214730 seconds

Successive Over Relaxation

Successive Over Relaxation for solving the 2D heat conduction is as shown in the above figure. Here we can observe that the Successive Over Relaxation method has taken 31 iterations to reach the convergence in the solution and we see that the convergence is reached by 0.072217 seconds.

TWO DIMENSIONAL HEAT CONDUCTION (UNSTEADY/TRANSCIENT STATE)

When there are any changes in temperature with respect to time within the body then mode of heat conduction is known as unsteady or transient heat conduction. The governing equation for transient conduction is 

`frac[del T][delt]=alpha[frac[del^2T][delx^2]+frac[del^2T][dely^2]]`

In order to solve the above equation we can write the diiference equation as

`frac[T_p^(n+1)-T_p^n][trianglet]=alpha[frac[T_l-2T_p+T_r][trianglex^2]+frac[T_t-2T_p+T_b][triangley^2]]`

Since we are first going to solve the equation using the explicit scheme we need to evaluate Tp at the next time step. Re arraging the terms in the above equation we get

`T_p^(n+1)=T_p+k_1(T_l-2T_p+T_r)^n+k_2(T_t-2T_P+t_b)^n`

where k1 and k2 are constants `k_1=alphafrac[trianglet][trianglex^2] and k_2=alphafrac[trianglet][triangley^2]`

MATLAB CODE FOR EXPLICIT SCHEME

clear all;
close all;
clc;
   
%Input Parameters
L =1;
nx = 10;
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;
dt=0.01;
alpha = 0.1;
t=5;
   
%CFL number
CFL = alpha*dt/(dx^2);
  
%Boundary Conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T = 298*ones(nx,ny);

%Temperature at boundary
T(:,1)=T_left;
T(:,end)=T_right;
T(end,:)=T_bottom;
T(1,:)=T_top;
 
%Calculation of temperature distributon explicitly using iterative solvers
 k1 = (alpha*dt)/(dx^2);
 k2 = (alpha*dt)/(dy^2);
 
 %creating a copy of T
  [X,Y] = meshgrid(x,y);
  Told =T;
  nt=t/dt;
  % 2d transient unsteady state heat conduction equation in explicit method
  iterative_solver =1;
  tic;
  iterations =1;
  for v =1:nt
  if iterative_solver ==1
    
    error = 9e9;
    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;
    iterations = iterations+1;
  end
end
time=toc;
% Transcient heat conduction plot  
figure(1);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Explicit Scheme,Iteration number=%d,computation time=%fs,CFL=%f ',iterations,time,CFL);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
axis("ij");
pause(0.03);
end

We have to define the square plate of length 1m with 10 nodes so we define some variables in the start of the program.

• Length of the plate (l) = 1m
• number of nodes on x-axis (nx) = 10
• number of nodes on y-axis (ny) = 10
• thermal Diffusivity (`alpha`) = 1e-3 `m^2/s`
• X is the array of nodes in the x axis
• Y is the array of nodes in y axis.
• error term for temperature= 9e9
• tolerance for the temperature criteria = 1e-4

After defining the square plate we must assign the boundary conditions at which the temperatures at the boundaries is assigned as

• Left: 400K
• Top: 600K
• Right: 800K
• Bottom: 900K

The above condition should be applied to all the nodes in the outer boundary so the next 4 lines indicates that we are assigning the same boundary conditions for all the outer nodes. We the define the CFL number for the explicit scheme and then start the iterations. Now we are saying to enter the space and time loops using if command statement by defining the explicit scheme formula in the code. The "while" statement is also made use by defining the error and the tolerance to run the number of iterations till we get convergent solution. In the iterative process we keep updating the new values of temperature into the domain. The transient heat plot for explicit scheme is then being plotted using the contourf statement by adding suitable title and labels to the plot.

RESULT AND DISCUSSION

Explicit scheme plot for unsteady heat conduction can be viewed and analysed from this plot as we can observe that it took 1081 iterations to reach a stable solution. The CFL number for the above scheme shows a value of 0.00081 which is well less than 0.5 so it can be infered from the plot that the solution that was obtained from this code is in stable condition. It has taken approximately 253 seconds for the simulation to run and come to a stable condition.

 TWO DIMENSIONAL HEAT CONDUCTION (UNSTEADY/TRANSCIENT) USING IMPLICIT SCHEME

The equation for 2D unsteady state heat conduction is given by

 `frac[del T][delt]=alpha[frac[del^2T][delx^2]+frac[del^2T][dely^2]]`

The diffrerence equation for implicit scheme is 

`T_p^(n+1)=frac[1][1+2k_1+2k_2][T_p^n+k1(T_l+T_r)^(n+1)+k_2(T_t+T_b)^(n+1)]`

 where k1 and k2 are constants

MATLAB CODE FOR IMPLICIT SCHEME

clear all;
close all;
clc;
   
%inputs
L =1;
nx = 10
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;
dt=0.1;
alpha = 1;
t=10;
   
%CFL number
CFL = alpha*dt/(dx^2);
  
%Boundary Conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T =ones(nx,ny);
T(:,1)=T_left;
T(:,end)=T_right;
T(end,:)=T_bottom;
T(1,:)=T_top;
 
 
% definining k1 and k2 terms
 k1 = (alpha*dt)/(dx^2);
 k2 = (alpha*dt)/(dy^2);
 nt=t/dt;
 %creating a copy of T
  [X,Y] = meshgrid(x,y);
  Told =T;
  T_prev_dt=Told;

  %jacobi iterative solver
  nt=t/dt;
    iterative_solver =1;
    tic;
     jacobi_iterations =1;
    %if iterative_solver ==1
    term1=(1+2*k1+2*k2)^(-1);
    term2=k1*term1;
    term3=k2*term1;
  for v=1:nt  

    error=100;
    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)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
    end
  end
    error = max(max(abs(Told-T)));
    Told = T;
    jacobi_iterations = jacobi_iterations+1;
end
T_prev_dt=T;
 %end
% Plot for jacobi unsteady heat conduction
figure(1);
[P,Q]=contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
F = toc;
title_text=sprintf('Jacobi iterative solver,Iteration number=%d,computation time F=%fs,CFL=%f',jacobi_iterations,F,CFL);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
axis("ij");
pause(0.03);
end
%Gauss seidel method
L =1; 
nx = 10;
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;
dt=0.1;
t=10;
alpha = 1;

%CFL number
CFL = alpha*dt/(dx^2);
  
% Boundary Conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T = ones(nx,ny);
T(:,1)=T_left;
T(:,end)=T_right;
T(end,:)=T_bottom;
T(1,:)=T_top;
 

% definining k1 and k2 terms
 k1 = (alpha*dt)/(dx^2);
 k2 = (alpha*dt)/(dy^2);
 nt=t/dt;
 %creating a copy of T
  [X,Y] = meshgrid(x,y);
  Told =T;
  T_prev_dt=Told;
  iterative_solver =2;
    tic;
     gs_iterations =1;
    for w = 1:nt
  %if iterative_solver ==2
    term1=(1+2*k1+2*k2)^(-1);
  term2=k1*term1;
  term3=k2*term1;
   
    error=100;
    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)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
    end
  end
    error = max(max(abs(Told-T)));
    Told = T;
    gs_iterations = gs_iterations+1;
end
T_prev_dt=T;
%end
F = toc;
% Plot for Gauss seidel iterative solver
figure(2);
[P,Q] = contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('Gauss Seidel Iterative solver ,Iteration number=%d,computation time F=%fs,CFL=%f',gs_iterations,F,CFL);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
axis("ij");
pause(0.03);
end


%SOR method
iterative_solver =3;
L =1; 
nx = 10;
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;
dt=0.1;
t=10;
alpha = 1;
   
%CFL number
CFL = alpha*dt/(dx^2);
% calculation of relaxation factor
omega = 2/(1+sin(pi*dx));
  
%Boundary Conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom= 900;
T = ones(nx,ny);
T(:,1)=T_left;
T(:,end)=T_right;
T(end,:)=T_bottom;
T(1,:)=T_top;
 
% definining k1 and k2 terms
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
 
%creating a copy of T
[X,Y] = meshgrid(x,y);
Told =T;
nt=t/dt;
T_prev_dt=Told
tic;
SOR_iterations =1;
for z=1:nt
%if iterative_solver ==3
   term1=(1+2*k1+2*k2)^(-1);
   term2=k1*term1;
   term3=k2*term1;
    
    error=100;
    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_gs(i,j)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
         T(i,j) = (1-omega)*Told(i,j) + omega*(T_gs(i,j));
    end
  end
    error = max(max(abs(Told-T)));
    Told = T;
    SOR_iterations = SOR_iterations+1;
end
T_prev_dt=T;
end
F = toc;
%plot for sor iterative solver
figure(3);
[P,Q] = contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('SOR iterative solver ,Iteration number=%d,computation time F=%fs,CFL=%f',SOR_iterations,F,CFL);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
axis("ij");
pause(0.03);
end

We have to define the square plate of length 1m with 20 nodes so we define some variables in the start of the program.

• Length of the plate (l) = 1m
• number of nodes on x-axis (nx) = 20
• number of nodes on y-axis (ny) = 20
• thermal Diffusivity (`alpha`) = 1.88E-5 `m^2/s`
• X is the array of nodes in the x axis
• Y is the array of nodes in y axis.
• error term for temperature= 9e9
• tolerance for the temperature criteria = 1e-4

After defining the square plate we must assign the boundary conditions at which the temperatures at the boundaries is assigned as

• Left: 400K
• Top: 600K
• Right: 800K
• Bottom: 900K

The above condition should be applied to all the nodes in the outer boundary so the next 4 lines indicates that we are assigning the same boundary conditions for all the outer nodes. We the define the CFL number for the implicit scheme and then start the iterations. Now we are saying to enter the space and time loops using if command statement by defining the explicit scheme formula in the code. The "while" statement is also made use by defining the error and the tolerance to run the number of iterations till we get convergent solution. In the iterative process we keep updating the new values of temperature into the domain. The transient heat plot for implicit scheme is then being plotted using the contourf statement by adding suitable title and labels to the plot.

RESULTS AND DISCUSSIONS

In the above three figures we observe that the jacobi method takes 1078 iterations, gauss seidel takes 648 iterations and SOR method takes 310 iterations for the convergent solutions. The implicit method leads to a coupled linear system of equations. 

CONCLUSION

Thus we have evaluated the 2D heat conduction in a square plate using the iterative solvers and we see the schemes which is is used in evaluating the equations.