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

MMohammed Fares Batesh4 Jun 2021

First defining the boundary conditions, the space nodes size and important study parameters

Then define “if else” to switch the method regarding iterative solver identifier, after that a “while” loop to get the error to the required accuracy. The space will be covered by two “for” loop for x and y axis. Inside this loop the governing equation will be handled. The error, the old temperature and the plotting is updated inside the “while” loop after the “for” loop. The plot title is the method state status iteration count number.

The governing equation for each method and the plotted result as following.

dx = 1/nx;

dy = 1/nx;

alpha = 1.2;

k = 2*((dx^2+dy^2)/(dx^2*dy^2));

k1 = alpha*dt/dx^2;

k2 = alpha*dt/dy^2;

w = 1.2;

Steady state Jacobi

term1 = (Told(i-1,j)+Told(i+1,j));

term2 = (Told(i,j-1)+Told(i,j+1));

T(i,j) = (term1/dx^2 + term2/dy^2)/k;

Steady state gauss-seidel

term1 = (T(i-1,j)+Told(i+1,j));

term2 = (T(i,j-1)+Told(i,j+1));

T(i,j) = (term1/dx^2 + term2/dy^2)/k;

Steady state successive over-relaxation

term1 = (T(i-1,j)+Told(i+1,j));

term2 = (T(i,j-1)+Told(i,j+1));

T(i,j) = (1-w)*Told(i,j) + w*(term1/dx^2 + term2/dy^2)/k;

Transient state explicit Jacobi

term1 = (T_prev_dt(i-1,j)-2*T_prev_dt(i,j)+T_prev_dt(i+1,j));

term2 = (T_prev_dt(i,j-1)-2*T_prev_dt(i,j)+T_prev_dt(i,j+1));

T(i,j) = T_prev_dt(i,j)+k1*term1 + k2*term2;

Transient state explicit gauss-seidel

term1 = (T(i-1,j)-2*T_prev_dt(i,j)+T_prev_dt(i+1,j));

term2 = (T(i,j-1)-2*T_prev_dt(i,j)+T_prev_dt(i,j+1));

T(i,j) = T_prev_dt(i,j)+k1*term1 + k2*term2;

Transient state explicit successive over-relaxation

term1 = (T(i-1,j)-2*T_prev_dt(i,j)+T_prev_dt(i+1,j));

term2 = (T(i,j-1)-2*T_prev_dt(i,j)+T_prev_dt(i,j+1));

T(i,j) = (1-w)*T_prev_dt(i,j) + w*(T_prev_dt(i,j)+k1*term1 + k2*term2);

Transient state implicit Jacobi

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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);

Transient state implicit gauss-seidel

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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);

Transient state implicit successive over-relaxation

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) = (1-w)*Told(i,j) + w*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));

The conclusion: the successive over relaxation is the fastest method to get the final solution then the gauss-seidel and finally the jacobi method is the slowest method to get the required solution.

clear all
clear all
clc

% Inputs
%boundary conditions
nx = 20;
T = 650*ones(nx,nx);
Ttop = 600;
Tbottom = 900;
Tleft = 400;
Tright = 800;
T(1,:) = Ttop;
T(end,:) = Tbottom;
T(:,1) = Tleft;
T(:,end) = Tright;
T(1,1) = 500;
T(1,end) = 700;
T(end,1) = 650;
T(end,end) = 850;
Told = T;
T_prev_dt = T;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
[X,Y] = meshgrid(x,y);
% Parameter
dx = 1/nx;
dy = 1/nx;
dt = 5e-4; 
alpha = 1.2;
k = 2*((dx^2+dy^2)/(dx^2*dy^2));
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
CFL = k1 +k2;
iterative_solver = 9;
error = 1e-1;
tol = 1e-4;
w_best = 1.7;
w_best1 = 1.8;
w = 1.2;
% videofile = VideoWriter('forward_kinematica.avi','Uncompressed AVI');
% open(videofile)
if iterative_solver == 1 % Jacobi method
    jacobi_iter = 1;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:nx-1
                term1 = (Told(i-1,j)+Told(i+1,j));
                term2 = (Told(i,j-1)+Told(i,j+1));
                T(i,j) =  (term1/dx^2 + term2/dy^2)/k;
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        jacobi_iter = jacobi_iter +1;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Jacobi Iteration number = %d', jacobi_iter);
        title(title_text)
        pause(0.03)
    end
elseif iterative_solver == 2 % Gauss-seidel method
    gs_iter = 1;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:nx-1
                term1 = (T(i-1,j)+Told(i+1,j));
                term2 = (T(i,j-1)+Told(i,j+1));
                T(i,j) =  (term1/dx^2 + term2/dy^2)/k;
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        gs_iter = gs_iter +1;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Gauss-seidel Iteration number = %d', gs_iter);
        title(title_text)
        pause(0.03)
    end
    elseif iterative_solver == 3 % Successive over-relaxation
    sor_iter = 1;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:nx-1
                term1 = (T(i-1,j)+Told(i+1,j));
                term2 = (T(i,j-1)+Told(i,j+1));
                T(i,j) = (1-w)*Told(i,j) + w*(term1/dx^2 + term2/dy^2)/k;
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        sor_iter = sor_iter +1;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Successive over-relaxation Iteration number = %d', sor_iter);
        title(title_text)
        pause(0.03)
    end
elseif iterative_solver == 4  %Jacobi
    jacobi_explicit_iter = 1;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:nx-1
                term1 = (T_prev_dt(i-1,j)-2*T_prev_dt(i,j)+T_prev_dt(i+1,j));
                term2 = (T_prev_dt(i,j-1)-2*T_prev_dt(i,j)+T_prev_dt(i,j+1));
                T(i,j) = T_prev_dt(i,j)+k1*term1 + k2*term2;                
            end
        end
        error = max(max(abs(T_prev_dt- T)));
        T_prev_dt = T;
        jacobi_explicit_iter = jacobi_explicit_iter +1;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Jacobi Explicit Iteration number = %d', jacobi_explicit_iter);
        title(title_text)
        pause(0.03)
    end
    elseif iterative_solver == 5  % Gauss-seidel method
    gs_explicit_iter = 1;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:nx-1
                term1 = (T(i-1,j)-2*T_prev_dt(i,j)+T_prev_dt(i+1,j));
                term2 = (T(i,j-1)-2*T_prev_dt(i,j)+T_prev_dt(i,j+1));
                T(i,j) = T_prev_dt(i,j)+k1*term1 + k2*term2;                
            end
        end
        error = max(max(abs(T_prev_dt- T)));
        T_prev_dt = T;
        gs_explicit_iter = gs_explicit_iter +1;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Gauss-seidel Explicit Iteration number = %d', gs_explicit_iter);
        title(title_text)
        pause(0.03)
    end
    elseif iterative_solver == 6  % Successive over-relaxation
    sor_explicit_iter = 1;
    while (error > tol)
        for i = 2:nx-1
            for j = 2:nx-1
                term1 = (T(i-1,j)-2*T_prev_dt(i,j)+T_prev_dt(i+1,j));
                term2 = (T(i,j-1)-2*T_prev_dt(i,j)+T_prev_dt(i,j+1));
                T(i,j) = (1-w)*T_prev_dt(i,j) + w*(T_prev_dt(i,j)+k1*term1 + k2*term2);                
            end
        end
        error = max(max(abs(T_prev_dt- T)));
        T_prev_dt = T;
        sor_explicit_iter = sor_explicit_iter +1;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Successive over-relaxation Explicit Iteration number = %d', sor_explicit_iter);
        title(title_text)
        pause(0.03)
    end
    elseif iterative_solver == 7  %Jacobi
    jacobi_implicit_iter = 1;
    for tn = 1:400
        error = 1e-1;
        while (error > tol)
            for i = 2:nx-1
                for j = 2:nx-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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);                
                end
            end
            error = max(max(abs(Told- T)));
            Told = T;
            jacobi_implicit_iter = jacobi_implicit_iter +1;
        end
        T_prev_dt = T;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Jacobi Implicit Iteration number = %d  Time = %d', jacobi_implicit_iter, tn);
        title(title_text)
        pause(0.03)
    end
        elseif iterative_solver == 8  % Gauss-seidel
    jacobi_implicit_iter = 1;
    for tn = 1:400
        error = 1e-1;
        while (error > tol)
            for i = 2:nx-1
                for j = 2:nx-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) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);                
                end
            end
            error = max(max(abs(Told- T)));
            Told = T;
            jacobi_implicit_iter = jacobi_implicit_iter +1;
        end
        T_prev_dt = T;
        contourf(X,Y,T);
        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Gauss-seidel Implicit Iteration number = %d  Time = %d', jacobi_implicit_iter, tn);
        title(title_text)
        pause(0.03)
    end
        elseif iterative_solver == 9  %Jacobi
    jacobi_implicit_iter = 1;
    for tn = 1:400
        error = 1e-1;
        while (error > tol)
            for i = 2:nx-1
                for j = 2:nx-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) = (1-w)*Told(i,j) + w*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));                
                end
            end
            error = max(max(abs(Told- T)));
            Told = T;
            jacobi_implicit_iter = jacobi_implicit_iter +1;
        end
        T_prev_dt = T;
        contourf(X,Y,T);

        xlabel('X axis')
        ylabel('Y axis')
        title_text = sprintf('Successive over-relaxation Implicit Iteration number = %d  Time = %d', jacobi_implicit_iter, tn);
        title(title_text)
        pause(0.03)
    end
end
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.