Menu

Executive Programs

Workshops

Projects

Blogs

Careers

Placements

Student Reviews


For Business


More

Academic Training

Informative Articles

Find Jobs

We are Hiring!


All Courses

Choose a category

Loading...

All Courses

All Courses

logo

Loading...
Executive Programs
Workshops
For Business

Success Stories

Placements

Student Reviews

More

Projects

Blogs

Academic Training

Find Jobs

Informative Articles

We're Hiring!

phone+91 9342691281Log in
  1. Home/
  2. Dushyanth Srinivasan/
  3. Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

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

Steady state is referred to when the temperature with respect to time does not change. Unsteady/Transient State is shen the temperature with respect to time does change. The 2nd order partial diffential equation for Unsteady State is: δTδt-α⋅(δ2Tδx2+δ2Tδy2)=0δTδt−α⋅(δ2Tδx2+δ2Tδy2)=0…

    • Dushyanth Srinivasan

      updated on 20 Oct 2021

    Steady state is referred to when the temperature with respect to time does not change.

    Unsteady/Transient State is shen the temperature with respect to time does change.

    The 2nd order partial diffential equation for Unsteady State is:

    δTδt-α⋅(δ2Tδx2+δ2Tδy2)=0

    During steady state, δTδt=0

    hence the equation becomes,

    δ2Tδx2+δ2Tδy2=0

    We solve the equations in a total of 7 ways:

    1. Steady State - Jacobi Method

    this is the PDE to be solved.

    δ2Tδx2+δ2Tδy2=0

    discretising this and rearranging the terms, we get:

     TP=1k(TL+TRΔx2)+1k(TT+TBΔy2)

    where k=2(Δx2+Δy2Δx2⋅Δy2)

    here, the values of T on the right hand side are the values of the previous iteration and not the current iteration

    Function code which takes dx, dy, nx, ny, T, error and tolerance as arguments

    filename: steady_jacobi_solver.m
    function [T,jacobi_iterations] = steady_jacobi_solver (dx, dy,nx,ny,T,error,tol)
        % creating copy of solution array
        Told = T;
        k = 2 * ( (dx^2 + dx^2)/(dx*dy)^2);
        jacobi_iterations = 1;
        % convergence loop
        while(error>tol)
                for i = 2:nx-1
                    for j = 2:ny-1
                        T(i,j)=((Told(i-1,j) + Told(i+1,j))/(k*dx^2)) + ((Told(i,j-1) + Told(i,j+1))/(k*dy^2));
                    end
                end
            % calculating error and updating Told
            error = max(max(abs(Told- T)));
            Told = T;
            jacobi_iterations = jacobi_iterations + 1;
        end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by a Jacobi solver for steady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied

    2. Steady State - Gauss Seidel Method

    this is the PDE to be solved.

    δ2Tδx2+δ2Tδy2=0

    discretising this and rearranging the terms, we get:

     TP=1k(TL+TRΔx2)+1k(TT+TBΔy2)

    where k=2(Δx2+Δy2Δx2⋅Δy2)

    here, the values of T on the right hand side are the values of the current iteration and not the previous iteration

    Function code which takes dx, dy, nx, ny, T, error and tolerance as arguments

    filename:steady_gauss_seidel_solver.m
    function [T,gs_iterations] = steady_gauss_seidel_solver (dx, dy,nx,ny,T,error,tol)
        % creating copy of solution array
        Told = T;
        k = 2 * ( (dx^2 + dx^2)/(dx*dy)^2);
        gs_iterations = 1;
        % convergence loop
        while(error>tol)
                for i = 2:nx-1
                    for j = 2:ny-1
                        T(i,j)=((T(i-1,j) + T(i+1,j))/(k*dx^2)) + ((T(i,j-1) + T(i,j+1))/(k*dy^2));
                    end
                end
            % calculating error and updating Told
            error = max(max(abs(Told- T)));
            Told = T;
            gs_iterations = gs_iterations + 1;
         end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by a Gauss Seidel solver for steady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied

    3. Steady State - Succesive over relaxation Method

    this is the PDE to be solved.

    δ2Tδx2+δ2Tδy2=0

    discretising this, rearranging the terms and applying succiesive over relaxation for gauss seidel, we get:

    Tn+1P=(1-ω)nT+ω⋅(1k(TL+TRΔx2)+1k(TT+TBΔy2))

    where ω is the relaxation factor

    w > 1 overrelaxed

    w < 1 underrelaxed

    and k=2(Δx2+Δy2Δx2⋅Δy2)

    here, the values of T on the right hand side are the values of the current iteration and not the previous iteration

    Function code which takes dx, dy, nx, ny, T, error and tolerance as arguments

    filename:steady_sor_solver.m
    function [T,sor_iterations] = steady_sor_solver (dx, dy,nx,ny,T,error,tol)
        % creating copy of solution array
        Told = T;
        k = 2 * ( (dx^2 + dx^2)/(dx*dy)^2);
        w = 1.2; % relaxation factor or 1.2
        sor_iterations = 1;
        % convergence loop
        while(error>tol)
                for i = 2:nx-1
                    for j = 2:ny-1
                        T(i,j)=(1-w)*(Told(i,j))+w*(((T(i-1,j) + T(i+1,j))/(k*dx^2)) + ((T(i,j-1) + T(i,j+1))/(k*dy^2)));
                    end
                end
            % calculating error and updating Told
            error = max(max(abs(Told- T)));
            Told = T;
            sor_iterations = sor_iterations + 1;
        end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by a succesive relaxation solver for gauss seidel steady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied

    4. Unsteady State - Explicit

    this is the PDE to be solved.

    δTδt-α⋅(δ2Tδx2+δ2Tδy2)=0

    discretising this and rearranging the terms, we get:

    Tn+1P=TnP+k1â‹…(TL-2â‹…TP+TR)n+k2â‹…(TT-2â‹…TP+TB)n

    where α is the heat transfer coefficient, k1=αΔtΔx2 and k2=αΔtΔy2

    Note: the exponent n denotes that this takes place in the n th timestep

    here, the values of T on the right hand side are the values of the previous timestep and not the current timestep

    Function code which takes dx, dy, nx, ny and T as arguments

    filename: unsteady_explicit_solver.m
    function [T,explicit_iterations] = unsteady_explicit_solver (dx, dy,nx,ny,T)
        % initial conditions
        dt = 1e-3;
        alpha = 1.4;
        k1 = alpha * dt / dx^2;
        k2 = alpha * dt / dy^2;
        
        explicit_iterations = 1;
        nt = 1400;
        % time loop
        for k = 1:nt
            % convergence loop
            for i = 2:nx-1
                for j = 2:ny-1
                    T(i,j)= T(i,j) + k1 *(T(i-1,j) - 2*T(i,j) + T(i+1,j)) + k2 * (T(i,j-1) - 2*T(i,j) + T(i,j+1));
                end
            end
            explicit_iterations = explicit_iterations + 1;
        end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by an explicit solver for unsteady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied

    5. Unsteady State - Implicit - Jacobi Solver

    this is the PDE to be solved.

    δTδt-α⋅(δ2Tδx2+δ2Tδy2)=0

    discretising this and rearranging the terms, we get:

    Tn+1P=(11+2â‹…k1+2â‹…k2)â‹…(TnP+k1â‹…(TL+TR)n+1+k2â‹…(TT+TB)n+1)

    where α is the heat transfer coefficient, k1=αΔtΔx2 and k2=αΔtΔy2

    Note: the exponent n+1 denotes that this takes place in the n+1 th timestep

    here, the values of T on the right hand side are the values of the previous iteration and not the current iteration

    Function code which takes dx, dy, nx, ny, T, error and tolerance as arguments

    filename: unsteady_jacobi_solver.m
    function [T,unsteady_jacobi_iterations] = unsteady_jacobi_solver (dx, dy,nx,ny,T,error,tol)
        % creating copies of solution array
        T_prev_dt = T;
        Told = T;
    
        % initial conditions
        dt = 1e-3;
        alpha = 1.4;
    
        k1 = alpha * dt / dx^2;
        k2 = alpha * dt / dy^2;
        unsteady_jacobi_iterations = 1;
        nt = 1400;
        % Time loop
        for k = 1:nt
            % convergence loop
            while(error>tol)
                for i = 2:nx-1
                    for j = 2:ny-1
                        T(i,j)= (1/(1+2*k1+2*k2)) * ( T_prev_dt(i,j) +( k1 * (Told(i-1,j) + Told (i+1,j))) +( k2 * (Told (i,j-1) + Told (i,j+1))));
                    end
                end
                % calculating error and updating Told
                error = max(max(abs(Told- T)));
                Told = T;
                unsteady_jacobi_iterations = unsteady_jacobi_iterations + 1;
            end
            % resetting error and updating T of previous time step
            T_prev_dt = T;
            error = 1;
        end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by a jacobi solver for unsteady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied

    6. Unsteady State - Implicit - Gauss Seidel Solver

    this is the PDE to be solved.

    δTδt-α⋅(δ2Tδx2+δ2Tδy2)=0

    discretising this and rearranging the terms, we get:

    Tn+1P=(11+2â‹…k1+2â‹…k2)â‹…(TnP+k1â‹…(TL+TR)n+1+k2â‹…(TT+TB)n+1)

    where α is the heat transfer coefficient, k1=αΔtΔx2 and k2=αΔtΔy2

    Note: the exponent n+1 denotes that this takes place in the n+1 th timestep

    here, the values of T on the right hand side are the values of the current iteration and not the previous iteration

    Function code which takes dx, dy, nx, ny, T, error and tolerance as argument

    filename: unsteady_sor_solver.m
    function [T,unsteady_sor_iterations] = unsteady_sor_solver (dx, dy,nx,ny,T,error,tol)
        % creating copies of solution array
        T_prev_dt = T;
        Told = T;
        
        % initial conditions
        dt = 1e-3;
        alpha = 1.4;
        w = 1.1; % relaxation factor or 1.2
        k1 = alpha * dt / dx^2;
        k2 = alpha * dt / dy^2;
        unsteady_sor_iterations = 1;
        nt = 1400;
        % time loop
        for k = 1:nt
            % convergence loop
            while(error>tol)
                for i = 2:nx-1
                    for j = 2:ny-1
                        T(i,j)= (1-w)*(Told(i,j))+w*((1/(1+2*k1+2*k2)) * ( T_prev_dt(i,j) +( k1 * (T(i-1,j) + T (i+1,j))) +( k2 * (T (i,j-1) + T (i,j+1)))));
                    end
                end
                % calculating error and updating Told
                error = max(max(abs(Told- T)));
                Told = T;
                unsteady_sor_iterations = unsteady_sor_iterations + 1;
            end
            % resetting error and updating T of previous time step
            T_prev_dt = T;
            error = 1;
        end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by a gauss seidel solver for unsteady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied

    7. Unsteady State - Implicit - Succesive over relaxation solver

    this is the PDE to be solved.

    δTδt-α⋅(δ2Tδx2+δ2Tδy2)=0

    discretising this, rearranging the terms and applying succiesive over relaxation for gauss seidel, we get:

    Tn+1P=(1-ω)nT+ω⋅((11+2⋅k1+2⋅k2)⋅(TnP+k1⋅(TL+TR)n+1+k2⋅(TT+TB)n+1))

    where α is the heat transfer coefficient, ω is the relaxation factork1=αΔtΔx2 and k2=αΔtΔy2

    w > 1 overrelaxed

    w < 1 underrelaxed

    Note: the exponent n+1 denotes that this takes place in the n+1 th timestep

    here, the values of T on the right hand side are the values of the current iteration and not the previous iteration

    Function code which takes dx, dy, nx, ny, T, error and tolerance as argument

     

    filename: unsteady_sor_solver.m
    function [T,unsteady_sor_iterations] = unsteady_sor_solver (dx, dy,nx,ny,T,error,tol)
        % creating copies of solution array
        T_prev_dt = T;
        Told = T;
        
        % initial conditions
        dt = 1e-3;
        alpha = 1.4;
        w = 1.1; % relaxation factor or 1.2
        k1 = alpha * dt / dx^2;
        k2 = alpha * dt / dy^2;
        unsteady_sor_iterations = 1;
        nt = 1400;
        % time loop
        for k = 1:nt
            % convergence loop
            while(error>tol)
                for i = 2:nx-1
                    for j = 2:ny-1
                        T(i,j)= (1-w)*(Told(i,j))+w*((1/(1+2*k1+2*k2)) * ( T_prev_dt(i,j) +( k1 * (T(i-1,j) + T (i+1,j))) +( k2 * (T (i,j-1) + T (i,j+1)))));
                    end
                end
                % calculating error and updating Told
                error = max(max(abs(Told- T)));
                Told = T;
                unsteady_sor_iterations = unsteady_sor_iterations + 1;
            end
            % resetting error and updating T of previous time step
            T_prev_dt = T;
            error = 1;
        end
    end

    Figure shows a filled contour plot of the temperature distribution in the domain calculated by succesive overrelaxation for gauss seidel for unsteady state. As expected, the boundary conditions (bottom: 900K, left: 400K, right: 800K and top: 600K) are satisfied.

    Conclusions and Observations

    State of Domain Method Used Number of Iterations
    Steady State Jacobi 208
    Steady State Gauss Seidel 112
    Steady State Succesive Overrelaxation 76
    Unsteady State Explicit 1401
    Unsteady State Jacobi 3504
    Unsteady State Gauss Seidel 2964
    Unsteady State Succesive Overrelaxation 1401

     

     

     

     

     

     

     

     

     

    We can notice that steady state solutions convergence significantly faster than unsteady state equations, this is because steady state equations are more simple than unsteady equations.

     

    We can also notice a pattern between the 3 solvers, both in steady and unsteady states. Jacobi solvers take the most number of iterations, followed by gauss seidel and then succesive overrelaxation.

    The decrease in iterations is also not linear, it seems to follow an exponential decrease

     

    Main Function

     

    clear all
    close all
    clc
    
    nx = 10; % gridpoints alon x more gridpoints makes more iterations
    ny = nx;
    L = 1;
    
    % initial conditions top = 600K bottom = 900K left = 400K right = 800K
    
    T = 300*ones(nx);
    T(:,1) = 400;
    T(:,nx) = 800;
    T(1,:) = 600;
    T(nx,:) = 900;
    
    % generating nodes and grid spacing
    x = linspace (0,1,nx+1);
    y = linspace(0,1,ny+1);
    dx = L/(nx);
    dy = L/(ny);
    
    % predetermined error and tolerance
    error = 1;
    tol = 1e-4;
    
    % these functions calculate the solution array and number of iterations
    % used. they take dx, dy, nx, ny, T, error and tolerance as arguments
    [T_solution,i] = steady_jacobi_solver (dx, dy,nx,ny,T,error,tol);
    figure(1)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Steady State - Jacobi Method. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')
    
    [T_solution,i] = steady_gauss_seidel_solver (dx, dy,nx,ny,T,error,tol);
    figure(2)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Steady State - Gauss Seidel Method. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')
    
    [T_solution,i] = steady_sor_solver (dx, dy,nx,ny,T,error,tol);
    figure(3)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Steady State - SOR Method. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')
    
    [T_solution,i] = unsteady_explicit_solver (dx, dy,nx,ny,T);
    figure(4)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Unsteady State - Explicit. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')
    
    [T_solution,i] = unsteady_jacobi_solver (dx, dy,nx,ny,T,error,tol);
    figure(5)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Unsteady State - Jacobi. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')
    
    [T,i] = unsteady_gauss_seidel_solver (dx, dy,nx,ny,T,error,tol);
    figure(6)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Unsteady State - Gauss Seidel. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')
    
    [T,i] = unsteady_sor_solver (dx, dy,nx,ny,T,error,tol);
    figure(7)
    contourf(T_solution,[450:50:900],'ShowText','on')
    title(['Unsteady State - SOR. Number of Iterations :',num2str(i)])
    set(gca,'YDIR','reverse')

     

     

    Leave a comment

    Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.

    Please  login to add a comment

    Other comments...

    No comments yet!
    Be the first to add a comment

    Read more Projects by Dushyanth Srinivasan (45)

    Project 2 - Rankine cycle Simulator

    Objective:

      In this project, I will be writing code in MATLAB to simulate a Rankine Cycle for the given parameters. A Rankine Cycle is an ideal thermodynamic heat cycle where mechanical work is extracted from the working fluid as it passes between a heat source and heat sink. This cycle or its derivatives is used in steam engines…

    calendar

    04 Sep 2022 12:52 PM IST

    • MATLAB
    Read more

    Project 1 - Parsing NASA thermodynamic data

    Objective:

    In this project, I will be parsing a data file prepared by NASA. The contents of the data file can be used to generated thermodynamic properties such as Specific Heat at Constant Pressure 'C_p' (J/(kg.K)), Enthalpy HorQ (J) and Entropy S (J/(kg.mol)) at various temperatures. The files will be parsed in MATLAB…

    calendar

    31 Aug 2022 01:07 PM IST

    • MATLAB
    Read more

    Week 5 - Genetic Algorithm

    Objective:

    In this project, I will be generating a stalagmite function in MATLAB and find the global maxima of the function using Genetic Algorithm. A stalagmite function is a function which generates Stalactites, which are named after a natural phenomenon where rocks rise up from the floor of due to accumulation of droppings of…

    calendar

    29 Aug 2022 07:55 AM IST

    • MATLAB
    Read more

    Week 4.1 - Solving second order ODEs

    Objective:

    In this project, I will be writing code in MATLAB to solve the motion of a simple pendulum. A simple pendulum motion's depends on Newton's Second Law. The equation which governs the motion of a simple pendulum is (with damping) d2θdt2+bmdθdt+gLsinθ=0 Where, θ is the angular displacement…

    calendar

    23 Aug 2022 08:06 AM IST

    • MATLAB
    Read more

    Schedule a counselling session

    Please enter your name
    Please enter a valid email
    Please enter a valid number

    Related Courses

    coursecardcoursetype

    Accelerated Career Program in Embedded Systems (On-Campus) - Powered by NASSCOM

    Recently launched

    0 Hours of Content

    coursecard

    5G Protocol and Testing

    Recently launched

    4 Hours of Content

    coursecard

    Automotive Cybersecurity

    Recently launched

    9 Hours of Content

    coursecardcoursetype

    Pre-Graduate Program in Bioengineering and Medical Devices

    Recently launched

    90 Hours of Content

    coursecardcoursetype

    Pre-Graduate Program in 5G Design and Development

    Recently launched

    49 Hours of Content

    Schedule a counselling session

    Please enter your name
    Please enter a valid email
    Please enter a valid number

                Do You Want To Showcase Your Technical Skills?
                Sign-Up for our projects.