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

AIM:

To solve the 2D heat conduction equation by using the point iterative technique using the methods like Jacobi, Gauss-Seidel, Successive over relaxation for both implicit and explicit schemes.

Theory:

Inputs and Boundary conditions:

• Domain is assumed to be a shape of a square of length 1m.
• of nodes in x and y direction be 30. (n = 30)
• Temperature at nodes be
  • Top Boundary = 600 K
  • Bottom Boundary = 900 K
  • Left Boundary = 400 K
  • Right Boundary = 800 K
  • Initial temperature inside the domain = 298 K
  • Top - Left corner = 500 K
  • Top - Right corner = 700 K
  • Bottom - Left corner = 650 K
  • Bottom - Right corner = 850 K
  • Absolute error criteria = 1e-4

 

Steady State heat conduction equation is given as,

Discretizing the above equation using the Central difference scheme (CDS), we get

Further simplifying the above equation, we get

Jacobi Method:

Gauss-Seidel Method:

Successive over Relaxation Method:

 

Unsteady State heat conduction equation is given as,

Explicit Approach:

 

Discretizing the above time term using the forward differencing scheme (FDS) and for space term using central difference scheme (CDS).

For explicit diffusion problem, CFL_x + CFL_y  <= 0.5 for solution to be stable.

Implicit Approach:

Jacobi Method:

Gauss-Seidel Method:

Successive over Relaxation Method:

[MATLAB CODE]

STEADY STATE:

clc
clear all
close all


n = 30;
x = linspace(0,1,n);
y = linspace(0,1,n);

dx = x(2) - x(1);
dy = dx;

T = 298*ones(n,n);
T(1,:) = 600;       % Top Boundary Temp. in Kelvin
T(end,:) = 900;     % Bottom Boundary Temp. in Kelvin
T(:,1) = 400;       % Left Boundary Temp. in Kelvin
T(:,end) = 800;     % Right Boundary Temp. in Kelvin

T(1,1) = 500;       % Top - Left corner Temp. in Kelvin
T(1,end) = 700;     % Top - Right corner Temp. in Kelvin
T(end,1) = 650;     % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850;   % Bottom - Right corner Temp. in Kelvin

Told1 = T; 
Told2 = T;
Told3 = T;

T1 = T;
T2 = T;
T3 = T;

jacobi_error = 1;
gauss_sidel_error = 1;
SOR_error = 1;
error_req = 1e-4;

jacobi_iter = 0;
gauss_sidel_iter = 0;
SOR_iter = 0;

w = 2/(1+sin(pi*dx)); % Relaxation parameter for SOR Method

for method = 1:3
    tic
    if method == 1
        while jacobi_error > error_req
            for i = 2:(n-1)
                for j = 2:(n-1)
                    T1(i,j) = 0.25*(Told1(i+1,j)+Told1(i-1,j)+Told1(i,j+1)+Told1(i,j-1));
                end
            end
            jacobi_iter = jacobi_iter + 1;
            jacobi_error = max(max(abs(Told1-T1)));
            Told1 = T1;
            time_jc = toc;
        end
    end
    
    figure(1)
    contourf(x,y,T1,'Showtext','on')
    colormap(jet)
    colorbar
    xlabel('X Axis')
    ylabel('Y Axis')
    title(['After' num2str(jacobi_iter) 'Jacobi'])
    pause(0.3)
    
    tic
    if method == 2
        while gauss_sidel_error > error_req
            for i = 2:(n-1)
                for j = 2:(n-1)
                    T2(i,j) = 0.25*(Told2(i+1,j)+T2(i-1,j)+Told2(i,j+1)+T2(i,j-1));
                end
            end
            gauss_sidel_iter = gauss_sidel_iter + 1;
            gauss_sidel_error = max(max(abs(Told2-T2)));
            Told2 = T2;
            time_gs = toc;
        end
    end
    
    figure(2)
    contourf(x,y,T2,'Showtext','on')
    colormap(jet)
    colorbar
    xlabel('X Axis')
    ylabel('Y Axis')
    title(['After'  num2str(gauss_sidel_iter)  'Gauss Sidel'])
    pause(0.3)
    
    tic
    if method == 3
        while SOR_error > error_req
            for i = 2:(n-1)
                for j = 2:(n-1)
                    T3(i,j) = (1-w)*Told3(i,j) + w*(0.25*(Told3(i+1,j)+T3(i-1,j)+Told3(i,j+1)+T3(i,j-1)));
                end
            end
            SOR_iter = SOR_iter + 1;
            SOR_error = max(max(abs(Told3-T3)));
            Told3 = T3;
            time_sor = toc;
        end
    end
    
    figure(3)
    contourf(x,y,T3,'Showtext','on')
    colormap(jet)
    colorbar
    xlabel('X Axis')
    ylabel('Y Axis')
    title(['After' num2str(SOR_iter) 'SOR'])
    pause(0.3)
    
end

 

UNSTEADY STATE EXPLICIT:

[MATLAB CODE]

clc
clear all
close all

%% 
n = 30;
L = 1;
x = linspace(0,L,n);
y = linspace(0,L,n);
%% 
dx = x(2) - x(1);
dy = dx;
dt = 2e-4;
alpha = 1.2;
%%
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
CFL = k1+k2;

%% 
T = 298*ones(n,n);
T(1,:) = 600;       % Top Boundary Temp. in Kelvin
T(end,:) = 900;     % Bottom Boundary Temp. in Kelvin
T(:,1) = 400;       % Left Boundary Temp. in Kelvin
T(:,end) = 800;     % Right Boundary Temp. in Kelvin

T(1,1) = 500;       % Top - Left corner Temp. in Kelvin
T(1,end) = 700;     % Top - Right corner Temp. in Kelvin
T(end,1) = 650;     % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850;   % Bottom - Right corner Temp. in Kelvin

Told = T;
T1 = T;

error_req = 1e-4;
error_mag = 1;
iter = 0;

tic
for nt = 1:1000
        for i = 2:(n-1)
            for j = 2:(n-1)
                term1 = Told(i,j);
                term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
                term3 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
                T1(i,j) = term1 + term2 + term3;
            end
        end
        iter = iter+1;
        error_mag = max(max(abs(Told - T1)));
        Told = T1;
        time_req = toc;
end
  
figure(1)
contourf(x,y,T1,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['Explicit Method at the end of total'  num2str( iter) 'Iteration'])

 

UNSTEADY STATE IMPLICIT JACOBI:

[MATLAB CODE]

clc
clear all
close all

n = 30;
l = 1;
x = linspace(0,l,n);
y = linspace(0,l,n);

dx = x(2) - x(1);
dy = dx;
nt = 1000;
dt = 0.02;
alpha = 1.2e3;
%%
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
CFL = k1+k2;

T = 298*ones(n,n);
T(1,:) = 600;       % Top Boundary Temp. in Kelvin
T(end,:) = 900;     % Bottom Boundary Temp. in Kelvin
T(:,1) = 400;       % Left Boundary Temp. in Kelvin
T(:,end) = 800;     % Right Boundary Temp. in Kelvin

T(1,1) = 500;       % Top - Left corner Temp. in Kelvin
T(1,end) = 700;     % Top - Right corner Temp. in Kelvin
T(end,1) = 650;     % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850;   % Bottom - Right corner Temp. in Kelvin

jacobi_iter = 0;
jacobi_error = 9e9;
tol = 1e-4;

Told = T;
T_previous = Told;

term1 = (1+2*k1+2*k2)^(-1);
term2 = (k1*term1);
term3 = (k2*term1);

tic
for k = 1:nt
    while jacobi_error > tol
        for i = 2:(n-1)
            for j = 2:(n-1)
                H = Told(i-1,j)+Told(i+1,j);
                V = Told(i,j-1)+Told(i,j+1);
                T(i,j) = (T_previous(i,j)*term1)+(H*term2)+(V*term3);                    
            end
        end
        jacobi_iter = jacobi_iter + 1;
        jacobi_error = max(max(abs(Told-T)));
        Told = T;
        time_jc = toc;
    end
    T_previous = T;
end

figure(1)
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['Implicit Method at the end of total'  num2str( jacobi_iter) 'Iteration'])

 

UNSTEADY STATE IMPLICIT GAUSS SEIDEL:

[MATLAB CODE]

clc
clear all
close all

n = 30;
l = 1;
x = linspace(0,l,n);
y = linspace(0,l,n);

dx = x(2) - x(1);
dy = dx;
nt = 1000;
dt = 0.02;
alpha = 1.2e3;
%%
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
CFL = k1+k2;

T = 298*ones(n,n);
T(1,:) = 600;       % Top Boundary Temp. in Kelvin
T(end,:) = 900;     % Bottom Boundary Temp. in Kelvin
T(:,1) = 400;       % Left Boundary Temp. in Kelvin
T(:,end) = 800;     % Right Boundary Temp. in Kelvin

T(1,1) = 500;       % Top - Left corner Temp. in Kelvin
T(1,end) = 700;     % Top - Right corner Temp. in Kelvin
T(end,1) = 650;     % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850;   % Bottom - Right corner Temp. in Kelvin

gauss_iter = 0;
gauss_error = 9e9;
tol = 1e-4;

Told = T;
T_previous = Told;

term1 = (1+2*k1+2*k2)^(-1);
term2 = (k1*term1);
term3 = (k2*term1);

tic
for k = 1:nt
    while gauss_error > tol
        for i = 2:(n-1)
            for j = 2:(n-1)
                H = T(i-1,j)+T(i+1,j);
                V = T(i,j-1)+T(i,j+1);
                T(i,j) = (T_previous(i,j)*term1)+(H*term2)+(V*term3);                    
            end
        end
        gauss_iter = gauss_iter + 1;
        gauss_error = max(max(abs(Told-T)));
        Told = T;
        time_gs = toc;
    end
    T_previous = Told;
end

figure(1)
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['Implicit Method at the end of total'  num2str( gauss_iter) 'Iteration'])

 

UNSTEADY STATE IMPLICIT SOR:

[MATLAB CODE]

clc
clear all
close all

n = 30;
l = 1;
x = linspace(0,l,n);
y = linspace(0,l,n);

dx = x(2) - x(1);
dy = dx;
nt = 1000;
dt = 0.02;
alpha = 1.2e3;
%%
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
CFL = k1+k2;

T = 298*ones(n,n);
T(1,:) = 600;       % Top Boundary Temp. in Kelvin
T(end,:) = 900;     % Bottom Boundary Temp. in Kelvin
T(:,1) = 400;       % Left Boundary Temp. in Kelvin
T(:,end) = 800;     % Right Boundary Temp. in Kelvin

T(1,1) = 500;       % Top - Left corner Temp. in Kelvin
T(1,end) = 700;     % Top - Right corner Temp. in Kelvin
T(end,1) = 650;     % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850;   % Bottom - Right corner Temp. in Kelvin

SOR_iter = 0;
SOR_error = 9e9;
tol = 1e-4;

Told = T;
T_previous = Told;

term1 = (1+2*k1+2*k2)^(-1);
term2 = (k1*term1);
term3 = (k2*term1);

w = 2/(1+sin(pi*dx)); % Relaxation parameter for SOR Method
%w = 1.2;
%%

tic
for k = 1:nt
    while SOR_error > tol
        for i = 2:(n-1)
            for j = 2:(n-1)
                H = T(i-1,j)+T(i+1,j);
                V = T(i,j-1)+T(i,j+1);
                T(i,j) = w*((T_previous(i,j)*term1)+(H*term2)+(V*term3)) - (w-1)*T(i,j);    
            end
        end
        SOR_iter = SOR_iter + 1;
        SOR_error = max(max(abs(Told-T)));
        Told = T;
        time_sor = toc;
    end
    T_previous = T;
end

figure(1)
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['Implicit Method at the end of total'  num2str( SOR_iter) 'Iteration'])

 

EXPLANATION:

• Initially, the number of nodes, boundary conditions (the temperature at side and edges) are assigned.
• x, dx, y, dy, and dt values are assigned.
• The constants ‘k’ (for steady-state), alpha, k1, k2 (CFL number for transient analysis) are assigned.
•  For Steady-state:
• 1. The steady-state solution is performed by using the iterative solvers (Jacobi, Gauss-Seidel, SOR) for which the initial errors and the tolerance error value is assigned and the value of omega ‘w’ is calculated for the SOR technique.
  2. The Jacobi, Gauss seidel, and SOR techniques are implemented under a for loop which uses the ‘for method = 1:3’ for Jacobi, ‘for method = 2:3’ for Gauss seidel, and ‘for method = 3:3’ for SOR.
  3. Under this convergence loop (while loop) is executed, in which T(i,j) is calculated under for loop, and after the for loop, the error value and iteration values are updated and the value of error is calculated.
  4. A ‘Tic’ and ‘Toc’ command is used to calculate the simulation run time.
  5. While loop comes to an end when error falls below the tolerance value and T(i,j) is plotted using ‘contour’ and ‘colormap’
  6. The procedure from 3 to 5 is repeated for Gauss-Seidel and SOR method.

• For Unsteady state Explicit method:
  • In the explicit method iterative solver (convergence loop) is not used.
  • Time loop (for nt = 1:1000) is used. After the loop ends, the temperature T(i,j) value is found after which the error and old temperature value is updated.
  • The T(i,j) is plotted at the end of the loop.

• For Unsteady state Implicit method:
  1. In the implicit method, the Told, T_previous, and iteration numbers are assigned.
  2. Under this convergence loop (while loop) is executed, in which T(i,j) is calculated under for loop, and after the for loop, the error value and iteration values are updated and the value of error is calculated.
  3. A ‘Tic’ and ‘Toc’ command is used to calculate the simulation run time.
  4. After each time loop, the T_previous value is updated.
  5. While loop comes to an end when error falls below the tolerance value and T(i,j) is plotted using ‘contour’ and ‘colormap’
  6. The procedure from 1 to 5 is repeated for Gauss-Seidel and SOR method.

 

OUTPUT:

PLOTS:

STEADY STATE:

STEADY STATE SIMULATION TIME:

 

UNSTEADY STATE EXPLICIT:

UNSTEADY STATE EXPLICIT SIMULATION TIME:

 

UNSTEADY STATE IMPLICIT:

• Unsteady state implicit Jacobi:

UNSTEADY STATE IMPLICIT SIMULATION TIME JACOBI:

 

• Unsteady state implicit Gauss – seidel:

UNSTEADY STATE IMPLICIT SIMULATION TIME GAUSS-SEIDEL:

• Unsteady state implicit SOR:

UNSTEADY STATE IMPLICIT SIMULATION TIME SOR:

 

 

EXPLANATION:

• In each plot, the domain is considered as the square of side length 1m.
• Each color in the plot denotes the specific temperature values.

Sr. No.	Method	Iterations	Errors
1.	Steady State (Jacobi)	1784	9.987e-5
2.	Steady State (Gauss-Seidel)	956	9.911e-5
3.	Steady State (SOR)	90	8.80e-5
4.	Unsteady state Explicit	-	0
5.	Unsteady state Implicit (Jacobi)	1781	9.943e-5
6.	Unsteady state Implicit (Gauss-Seidel)	954	9.911e-5
7.	Unsteady state Implicit (SOR)	89	9.435e-5

• From the table, we can conclude that the Jacobi method took lots of iteration to converge and the Successive over-relaxation method took the least, whereas the Gauss – seidel iteration lies between the Jacobi and SOR.

 

CONCLUSION:

• Thus, in MATLAB, the 2D linear heat conduction equation is been solved by using the point iterative technique using the methods Jacobi, Gauss – seidel, Successive over Relaxation (SOR) for both implicit and explicit scheme.
• The temperature distribution plots are plotted for each method.
• The number of iteration and error values for each method is calculated and displayed in the above table.