Skill-Lync Launch Pad – Your Gateway to a Core Engineering Job by 2025! Only 1̶0̶0̶ 13 Seats Available.

00D 10H 32M 13S

Executive Programs

Workshops

Projects

Blogs

Careers

Placements

Student Reviews

For Business

Academic Training

Informative Articles

Find Jobs

We are Hiring!

All Courses

Choose a category

Mechanical

Electrical

Civil

Computer Science

Electronics

Offline Program

All Courses

CHOOSE A CATEGORY

Mechanical

Electrical

Civil

Computer Science

Electronics

Offline Program

Top Job Leading Courses

Automotive

CFD

FEA

Design

MBD

Med Tech

Courses by Software

Design

Solver

Automation

Vehicle Dynamics

CFD Solver

Preprocessor

Courses by Semester

First Year

Second Year

Third Year

Fourth Year

Courses by Domain

Automotive

CFD

Design

FEA

Tool-focused Courses

Design

Solver

Automation

Preprocessor

CFD Solver

Vehicle Dynamics

Machine learning

Machine Learning and AI

POPULAR COURSES

Post Graduate Program in Hybrid Electric Vehicle Design and Analysis

Post Graduate Program in Computational Fluid Dynamics

Post Graduate Program in CAD

Post Graduate Program in CAE

Post Graduate Program in Manufacturing Design

Post Graduate Program in Computational Design and Pre-processing

Post Graduate Program in Complete Passenger Car Design & Product Development

Executive Programs

Workshops

For Business

Success Stories

Placements

Student Reviews

Projects

Blogs

Academic Training

Find Jobs

Informative Articles

We're Hiring!

+91 9342691281 Log in

MATLAB - Mid term project - Solving the steady and unsteady 2D heat conduction

Problem Statement: - For A 2D surface maintained at a certain temperature at its edges. Solve using Iteration Methods(Jacobi, Gauss-Seidel, and S.O.R. using Gauss-Seidel) for Temperature distribution due to Heat Conduction that can be observed for the following States: - 1. Transient State until Steady State is Achieved…

Aditya Aanand
updated on 09 Sep 2022

Problem Statement: -

For A 2D surface maintained at a certain temperature at its edges. Solve using Iteration Methods(Jacobi, Gauss-Seidel, and S.O.R. using Gauss-Seidel) for Temperature distribution due to Heat Conduction that can be observed for the following States: -

1. Transient State until Steady State is Achieved

2. Steady State

The Given Temperature Maintained at the edges are as follows

$T_{Top} = 600 K, T_{Bottom} = 900 K, T_{Left} = 400 K, and T_{Right} = 800 K$ $T_("Top") = 600 K, T_("Bottom") = 900 K, T_("Left") = 400K, and T_("Right") = 800K$

Solution: -

1. AIM: -

A. To Solve 2D heat Conduction Equation for Steady and Transient State using Iterative Solver( Jacobi, Gauss-Seidel and Successive Over-Relaxation)

2. Given: -

A. $T_{top} = 600 K$ $T_("top") = 600 K$

B. $T_{Bottom} = 900 K$ $T_("Bottom") = 900 K$

C. $T_{Left} = 400 K$ $T_("Left") = 400 K$

D. $T_{Right} = 800 K$ $T_("Right") = 800K$

3. Assumed: -

A. $2D Transient Heat Conduction \to \frac{\partial T}{\partial t} - α (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}}) = 0$ $"2D Transient Heat Conduction" rarr (delT)/(delt) - alpha((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0$

B. $2D Steady Heat Conduction \to \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} = 0$ $"2D Steady Heat Conduction" rarr (del^2T)/(delx^2) + (del^2T)/(dely^2) = 0$

C. Dimensions of Surface $\Rightarrow length = 1 m and Width = 1 m$ $=> "length" = 1 m and "Width" = 1 m$

D. Number of Nodes in each Direction = 51

E. Using the Following notation for $T_{i, j}^{n}$ $T_(i,j)^n$

where $i, j \to x and y positions of Node$ $i,j rarr x and y " positions of Node"$

$n \to time Step iteration$ $n rarr "time Step iteration"$

4. Deriving Numerical Methods: -

A. 2D Steady Heat Conduction: -

For 2D Steady Heat Conduction, $\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} = 0$ $(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0$

First Order Numerical Method,

$\Rightarrow (\frac{T_{i + 1, j}^{n} - 2 T_{i, j}^{n} + T_{i - 1, j}^{n}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n} - 2 T_{i, j}^{n} + T_{i, j - 1}^{n}}{{d y}^{2}}) = 0$ $=> ((T_(i+1,j)^n - 2T_(i,j)^n + T_(i-1,j)^n)/dx^2) + ((T_(i,j+1)^n - 2T_(i,j)^n + T_(i,j-1)^n)/dy^2) = 0$

$\Rightarrow 2 T_{i, j}^{n} (\frac{1}{{d x}^{2}} + \frac{1}{{d y}^{2}}) = (\frac{T_{i + 1, j}^{n} + T_{i - 1, j}^{n}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n} + T_{i, j - 1}^{n}}{{d y}^{2}})$ $=>2T_(i,j)^n(1/dx^2 + 1/dy^2) = ((T_(i+1,j)^n + T_(i-1,j)^n)/dx^2) + ((T_(i,j+1)^n + T_(i,j-1)^n)/dy^2)$

$\Rightarrow 2 T_{i, j}^{n} (\frac{{d x}^{2} + {d y}^{2}}{{d x}^{2} {d y}^{2}}) = (\frac{T_{i + 1, j}^{n} + T_{i - 1, j}^{n}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n} + T_{i, j - 1}^{n}}{{d y}^{2}})$ $=>2T_(i,j)^n((dx^2 + dy^2)/(dx^2dy^2)) = ((T_(i+1,j)^n + T_(i-1,j)^n)/dx^2) + ((T_(i,j+1)^n + T_(i,j-1)^n)/dy^2)$

$\Rightarrow T_{i, j}^{n} = \frac{1}{2} [(\frac{T_{i + 1, j}^{n} + T_{i - 1, j}^{n}}{{d x}^{2} + {d y}^{2}}) {d y}^{2} + (\frac{T_{i, j + 1}^{n} + T_{i, j - 1}^{n}}{{d x}^{2} + {d y}^{2}}) {d x}^{2}]$ $=>T_(i,j)^n = 1/2[((T_(i+1,j)^n + T_(i-1,j)^n)/(dx^2 + dy^2))dy^2 + ((T_(i,j+1)^n + T_(i,j-1)^n)/(dx^2 + dy^2))dx^2]$

$Let \frac{{d y}^{2}}{{d x}^{2} + {d y}^{2}} be k_{1} and Let \frac{{d x}^{2}}{{d x}^{2} + {d y}^{2}} be k_{2}$ $"Let "(dy^2)/(dx^2 + dy^2) " be " k_1 and "Let "(dx^2)/(dx^2 + dy^2) " be " k_2$

$\Rightarrow T_{i, j}^{n} = \frac{1}{2} [(T_{i + 1, j}^{n} + T_{i - 1, j}^{n}) k_{1} + (T_{i, j + 1}^{n} + T_{i, j - 1}^{n}) k_{2}]$ $=>T_(i,j)^n = 1/2[(T_(i+1,j)^n + T_(i-1,j)^n)k_1 + (T_(i,j+1)^n + T_(i,j-1)^n)k_2]$

Now,

$For d x = \frac{1}{51 - 1} and d y = \frac{1}{51 - 1}$ $"For " dx = 1/(51-1) and dy = 1/(51-1)$

$k_{1} = \frac{{(d y = \frac{1}{50})}^{2}}{{(d x = \frac{1}{50})}^{2} + {(d y = \frac{1}{50})}^{2}} = \frac{1}{2}$ $k_1 = (dy = 1/50)^2/((dx = 1/50)^2 + (dy = 1/50)^2) = 1/2$

$k_{2} = \frac{{(d x = \frac{1}{50})}^{2}}{{(d x = \frac{1}{50})}^{2} + {(d y = \frac{1}{50})}^{2}} = \frac{1}{2}$ $k_2 = (dx = 1/50)^2/((dx = 1/50)^2 + (dy = 1/50)^2) = 1/2$

$\Rightarrow T_{i, j}^{n} = \frac{1}{2} [(T_{i + 1, j}^{n} + T_{i - 1, j}^{n}) (k_{1} = \frac{1}{2}) + (T_{i, j + 1}^{n} + T_{i, j - 1}^{n}) (k_{2} = \frac{1}{2})]$ $=>T_(i,j)^n = 1/2[(T_(i+1,j)^n + T_(i-1,j)^n)(k_1 = 1/2) + (T_(i,j+1)^n + T_(i,j-1)^n)(k_2 = 1/2)]$

Hence,

$T_{i, j}^{n} = \frac{T_{i + 1, j}^{n} + T_{i - 1, j}^{n}}{4} + \frac{T_{i, j + 1}^{n} + T_{i, j - 1}^{n}}{4}$ $T_(i,j)^n = (T_(i+1,j)^n + T_(i-1,j)^n)/4 + (T_(i,j+1)^n + T_(i,j-1)^n)/4$

2D Steady State Heat Conduction Equation

T_{i, j}^{n} = \frac{T_{i + 1, j}^{n} + T_{i - 1, j}^{n}}{4} + \frac{T_{i, j + 1}^{n} + T_{i, j - 1}^{n}}{4}

$T_(i,j)^n = (T_(i+1,j)^n + T_(i-1,j)^n)/4 + (T_(i,j+1)^n + T_(i,j-1)^n)/4$

B. 2D Transient Heat Conduction: -

For 2D Steady Heat Conduction, $\frac{\partial T}{\partial t} - α (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}}) = 0$ $(delT)/(delt)-alpha((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0$

First Order Numerical Method,

$\Rightarrow (\frac{T_{i, j}^{n + 1} - T_{i, j}^{n}}{d t}) - α [(\frac{T_{i + 1, j}^{n + 1} - 2 T_{i, j}^{n + 1} + T_{i - 1, j}^{n + 1}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n + 1} - 2 T_{i, j}^{n + 1} + T_{i, j - 1}^{n + 1}}{{d y}^{2}})] = 0$ $=> ((T_(i,j)^(n+1) - T_(i,j)^n)/dt)-alpha[((T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1))/dy^2)] = 0$

$\Rightarrow T_{i, j}^{n + 1} - T_{i, j}^{n} = (d t) α [(\frac{T_{i + 1, j}^{n + 1} - 2 T_{i, j}^{n + 1} + T_{i - 1, j}^{n + 1}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n + 1} - 2 T_{i, j}^{n + 1} + T_{i, j - 1}^{n + 1}}{{d y}^{2}})]$ $=> T_(i,j)^(n+1) - T_(i,j)^n = (dt)alpha[((T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1))/dy^2)]$

$\Rightarrow T_{i, j}^{n + 1} = T_{i, j}^{n} + (d t) α [(\frac{T_{i + 1, j}^{n + 1} - 2 T_{i, j}^{n + 1} + T_{i - 1, j}^{n + 1}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n + 1} - 2 T_{i, j}^{n + 1} + T_{i, j - 1}^{n + 1}}{{d y}^{2}})]$ $=> T_(i,j)^(n+1) = T_(i,j)^n + (dt)alpha[((T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1))/dy^2)]$

$\Rightarrow T_{i, j}^{n + 1} (1 + 2 (d t) α (\frac{1}{{d x}^{2}} + \frac{1}{{d y}^{2}})) = T_{i, j}^{n} + (d t) α [(\frac{T_{i + 1, j}^{n + 1} + T_{i - 1, j}^{n + 1}}{{d x}^{2}}) + (\frac{T_{i, j + 1}^{n + 1} + T_{i, j - 1}^{n + 1}}{{d y}^{2}})]$ $=> T_(i,j)^(n+1)(1 + 2(dt)alpha(1/dx^2+1/dy^2)) = T_(i,j)^n + (dt)alpha[((T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1))/dx^2) + ((T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))/dy^2)]$

As, $C F L_{x} = \frac{α (d t)}{{d x}^{2}}$ $CFL_x = (alpha(dt))/(dx^2)$ and $C F L_{y} = \frac{α (d t)}{{d y}^{2}}$ $CFL_y = (alpha(dt))/(dy^2)$

For Stable Numerical Method $C F L < 1$ $CFL < 1$

Now,

Let $\frac{α (d t)}{{d x}^{2}} = k_{1}$ $(alpha(dt))/dx^2 = k_1$

Let $\frac{α (d t)}{{d y}^{2}} = k_{2}$ $(alpha(dt))/dy^2 = k_2$

Let $1 + 2 (d t) α (\frac{1}{{d x}^{2}} + \frac{1}{{d y}^{2}}) = k_{3}$ $1 + 2(dt)alpha(1/dx^2+1/dy^2) = k_3$

$\Rightarrow k_{3} = 1 + 2 (k_{1} + k_{2})$ $=> k_3 = 1 + 2(k_1 + k_2)$

Hence,

$\Rightarrow T_{i, j}^{n + 1} (k_{3}) = T_{i, j}^{n} + k_{1} (T_{i + 1, j}^{n + 1} + T_{i - 1, j}^{n + 1}) + k_{2} (T_{i, j + 1}^{n + 1} + T_{i, j - 1}^{n + 1})$ $=> T_(i,j)^(n+1)(k_3) = T_(i,j)^n + k_1(T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1)) + k_2(T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))$

$\Rightarrow T_{i, j}^{n + 1} = \frac{1}{k_{3}} [T_{i, j}^{n} + k_{1} (T_{i + 1, j}^{n + 1} + T_{i - 1, j}^{n + 1}) + k_{2} (T_{i, j + 1}^{n + 1} + T_{i, j - 1}^{n + 1})]$ $=> T_(i,j)^(n+1) = 1/k_3[T_(i,j)^n + k_1(T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1)) + k_2(T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))]$

2D Unsteady State Heat Conduction Equation

T_{i, j}^{n + 1} = \frac{1}{k_{3}} [T_{i, j}^{n} + k_{1} (T_{i + 1, j}^{n + 1} + T_{i - 1, j}^{n + 1}) + k_{2} (T_{i, j + 1}^{n + 1} + T_{i, j - 1}^{n + 1})]

$T_(i,j)^(n+1) = 1/k_3[T_(i,j)^n + k_1(T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1)) + k_2(T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1))]$

5. Explanation Code: -

A. 2D Steady-State Heat Conduction: -

i. Section 1(Defining Variables)

a. Defining Geometry of Surface

$\to$ $rarr$ Defined Length and Width

$\to$ $rarr$ Defined nodes along x and y

$\to$ $rarr$ Defined dx and dy based on nodes and Dimension

$\to$ $rarr$ Created X and Y Vectors with the geometric distance of each node

b. Defining Boundary Condition and Initial Surface Temperature

$\to$ $rarr$ Defined each node of Surface to be at 500 K

$\to$ $rarr$ Defined Top, Bottom, Left, and Right Nodes at a Certain Temperature and the Corners as the average of adjacent nodes.

c. Variables for Iteration

$\to$ $rarr$ Defined Set the tolerance for the Temperature change to be $10^{- 5}$ $10^(-5)$

$\to$ $rarr$ Defined initial error = 1

$\to$ $rarr$ Defined iteration counter and set value = 0

d. Method Selection

$\to$ $rarr$ Defined a Variable "Method" whose value is used to select the method for iteration

1 - Jacobi Method, 2 - Gauss-Siedel, and 3 - S.O.R. with Gauss-Siedel( $ω = 1.95$ $omega = 1.95$ )

ii. Section 2(Computation)

a. Using Jacobi

$\to$ $rarr$ Checked if Method Selected is Jacobi using the "Method" Variable

$\to$ $rarr$ displayed that the iteration in use is Jacobi

$\to$ $rarr$ Created a While loop with condition that the loop terminates after (error for Temperature < Tolerance of Temperature)

$\to$ $rarr$ Applied Jacobi Iteration on each Node for Value of Temperature based on the derived formula.

$\to$ $rarr$ Updated Value of error of Temperature that is maximum error Node calculated from $T_{new} - T_{old}$ $T_("new") - T_("old")$

b. Using Gauss-Siedel

$\to$ $rarr$ Checked if Method Selected is Gauss-Siedel using the "Method" Variable

$\to$ $rarr$ displayed that the iteration in use is Gauss-Siedel

$\to$ $rarr$ Created a While loop with condition that the loop terminates after (error for Temperature < Tolerance of Temperature)

$\to$ $rarr$ Applied Gauss-Siedel Iteration on each Node for Value of Temperature based on the derived formula.

$\to$ $rarr$ Updated Value of error of Temperature that is maximum error Node calculated from $T_{new} - T_{old}$ $T_("new") - T_("old")$

c. Using S.O.R. using Gauss-Siedel

$\to$ $rarr$ Checked if Method Selected is S.O.R. using Gauss-Siedel using the "Method" Variable

$\to$ $rarr$ displayed that the iteration in use is S.O.R. using Gauss-Siedel

$\to$ $rarr$ Created a While loop with condition that the loop terminates after (error for Temperature < Tolerance of Temperature)

$\to$ $rarr$ Applied S.O.R. using Gauss-Siedel Iteration on each Node for Value of Temperature based on the derived formula.

$\to$ $rarr$ Updated Value of error of Temperature that is maximum error Node calculated from $T_{new} - T_{old}$ $T_("new") - T_("old")$

iii. Section 3(OUTPUT and Plotting)

$\to$ $rarr$ displaying the Total number of iterations

$\to$ $rarr$ Displaying the latest updated error

$\to$ $rarr$ Plotted the Graph of Temperature Distribution

B. 2D UnSteady-State Heat Conduction: -

i. Section 1(Defining Variables)

a. Defining Geometry of Surface

$\to$ $rarr$ Defined Length and Width

$\to$ $rarr$ Defined nodes along x and y

$\to$ $rarr$ Defined dx and dy based on nodes and Dimension

$\to$ $rarr$ Created X and Y Vectors with the geometric distance of each node

b. Defining Boundary Condition and Initial Surface Temperature

$\to$ $rarr$ Defined each node of Surface to be at 500 K

$\to$ $rarr$ Defined Top, Bottom, Left, and Right Nodes at a Certain Temperature and the Corners as the average of adjacent nodes.

c. Variables for Iteration

$\to$ $rarr$ Defined Set the tolerance for the Temperature change to be $10^{- 4}$ $10^(-4)$

$\to$ $rarr$ Defined initial error = 1

$\to$ $rarr$ Defined iteration counter and set value = 0

d. Time Step Variable

$\to$ $rarr$ Defined value of CFL = 0.9

$\to$ $rarr$ Defined $α = 1.2$ $alpha = 1.2$

$\to$ $rarr$ Calculated dt

$\to$ $rarr$ Defined k1, k2 and k3

$\to$ $rarr$ Defined total time Step and a Vector to define each time step

e. Method Selection

$\to$ $rarr$ Defined a Variable "Method" whose value is used to select the method for iteration

1 - Jacobi Method, 2 - Gauss-Siedel, and 3 - S.O.R. with Gauss-Siedel( $ω = 1.44$ $omega = 1.44$ )

ii. Section 2(Computation)

a. Time Marching Loop

$\to$ $rarr$ Defined a for loop from 1 to total time steps

$\to$ $rarr$ Recorded current value of T in $T_{t}$ $T_t$ Matrix

$\to$ $rarr$ Created a While Loop based on error value and Tolerance Value

$\to$ $rarr$ Selected Each node and based on the value of "Method" Selected each iteration and applied formula based on derivation.

$\to$ $rarr$ Within each time Step Plotted the Change in temperature

$\to$ $rarr$ Placed a Break Command which is based on Tolerance of Temperature and breaks the loop when Stability is Reached.

iii. Section 3(OUTPUT)

$\to$ $rarr$ Displayed The Type of Iteration under use

$\to$ $rarr$ Displayed the Total number of Iterations and Total Time taken to reach Stability.

6. Code: -

A. 2D Steady-State Heat Conduction: -

% 2D Steady State Heat Conduction using Iterative Methods
  clear all
  close all
  clc

%% Defining Variables
 % Geometry of Surface
   L = 1;
   W = 1;
   Node_x = 51;
   Node_y = 51;
   dx = L/(Node_x-1);
   dy = W/(Node_y-1);
   X = [0:dx:L];
   Y = [0:dy:W];

 % Boundary Conditions and Initial Surface Temperature
   T = 500*ones(Node_x,Node_y);
   T_top = 600;
   T_bottom = 900;
   T_left = 400;
   T_right = 800;

   T(2:(Node_x -1),Node_y) = T_top;
   T(2:(Node_x - 1),1) = T_bottom;
   T(1,2:(Node_y - 1)) = T_left;
   T(Node_x,2:(Node_y -1)) = T_right;

   T(1,1) = 0.5*(T(1,2)+T(2,1));
   T(Node_x,Node_y) = 0.5*(T(Node_x,(Node_y-1))+T((Node_x-1),Node_y));
   T(Node_x,1) = 0.5*(T(Node_x,2)+T((Node_x-1),1));
   T(1,Node_y) = 0.5*(T(1,(Node_y-1))+T(2,Node_y));

 % Variables for Iterations
   Tol_T = 1e-5;
   err_T = 1;
   iter = 0;
   
 % Method Selection 
 % 1 - Jacobi, 
 % 2 - Gauss-Siedel
 % 3 - S.O.R. with Gauss-Siedel
   Method = 2;

%% Computation
 % Using Jacobi
   if Method == 1
       disp('Iteration Method in Use is Jacobi')
       while (err_T > Tol_T)
           iter = iter+1;
           T_old = T;
           for i = 2:(Node_x - 1)
               for j = 2:(Node_y - 1)
                   T(i,j) = 0.25*(T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1));
               end
           end
           err_T = max(abs(T - T_old),[],'all');
       end
   end

  % Using Gauss-Siedel
    if Method == 2
        disp('Iteration Method in Use is Gauss-Siedel')
        while (err_T > Tol_T)
            iter = iter+1;
            T_old = T;
            for i = 2:(Node_x - 1)
                for j = 2:(Node_y - 1)
                    T(i,j) = 0.25*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1));
                end
            end
            err_T = max(abs(T - T_old),[],'all');
        end
    end

  % Using S.O.R. and Gauss-Siedel
    if Method == 3
        disp('Iteration Method in Use is S.O.R.')
        omega_val = 1.95;
        disp(['Value of Omega = ',num2str(omega_val)])
        while (err_T > Tol_T)
            iter = iter+1;
            T_old = T;
            for i = 2:(Node_x - 1)
                for j = 2:(Node_y - 1)
                    T(i,j) = T(i,j)*(1-omega_val) + (omega_val)*0.25*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1));
                end
            end
            err_T = max(abs(T - T_old),[],'all');
        end
    end

%% OUTPUT and Plotting
   disp(['Total Iteration = ',num2str(iter)])
   disp(['Error = ',num2str(err_T)])
   colormap(jet);
   contourf(X,Y,T');
   colorbar
   title('Temperature Distribution','FontSize',16)
   xlabel('Length of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)
   ylabel('Width of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)

B. 2D UnSteady-State Heat Conduction: -

% 2D Un-Steady State Heat Conduction using Iterative Methods
  clear all
  close all
  clc

%% Defining Variables
 % Geometry of Surface
   L = 1;
   W = 1;
   Node_x = 51;
   Node_y = 51;
   dx = L/(Node_x-1);
   dy = W/(Node_y-1);
   X = [0:dx:L];
   Y = [0:dy:W];

 % Boundary Conditions and Initial Guess
   T = 500*ones(Node_x,Node_y);
   T_top = 600;
   T_bottom = 900;
   T_left = 400;
   T_right = 800;

   T(2:(Node_x -1),Node_y) = T_top;
   T(2:(Node_x - 1),1) = T_bottom;
   T(1,2:(Node_y - 1)) = T_left;
   T(Node_x,2:(Node_y -1)) = T_right;

   T(1,1) = 0.5*(T(1,2)+T(2,1));
   T(Node_x,Node_y) = 0.5*(T(Node_x,(Node_y-1))+T((Node_x-1),Node_y));
   T(Node_x,1) = 0.5*(T(Node_x,2)+T((Node_x-1),1));
   T(1,Node_y) = 0.5*(T(1,(Node_y-1))+T(2,Node_y));

 % Variables for Iterations
   Tol_T = 1e-4;
   err_T = 1;
   iter = 0;

 % Time Step Variables
   CFL = 0.9;
   alp = 1.2;
   dt = CFL*(dx^2)/alp;
   k1 = CFL;
   k2 = CFL;
   k3 = 1 + 2*k1 + 2*k2;
   t_steps = 1500;
   total_t = [0:dt:(t_steps*dt)];

 % Method Selection 
 % 1 - Jacobi, 
 % 2 - Gauss-Siedel
 % 3 - S.O.R. with Gauss-Siedel
   Method = 3;
   omega_val = 1.44;

%% Computation
 % Time Marching Loop
   for t_count = 1:t_steps
       T_t = T;

     % Iteration using Gauss-Siedel
       err_T = 1;
       while (err_T > Tol_T)
            iter = iter+1;
            T_old = T;
            for i = 2:(Node_x - 1)
                for j = 2:(Node_y - 1)
                    if Method == 1
                        T(i,j) = (T_t(i,j) +  k1*(T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1)))/k3;
                    end
                    if Method == 2
                        T(i,j) = (T_t(i,j) +  k1*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1)))/k3;
                    end
                    if Method == 3
                        T(i,j) = T(i,j)*(1 - omega_val) + (omega_val)*(T_t(i,j) +  k1*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1)))/k3;
                    end
                end
            end
            err_T = max(abs(T - T_old),[],'all');
       end

     % Plotting after each Iteration
       colormap(jet);
       contourf(X,Y,T');
       colorbar
       title(['Temperature Distribution at time = ',num2str(total_t(t_count)),'sec'],'FontSize',12)
       xlabel('Length of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)
       ylabel('Width of Surface(m)','Color',[0.15 0.5 0.5],'FontSize',14)
       pause(0.05)
       err_T_t = max(abs(T - T_t),[],'all');

     % Breaking Time Marching if Change become less than Tolerance
       if err_T_t < Tol_T
           break
       end
   end

%% OUTPUTS
   if Method == 1
       disp('The Iteration method in use was Jacobi')
   end
   if Method == 2
       disp('The Iteration method in use was Gauss-Siedel')
   end
   if Method == 3
       disp('The Iteration method in use was S.O.R. with Gauss-Siedel')
       disp(['Value of Omega for S.O.R. = ',num2str(omega_val)])
   end
   disp(['Total Iteration = ',num2str(iter)])
   disp(['Total Time for reaching Stability = ',num2str(total_t(t_count)),'sec'])

7. OUTPUTS and Plots Generated: -

A. 2D Steady-State Heat Conduction: -

i. Plot Generated: -

Temp Distribution in 2D Steady Heat Conduction