Solving the steady and unsteady 2D heat conduction problem

For the Explicit method for solving the Transient case of current work, the code are as follows:

%% MATLAB Program to generate temperature profile of 2D Transient State Heat Conduction Equation explicitly:
clear all;
close all;
clc;

%% Initializing Variables:
Nx = 51;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
residue = 1e-4;
t_start = 0;
t_end = 1/residue;
dt = 0.01;
error = 5e8;
tol = 1e-4;

%% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature


%% CFL Number calculation and adaptive Time-step control:
CFL_Number = residue*dt/(dx^2);

%% Calculation of Temperature Distribution explicitly using iterative solvers:

k1 = (residue*dt)/(dx^2);
k2 = (residue*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
initi = T;
T_prev_dt = initi;
tic;
counter = 1;
for z = 1:20000
    %% Gauss Jacobi explicit Iteration Scheme
    for i = 2:Nx-1
        for j = 2:Ny-1
            term1 = initi(i,j);
            term2 = k1*(initi(i-1,j)-2*initi(i,j)+initi(i+1,j));
            term3 = k2*(initi(i,j-1)-2*initi(i,j)+initi(i,j+1));
            T(i,j) = term1+term2+term3;
        end
    end
    initi = T;
    counter = counter+1;
end
time_counter = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf(\'2D Transient state heat conduction Explicit Method n Number of iterations: %d (Time: %0.2fs) n Computation Time: %0.4f n CFL Number: %0.4f\',counter-1,z*dt,time_counter,CFL_Number));
xlabel(\'xx\');
ylabel(\'yy\');

The output for the current Transient case is as follows:

This above method is for explicite approach.

In the Implicite approach, We go for various methods for sovling the equations as shown below:

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

The code for implicite approach for solving the approach for transient case.

%% MATLAB Program to solve the transient 2D Heat Conduction equation Implictly:
clear all;
close all;
clc;

%% Initializing Variables:
Nx = 51;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
Nt = [t_start:dt*50:t_end];
error = 5e8;
residue = 1e-4;

%% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% Calculation of optimal value of over-relaxation parameter (omega):
omega = 2/(1+sin(pi*dx));

%% Calculation of Temperature Distribution explicitly using iterative solvers:

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
counter = 1;
solver = input(\'Enter the number for the dedicated solver type: n 1 for Gauss Jacobi n 2 for Gauss Siedel n 3 for SOR n\');
tic;
for k = 1:20000
    %% Gauss Jacobi implicit Iteration Scheme:
    if solver == 1
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > residue)
            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) = term1*T_prev_dt(i,j)+term2*H+term3*V;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    if solver == 2
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > residue)
            for i = 2:Nx-1
                for j = 2:Ny-1
                    H = (T(i-1,j)+Told(i+1,j));
                    V = (T(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    if solver == 3
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > residue)
            for i = 2:Nx-1
                for j = 2:Ny-1
                    %H = (T(i-1,j)+Told(i+1,j));
                    %V = (T(i,j-1)+Told(i,j+1));
                    %T(i,j) = (omega)*(term1*T_prev_dt(i,j)+term2*H+term3*V) + (1-omega)*T_prev_dt(i,j);
                    H=(omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
                    V=(omega*T(i,j-1)+(1-omega)*Told(i,j-1)+Told(i,j+1));
                    T(i,j)=T_prev_dt(i,j)*term1+(term2*H+term3*V);
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    T_prev_dt = T;
end
execution_time = toc;
figure(1);
str_array = [\"Gauss Jacobi approach\",\"Gauss Siedel approach\",\"SOR approach\"];
formatSpec = \"Transient state 2D Heat Transfer by implicit: (%s Method) n Time: %0.3fs n Iteration Counter: %d n Computation time: %0.3fs\";

[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel(\'xx\');
ylabel(\'yy\');
title(sprintf(formatSpec,str_array(solver),k*dt,counter,execution_time));

The Code for Transient case is shown above. when runnning the code we get the input option for{ n1 for Gauss Jacobi Method, n2 for Gauss Siedel method and n3 for Successive over-relaxation method.

for n1 which is Gauss Jacobi method the output is as follows:

For n2 which is Gauss Siedel Method the output is as follows:

Similary for Successive over-relaxation method we have to give the input for n3:

Till now the code was for the the Transient State by implicite method.

For Steady State we solve the current work by Implicite method, which is also the further divided by following method:

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

The code for Steady State implicite method as shown below:

%% MATLAB Program to generate the 2D temperature distribution at steady state:
clear all;
close all;
clc;

%% Initial Conditions:
Nx = 81;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);

k = 2*(1/(dx^2)+1/(dy^2));
error = 5e8;
residue = 1e-4;

%% Initializing the temperature matrix:

T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% Calculation of over-relaxation parameter Omega:

omega = 2/(1+sin(pi*dx)); 
%% Solving the 2D steady state system iteratively:
[xx,yy] = meshgrid(x,y);
Initia = T;

solver = input(\'Enter the desired number: \')

tic;
%% Gauss Jacobi Iteration Scheme:
if solver == 1
    iter_count_jac = 1;
    while (error > residue)
        for i = 2:Nx-1
            for j = 2:Ny-1
                term1 = (Initia(i-1,j)+Initia(i+1,j))/(k*(dx^2));
                term2 = (Initia(i,j-1)+Initia(i,j+1))/(k*(dy^2));
                T(i,j) = term1+term2;
            end
        end
        error = max(max(abs(Initia-T)));
        Initia = T;
        iter_count_jac = iter_count_jac + 1;
    end
end
%% Gauss Siedel Iteration Scheme:
if solver == 2
    iter_count_gs = 1;
    while (error > residue)
        for i = 2:Nx-1
            for j = 2:Ny-1
                term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
                term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
                T(i,j) = term1+term2;
            end
        end
        error = max(max(abs(Initia-T)));
        Initia = T;
        iter_count_gs = iter_count_gs + 1;
    end
end
%% Successive Over Relaxation (SOR) Iteration Scheme:
if solver == 3
    iter_count_sor = 1;
    while (error > residue)
        for i = 2:Nx-1
            for j = 2:Ny-1
                term1 = (T(i-1,j)+Initia(i+1,j))/(k*(dx^2));
                term2 = (T(i,j-1)+Initia(i,j+1))/(k*(dy^2));
                 T(i,j) = omega*(term1+term2)-(omega-1)*Initia(i,j);
            end
        end
        error = max(max(abs(Initia-T)));
        Initia = T;
        iter_count_sor = iter_count_sor + 1;
    end
end
time_required = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if solver == 1
    title(sprintf(\'Temperature profile of steady state heat conduction using Gauss-Jacobi approach n No. of Iterations: %d n Computation Time: %0.3f s\',iter_count_jac,time_required));
else if solver == 2
        title(sprintf(\'Temperature profile of steady state heat conduction using Gauss-Siedel aproach n No. of Iterations: %d n Computation Time: %0.3f s\',iter_count_gs,time_required));
    else if solver == 3
            title(sprintf(\'Temperature profile of steady state heat conduction using SOR appraoch n No. of Iterations: %d n Computation Time: %0.3f s n Relaxation Parameter: %f\',iter_count_sor,time_required,real(omega)));
        end
    end
end
xlabel(\'xx\');
ylabel(\'yy\');

In the above code, if we type 1 in the commmand window we will get the the following countours for Steady state Solution solved by Gauss Jacobis method:

Similary for Steady State approach for Gauss Siedel method, we get the following output by typing the 2 in the command window:

At the end for Successive over-relaxation method for steady state the output will be following by typing 3 in command window:

This is the required project for steady and unsteady problem solved implicitly and explicitly. We have found the Explicite method is quite easy as compaired to implicite method. Because there is more number os steps for solving the problem in the implicite method. But at the same time at implicite method more accurate answer is found.