Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

Steady and Unsteady State 2D Heat Conduction

We will write a code for the steady and unsteay state 2D heat conduction 

First consider the steady state heat conduction using implicit method :

The code for the steady state implicit method is given below in the form of a function

Here in this function we have three inputs the first is the length of the plate which we have assumed to be square, the second input is the number of nodes in x direction also the number of nodes in the y direction will be the same and our third input will be the iteration method by which we want to solve our problem. Here we have iterative three methods so we will use number 1 , 2 and 3 to give as an input for jacobian iteration method, Gauss_seidel iteration method and Successive over relaxation method respectively.

function [out] = steady_state_2d_heat_conduction_implicit_methods(L,nx,iterative_solver)
%this function has three inputs
% 1. length of the square plate
% 2. numbers of nodes inx direction ( the number of nodes in x and y directions is assumed to be equal)
% 3. the iterative solver we have given numbers to the solver 1-jacobi, 2-gauss seidel, 3-successive over relaxation

% inputs
% L = 1 ;
% nx = 10;
ny = nx; % number of nodes in the x and y direction is same
dx = L/nx;
dy = L/ny;

% creating 1D srrays for x and y directions
x = linspace(0,L,nx);
y = linspace(0,L,ny);

% creating a 2D grid
[xx yy] = meshgrid(x,y);

%assigning temperature over the grid
T = 300*ones(nx,ny);

%assigning temperature at the edges
T(1,1:end) = 600;
T(end,1:end) = 900;
T(1:end,1) = 400;
T(1:end,end) =800;

% assigning temperature at the corners
T(1,1) = (T(1,1+1)+T(1+1,1))/2;
T(1,end) = (T(1,end-1)+T(1+1,end))/2;
T(end,1) = (T(end-1,1)+T(end,1+1))/2;
T(end,end) = (T(end-1,end)+T(end,end-1))/2;

% creating a copy of temp
Told = T;
error = 9e9;
tol = 1e-4;
% iterative_solver = 1;

%jacobi method
if iterative_solver == 1 % if the value of the iterative solver is 1 than the function will use this while loop
    jacobi_iter = 1; % initiating the number of jacobi iterations
    while(error > tol) % to check whether the error is less than the tolerence
                
        k = 2*((dx^2+dy^2)/(dx^2*dy^2)); 
        
        for i = 2:nx-1 % saptial loop in x diration
            for j = 2:ny-1 % saptial loop in y direction
                
                % equation for the calculation of temperature for the jacobi method
                T(i,j) = 1/k*((Told(i-1,j)+Told(i+1,j))/dx^2) + 1/k*((Told(i,j-1)+Told(i,j+1))/dy^2);
                
            end
        end
        error = max(max(abs(Told - T))); % finding the maximum value of error
        Told = T; % updating the temperature values for the next iteration
        jacobi_iter = jacobi_iter + 1; % to calculate the total number of iterations
        %plotting the results in the form of contours 
        figure(1)
        contourf(x,y,T,'showtext','on')
        colorbar % shows the colorbar
        colormap(jet)  
        set(gca,'YDIR','reverse') % to see the proper allignment of the square plate
        title_text= sprintf('countour of temperature across a square plate at %d iteration',jacobi_iter);
        title(title_text)
        xlabel('X')
        ylabel('Y')
    end
end



% Gauss-seidel method
if iterative_solver == 2 % if the value of the iterative solver is 2 than teh function will use this loop 
    gs_iter = 1; % initiating the number of gauss seidel iterations
    while (error > tol) % to check whether the error is less than the tolerence
        % calculcating the value of k 
        k = 2*((dx^2+dy^2)/(dx^2*dy^2));
        
        for i = 2:nx-1 % spatial loop in x direction
            for j = 2:ny-1 % spatial loop in Y direction
                % equation to calculate the value of temperature for the gauss-seidel method
                T(i,j) = 1/k*((T(i-1,j)+Told(i+1,j))/dx^2) + 1/k*((T(i,j-1)+Told(i,j+1))/dy^2);
                
            end
        end
        error = max(max(abs(Told - T))); % finding the maximum value of error
        Told = T; % updating the values of temerature for teh next iterations
        gs_iter = gs_iter + 1; % to calculate the total number of iterations
        %plotting results in the form of contours
        figure(1) 
        contourf(x,y,T,'showtext','on')
        colorbar % to show the colorbar 
        colormap(jet) 
        set(gca,'YDIR','reverse') % to see the square plate in proper allignment
        title_text= sprintf('countour of temperature across a square plate at %d iteration',gs_iter);
        title(title_text)
        xlabel('X')
        ylabel('Y')
    end
end



% successive over relaxation method
if iterative_solver == 3 % if the value of iterative_solver is 3 than the function will use this loop
    sor_iter = 1; % initiating the number of SOR iterations
    while (error > tol) % to check whether the error is less than the tolerence
        % calculating the value of k
        k = 2*((dx^2+dy^2)/(dx^2*dy^2));
        omega = 1.5; % omega is the relaxation factor 
        
        for i = 2:nx-1 % spatial loop in X direction
            for j = 2:ny-1 % spatial loop in y direction
                % equation to calculate the value of temperature with SOR method
                T(i,j) = (1- omega)*Told(i,j) + omega*(1/k*((T(i-1,j)+T(i+1,j))/dx^2) + 1/k*((T(i,j-1)+T(i,j+1))/dy^2));
                
            end
        end
        error = max(max(abs(Told - T))); % to find the maximum value of error
        Told = T; % updating the values of temperature for the next loop
        sor_iter = sor_iter + 1; % to calculate teh total number of iterations
        % plotting the output in the form of contours
        figure(1)
        contourf(x,y,T,'showtext','on')
        colorbar % to show the colorbar
        colormap(jet)
        set(gca,'YDIR','reverse')% to see the square plate in proper allignment
        title_text= sprintf('countour of temperature across a square plate at %d iteration',sor_iter);
        title(title_text)
        xlabel('X')
        ylabel('Y')
    end
end

Now to run this code we will write the function name in our command window along with the inputs. 

let us first consider the inputs as 1m length , 10 nodes in x direction and Jacobi iteration method for that our input command will be "steady_state_2d_heat_conduction_implicit_methods(1,10,1)" 

So after enterring our command in the command window as shown

>> steady_state_2d_heat_conduction_implicit_methods(1,10,1)

we get the contour as the output as shown below

from the image we can see the contour of temperature for a 2D surface with length 1m.

also we can see that it takes 208 iteration to reach the steady state using jacobi method

Now to implement the Gauss-Seidel iterative method

The input command for the same length and number of nodes will be "steady_state_2d_heat_conduction_implicit_methods(1,10,2)"

So after entering our command in the command window as shown below

>> steady_state_2d_heat_conduction_implicit_methods(1,10,2)

We get the following contour for temperature as shown below

we can see the steady state temperature contour and also here it take 112 iteration to reach the steady state. so we can say the gauss seidel method is faster than jacobi method and gives results in less number of iterations.

Now we will implement the Successive over-relaxation iterative method 

for this method our input command will be as shown below 

>> steady_state_2d_heat_conduction_implicit_methods(1,10,3)

we get the following contour for the temperature over the 2D surface

here we can see the steady state for the 2D heat conduction using the successive over relaxation iterative method. the relaxation used here is 1.5 and as we can see we reach the steady state in very less number of iterations, here it takes only 29 iterations. It is the fastest method but only when the relaxation factor has a specific value. The number of iteration may increase or decrease depending on the relaxation factor.

Now we will write a code for the Unsteady State:

The code for the Unsteady state using explicit method is given below

we have used the thermal conductivity of the material as 1.4

% transient state 2d heat conduction using explicit method 

clear all
close all
clc

% inputs
L = 1; % length of square plate
nx = 10; % number of nodes in x direction
ny = nx; % number of nodes in Y direction
dx = L/nx; % size of element in x direction
dy = L/ny; % size of element in y direction
% creating arrays in x and y
x = linspace(0,L,nx);
y = linspace(0,L,ny);

%time step
dt = 1e-4;

% creating a 2D grid
[xx yy] = meshgrid(x,y);

%assigning initial condition of temperature over the grid
T = 300*ones(nx,ny); 

%assigning initial temperature at the edges
T(1,1:end) = 600;
T(end,1:end) = 900;
T(1:end,1) = 400;
T(1:end,end) =800;

% assigning initial temperature at the corners
T(1,1) = (T(1,1+1)+T(1+1,1))/2;
T(1,end) = (T(1,end-1)+T(1+1,end))/2;
T(end,1) = (T(end-1,1)+T(end,1+1))/2;
T(end,end) = (T(end-1,end)+T(end,end-1))/2;

% creating a copy of temp
Told = T;

alpha = 1.4; % thermal conductivity

nt = 1500; % number of time steps

for k = 1:nt % time loop
    
    for i = 2:nx-1 % spatial loop in x direction
        for j =2:ny-1 % spatial loop in y direction
            % equation for the calculation of temperature at the new time step
            T(i,j) = Told(i,j) + alpha*dt*((Told(i-1,j)-2*Told(i,j)+Told(i+1,j))/dx^2 + (Told(i,j-1)-2*Told(i,j)+Told(i,j+1))/dy^2 );
            
        end
    end
    Told = T; % updating the value of temerature(Told) for the new time loop
    % plotting the result in the form of contours
    figure(1)
    contourf(x,y,T,'showtext','on')
    %axis([0 1 0 1])
    colorbar
    colormap(jet)
    set(gca,'YDIR','reverse')
    title_text= sprintf('countour of temperature across a square plate at %d time step',k);
    title(title_text)
    xlabel('X')
    ylabel('Y')
    
    
    
end

After running the code we can see the contour of temperature for the 2D surface at every time step here i have attached the screen shot of the last time step i.e at time step = 1500 , we can see that almost steady state is reached

Now we will write the code for unsteady state using the implicit method.

The code the unsteady state implicit method is given below, this code is in the form of function 

Here again for this function we have there inputs which are the same as we used in the steady state. Just the name of the function is different. 

function [out] = transient_implicit_2d_heat_conduction(L,nx,iterative_solver)
    %this function has three inputs
    % 1. length of the square plate
    % 2. numbers of nodes inx direction ( the number of nodes in x and y directions is assumed to be equal)
    % 3. the iterative solver we have given numbers to the solver 1-jacobi, 2-gauss seidel, 3-successive over relaxation

    %inputs
    % L = 1; %length = 1m
    % nx = 10; % number of nodes along x
    ny = nx; % number of nodes along y
    dx = L/nx; % element size along x
    dy = L/ny; % element size along y

    % creating arrays in x and y
    x = linspace(0,L,nx);
    y = linspace(0,L,ny);

    %time step
    dt = 1e-4;

    % creating a 2D grid
    [xx yy] = meshgrid(x,y);

    %assigning initial condition of temperature over the grid
    T = 300*ones(nx,ny); 

    %assigning initial temperature at the edges
    T(1,1:end) = 600;
    T(end,1:end) = 900;
    T(1:end,1) = 400;
    T(1:end,end) =800;

    % assigning initial temperature at the corners
    T(1,1) = (T(1,1+1)+T(1+1,1))/2;
    T(1,end) = (T(1,end-1)+T(1+1,end))/2;
    T(end,1) = (T(end-1,1)+T(end,1+1))/2;
    T(end,end) = (T(end-1,end)+T(end,end-1))/2;

    % creating a copy of temp
    Told = T; % for the iterations
    T_prev_dt = T % for time loop

    % iterative_solver = 1
    tol = 1e-4;
    alpha = 1.4;
    
    nt = 1500; % number of time steps
    
    % jacobi method
    jacobi_iter = 1; % initiation number of iterations for jacobi method  
   
    if iterative_solver == 1 % if the value of iterative_solver is 1 than the function will use this command and it will follow jacobi mehtod
            
        for k = 1:nt % time loop

            error = 9e9; % assigning value of error
            while (error > tol) % to check whether the error is less than tolerence

                for  i = 2:nx-1 % saptial loop in X direcction
                    for j = 2:ny-1 % spatial loop in Y direction
                        % braking the large equation into small terms
                        k1 = alpha*dt/dx^2;
                        k2 = alpha*dt/dy^2;

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

                        H = Told(i-1,j) + Told(i+1,j);
                        V = Told(i,j-1) + Told(i,j+1);
                        % equation for the temperature using jacobi method
                        T(i,j) = term1*T_prev_dt(i,j) + term2*H + term3*V;

                    end
                end

                error = max(max(abs(Told - T))); % finding the maximum value of error
                Told = T; % updating the temperature of Told for teh next iteration
                jacobi_iter = jacobi_iter + 1 % to find the total number of iteration 
   
            end
            T_prev_dt = T; % updating the value of T_prev_dt for the next time step

            
        end
        % plotting the output in the form of contours
        figure(1)
        contourf(x,y,T,'showtext','on')
        colorbar
        colormap(jet)
        title_text = sprintf('contour at time step = %d',k);
        title(title_text)
        set(gca,'YDIR','reverse')
        
    end
    
    %gauss seidel method
    gs_iter = 1; % initiating the number of iteration for Gauss Seidel method
    if iterative_solver == 2 % if the value of iterative_solver the functio will use this loop and will follow gauss seidel iterative method
            
        for k = 1:nt % time loop

            error = 9e9; % assigning value of error
            while (error > tol) % to check whether the error is less than the tolerence 

                for  i = 2:nx-1 % spatial loop in X direction
                    for j = 2:ny-1 % saptial loop in Y direction
                        
                        % breaking the large equation into small terms
                        k1 = alpha*dt/dx^2;
                        k2 = alpha*dt/dy^2;

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

                        H = T(i-1,j) + Told(i+1,j);
                        V = Told(i,j-1) + T(i,j+1);
                        % equation to find the temperature for the gauss-seidel method
                        T(i,j) = term1*T_prev_dt(i,j) + term2*H + term3*V;

                    end
                end
                error = max(max(abs(Told - T))); % finding the maximum value in error
                Told = T; % updating the value of Told for the next iteration 
                gs_iter = gs_iter + 1 % to find the total number of iterations
   
            end
            T_prev_dt = T; % updating the value of T_prev_dt for the next time step

               

        end 
        % plotting the output in the form of contours
        figure(2)
        contourf(x,y,T,'showtext','on')
        colorbar
        colormap(jet)
        title_text = sprintf('contour at time step = %d',k);
        title(title_text)
        set(gca,'YDIR','reverse')
        
    end
    
    % SOR method
    sor_iter = 1; % initiating the number of iteration for SOR method
    if iterative_solver == 3 % if the value of iterative_solver is 3 than the function will use this loop and will use SOR method
        omega = 1.02; % assigning the value to the relaxation factor 
        for k = 1:nt % time loop

            error = 9e9; % assigning value to the error
            while (error > tol) % convergence loop to check if the error id lesss than the tolerence

                for  i = 2:nx-1 % saptial loop in X direction
                    for j = 2:ny-1 % spatial loop in Y direction
                        % breaking the large equation into small terms
                        k1 = alpha*dt/dx^2;
                        k2 = alpha*dt/dy^2;

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

                        H = Told(i-1,j) + Told(i+1,j);
                        V = Told(i,j-1) + Told(i,j+1);
                        % equation for the temperature using the SOR method
                        T(i,j) = (1 - omega)*Told(i,j) + omega* (term1*T_prev_dt(i,j) + term2*H + term3*V);
                        
                    end
                end
                error = max(max(abs(Told - T)));% finding the maximum in the error
                Told = T; % updating the value of Told for the next iteration
                sor_iter = sor_iter + 1 % to calculate the total number of iterations
   
            end
            T_prev_dt = T;% updating the value of T_prev_dt for the next time step
            
        end
        % plotting the output in the form of contour
        figure(3)
        contourf(x,y,T,'showtext','on')
        colorbar
        colormap(jet)
        title_text = sprintf('contour at time step = %d',k);
        title(title_text)
        set(gca,'YDIR','reverse')
        
    end
end

 To implement jacobi iterative method we will use 1 as the input, so our input command will be  "transient_implicit_2d_heat_conduction(1,10,1)". We will give this command in our command window

>> transient_implicit_2d_heat_conduction(1,10,1)

the output we get is a contour as shown below

It takes 5887 iteration for the jacobi method , we have used 1500 time steps in our code so that a steady state is reached.

For gauss seidel iterative method we will use input command as "transient_implicit_2d_heat_conduction(1,10,2)"

>> transient_implicit_2d_heat_conduction(1,10,2)

 we get the temperature contour as given below 

It takes 5561 iteration for the gauss seidel iterative method for 1500 time steps.

we can see the contour is similar to that of the steady state

For Successive over relaxation method our input will be  "transient_implicit_2d_heat_conduction(1,10,3)"

for the sor method the relaxation factor used here is 1.02 

>> transient_implicit_2d_heat_conduction(1,10,3)

the output we get here is shown below

this is the contour for temperature using successive over relaxation method. here it takes about 5368 iterations if use different value of relaxation factor denoted by omega in the code it shows more number of iterations. we can see that after 1500 iterations a steady state had reached.