BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

BBalavijay V2 Dec 2021

Objective: The aim is to solve the 2D heat conduction equation by using the point iterative techniques.

Given:

Boundary conditions

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

Error criteria is 1e-4

Iteration methods to implement

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Steady state:

2D heat conduction equation for the steady state condition is

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

Apply central differencing scheme, we get

(Δx2 * Δy2)/2*(Δx2 + Δy2))*(((Ti-1,j + Ti+1,j)/Δx2) + (Ti,j-1 + Ti,j+1)/Δy2) = Ti,j

Using this equation, we can solve the problem by implement it in different iteration method.

Code for steady condition:

clear all
close all
clc

%Inputs
%Length of the domain
x = linspace(0,1,10);
y = linspace(0,1,10);

%Grid spacing
dx = x(2) - x(1);
dy = y(2) - y(1);

%Node division for the length of the domain
nx = 10;
ny = 10;

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

%Initial condition
T = 300*ones(nx,ny);

%Defining boundary condition
T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 600;
T(end,:) = 900;

absolute_error = 1e-4;
omega = 1.1;
K = (dx^2)*(dy^2)/(2*((dx^2)+(dy^2)));
Told = T;

% method = [1 2 3];
method=input('enter the number of method you want=');

if method==1 %Jacobi
    [T,ct,j_elapse] = J_steady(T,absolute_error,K,dx,dy,nx,ny,Told,xx,yy,method);
end

if method==2 %Gauss seidal
    [T,ct,GS_elapse] = GS_steady(Told,absolute_error,K,dx,dy,nx,ny,T,xx,yy,method);
end

if method==3 %SOR
    [T,ct,SOR_elapse] = SOR_steady(T,absolute_error,omega,nx,ny,xx,yy,Told,dx,dy,K,method);
end

%Jacobi 
function [T,ct,j_elapse] = J_steady(T,absolute_error,K,dx,dy,nx,ny,Told,xx,yy,method)
ct = 0;
error = 1;
tic
while error > absolute_error
    for i = 2:(nx-1)
        for j = 2:(ny-1)
            H=(Told(i-1,j) + Told(i+1,j))/(dx^2);
            V=(Told(i,j-1) + Told(i,j+1))/(dy^2);
            T(i,j) = (H+V)*K;
        end
    end
    error = max(max(abs(Told - T)));
    Told = T;
    ct = ct+1;
    j_elapse = toc;
end
    figure(1)
    contourf(xx,yy,T,'showText','on');
    colormap(jet)
    colorbar
    set(gca,'YDir','reverse')
    xlabel('length of x domain')
    ylabel('length of y domain')
    title({'steady state temperature',['method=' num2str(method)],['no. of iteration =' num2str(ct)]})
    pause(0.03)
end

%Gauss seidal
function [T,ct,GS_elapse] = GS_steady(Told,absolute_error,K,dx,dy,nx,ny,T,xx,yy,method)
error = 1;
ct = 0;
tic
while error > absolute_error
    for i = 2:(nx-1)
        for j = 2:(ny-1)
            H=(T(i-1,j) + Told(i+1,j))/(dx^2);
            V=(T(i,j-1) + Told(i,j+1))/(dy^2);
            T(i,j) = (H+V)*K;
        end
    end
    error = max(max((abs(Told - T))));
    Told = T;
    ct = ct+1;
    GS_elapse = toc;
end
    figure(2)
    contourf(xx,yy,T,'showText','on');
    colormap(jet)
    colorbar
    set(gca,'YDir','reverse')
    xlabel('length of x domain')
    ylabel('length of y domain')
    title({'steady state temperature',['method=' num2str(method)],['no. of iteration =' num2str(ct)]})
    pause(0.03)
end

%SOR
function [T,ct,SOR_elapse] = SOR_steady(T,absolute_error,omega,nx,ny,xx,yy,Told,dx,dy,K,method)
error = 1;
ct = 0;
tic
while error > absolute_error
    for i = 2:(nx-1)
        for j = 2:(ny-1)
            H=(T(i-1,j) + Told(i+1,j))/(dx^2);
            V=(T(i,j-1) + Told(i,j+1))/(dy^2);
            T(i,j) = (H+V)*K;
            T(i,j) = (Told(i,j)*(1-omega)) + (omega*T(i,j));
        end
    end
    error = max(max((abs(Told - T))));
    Told = T;
    ct = ct+1;
    SOR_elapse = toc;
end
    figure(3)
    contourf(xx,yy,T,'showText','on');
    colormap(jet)
    colorbar
    set(gca,'YDir','reverse')
    xlabel('length of x domain')
    ylabel('length of y domain')
    title({'steady state temperature',['method=' num2str(method)],['no. of iteration =' num2str(ct)]})
    pause(0.03)
end

The solution result shows how much the number of iteration differs for each methods to converge the temperature for the length of the domain x & y.

Simply, Jacobi > Gauss seidel > SOR.

Transient state:

Explicit solver:

`α*((∂^2T)/(∂x^2) + (∂^2T)/(∂y^2)) = (∂^2T)/(∂t^2)`

Apply central differencing scheme, we get

(Ti,j^n) + α*(Δ2t/Δx2)*( (Ti-1,j - 2 Ti,j + Ti+1,j)/Δx2)^n + α*(Δ2t/Δy2)*( (Ti,j-1 - 2 Ti,j + Ti,j+1)/Δy2)^n = (Ti,j^n+1)

Using this equation, we are solving for "n" term on LHS by time(t) to converge the temperature.

Code for transient explicit scheme:

clear all
close all
clc

%Inputs
%Length of the domain
x = linspace(0,1,10);
y = linspace(0,1,10);

%Grid spacing
dx = x(2) - x(1);
dy = y(2) - y(1);

%Node division for the length of the domain
nx = 10;
ny = 10;

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

%Initial condition
T = 300*ones(nx,ny);

%Defining boundary condition
T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 600;
T(end,:) = 900;

Told = T;
t = input('Time(in sec) =');
dt = 0.01;
nt = t/dt;
alpha = 1.4e-4;
CFL = alpha*dt*((1/dx^2)+(1/dy^2));
k1 = alpha*dt*(1/dx^2);
k2 = alpha*dt*(1/dy^2);

ct = 0;
tic
for k = 1:nt
    for i = 2:(nx-1)
        for j = 2:(ny-1)
            H = (Told(i-1,j))-(2*Told(i,j))+(Told(i+1,j));
            V = (Told(i,j-1))-(2*Told(i,j))+(Told(i,j+1));
            T(i,j) = Told(i,j) + (H*k1) + (V*k2);
        end
    end
    Told = T;
    ct = ct+1;
end
elapse_time = toc;

contourf(xx,yy,T,'showText','on');
colormap(jet)
colorbar
set(gca,'YDir','reverse')
xlabel('length of x domain')
ylabel('length of y domain')
title({'transient temperature(Explicit scheme)',['no. of iteration =' num2str(ct)],['time(in sec) =' num2str(t)]})

The image shows how temperature converge with increase of time

Implicit solver:

`α*((∂^2T)/(∂x^2) + (∂^2T)/(∂y^2)) = (∂^2T)/(∂t^2)`

Apply central differencing scheme, we get

(Ti,j^n) + α*(Δ2t/Δx2)*( (Ti-1,j - 2 Ti,j + Ti+1,j)/Δx2)^(n+1) + α*(Δ2t/Δy2)*( (Ti,j-1 - 2 Ti,j + Ti,j+1)/Δy2)^(n+1) = (Ti,j^n+1)

Using this equation, we are solving for "n+1" term on LHS for different iteration methods to converge the temperature.

Code for transient implicit scheme:

clear all
close all
clc

%Inputs
%Length of the domain
x = linspace(0,1,10);
y = linspace(0,1,10);

%Grid spacing
dx = x(2) - x(1);
dy = y(2) - y(1);

%Node division for the length of the domain
nx = 10;
ny = 10;

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

%Initial condition
T = 300*ones(nx,ny);

%Defining boundary condition
T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 600;
T(end,:) = 900;

Told = T;
Tprev = T;
t = input('Time(in sec) =');
dt = 10;
nt = t/dt;
alpha = 1.4;
absolute_error = 1e-4;
omega = 1.8;
CFL = alpha*dt*((1/dx^2)+(1/dy^2));
k1 = alpha*dt*(1/dx^2);
k2 = alpha*dt*(1/dy^2);

% method = [1 2 3];
method=input('enter the number of method you want=');

if method==1 %Jacobi
    [T,ct] = J_steady(T,absolute_error,k1,k2,nx,ny,Told,xx,yy,method,nt,Tprev,t);
end

if method==2 %Gauss seidal
    [T,ct] = GS_steady(Told,absolute_error,k1,k2,nt,nx,ny,T,xx,yy,method,Tprev,t);
end

if method==3 %SOR
    [T,ct] = SOR_steady(T,absolute_error,omega,nt,nx,ny,xx,yy,Told,k1,k2,method,Tprev,t);
end

%Jacobi 
function [T,ct] = J_steady(T,absolute_error,k1,k2,nx,ny,Told,xx,yy,method,nt,Tprev,t)
ct = 0;
error = 1;
for k = 1:nt
while error > absolute_error
    for i = 2:(nx-1)
        for j = 2:(ny-1)
            term1 = (1+(2*k1)+(2*k2))^-1;
            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) + (H*term2) +(V*term3);
        end
    end
    error = max(max(abs(Told - T)));
    Told = T;
    ct = ct+1;
end
Tprev = Told;
end
    contourf(xx,yy,T,'showText','on');
    colormap(jet)
    colorbar
    set(gca,'YDir','reverse')
    xlabel('length of x domain')
    ylabel('length of y domain')
    title({'transient temperature(Implicit scheme)',['method=' num2str(method)],['no. of iteration =' num2str(ct)],['time(in sec) =' num2str(t)]})
end

%Gauss seidal
function [T,ct] = GS_steady(Told,absolute_error,nt,k1,k2,nx,ny,T,xx,yy,method,Tprev,t)
error = 1;
ct = 0;
for k = 1:nt
while error > absolute_error
    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) + Told(i+1,j));
            V=(T(i,j-1) + Told(i,j+1));
            T(i,j) = (Tprev(i,j)*term1) + (H*term2) +(V*term3);
        end
    end
    error = max(max((abs(Told - T))));
    Told = T;
    ct = ct+1;
end
    Tprev = T;
end
    contourf(xx,yy,T,'showText','on');
    colormap(jet)
    colorbar
    set(gca,'YDir','reverse')
    xlabel('length of x domain')
    ylabel('length of y domain')
    title({'transient temperature(Implicit scheme)',['method=' num2str(method)],['no. of iteration =' num2str(ct)],['time(in sec) =' num2str(t)]})
end

%SOR
function [T,ct] = SOR_steady(T,absolute_error,omega,nt,nx,ny,xx,yy,Told,k1,k2,method,Tprev,t)
error = 1;
ct = 0;
for k = 1:nt
while error > absolute_error
    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) + Told(i+1,j));
            V=(T(i,j-1) + Told(i,j+1));
            T(i,j) = (Tprev(i,j)*term1) + (H*term2) +(V*term3);
            T(i,j) = (Told(i,j)*(1-omega)) + (omega*T(i,j));
        end
    end
    error = max(max((abs(Told - T))));
    Told = T;
    ct = ct+1;
end
    Tprev = T;
end
    contourf(xx,yy,T,'showText','on');
    colormap(jet)
    colorbar
    set(gca,'YDir','reverse')
    xlabel('length of x domain')
    ylabel('length of y domain')
    title({'transient temperature(Implicit scheme)',['method=' num2str(method)],['no. of iteration =' num2str(ct)],['time(in sec) =' num2str(t)]})
end

When CFL condition is less than or equal to 0.5, solution is steady. Otherwise, unsteady.

For t = 1000 sec, see how quickly temperature converges in SOR method than other two methods.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.