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

03D 18H 26M 20S

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

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

Objective: To perform steady and unsteady analysis of 2D heat conduction using different point iterative techniques. Given: boundary conditions of domain top 600K, bottom 900K, left 400K, and right boundary 800K Steady-state 2D 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$ Equation for 2D steady-state…

Bharath Gaddameedi
updated on 16 Nov 2021

Objective: To perform steady and unsteady analysis of 2D heat conduction using different point iterative techniques.

Given: boundary conditions of domain top 600K, bottom 900K, left 400K, and right boundary 800K

Steady-state 2D 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$ Equation for 2D steady-state conduction

Using finite differencing for discretizing the above equation with the central differencing scheme the equation would become

$\frac{T_{L} - 2 T_{P} + T_{R}}{Δ x^{2}} + \frac{T_{T} - 2 T_{P} + T_{B}}{Δ y^{2}} = 0$ $(T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2) = 0$

Rearranging the equation would give $T_{P} = (\frac{1}{k}) \cdot \frac{T_{L} + T_{R}}{Δ x^{2}} \cdot \frac{T_{T} + T_{R}}{Δ y^{2}}$ $T_P = (1/k)*(T_L+T_R)/(Deltax^2)*(T_T+T_R)/(Deltay^2)$ where k is defined as $2 \cdot \frac{Δ x^{2} + Δ y^{2}}{Δ x^{2} \cdot Δ y^{2}}$ $2*(Deltax^2+Deltay^2)/(Deltax^2*Deltay^2)$

Code

Assuming that the given domain is a square.

taking number of nodes in x and y direction same, which I have taken as 10.

Then we define domain dimensions, number of nodes and calculate mesh size.

As we are going to calculate the temperature, the temperature is initialized with a guess. Although the guess value does not affect the end state, it helps reach the solution faster. So an optimal guess would be beneficial. Here I have taken 300K as the initial value.

Later boundary conditions are assigned to the temperature matrix we defined earlier.

While solving the equation with more terms to avoid missing out on some terms by mistake we can split the actual equations into parts.

If loop is used to select the iteration methods.

The code structure is first we have convergence loop and spacial loops.

In the Jacobi method, we use old values to calculate the present values so a variable is defined to store the value as T_old.

In the case of the Gauss-Siedel method, we use the values that are already calculated in the iteration, unlike Jacobi which uses values from the previous iteration.

clear 
close all
clc

%code for solving 2d steady heat conduction equation
%domain we are considering is a square of 1 unit
nx = 10;
ny = nx; %number of nodes are same in the x and y directions.
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);%grid spacing in x and y directions

%boundary conditions 

%top = 600k bottom = 900k left = 400k and right = 800k
%initializing temperature in all the nodes except boundary nodes.
T = 300*ones(nx,ny);

T(1,:) = 600;%top face
T(nx,:) = 900;%bottom face
T(:,1) = 400;%left face
T(:,ny) = 800;%right face

%At the corner nodes we are taking average of two adjacent nodes
T(1,1) = (T(1,2)+T(2,1))/2;%top left corner
T(10,1) = (T(10,2)+T(9,1))/2;%bottom left corner
T(10,10) = (T(10,9)+T(9,10))/2;%bottom right corner
T(1,10) = (T(1,9)+T(2,10))/2;%top right corner

term = (dx^2*dy^2)/(2*(dx^2+dy^2));
tol = 1e-4;
error = 9e9;
T_old = T;

iteration_solver = 3;%Iteration_solver = 1-Jacobi ; 2-Gauss seidel ; 3-SOR

if iteration_solver == 1 %Jacobi method
    tic
    iter = 0;
    while (error > tol)
        for j = 2:ny-1
            for i = 2:nx-1
                T(i,j) = term*((T_old(i-1,j)+T_old(i+1,j))/(dx^2)+(T_old(i,j-1)+T_old(i,j+1))/(dy^2));
            end
        end
        error = max(max(abs(T-T_old)));
        T_old = T;
        iter = iter+1;
    end
    time = toc;
    
elseif iteration_solver == 2 %Gauss Seidel method
    tic
    iter = 0;
    while (error > tol)
        for j = 2:ny-1
            for i = 2:nx-1
                T(i,j) = term*((T(i-1,j)+T_old(i+1,j))/dx^2 + (T(i,j-1)+T_old(i,j+1))/dy^2);
            end
        end
        error = max(max(abs(T-T_old)));
        T_old = T;
        iter = iter+1;
    end
    time = toc;
    
elseif iteration_solver == 3 %SOR method
    omega = 1.5;
    tic
    iter = 0;
    while (error > tol)
        for j = 2:ny-1
            for i = 2:nx-1
                T(i,j) = term*((T(i-1,j)+T_old(i+1,j))/dx^2 + (T(i,j-1)+T_old(i,j+1))/dy^2);
                T(i,j) = T_old(i,j)*(1-omega)+(omega*T(i,j));
            end
        end
        error = max(max(abs(T-T_old)));
        T_old = T;
        iter = iter+1;
    end
    time = toc;
end

%Plotting temperature contour on X and Y axes
[c,h] = contourf(x,y,T);
clabel(c,h)
colorbar
colormap(jet)
set(gca,'YDir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
grid on
title({'Steady state temperature field';['Iteration Method = ',num2str(iteration_solver)];['No. of iterations = ',num2str(iter)];['Simulation time = ',num2str(time),' s']})

Graphs

By selecting different iterative methods we have obtained the plots.

From the graphs of temperature profiles, we can see that the number of iterations in the Jacobi method is higher than the other two.

Iacobi iterations > Gauss-Seidel iterations > SOR iterations

Unsteady 2D heat conduction

$\frac{\partial T}{\partial t} = α \cdot (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}})$ $(delT)/(delt) = alpha*((del^2T)/(delx^2)+(del^2T)/(dely^2))$ - equation fo 2D Transient heat conduction

$α =$ $alpha =$ thermal diffusivity.

this is a parabolic differential equation. Then using time marching solution and finite difference scheme we obtain

$\frac{T_{p}^{n + 1} - T_{p}}{Δ t} = α \cdot {(\frac{T_{L} - 2 T_{P} + T_{R}}{Δ x^{2}} + \frac{T_{T} - 2 T_{P} + T_{B}}{Δ y^{2}})}^{n or n + 1}$ $(T_p^(n+1)-T_p)/(Deltat) = alpha*((T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2))^(n or n+1)$

when the RHS is calculated at nth step it is called explicit solution and at (n+1)th step, it is called implicit method.

Explicit method

At the the nth level every value is known except $T_{p}^{n + 1}$ $T_p^(n+1)$ , so

$T_{p}^{n + 1} = T_{p}^{n} + k_{1} \cdot {(T_{L} - 2 \cdot T_{P} + T_{R})}^{n} + k_{2} \cdot {(T_{T} - 2 \cdot T_{P} + T_{B})}^{n}$ $T_p^(n+1) = T_p^(n) + k_1*(T_L-2*T_P+T_R)^n+k_2*(T_T-2*T_P+T_B)^n$

$k_{2} = \frac{α \cdot Δ t}{Δ y^{2}}$ $k_2 = (alpha*Deltat)/(Deltay^2)$ $k_{2} = \frac{α \cdot Δ t}{Δ y^{2}}$ $k_2 = (alpha*Deltat)/(Deltay^2$

Explicit methods are not stable inherently. CFL criteria helps us setup for stable solution.

the CFL number should be less than 0.5 so, $C F L_{x} + C F L_{y} \leq 0.5$ $CFL_x+CFL_y <=0.5$

In the code we use time loop then convergence loop and spatial loop.

clear
close all
clc

%code for solving 2d unsteady heat conduction equation
%domain we are considering is a square of 1 unit
nx = 10;
ny = nx; %number of nodes are same in the x and y directions.
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);%grid spacing in x and y directions
A = 1.2;
dt = 0.002;

%boundary conditions 

%top = 600k bottom = 900k left = 400k and right = 800k
%initializing temperature in all the nodes except boundary nodes.
T = 300*ones(nx,ny);

T(1,:) = 600;%top face
T(nx,:) = 900;%bottom face
T(:,1) = 400;%left face
T(:,ny) = 800;%right face

%At the corner nodes we are taking average of two adjacent nodes
T(1,1) = (T(1,2)+T(2,1))/2;%top left corner
T(10,1) = (T(10,2)+T(9,1))/2;%bottom left corner
T(10,10) = (T(10,9)+T(9,10))/2;%bottom right corner
T(1,10) = (T(1,9)+T(2,10))/2;%top right corner


k1 = A*dt/(dx^2);
k2 = A*dt/(dy^2);
tol = 1e-4;
error = 9e9;

%Assigning orginal values to T
T_initial = T;
T_old = T;


%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

iteration = 1;
    for k = 1:1400    
        error = 9e9;
        while(error > tol)
             for i = 2:nx - 1
                for j = 2:ny - 1
                    Term_1 = T_old(i,j);
                    Term_2 = k1*(T_old(i+1,j) - 2*T_old(i,j) + (T_old(i-1,j)));
                    Term_3 = k2*(T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1));
                    T(i,j) = Term_1 + Term_2 + Term_3;
                end
            
            end
        
            error = max(max(abs(T_old - T)));
            T_old = T;
        end
           iteration = iteration + 1;
    end
    
%Plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('unsteady explicit \n number of iterations = %d',iteration))
grid on

First, we have defined the mesh and initialized the Temperature with a guess value. Boundary conditions are defined at faces and the corner nodal value is calculated by averaging the adjacent values.

Then we defined k1 and k2 terms. error and tolerance values are assigned.

T_old variable is used for storing the previous iteration values.

later we have used the time loop and convergence loop inside it.

Plots are made using the contour command.

Graphs

It took 1401 iterations for solving the equation.

Implicit method

$\frac{T_{p}^{n + 1} - T_{p}}{Δ t} = α \cdot {(\frac{T_{L} - 2 T_{P} + T_{R}}{Δ x^{2}} + \frac{T_{T} - 2 T_{P} + T_{B}}{Δ y^{2}})}^{n + 1}$ $(T_p^(n+1)-T_p)/(Deltat) = alpha*((T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2))^( n+1)$

We only know $T_{p}^{n}$ $T_p^n$ and the rest are unknowns. So we use iterative solvers.

Implicit methods are unconditionally stable so even the larger time steps don't affect the solution much.

$T_{p}^{n + 1} = \frac{1}{1 + 2 \cdot k_{1} + 2 \cdot k_{2}} (T_{p}^{n} + k_{1} {(T_{L} + T_{R})}^{n + 1} + k_{2} {(T_{T} + T_{B})}^{n + 1})$ $T_p^(n+1) = 1/(1+2*k_1+2*k_2)(T_p^n+k_1(T_L+T_R)^(n+1)+k_2(T_T+T_B)^(n+1))$

In the Jacobi method, we have used old values to calculate the temperature at a particular time step.

In the Gauss Siedel method, we used the updates solution. Similarly in the SOR method, we have used updated values and relaxation factors.

$T_{S O R} = (1 - ω) T_{o l d} + ω (T_{G S})$ $T_(SOR) = (1-omega)T_(old)+omega(T_(GS))$

the code structure is similar to the Explicit method except for the formula for calculating temperature is iterative solvers.

T_prev variable is used for storing the values at the previous time step, unlike T_old which is for iteration.

By changing the variable iterative_method from 1,2 and 3 we can select different methods.

For omega between 1.1 and 1.9, the optimum value of 1.3 is chosen.

clear 
close all
clc

%code for solving 2d unsteady heat conduction equation
%domain we are considering is a square of 1 unit
nx = 10;
ny = nx; %number of nodes are same in the x and y directions.
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);%grid spacing in x and y directions
dt = 1e-2;
nt = 1400;
omega = 1.3; %relaxation factor

%boundary conditions 

%top = 600k bottom = 900k left = 400k and right = 800k
%initializing temperature in all the nodes except boundary nodes.
T = 300*ones(nx,ny);

T(1,:) = 600;%top face
T(nx,:) = 900;%bottom face
T(:,1) = 400;%left face
T(:,ny) = 800;%right face

%At the corner nodes we are taking average of two adjacent nodes
T(1,1) = (T(1,2)+T(2,1))/2;%top left corner
T(10,1) = (T(10,2)+T(9,1))/2;%bottom left corner
T(10,10) = (T(10,9)+T(9,10))/2;%bottom right corner
T(1,10) = (T(1,9)+T(2,10))/2;%top right corner

tol = 1e-4;
error = 9e9;
T_old = T;
T_prev = T;

k1 = 1.4 *dt/dx^2;
k2 = 1.4 *dt/dx^2;
   iter = 1;
   iteration_method = 2;
if iteration_method == 1 %1-jacobi ; 2-GS ; 3-SOR
for k = 1:nt
    error = 9e9; % define error inside time loop to reset the error after each space loop calculation
   while (error > tol)
               for i = 2:nx-1
                 for j = 2:ny-1
                    term1 = (1+(2*k1)+(2*k2))^-1;
                    term2 = k1*term1;
                    term3 = k2*term1;
                    H = (T_old(i-1,j)+T_old(i+1,j));
                    V = (T_old(i,j-1)+T_old(i,j+1));
                    T(i,j) = (T_prev(i,j)*term1)+(H*term2)+(V*term3);
                end
               end
                    error = max(max(abs(T_old-T)));
                    T_old = T;
                    iter = iter+1
   end
    T_prev = T;
    
end 
%Plotting
    figure(1)
    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca, 'YDir', 'reverse')
    title(sprintf('unsteady explicit jacobi iterations = %d',iter))
    xlabel('X-Axis')
    ylabel('Y-Axis')
elseif iteration_method == 2
    for k = 1:nt
    error = 9e9; % define error inside time loop to reset the error after each space loop calculation
   while (error > tol)
               for i = 2:nx-1
                 for j = 2:ny-1
                    term1 = (1+(2*k1)+(2*k2))^-1;
                    term2 = k1*term1;
                    term3 = k2*term1;
                    H = (T(i-1,j)+T_old(i+1,j));
                    V = (T(i,j-1)+T_old(i,j+1));
                    T(i,j) = (T_prev(i,j)*term1)+(H*term2)+(V*term3);
                end
               end
                    error = max(max(abs(T_old-T)));
                    T_old = T;
                    iter = iter+1;
   end
    T_prev = T;
    
    end
    %Plotting
    figure(1)
    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca, 'YDir', 'reverse')
    title('unsteday explicit GS')
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title(sprintf('Gauss sidel iteration number = %d',iter))
elseif iteration_method == 3
for k = 1:nt
     error = 9e9; % define error inside time loop to reset the error after each space loop calculation
   while (error > tol)
               for i = 2:nx-1
                 for j = 2:ny-1
                    term1 = (1+(2*k1)+(2*k2))^-1;
                    term2 = k1*term1;
                    term3 = k2*term1;
                    H = (T(i-1,j)+T_old(i+1,j));
                    V = (T(i,j-1)+T_old(i,j+1));
                    T(i,j) = ((1-omega)*T_old(i,j))+ omega*((T_prev(i,j)*term1)+(H*term2)+(V*term3));
                end
               end
                    error = max(max(abs(T_old-T)));
                    T_old = T;
                    iter = iter+1;
   end
    T_prev = T;
    
   
end
%Plotting
    figure(1)
    contourf(x,y,T)
    clabel(contourf(x,y,T))
    colorbar
    colormap(jet)
    set(gca, 'YDir', 'reverse')
    title(sprintf('unsteday explicit SOR iterations = %d',iter))
    xlabel('X-Axis')
    ylabel('Y-Axis')
end

The plotting of graphs is done outside the time loop for faster output.

Results

In both unsteady and steady solutions Jacobi iterations > gauss Sidel > SOR iterations.

SOR method has been efficient in obtaining the solution.

Steady-state solutions are obtained faster than the transient state, as the transient state has a time term in the equation.

We have solved the Transient state both explicitly and implicitly. Coding explicit was easier compared to implicit solution as the equation is simpler, but the CFL constraint has to be satisfied to obtain a stable solution.

In the implicit method, the scheme is stable unconditionally so the time step size doesn't affect the solution significantly.

Errors made

Used the old values instead of updated values in the GS method which I later corrected.

Defining error term outside time loop terminated iterations after just 12 iterations. The error term gets small after each iteration so when the tolerance limit is reached the convergence loop gets terminated. so Defined the error term inside the time loop.