Talk to a counsellor
CoursesIIT ProgramsProjectsPlacementsFor BusinessBlogLog inTalk to a counsellor
Challenge

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

SSheik Allabakash28 Feb 2024
HTML

1. Steady State solver:

We know that the steady state equation for 2D heat conduction is

∂2T∂x2+∂2T∂y2=0">∂2T∂x2+∂2T∂y2=0∂2T∂x2+∂2T∂y2=0

By applying central differencing for the space derivatives and solving, we get

Ti,j=∆x2∆y22(∆x2+∆y2)⋅[Ti-1,j+Ti+1,j∆x2+Ti,j-1+Ti,j+1∆y2]">Ti,j=Δx2Δy22(Δx2+Δy2)⋅[Ti−1,j+Ti+1,jΔx2+Ti,j−1+Ti,j+1Δy2]Ti,j=∆x2∆y22(∆x2+∆y2)⋅[Ti-1,j+Ti+1,j∆x2+Ti,j-1+Ti,j+1∆y2]

The code solves the above equation using Jacobi, Gauss seidel and SOR methods and plots the temperature contour on X and Y axes.

clear all
close all
clc

%Parameters
L = 1;
nx = 100;
ny = nx;

x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);

%Guess temperature matrix
T = 300*ones(ny,nx);
T(:,1) = 400;
T(1,:) = 600;
T(:,end) = 800;
T(end,:) = 900;

T_old = T;

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

tol = 1e-4;
error = 1;

iteration_method = 1;%Iteration_method = 1-Jacobi ; 2-Gauss seidel ; 3-SOR

if iteration_method == 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_method == 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_method == 3 %SOR method
    omega = 1.8;
    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)
ax = gca;
ax.YDir = 'reverse';
xlabel('X-axis')
ylabel('Y-axis')
title({'Steady state temperature profile';['Iteration Method = ',num2str(iteration_method)];['No. of iterations = ',num2str(iter)];['Simulation time = ',num2str(time),' s']})
grid on

The plots show the temperature profiles of steady state condition using each of the methods and also additional information such as simulation time and no. of iterations performed.

Here, Iteration method = 1 =>Jacobi method; 2 =>Gauss seidel method; 3 =>SOR method

2. Transient state solvers:

a. Explicit solver

The transient equation of 2D conduction is

∂2T∂t2=α(∂2T∂x2+∂2T∂y2)">∂2T∂t2=α(∂2T∂x2+∂2T∂y2)∂2T∂t2=α(∂2T∂x2+∂2T∂y2)

Applying forward differencing to time derivative and central differencing for the space derivatives and solving, we get

Ti,jn+1=Ti,jn+α∆t∆x2(Ti-1,j-2Ti,j+Ti+1,j)n+α∆t∆y2(Ti,j-1-2Ti,j+Ti,j+1)n">Tn+1i,j=Tni,j+αΔtΔx2(Ti−1,j−2Ti,j+Ti+1,j)n+αΔtΔy2(Ti,j−1−2Ti,j+Ti,j+1)nTi,jn+1=Ti,jn+α∆t∆x2(Ti-1,j-2Ti,j+Ti+1,j)n+α∆t∆y2(Ti,j-1-2Ti,j+Ti,j+1)n

Since, we are solving explicitly, there are terms calculate at time 'n' in the RHS.

The below code uses a time marching solution to the above equation to solve it and lso plots the temperature.

clear all
close all
clc

%Parameters
L = 1;
nx = 100;
ny = nx;
dt = 0.14;
t = 1000;
nt = t/dt;
alpha = 1.75e-4;

x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);
CFL = alpha*dt*((1/dx^2)+(1/dy^2));

%Guess temperature matrix
T = 300*ones(ny,nx);
T(:,1) = 400;
T(1,:) = 600;
T(:,end) = 800;
T(end,:) = 900;
T_old = T;

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

iter = 0;
tic
%Time loop
for k = 1:nt
    %Y-Space loop
    for j = 2:ny-1
        %X-Space loop
        for i = 2:nx-1
            H = (T_old(i-1,j)-2*T_old(i,j)+T_old(i+1,j));
            V = (T_old(i,j-1)-2*T_old(i,j)+T_old(i,j+1));
            T(i,j) = T_old(i,j)+(H*k1)+(V*k2);
        end
    end
    T_old = T;
    iter = iter+1;
end
time = toc;

%Plotting temperature contour on X and Y axes
[c,h] = contourf(x,y,T);
clabel(c,h)
colorbar
colormap(jet)
ax = gca;
ax.YDir = 'reverse';
xlabel('X')
ylabel('Y')
title({'Transient state-Explicit temperature profile';['No. of iterations = ',num2str(iter),'    Time = ',num2str(t),' s'];['CFL = ',num2str(CFL),'     Simulation time = ',num2str(time)]})
grid on

The plots shows temperature profiles at time t=300s and 1000s :

  • The temperature profile changes as the time changes
  • As the temperature reaches 1000 s, the solution seems to achieve steady state conditions

b. Implicit solver

The transient equation of 2D conduction is

∂2T∂t2=α(∂2T∂x2+∂2T∂y2)">∂2T∂t2=α(∂2T∂x2+∂2T∂y2)∂2T∂t2=α(∂2T∂x2+∂2T∂y2)

Applying forward differencing to time derivative and central differencing for the space derivatives and solving, we get

Ti,jn+1=Ti,jn+α∆t(Ti-1,j-2Ti,j+Ti+1,j∆x2+Ti,j-1-2Ti,j+Ti,j+1∆y2)n+1">Tn+1i,j=Tni,j+αΔt(Ti−1,j−2Ti,j+Ti+1,jΔx2+Ti,j−1−2Ti,j+Ti,j+1Δy2)n+1Ti,jn+1=Ti,jn+α∆t(Ti-1,j-2Ti,j+Ti+1,j∆x2+Ti,j-1-2Ti,j+Ti,j+1∆y2)n+1

Let k1=α∆t∆x2">k1=αΔtΔx2k1=α∆t∆x2 & k2=α∆t∆y2">k2=αΔtΔy2k2=α∆t∆y2

We get `T_(i,j)^(n+1)=(T_(i,j)^n+k_1 (T_(i-1,j)+T_(i+1,j) )^(n+1)+k_2 (T_(i,j-1)+T_(i,j+1) )^(n+1))/(1+2k_1+2k_2 )

There are (n+1) terms in the RHS

The code solves the implicit equation by using Jacobi, Gauss seidel and SOR methods and plots the temperature.

clear all
close all
clc

%Parameters
L = 1;
nx = 50;
ny = nx;
dt = 10;
t = 1000;
nt = t/dt;
alpha = 1.75e-4;

x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);
CFL = alpha*dt*((1/dx^2)+(1/dy^2));

%Guess temperature matrix
T = 300*ones(ny,nx);
T(:,1) = 400;
T(1,:) = 600;
T(:,end) = 800;
T(end,:) = 900;
T_old = T;
T_prev = T;

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

tol = 1e-4;

iteration_method = 3;%Iteration_method = 1-Jacobi ; 2-Gauss seidel ; 3-SOR

%Jacobi method
if iteration_method == 1
    tic
    iter = 1;
    for k = 1:nt
        error = 1;
        while error > tol
            for j = 2:ny-1
                for i = 2:nx-1
                    term1 = 1/(1+2*k1+2*k2);
                    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-T_old)));
            T_old = T;
            iter = iter+1;
        end
        T_prev = T;
    end
    time = toc;
    
%Gauss Seidel method
elseif iteration_method == 2
    tic
    iter = 1;
    for k = 1:nt
        error = 1;
        while error > tol
            for j = 2:ny-1
                for i = 2:nx-1
                    term1 = 1/(1+2*k1+2*k2);
                    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-T_old)));
            T_old = T;
            iter = iter+1;
        end
        T_prev = T;
    end
    time = toc;
    
%SOR method
elseif iteration_method == 3
    omega = 1.1;
    tic
    iter = 1;
    for k = 1:nt
        error = 1;
        while error > tol
            for j = 2:ny-1
                for i = 2:nx-1
                    term1 = 1/(1+2*k1+2*k2);
                    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);
                    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
        T_prev = T; 
    end 
    time = toc;
end

%Plotting temperature contour on X and Y axes
[c,h] = contourf(x,y,T);
clabel(c,h)
colorbar
colormap(jet)
ax = gca;
ax.YDir = 'reverse';
xlabel('X')
ylabel('Y')
title({'Transient state-Implicit temperature profile';['Iteration method = ',num2str(iteration_method)];['No. of iterations = ',num2str(iter),'    Time = ',num2str(t),' s'];['CFL = ',num2str(CFL),'     Simulation time = ',num2str(time)]})
grid on

The plots obtained using each of the methods are shown below at time, t = 200s :

  • The solution obtained is similar to that of the explicit method
  • The no. of iterations decreased from jacobi to gauss seidel to SOR 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.