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

Steady and Unsteady 2D heat conduction

KKamal Patel22 Sept 2019
MATLAB

Objective: Implement iterative techinque to solve 2D heat conduction equation.

2D Heat Equation is of the form

`\\frac{\\partial^2T}{\\partialx^2} + \\frac{\\partial^2T}{\\partialy^2} = 0`

Discretizing the 2D PDE by central differencing method as follows

First with respect to x

`\\frac{\\partial^2T}{\\partialx^2} = \\frac{T_{(i-1, j)} - 2T_{(i, j)} + T_{(i+1, j)}}{\\Deltax^2}`

With respect to y

`\\frac{\\partial^2T}{\\partialx^2} = \\frac{T_{(i, j-1)} - 2T_{(i, j)} + T_{(i, j+1)}}{\\Deltay^2}`

Since nx = ny we can take `\\Deltax^2 = \\Deltay^2`

So the equation we get is

`T(i, j) = \\frac{1}{4}(T(i-1, j) + T(i+1, j) + T(i, j-1) + T(i, j+1))`

steady state analysis:

MATLAB CODE for

main.m ( main code)

clear all
close all
clc

L = 1;
n_x = 20;
n_y = 20;
x = linspace(0, L, n_x);
y = linspace(0, L, n_y);
dx = x(2) - x(1);
dy = y(2) - y(1);
error = 9e9;
tol = 1e-4;

%Boundary conditions
T_left = 400;
T_top = 600;
T_right = 800;
T_bottom = 900;
T = ones(n_x, n_y);

% temperature at boundary

T(1,1:n_x) = T_top;
T(n_y, 1:n_x) = T_bottom;
T(2:n_y-1, 1) = T_left;
T(2:n_y-1, n_x) = T_right;

told = T;

k1 = 1/4;
k2 = 1/4;

iterative_method = 1;
[X, Y] = meshgrid(x, y);

if iterative_method == 1
    jacobi_iter(T, x, y, n_x, n_y, tol, error, k1, k2, told);
end

if iterative_method == 2
    gauss_iter(T, x, y, n_x, n_y, tol, error, k1, k2, told);
end

if iterative_method == 3
    SOR_iter(T, x, y, n_x, n_y, tol, error, k1, k2, told);
end

jacobi_iter.m

function jacobi_iter(T, x, y, nx, ny, tol, error, k1, k2, told)
    
    %jacobi iteration
    jacobi_iter = 1;
    while(error > tol)
        
        for i=2:nx-1
            for j=2:ny-1
                T(i, j) = k1*(told(i-1, j)+told(i+1, j))+k2*(told(i, j-1)+told(i, j+1));
            end
        end
        error = max(max(abs(told - T)));
        told = T;
        jacobi_iter = jacobi_iter + 1;

        figure(1)
        contour(x, y, T);
        [a, b] = contourf(x, y, T);
        clabel(a, b);
        xlabel(\'X\')
        ylabel(\'Y\')
        colorbar

        title_text = sprintf(\'jacobi iteration = %d\', jacobi_iter);
        title(title_text)
    end
end

gauss_iter.m

function gauss_iter(T, x, y, nx, ny, tol, error, k1, k2, told)
    %Gauss Siedel Iteration
    gs_iter = 1;

    while(error > tol)

        for i=2:nx-1
            for j=2:ny-1
                T(i, j) = k1*(T(i, j-1)+told(i, j+1))+k2*(T(i-1, j)+told(i+1, j));
            end
        end
        error = max(max(abs(told - T)));
        told = T;
        gs_iter = gs_iter + 1;

        figure(2)
        contour(x, y, T);
        [a, b] = contourf(x, y, T);
        clabel(a, b);
        xlabel(\'X\')
        ylabel(\'Y\')
        colorbar

        title_text = sprintf(\'gauss iteration = %d\', gs_iter);
        title(title_text)
    end
end

SOR_iter.m

function SOR_iter(T, x, y, nx, ny, tol, error,k1, k2, told)
    
    %Successive over-relaxtion iteration
    sor_iter = 1;
    
    omega=2/(1+sin(2*pi/nx));
    while(error > tol)
        for i=2:nx-1
            for j=2:ny-1
                T(i, j) = omega*(k1*(T(i,j-1)+T(i,j+1))+k2*(T(i-1,j)+T(i+1,j)))+(1-omega)*told(i, j);
            end
        end
        error = max(max(abs(told-T)));
        told = T;
        sor_iter = sor_iter + 1;
        
        figure(3)
        contour(x, y, T);
        [a, b] = contourf(x, y, T);
        clabel(a, b);
        xlabel(\'X\')
        ylabel(\'Y\')
        colorbar
        title_text = sprintf(\'SOR iteration = %d\', sor_iter);
        title(title_text)
    end
    
end

Result

Unsteady State Analysis

General 2D equation is

`\\frac{\\partialT}{\\partial t} - \\alpha(\\frac{\\partial^2T}{\\partialx^2} + \\frac{\\partial^2T}{\\partialy^2}) = 0`

We can discretize the as follows

`\\frac{T_{i, j}^(n+1) - T_{i,j}^n}{\\Deltat} - \\alpha{\\frac{T_{(i-1, j)} - 2T_{(i, j)} + T_{(i+1, j)}}{\\Deltax^2}+\\frac{T_{(i, j-1)} - 2T_{(i, j)} + T_{(i, j+1)}}{\\Deltay^2}}`

On simplification we get

`T_{i,j}^(n+1) = T_{i,j}^n + \\alpha *Deltat*{\\frac{T_{(i-1, j)} - 2T_{(i, j)} + T_{(i+1, j)}}{\\Deltax^2}+\\frac{T_{(i, j-1)} - 2T_{(i, j)} + T_{(i, j+1)}}{\\Deltay^2}}`

`k1 = \\frac{\\alpha*Deltat}{\\Deltax^2}, k2 = \\frac{\\alpha*Deltat}{\\Deltay^2}`

MATLAB CODE:

Explicit method

clear all
close all
clc

L = 1;
n_x = 20;
n_y = 20;
x = linspace(0, L, n_x);
y = linspace(0, L, n_y);
dx = x(2) - x(1);
dy = y(2) - y(1);
error = 9e9;
tol = 1e-4;
dt = 30;
alpha = 1.88e-5;
%boundary condition
t_l = 400;
t_t = 600;
t_r = 800;
t_b = 900;

t = ones(n_x, n_y);
t(1,2:n_x) = t_t;
t(n_x, 2:n_x) = t_b;
t(2:n_y-1, 1) = t_l;
t(2:n_y-1, n_x) = t_r;

t(1, 1) = (t_t + t_l)/2;
t(n_x, n_y) = (t_r + t_b)/2;
t(1, n_y) = (t_t + t_r)/2;
t(n_x, 1) = (t_l + t_b)/2;


k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

[X, Y] = meshgrid(x, y);
told = t;
nt = 1000;
ctr = 1; %counter for number of iteration
for k = 1:nt
    for i=2:n_x-1
        for j=2:n_y-1
            t1 = told(i, j);
            t2 = k1*(told(i+1,j) - 2*told(i,j) + told(i-1, j));
            t3 = k1*(told(i,j+1) - 2*told(i,j) + told(i, j-1));
            t(i, j) = t1 + t2 + t3;
        end
    end
    error = max(max(abs(told - t)));
    told = t;
    ctr = ctr + 1;
end
figure(1)
[a, b] = contourf(X, Y, t);
clabel(a, b);
colormap(jet)
%contour(x, y, t);
xlabel(\'X\')
ylabel(\'Y\')
colorbar
title_text = sprintf(\'no of iteration = %d\', ctr);
title(title_text)

Result:

Implicit Method:

clear all
close all
clc

L = 1;
n_x = 20;
n_y = 20;
x = linspace(0, L, n_x);
y = linspace(0, L, n_y);
dx = x(2) - x(1);
dy = y(2) - y(1);
error = 9e9;
tol = 1e-4;
dt = 30;
alpha = 1.88e-5;

%boundary condition
t_l = 400;
t_t = 600;
t_r = 800;
t_b = 900;

t = ones(n_x, n_y);
t(1,2:n_x) = t_t;
t(n_x, 2:n_y-1) = t_b;
t(2:n_y-1, 1) = t_l;
t(2:n_y-1, n_x) = t_r;

t(1, 1) = (t_t + t_l)/2;
t(n_x, n_y) = (t_r + t_b)/2;
t(1, n_y) = (t_t + t_r)/2;
t(n_x, 1) = (t_l + t_b)/2;


k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

[X, Y] = meshgrid(x, y);
nt = 1000;

for i=1:3
    if i==1
        jacobi(t, x, y, n_x, n_y, tol, k1, k2, nt);
    end
    if i ==2
        gs_iteration(t, x, y, n_x, n_y, tol, k1, k2, nt);
    end
    if i==3
        SOR_iteration(t, x, y, n_x, n_y, tol, k1, k2, nt);
    end
end

%{
iterative_method =3;

if iterative_method==1
    jacobi(t, x, y, n_x, n_y, tol, k1, k2, nt);
end

if iterative_method==2
    gs_iteration(t, x, y, n_x, n_y, tol, k1, k2, nt);
end

if iterative_method==3
    SOR_iteration(t, x, y, n_x, n_y, tol, k1, k2, nt);
end
%}

Jacobi_iter.m

function jacobi(t, x, y, nx, ny, tol, k1, k2, nt)
    
    %jacobi iteration
    jacobi_iter = 1;
    told = t;
    t_prev_dt = told;
    term1 = (1+2*k1+2*k2)^(-1);
    term2 = k1*term1;
    term3 = k2*term1;
    
    for k = 1:nt %time-loop
        error = 9e9;
        while(error > tol) %convergence loop
            %nodal loop
            for i=2:nx-1
                for j=2:ny-1
                    H = (told(i-1,j)+told(i+1, j));
                    V = (told(i, j-1)+told(i, j+1));
                    t(i, j) = ((t_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
                end
            end
            error = max(max(abs(told - t)));
            told = t;
            jacobi_iter = jacobi_iter + 1;

            %figure(1)
            %contour(x, y, t);
            %[a, b] = contourf(x, y, t);
            %clabel(a, b);
            %xlabel(\'X\')
            %ylabel(\'Y\')
            %colorbar

            %title_text = sprintf(\'jacobi iteration = %d\', jacobi_iter);
            %title(title_text)
        end
        t_prev_dt = t;
    end
    figure(1)
    contour(x, y, t);
    [a, b] = contourf(x, y, t);
    clabel(a, b);
    xlabel(\'X\')
    ylabel(\'Y\')
    colorbar

    title_text = sprintf(\'jacobi iteration = %d\', jacobi_iter);
    title(title_text)
end

gauss_iter.m

function gs_iteration(t, x, y, nx, ny, tol, k1, k2, nt)
    
    %gs iteration
    gs = 1;
    told = t;
    t_prev_dt = told;
    term1 = (1+2*k1+2*k2)^(-1);
    term2 = k1*term1;
    term3 = k2*term1;
    
    for k = 1:nt %time-loop
        error = 9e9;
        while(error > tol) %convergence loop
            %nodal loop
            for i=2:nx-1
                for j=2:ny-1
                    H = (t(i-1,j)+t(i+1, j));
                    V = (t(i, j-1)+t(i, j+1));
                    t(i, j) = ((t_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
                end
            end
            error = max(max(abs(told - t)));
            told = t;
            gs = gs + 1;

            %figure(1)
            %contour(x, y, t);
            %[a, b] = contourf(x, y, t);
            %clabel(a, b);
            %xlabel(\'X\')
            %ylabel(\'Y\')
            %colorbar

            %title_text = sprintf(\'jacobi iteration = %d\', jacobi_iter);
            %title(title_text)
        end
        t_prev_dt = t;
    end
    figure(2)
    contour(x, y, t);
    [a, b] = contourf(x, y, t);
    clabel(a, b);
    xlabel(\'X\')
    ylabel(\'Y\')
    colorbar

    title_text = sprintf(\'gs iteration = %d\', gs);
    title(title_text)
end

SOR.m

function SOR_iteration(t, x, y, nx, ny, tol, k1, k2, nt)
    
    %SOR iteration
    SOR = 1;
    told = t;
    t_prev_dt = told;
    term1 = (1+2*k1+2*k2)^(-1);
    term2 = k1*term1;
    term3 = k2*term1;
    omega=2/(1+sin(2*pi/nx));
    for k = 1:nt %time-loop
        error = 9e9;
        while(error > tol) %convergence loop
            %nodal loop
            for i=2:nx-1
                for j=2:ny-1
                    H = (t(i-1,j)+t(i+1, j));
                    V = (t(i, j-1)+t(i, j+1));
                    t(i, j) = omega * ((t_prev_dt(i,j)*term1)+(H*term2)+(V*term3))-(omega-1)*t_prev_dt(i, j);
                end
            end
            error = max(max(abs(told - t)));
            told = t;
            SOR = SOR + 1;

            %figure(1)
            %contour(x, y, t);
            %[a, b] = contourf(x, y, t);
            %clabel(a, b);
            %xlabel(\'X\')
            %ylabel(\'Y\')
            %colorbar

            %title_text = sprintf(\'jacobi iteration = %d\', jacobi_iter);
            %title(title_text)
        end
        t_prev_dt = t;
    end
    figure(3)
    contour(x, y, t);
    [a, b] = contourf(x, y, t);
    clabel(a, b);
    xlabel(\'X\')
    ylabel(\'Y\')
    colorbar

    title_text = sprintf(\'SOR iteration = %d\', SOR);
    title(title_text)
end

Result:

From the plot we can observe that iterative solver showed similar trend when solved implicitly. Compared to all the solvers successive over relaxation solver converged faster compared to other, with jacobi converging the slowest.

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.