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

AIM:- Mid term project - Solving the steady and unsteady 2D heat conduction problem

OBJECTIVE:-

Steady-state analysis & Transient State Analysis

Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class.

The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

You will implement the following methods for solving implicit equations.

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Your absolute error criteria are 1e-4

Make sure to use the "INSERT CODE SAMPLE" option to paste the code here.

Post screenshots of the plots and use the projects.skill-lync.com page to Publish your findings.

THEORY:-

Numerical Solution schemes are different in being explicit or implicit.

When a  direct computation of the dependent variable can be made in terms of known quantities the computation is said to be explicit.

When the dependent variable is defined by sets of equations and either a matrix or iterative.

The technique is n added to obtain the solution to the numerical method is said to be implicit.

EXPLICIT METHOD:-

When the dependent variables are defined by coupled sets o equations, and either a matrix or iterative technique is needed to obtain the solution the numerical method is said to be explicit.

In an explicit numerical method, s would be evaluated in terms of known quantities at the previous time step n.

As a general rule, it can be shown that the condition cd<1 is very nearly equivalent to the stability condition for an explicit approximation.

Another general rule is that the time step sizes for explicit stability and accuracy are usually equivalent 

it is this increase damping with the increase in time step size which produces inaccuracies in transient behavior

IMPLICIT METHOD:-

when the dependent variables are defined by coupled sets of equation sets of equations and either a matrix or iterative technique is needed to obtain the solution the numerical method is said to be implicit.

in Computational fluid dynamics (CFD) the governing equations are nonlinear and the number of unknown variables is typically very large.

iteration is used to advance a solution through a sequence of steps from a starting state to a final converged state.

The principal reason for using implicit solution methods which are more complex to program and require more computational either in each solution step is to allow for large time step sizes.

A simple qualitative model will help to illustrate how this works.

In order to solve for dependent variables iterative solvers used are

1. jacobian iterative solver

2. Gauss seidel iterative solver

3. Successive over-relaxation iterative solver.

* jacobian iterative solver

in numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of numerical linear algebra.

Each diagonal element is solved for and an approximate value is plugged in

This process is then iterative until it converges 

The algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization.

*Gauss  seidel iterative method:- is a technique used to solve a linear equations

The method is named other the German mathematics curl Friedrid gauss and Philip laud was for seidel 

The method is similar to the Jacobi method and in the same, a strict or irreducible diagonal dominance of the system is sufficient to ensure coverage meaning the method will work.

 * Successive over-relaxation iterative solver.

In this type of iterative solver, the relaxation factor is used in increasing the coverage rate. it goes by the formula 

where`ω`>(`∂t/`∂x^2)+(∂t/`∂y^2) is the relaxation factor.

Solving 2D heat conduction equation for steady-state:-

Governing equation for the 2D heat conduction equation for steady-state is (∂t/`∂x^2)+(∂t/`∂y^2) =0

Now in order to solve the above governing equation, we first discretize the equation by using a central differencing scheme are as follows:-

Equation 1 will be solved by using the point iterative method.

The point iterative method which will be used to solve the above equation 1 is the above iteration method which are Jacobi, gauss Siedel, and the SOR method 

MATLAB CODE:-

. In the following program, we have to solve equation 1 by using different iterative methods.

Matlab code for steady-state 2D heat conduction profile :

clear all
close all
clc
%specifying length of domain and no of grid points
l=1;
nx=10;
ny=nx;
% Defining space domain
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
%creating tempeatures across given region
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,end) = 800;
t(end,:) = 900;
% Tempeature at edges
t(1,1) = (600+400)/2;
t(1,end)=(600+800)/2;
t(end,1)= (400+900)/2;
t(end,end) = (800+900)/2;

% updating old values
told = t;
% solving 2D steady state system
[xx,yy] = meshgrid(x,y);
% Assigning tolerence and error values 
tolerance=1e-4;
error = 1e8;
% Defining k factor
k = 2*((1/(dx^2))+(1/(dy^2)));
%iteration method input
iter_method = input('Enter the iteration method number= ')
% jacobi method
%Defining iteration varible to store and display no of iterations computated
 if iter_method == 1
tic 
jacobi_iter =1;
while(error>tolerance)
    for i = 2:nx-1
        for j  = 2:ny-1
            term1 = (told(i-1,j)+ told(i+1,j))/(k*(dx^2));
            term2 = (told(i,j-1)+told(i,j+1))/(k*(dy^2));
            t(i,j) = term1+term2;
        end
    end
    error = max(max(abs(told-t)));
    told = t;
    jacobi_iter = jacobi_iter+1;
end
time1=toc;
 end
%Gauss siedal method 
 if iter_method==2
    tic
    gs_iter=1;
    while(error>tolerance)
        for i=2:nx-1
            for j=2:ny-1
                % The brackets we have inserted was the problem. The
                % equation would be incorrect if we put the extra bracket
                % at the starting. I am sending the difference in equations because of the brackets. 
               % term1 = ((t(i-1,j)+ told(i+1,j))/k*(dx^2));
                %term2 = ((t(i,j-1)+ told(i,j+1))/k*(dy^2));
                term1 = (t(i-1,j)+ told(i+1,j))/(k*(dx^2));
                term2 = (t(i,j-1)+ told(i,j+1))/(k*(dy^2));
                t(i,j) = term1+term2;
            end
        end
        error = max(max(abs(told - t)));
        told = t;
        gs_iter=gs_iter+1;
    end
    time2 =toc;
 end
% SOR Method
if iter_method ==3
    omega =1.8;
    tic
    sor_iter=1;
    while(error>tolerance)
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (t(i-1,j)+told(i+1,j))/(k*(dx^2));
                term2 =(t(i,j-1)+told(i,j+1))/(k*(dy^2));
                t(i,j) = (1-omega)*told(i,j)+omega*(term1+term2);
            end
        end
        error = max(max(abs(told-t)));
        told=t;
        sor_iter=sor_iter+1;
    end
    time3=toc;
end
%plotting the results tempeature
[c,h] = contourf(xx,yy,t);
clabel(c,h);
colorbar;
set(gca,'YDir','reveres')
colormap(jet);
if iter_method ==1
title({'study state tempeature profile';['iteration Method = ',num2str(iter_method)];['No iterations=',num2str(jacobi_iter)];['Simulation time= ',num2str(time1),'S']})
elseif iter_method ==2
    title({'Steady state tempeature profile';['iteration Method=',num2str(iter_method)];['No of iterations = ',num2str(gs_iter)];['Simulation time= ',num2str(time2),'S']})
elseif iter_method ==3
        title({'Steady state tempeature profile';['iterations Method= ',num2str(iter_method)];['No of iterations = ',num2str(sor_iter)];['simulation time = ',num2str(time3),'S']})
end
    xlabel('x_length Domain')
    ylabel('y_length Domain')
    
        
    

RESULT:-

Plot for 2D heat conduction for steady-state using Jacobi iterative method:-

Steady-state temperature profile iteration Method = 1

No of iterations =  15940

simulation  time = 3.4096s

plot for 2D heat conduction for steady-state using gauss Siedel iterative method:

Steady-state temperature profile 

iteration Method =2

No of iterations = 2

Simulation time = 1.9161s

plot for 2D heat conduction for steady-state using SOR iterative method:-

Steady-state temperature profile 

iteration Method  = 3

Simulation time = 0.34395s 

Explicit schemes:-

After applying the first-order forward differencing scheme on L.H.S. and central differencing scheme on R.H.S. on equation 2 we get 

Matlab code for unsteady-state 2D heat conduction profile in explicit method:

clear all
close all
clc
l = 1
nx = 60;
ny = 60;

%Defining space domain
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);

%calculation of time step 
%dt = (0.5)/((alpha/dx^2)+(alpha/dy^2));
dt = 0.00001;
alpha = 1.4;
% creating tempeatures across given region
t = 300*ones(nx,ny);
t(:,1)= 400;
t(1,:) = 600;
t(:,end) = 800;
t(end,:) = 900;
% Tempeatures at edges
%t(1,1)= (600+400)/2;
%t(1,end) = (600+800)/2;
%t(end,1) = (400+900)/2;
%t(end,end)= (800+900)/2;
%updating at t values
told = t;
CFL_num =(alpha*dt)/((1/dx^2)+(1/dy^2));
% solving 2D unsteady state system
%[xx,yy] = meshgrid(x,y);
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

iter_method=1;
tic;
for k = 1:7000
    for i = 2:nx-1
        for j = 2:ny-1
            term1=(1-2*k1-2*k2);
            term2 = k1;
            term3 = k2;
            H = (told(i-1,j)+told(i+1,j));
            V = (told(i,j-1)+told(i,j+1));
            t(i,j) = (told(i,j)*term1)+(term2*H)+(term3*V);
            
        end
    end
    told = t;
end
time = toc;
figure(1)
colormap(jet)
contourf(x,y,t,'showtext','on');
colorbar
%clabel(c,h);
title({'unsteady tempeature profile';['simulation time = ',num2str(time),'s']})
xlabel('xlength Domain');
ylabel('ylength Domain');

            
        

RESULTS:-

 implicit scheme:-

governing equation for 2D heat conduction for the unsteady state is  

∂T/∂t=(∂t/`∂x^2)+(∂t/`∂y^2) =0..........(2)

After applying the first-order forward differencing scheme on L.H.S. and central differencing scheme on R.H.S. on equation 2D, we get.

Matlab code for unsteady-state 2D heat conduction profile for implicit method:-

clear all
close all
clc
nx = 10;
ny = 10;
x = linspace(0,1,nx);
dx = 1/(nx-1);
y = linspace(0,1,ny);
dy = 1/(ny -1);
alpha = 1.4;
%dt =(0.5)/((alpha/dx^2)+(alpha/dx^2));
dt = 0.001;
tolerance = 1e-4;
%error = 1e9;
t = 300*ones(nx,ny);
t(:,1) = 400;
t(1,:) = 600;
t(:,end) = 800;
t(end,:) = 900;
%Tempeature at edges
t(1,1) = (600+400)/2;
t(1,end) = (400+800)/2;
t(end,1) = (400+900)/2;
t(end,end) = (800+900)/2;
%upadting old t values
told = t;
tprev =told;
%solving 2D steady state system
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
cFL_num = k1+k2;
nt = 300; %number of time step

iter_method =input('Enter the iteration method number=');

if iter_method == 1
    %computing Time loop
    tic;
    jacobi_iter = 1;
    for k =1:nt 
        %jacobian iterative solve
        error =1e3;
        while(error>tolerance)
            for i = 2:nx-1
                for j=2:ny-1
                    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));
                   
                    t(i,j) = (tprev(i,j)*term1)+(term2*H)+(term3*V);
                end
            end
           % error = mx(max(abs(told - t)));
             error = max(max(abs(t - told)));
            told = t;
            jacobi_iter = jacobi_iter+1;
        end
        tprev = t;
    end
     time1 = toc;
     
%     figure(1)
%     [c,h] = contourf(x,y,t);
%     set(gca, 'YDir','reverse')
%     clabel(c,h);
%     colormap(jet);
end

%Gauss siedal iteration method
if iter_method == 2
    gs_iter =1;
    tic
    for nt = 1:600
        term1 = 1/(1+2*k1+ 2*k2);
        term2 = k1*term1;
        term3 =  k2*term1;
        error = 1e9;
        while(error>tolerance)
            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) = (told(i,j)*term1)+(term2*H)+(term3*V);
                end
            end
            error = max(max(abs(told-t)));
            told = t;
            gs_iter=gs_iter+1;
        end
        time2 =toc;
    end
end
%SOR iteration method
omega = 1.2;
if iter_method ==3
    sor_iter= 1;
    tic
    for nt = 1:600
        term1 =  1/(1+2*k1+2*k2);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while(error>tolerance)
            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)= (1-omega)*told(i,j)+omega*(term1*told(i,j)+term2*H+term3*V);
                end
            end
            error = max(max(abs(told - t)));
            told = t;
            sor_iter = sor_iter+1;
        end
        time3 = toc;
    end
end


figure(1)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar;
colormap(jet);
if iter_method ==1
    title({'unsteady state tempeature profile';['iteration Method =jacobi_iter'];['No of iterations= ',num2str('jacobi_iter')]})
elseif iter_method ==2
    title({'unsteady state tempeature profile';['iterations method = gauss siedal']; ['No of iterations = ',num2str(gs_iter)]})
elseif iter_method ==3
        title({'unsteady state tempeature profile';['iterations Method = SOR'];['No of iterations = ',num2str(sor_iter)]})
end
                   

Result-

plot for unsteady-state 2D heat conduction in the implicit method for jacobian iteration:-

plot for unsyeady state 2D heat conduction in implicit method for gauss siedel iteratons:

plot for unsteady state 2D heat conduction in implicit method for SOR iterations:-

Conclusion and Explanation:-

Though we are getting the same result for different iterations and computation time taken by the iterative solver to get convergence result.

By observing the above plots it is clear that Jacobi iterations such as gauss Siedel and SOR.

The computation time is also less in jacobian iterations than in the other two iterative methods.

The implicit method requires guess values iteration for  convergence as it also based on the the future state values

Which are not known at the time f computations

it checks convergence for each time step.

This caused a delay in the simulation.

The implicit method can be used for problems when converged results are needed in less computation time as they are unconditionally stable.

with an increase in time size, the convergence rate is increased which leads to a steady state with less no of iterations.

The explicit method results in a time-accurate solution as it is governed by stability in criteria in the solution and when is unstable it requires more no iterations