Solving the steady and unsteady 2D heat conduction problem.

Mid term project

Solving the steady and unsteady 2D heat conduction problem.

Given Data:

• Assume that the domain is a unit square i.e. domail length = 1m.
• Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction].

• Boundary conditions for steady and transient case:

  1. Left:400K

  2. Top:600K

  3. Right:800K

  4. Bottom:900K

• Absolute error criteria is 1e-4.

To Find:

• Steady state analysis & Transient State Analysis.
• Solve the 2D heat conduction equation by using the point iterative techniques i.e.

  1. Jacobi

  2. Gauss-seidel

  3. Successive over-relaxation

• To write code for both implicit and explicit schemes (This is for the Transient Analysis only)
• To write a detailed explanation for your results, plots or codes.

SOLUTION:-

STEADY STATE ANALYSIS:-

The 2- Dimensional equationis given by,

On Decretizing the above equation, we get;

Decretizing w.r.t. 'x' -

Decretizing w.r.t. 'y' -

As it is given to assume that nx=ny i.e. Δ">ΔΔx = Δ">ΔΔy.

By further solving & rearranging the term we get,

Program to solve 2D heat conduction equation by using the point iterative techniques i.e. by Jacobi, Gauss-Sedial & Successive Over Relaxation Methods:

% Solving Steady state 2D heat condition equation.
clear all
close all
clc
%JACOBI METHOD.
%Inputs.
L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
error = 9e9;
tol = 1e-4;   % absolute error value.
%Boundary Conditions.
T_L = 400;    %Temperature in 'K'.
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at the matrix corners.
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;
k1 = (dx^2)/(2*(dx^2 + dy^2));
k2 = (dy^2)/(2*(dx^2 + dy^2));
[X,Y] = meshgrid(x,y);
Told = T;     
% Jacobi Iteration.
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) = k1*(Told(i-1,j) + Told(i+1,j)) + k2*(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
 toc;
 time_required = toc;
%Plotting 
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
hold on
title_text_1 = sprintf('Steady State - JACOBI METHOD ');
title_text_2 = sprintf('iteration number=%d', jacobi_iter);
title({title_text_1,title_text_2});
end
%GAUSS-SEIDAL METHOD.

%Inputs.
L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
error = 9e9;
tol = 1e-4;   % absolute error value.
%Boundary Conditions.
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at the matrix corners.
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;
k1 = (dx^2)/(2*(dx^2 + dy^2));
k2 = (dy^2)/(2*(dx^2 + dy^2));
[X,Y] = meshgrid(x,y);
Told = T;
% Gauss-seidal Iteration.
for iterative_solver = 2;
tic;
 if iterative_solver == 2
    gs_iter = 2;
    while(error > tol) 
        for i = 2:nx-1 
            for j = 2:ny-1
                T(i,j) = k1*(T(i,j-1) + Told(i,j+1)) + k2*(T(i-1,j) + Told(i+1,j));
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        gs_iter =   gs_iter + 1;
    end
 end
 toc;
 time_required = toc;
%Plotting 
figure(2)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
hold on
title_text_1 = sprintf('Steady State - GAUSS-SEIDAL METHOD');
title_text_2 = sprintf('iteration number=%d', gs_iter);
title({title_text_1,title_text_2});
end
%SUCCESSIVE OVER RELAXATION METHOD.

%inputs
L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
error = 9e9;
tol = 1e-4;   % absolute error value.
alpha = 1.3 
omega = 1.2
%Boundary Conditions.
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at the matrix corners.
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;
k1 = (dx^2)/(2*(dx^2 + dy^2));
k2 = (dy^2)/(2*(dx^2 + dy^2));
[X,Y] = meshgrid(x,y);
Told = T;
% SOR Iteration.
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
                 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_iter = sor_iter + 1;
    end
 end
 toc;
 time_required = toc;
%Plotting 
figure(3)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('Steady State - SOR METHOD');
title_text_2 = sprintf('iteration number=%d', sor_iter);
title({title_text_1,title_text_2});  
end

EXPLANATION -

>> Initially declair all the values that are necessary for the programe. i.e. mention the inputs required for solving the equations. Here consider error= 9e9.

>> After defining all the inputs mention the boundary conditions which are give in the question.

>> Calculate the temperatures at all the boundaries i.e at left, right, top & bottom of a matrix respectively.

>> Now by using the above calculated temperatures at all the sides calculate the temperatures at corner points or the boundary points of a matrix. This is bascially done by doing average of two corresponding boundaries.

>> Then K1 & k2 is also declaired as calculated.

>> Initiallise a command for meshgrid.

>> With reference to the given coding stages , the code is written in a correct sequence of order to end up with a correct solution.

>> After the end of the nodal loop it is mandatory that we should update the old values.

>> Similarly after the end of convergence loop , update the previous values of timestep.

>> Now after checking if all the loop are closed with 'end ' properly , then write to plot the results with necessary labels, axis notations etc.

>> During plotting, ' contourf(X,Y,T) ' command is used to specify the x & y coordinates for the values in T.

>> ' colorbar ' - this is used to display a vertical colorbar to the right of the current axes.

>> ' clabel (a,b) ' - this is used to label the current contour plot with rotated text inserted into each contour line.

>> Write a line which help to show the text on the plots like iteration number, simulation time and title.

>> Similarly follow the same steps for the next two methods by changing the iterative values i.e iterative solver ==2 for Gauss -Sedial method & iterativesolver==3 for SOR method. For SOR method alpha is initiallised as relaxation factor. Here omega = 1.2.

RESULTS:

Plot showing 2D heat conduction equation by using the point iterative techniques i.e. by Jacobi, Gauss-Sedial & Successive Over Relaxation Methods:-

JACOBI METHOD -

>> Iteration number = 869

GAUSS-SEIDAL METHOD -

>> Iteration number = 464

SUCCESSIVE OVER RELAXATION METHOD -

>> Iteration number = 318

OUTCOMES :

Jacobi - Iteration number = 869

Gauss-sedial - Iteration number = 464

SOR - Iteration number = 318

>> From the above three plots it is clearly observed that for successive over relaxation method we get the faster results with the lower iteration number when compared to Gauss-Sedial method.

>> Here for SOR solver we use relaxation factor/scaling factor as alpha, there are two perameters for the value of 'alpha' i.e. if omega <1. It is under relaxed. & omega >1. It is over relaxed.

>> In a steady state equation , SOR method it is very effecient because it takes very less convergence with a very less number of iterations to give a steady state solution.

TRANSIENT STATE ANALYSIS:-

This method can be solved in two ways i.e. #Explicit Method.

#Implicit Method.

#Explicit Method :

In this method we take information from all the surrounding temperature values and we use those values to compute the new values in the time step.

Program to solve 2D heat conduction equation by using the Explicit Method.

             % Transient Analysis (Explicit schemes).
clear all
close all
clc
%inputs
L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
dt = 30;   % time step
error = 9e9;
tol = 1e-4;   % absolute error value.
alpha =2.3e-5
%Boundary Conditions.
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at matrix corners.
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;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
% JACOBI METHOD 
iterative_solver = 1;
tic;
if iterative_solver == 1
    iteration = 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;
        iteration =   iteration + 1;
toc;
time_taken = toc;
    end        
end
%Plotting 
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('EXPLICIT METHOD')
title_text_2 = sprintf('iteration number=%d', iteration);
title({title_text_1,title_text_2});

EXPLANATION -

>> Initially declair all the values that are necessary for the programe. i.e. mention the inputs required for solving the equations. Here consider error= 9e9. alpha is thermal diffusity.

>> After defining all the inputs mention the boundary conditions which are give in the question.

>> Calculate the temperatures at all the boundaries i.e at left, right, top & bottom of a matrix respectively.

>> Now by using the above calculated temperatures at all the sides calculate the temperatures at corner points or the boundary points of a matrix. This is bascially done by doing average of two corresponding boundaries.

>> Then K1 & k2 is also declaired as calculated.

>> " [X,Y] = meshgrid(x,y); " used to create a creating meshgrid.

>> After the end of the nodal loop it is mandatory that we should update the old values.

>> Similarly after the end of convergence loop , update the previous values of timestep.

>> Now after checking if all the loop are closed with 'end ' properly , then write to plot the results with necessary labels, axis notations etc.

>> During plotting, ' contourf(X,Y,T) ' command is used to specify the x & y coordinates for the values in T.

>> ' colorbar ' - this is used to display a vertical colorbar to the right of the current axes.

>> ' clabel (a,b) ' - this is used to label the current contour plot with rotated text inserted into each contour line.

>> Write a line which help to show the text on the plots like iteration number, simulation time and title.

RESULT:

Plot showing 2D heat conduction equation by using the Explicit scheme:

JACOBI -

>> Iteration number = 872

#Implicit Method:

Program to solve 2D heat conduction equation Implicitily by using Jacobi, Gauss-Sedal, SOR methods.

% Transient Analysis ( implicit schemes). 
clear all
close all
clc
%inputs
L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
dt = 30;   % time step
error = 9e9;
tol = 1e-4;   % absolute error value.
alpha =1.3 % thermal diffusity.
%Boundary Conditions.
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at the matrix corners.
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;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
iterative_solver = 1;
% Jacobi 
term1 = (1+ 2*k1 + 2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
            end    
        end
        error = max(max(abs(Told - T)));
        Told = T;
        jacobi_iter =   jacobi_iter + 1;
    end
    T_prev_dt = T;
 end 
toc;
time_required = toc;
%Plotting 
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf(' IMPLICT - JACOBI METHOD');
title_text_2 = sprintf('iteration number=%d', jacobi_iter);
title({title_text_1,title_text_2});  

%GAUSS-SEDIAL METHOD
%inputs

L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
dt = 30;   % time step
error = 9e9;
tol = 1e-4;   % absolute error value.
alpha =1.3
%Boundary Conditions.
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at the matrix corners.
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;
% calculating of temperature distribution using iterative solvers.
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
% Gauss-sidal 
[X,Y] = meshgrid(x,y);
 Told = T;
 T_prev_dt = Told;
term1 = (1+ 2*k1 + 2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
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)+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);
            end    
        end
        error = max(max(abs(Told - T)));
        Told = T;
       gs_iter = gs_iter + 1;
    end
    T_prev_dt = Told;
 end 
toc;
time_required = toc;
%Plotting 
figure(2)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT - GAUSS SEDIAL METHOD');
title_text_2 = sprintf('iteration number=%d', gs_iter);
title({title_text_1,title_text_2});  

% SOR METHOD
%inputs
L = 1;   % domain length- assumiong domain is a unit square.
nx = 20;
ny = 20;  % given that nx = ny;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);  % dx = L/(n-1).
dy = y(2) - y(1);  % dy = L/(n-1).
dt = 30;   % time step
error = 9e9;
tol = 1e-4;   % absolute error value.
omega =1.2;
alpha = 1.3;
%Boundary Conditions.
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900
T=ones(nx,ny);
% calculating the temperature at boundary.
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 average temperatures at the matrix corners.
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;
% calculating of temperature distribution using iterative solvers.
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
% SOR 
[X,Y] = meshgrid(x,y);
 Told = T;
 T_prev_dt = Told;
term1 = (1+ 2*k1 + 2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
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 = (omega*T(i-1,j) + (1-omega)*Told(i-1,j) + Told(i+1,j));
               V = (omega*T(i,j-1) + (1-omega)* 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;
       sor_iter = sor_iter + 1;
    end
    T_prev_dt = Told;
 end 
toc;
time_required = toc;
%Plotting 
figure(3)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT - SOR METHOD');
title_text_2 = sprintf('iteration number=%d', sor_iter);
title({title_text_1,title_text_2});  

EXPLANATION -

>> Initially declair all the values that are necessary for the programe. i.e. mention the inputs required for solving the equations. Here consider error= 9e9. alpha is thermal diffusity.

>> After defining all the inputs mention the boundary conditions which are give in the question.

>> Calculate the temperatures at all the boundaries i.e at left, right, top & bottom of a matrix respectively.

>> Here in SOR method alpha is the relaxation factor which is assumed to be alpha = 1.3.