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

This study is looking at modelling a 2D heat conduction problem employing three FDM methods (Jacobian, Gauss-Seidel and Gauss-Seidel with SOR) within steady-state, time-marching explicit and time-marching implicit formulations.

Below is the time-independent steady-state solution, written in MATLAB. This also shows what the rest of the solutions will look like.

clear all;
close all;
clc;

% Defining the domain

T = zeros(11);

T(:,1) = 400;
T(:,length(T)) =800;
T(1,:) = 600;
T(length(T),:) =900;

T(1,1)=500;
T(1,length(T))=700;
T(length(T),1)=650;
T(length(T),length(T))=850;
dx=1;
dy=1;
x=linspace(0, dx*length(T),length(T));
y=linspace(0, dy*length(T),length(T));

k=2*(dx^2 + dy^2)/(dx^2 * dy^2);
tolerance = 0.0001;
error = 10000000;
% the solving
% Jacobian
Tjac = T;
T_old = Tjac;
ct = 0;


while error > tolerance
    T_old = Tjac;
    for i = 2:length(T)-1
        for j = 2:length(T)-1
            Tjac(i,j) = T_Jac(T_old,i,j,k,dx,dy);
        end
    end
    error = max(max(abs(T_old - Tjac)));
    disp(error);
    ct = ct + 1;
end
figure
contourf(x,y,Tjac,'ShowText','on');
set(gca,'YDIR','reverse');
title("Jacobian");
colorbar
colormap(jet)

% Gauss-Seidel Method 
Tgs = T;
error2 = 1e6;
ct2=0;

while error2 > tolerance
    T_old = Tgs;
    for i = 2:length(T)-1
        for j = 2:length(T)-1
            Tgs(i,j) = T_Jac(Tgs,i,j,k,dx,dy);
        end
    end
    error2 = max(max(abs(T_old - Tgs)));
    disp(error2);
    ct2 = ct2 + 1;
end

figure
contourf(x,y,Tgs,'ShowText','on');
set(gca,'YDIR','reverse');
title("Gauss-Seidel");
colorbar
colormap(jet)

%SOR
error3 = 1e8;
ct3 = 0;
Tsor = T;
omega = 1.66;
while error3 > tolerance
    T_old = Tsor;
    
    for i = 2:length(T)-1
        for j = 2:length(T)-1
            Tsor(i,j) = T_Sor(Tsor,T_old,i,j,k,dx,dy,omega);
        end
    end
    error3 = max(max(abs(T_old - Tsor)));
    disp(error3);
    ct3 = ct3 + 1;
end
figure
contourf(x,y,Tsor,'ShowText','on');
set(gca,'YDIR','reverse');
title("Gauss-Seidel with SOR");
colorbar
colormap(jet)

% T_Jac Function referenced above
% function out = T_Jac(T,i,j,k,dx,dy)
%     out = 1/k * ( ((T(i-1,j) + T(i+1,j))/(dx^2)) + ((T(i,j-1) + T(i,j+1))/(dy^2)) );
% end

Running the code above uses the steady-state criterion to solve for the temperature field in an iterative fashion. Convergence was decided upon when all temperature values saw a change of less than 0.0001. Here it can be seen how the different FDM methods achieve convergence at widely varying rates. Jacobian method shows the slowest convergence, at 264 iterations required, followed by Gauss-Seidel, with 144 iterations. This shows a clear advantage of the Gauss-Seidel method as a computational method, being of roughly equivalent computational cost. When combined with SOR, convergence was achieved in only 44 iterations. Albeit requiring a few additional operations, it is fully justified by the far superior convergence speed.

Below is the explicit transient solution of the same problem:

clear all;
close all;
clc;

%Explicit 2D Heat solutions

dx=1;
dy=1;
nt=1000;
tstep=0.05;
alpha=0.99999;
k1=alpha*tstep/dx^2;
k2=alpha*tstep/dy^2;

T = ones(11);

x=linspace(0, dx*length(T),length(T));
y=linspace(0, dy*length(T),length(T));
T(:,1) = 400;
T(:,length(T)) =800;
T(1,:) = 600;
T(length(T),:) =900;

T(1,1)=500;
T(1,length(T))=700;
T(length(T),1)=650;
T(length(T),length(T))=850;

T_init = T;

%Jacobian

T_next = T;
for iter = 1:2000
    T = T_next;
    for i = 2:length(T)-1
        for j = 2:length(T)-1
            T_next(i,j) = T(i,j)+k1*(T(i-1,j)+T(i+1,j)-2*T(i,j))+k2*(T(i,j-1)+T(i,j+1)-2*T(i,j));
        end
    end
end
figure;
contourf(T_next,'ShowText','on')
colorbar
colormap(jet)
set(gca,'YDIR','reverse')

%Gauss-Seidel

T = T_init;
for iter = 1:2000
    t_old = T;
    for i = 2:length(T)-1
        for j = 2:length(T)-1
            T(i,j) = T(i,j)+k1*(T(i-1,j)+T(i+1,j)-2*T(i,j))+k2*(T(i,j-1)+T(i,j+1)-2*T(i,j));
        end
    end
end
figure;
contourf(T,'ShowText','on')
colorbar
colormap(jet)
set(gca,'YDIR','reverse')

%SOR

T=T_init;
Tsor = T;
Tgs=T_init;
omega = 1.1;

T = T_init;
t_old = T;
for iter = 1:2000
    for i = 2:length(T)-1
        for j = 2:length(T)-1
            T(i,j) = T(i,j)+k1*(T(i-1,j)+T(i+1,j)-2*T(i,j))+k2*(T(i,j-1)+T(i,j+1)-2*T(i,j));
        end
    end
    Tsor = T*omega +(1-omega)*t_old;
    t_old = Tsor;
end
figure;
contourf(Tsor,'ShowText','on')
colorbar
colormap(jet)
set(gca,'YDIR','reverse')

Running the code above solves for the temperature field at an arbitrarily set time step ahead of its current instance. Solutions are achieved almost innediately, albeit Jacobian method takes aboit 30% longer than the other methods. The speed of the explicit method can be attributed to its simplicity. However, it can be subject to divergence and cumulative error if initial conditions are defined wrongly. For example, setting an over-relaxation factor above 2 will cause the solution to diverge completely.

Below is the implicit transient model solving the same problem:

clear all;
close all;
clc;

dx=1;
dy=1;
nt=1000;
tstep=0.05;
T = zeros(11);
x=linspace(0, dx*length(T),length(T));
y=linspace(0, dy*length(T),length(T));

alpha=0.99999 ;
k1=alpha*tstep/dx^2;
k2=alpha*tstep/dy^2;
term1=(1+2*k1+2*k2)^-1;
term2=k1*term1;
term3=k2*term1;

T(:,1) = 400;
T(:,length(T)) =800;
T(1,:) = 600;
T(length(T),:) =900;

T(1,1)=500;
T(1,length(T))=700;
T(length(T),1)=650;
T(length(T),length(T))=850;

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

iter=0;
tic;
for z = 1:nt
    error=1e9;
    while(error>tol)
        for i = 2:length(T)-1
            for j = 2:length(T)-1
                T(i,j) = term1*T_prev(i,j) + term2*(T_old(i-1,j)+T_old(i+1,j)) + term3*(T_old(i,j-1)+T_old(i,j+1));
            end
        end
        error=max((max(abs(T_old-T))));  
        T_old=T;
        iter=iter+1;
    end
    T_prev=T;
end
t1=toc;

%GS

T_next = T_init;
error=1e9;
tol=1e-4;
T_old=T_init;
T_prev = T_init;

iter1=iter;
iter=0;
tic;
for z = 1:nt
    error=1e9;
    while(error>tol)
        for i = 2:length(T)-1
            for j = 2:length(T)-1
                T(i,j) = term1*T_prev(i,j) + term2*(T(i-1,j)+T(i+1,j)) + term3*(T(i,j-1)+T(i,j+1));            
            end
        end
        error=max((max(abs(T_old-T))));  
        T_old=T;
        iter=iter+1;
    end
    T_prev=T;
end
t2=toc;

%SOR-GS

T=T_init;
Tsor=T;
T_next = T_init;
error=1e9;
tol=1e-4;
T_old=T_init;
T_prev = T_init;
w=1.07;

iter2=iter;
iter=0;
tic;
for z = 1:nt
    error=1e9;
    while(error>tol)
        for i = 2:length(T)-1
            for j = 2:length(T)-1
                T(i,j) = T_old(i,j)*(1-w)+ w*(term1*T_prev(i,j) + term2*(T(i-1,j)+T_old(i+1,j)) + term3*(T(i,j-1)+T_old(i,j+1)));
            end
        end
        error=max(max(abs(T_old-T)));  
        T_old=T;
        iter=iter+1;
    end
    T_prev=T;
end
t3=toc;
iter3=iter;

The implicit method is inherently more computationally expensive than the explicit method, by virtue of taking the initial temperature field and calculating multiple temperature values (in our problem it outputs the whole temperature field) before moving on to the next iteration. Although the times are inconsistent, it takes about triple the time to compute the implicit solutions as it does the explicit solutions. This is because of the additional loop which ensures convergence at every time-step of the solution. This process is intended to substitute for the usual practice of solving for multiple unknown temperature values at once, with converging towards these values whilst still calculating each of them one by one - as per an explicit solution. What this means for the solution is that at every time-step it will be within a user-specified margin of error, which is larger in an explicit solution.