Challenge

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

MMd Zeeshan Zafer21 Mar 2023

AIM:

To solve the 2D heat conduction problem for the steady and unsteady state using Jacobi, gauss seidel, and successive over-relaxation iterative techniques.

THEORY:

Conduction is a mode of heat transfer, generally takes place in solids due to temperature differences associated with free electron and lattice vibration.
A generalized conduction equation is given

Thermal diffusivity is a thermophysical property of a material defined as the ratio of thermal conductivity of the material and its thermal capacity.
Since the problem is 2D and there is no heat generation is considered, therefore those terms are neglected. So the equation represents the unsteady 2D heat conduction (parabolic).

For steady-state, time derivative is not considered, their fore below equation represents the 2D steady-state Heat conduction (elliptic).

Difference equations for the steady-state and unsteady-state 2D conduction equation are derived to obtain the numerical solution.

Steady State 2D heat conduction
Transient 2D heat conduction (explicit)
Transient 2D heat conduction (implicit)

Here iterative methods are used to solve the difference equation, Iterative methods are numerical techniques used to converge the solution(accurately) of a problem with some initial guess.
Here Jacobi, gauss-seidel, and successive over-relaxation techniques are used.

MATLAB CODE:

There is 2 Matlab code written separately for steady and unsteady conditions.
But for both conditions inputs are the same
while programming it is assumed that, the domain is the surface of area (1mX1m). Discretized into 10 nodes both horizontally and vertically.
Constant values like alpha and time step 'data, by keeping the explicit condition in mind.
relaxation factor omega is also assumed while programming only.

Inputs

Top boundary temperature	600 K
Bottom Boundary temperature	900 K
Left Boundary temperature	400 K
Right Boundary temperature	800 K
error	1e-4

UNSTEADY CONVECTION:

clear 
close all
clc

%Domain and discretization length
n = 10;  % Number of nodes
l = 1;   % domain length
tol = 1e-4;

% Value of x and y at nodes in a discretization plane

x = linspace(0,l,n);
y = linspace(0,l,n);
dx = l/(n-1);
dy = l/(n-1);

%Temperature at all nodes with initial guess 300k

T = 300*ones(n,n);

% Initial boundary condition

T(1,1:end) = 600;   % Top
T(end,1:end) = 900; % Bottom
T(1:end,1) = 400;   % Left
T(1:end,end) = 800; % Right

%Average Temperature at corner

T(1,1) = (600+400)/2 ;
T(end,end) = (900+800)/2;
T(1,end) = (600+800)/2;
T(end,1) = (900+400)/2;

% Copy of temperature to compute
T_old = T;
T_prev = T;

% Constant
alpha = 1.4;
total_time = 3;
dt = 1/500;
nt = total_time/dt;
Omega = 1.2;
K1 = (alpha*dt)/(dx^2);
K2 = (alpha*dt)/(dy^2);
cfl = K1+K2;
no_iteration = 1;


approach = input('Enter 1 for implicit or 2 for explicit ');

 if approach ==1
    method = input('Enter 1 for jacobi, 2 for gauss seidal, 3 for sor ');
 end


%IMPLICIT JACOBIEN METHOD

if approach == 1 && method == 1

    for k = 1:nt
         error = 10;
        while(error>tol)
            for i = 2:n-1
                for j = 2:n-1
                    term1 = 1/(1+(2*K1)+(2*K2));
                    term2 = K1*(T_old(i+1,j)+T_old(i-1,j));
                    term3 = K2*(T_old(i,j+1)+T_old(i,j-1));
                    T(i,j) = term1*(T_prev(i,j)+term2+term3);
                end
            end
            error = max(max(abs(T_old-T)));
            T_old = T;
            no_iteration = no_iteration+1;
        end
        T_prev = T;
    end

%PLOTTING

[X,Y] = meshgrid(x,y);
T = flipud(T);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT METHOD - JACOBI METHOD');
title_text_2 = sprintf('ITERATION NUMBER = %d',no_iteration);
title_text_3 = sprintf('CFL = %d',cfl);
title_text_4 = sprintf(' Time Step dt = %d', dt);
title({title_text_1,title_text_2,title_text_3,title_text_4});

end


%Gauss Method

if approach == 1  && method == 2

    for k = 1:nt
        error = 10;
        while(error>tol)
            for i = 2:n-1
                for j = 2:n-1
                    term1 = 1/(1+(2*K1)+(2*K2));
                    term2 = K1*(T_old(i+1,j)+T(i-1,j));
                    term3 = K2*(T_old(i,j+1)+T(i,j-1));
                    T(i,j) = term1*(T_prev(i,j)+term2+term3);
                end
            end
            error = max(max(abs(T_old-T)));
            T_old = T;
            no_iteration = no_iteration +1;
        end
        T_prev = T;
    end

 %PLOTTNG

[X,Y] = meshgrid(x,y);
T = flipud(T);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT METHOD - GAUSS SEIDEL METHOD');
title_text_2 = sprintf('ITERATION NUMBER = %d',no_iteration);
title_text_3 = sprintf('CFL = %d',cfl);
title_text_4 = sprintf(' Time Step dt = %d', dt);
title({title_text_1,title_text_2,title_text_3,title_text_4});

end

 %SOR METHOD

    if   approach == 1  && method == 3
   
    for k = 1:nt
        error = 10;
        while(error>tol)
            for i = 2:n-1
                for j = 2:n-1
                    term1 = 1/(1+(2*K1)+(2*K2));
                    term2 =K1*(T_old(i+1,j)+T(i-1,j));
                    term3 = K2* (T_old(i,j+1)+T(i,j-1));
                    term4 = term1*(T_prev(i,j)+term2+term3);
                    T(i,j)=(1-Omega)*T_old(i,j) +(Omega*term4);
                end
            end
            error = max(max(abs(T_old-T)));
            T_old = T;
            no_iteration = no_iteration+1;
         end
        T_prev = T;
   end

 %PLOTTING

[X,Y] = meshgrid(x,y);
T = flipud(T);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT METHOD - SOR');
title_text_2 = sprintf('ITERATION NUMBER = %d',no_iteration);
title_text_3 = sprintf('CFL  = %d',cfl);
title_text_4 = sprintf(' Time Step dt = %d', dt);
title({title_text_1,title_text_2,title_text_3,title_text_4});
    end

%EXPLICIT

if approach == 2
    for k = 1:nt
        for i = 2:n-1
            for j= 2:n-1
                term1 = K1*(T_old(i+1,j)-2*T_old(i,j)+T_old(i-1,j));
                term2 = K2*(T_old(i,j+1)-2*T_old(i,j)+T_old(i,j-1));
                T(i,j) = T_old(i,j)+term1+term2;
            end
        end
        T_old = T;
        no_iteration = no_iteration+1;
    end
      

 %PLOTTING   

[X,Y] = meshgrid(x,y);
T = flipud(T);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('EXPLICIT METHOD');
title_text_2 = sprintf('ITERATION NUMBER = %d',no_iteration);
title_text_3 = sprintf('CFL = %d',cfl);
title_text_4 = sprintf(' Time Step dt = %d', dt);
title({title_text_1,title_text_2,title_text_3,title_text_4});

end

Unsteady-State 2D Heat conduction

Inputs are defined into variables and the temperature array is created.
Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged.
Constants like alpha, omega, tolerance are defined to run the iteration.
Here a user must provide, one or two inputs.
The first input is to whether solve the problem using an explicit approach or an implicit approach.
If an implicit approach is selected, then the user must provide another input, to select the iterative technique.
The if-else loop is used to solve the equation using the selected(input from the user) approach and iterative technique.
For loop is used to run the iteration for a large number of times to reach steady-state conditions.
while loop is used to run the iterative technique till the required accuracy is achieved.
Two for loop is used to calculate the temperature values at every node at each iterative step.
Finally, a temperature distribution graph is obtained using 'conturf' command.

STEADY CONVECTION MAIN CODE

clear 
close all
% clc

% Domain and discretization length
n = 10;  % Number of nodes
l = 1;   % domain length

% Value of x and y at nodes in a discretization plane

x = linspace(0,l,n);
y = linspace(0,l,n);
dx = l/(n-1);
dy = l/(n-1);
Timetaken = 3;
dt = 1/10;
nt = Timetaken/dt;
%Temperature at all nodes with initial guess 300k

T = 300*ones(n,n);

% Initial boundary condition

T(1,1:end) = 600;  % Top
T(end,1:end) = 900; % Bottom
T(1:end,1) = 400; % Left
T(1:end,end) = 800; % Right

%Average Temperature at corner

T(1,1) = (600+400)/2 ;
T(end,end) = (900+800)/2;
T(1,end) = (600+800)/2;
T(end,1) = (900+400)/2;

% Copy of temperature to compute
T_old = T;

% Constant
Omega = 1.2;
K = (2*(dx^2+dy^2))/(dx^2*dy^2);

% To enter loop initialize error value randomly

error = 10;
tol = 1e-4;
no_iteration = 0;

method = input('Enter(1)for jacobi,(2) for gauss seidal,(3) for sor   : ');


%Jacobi method

if method == 1
   % k =1:nt;
    while(error>tol)
        for i = 2:n-1
            for j = 2:n-1
                term1 = (T_old(i+1,j)+T_old(i-1,j))/dx^2;
                term2 = (T_old(i,j+1)+T_old(i,j-1))/dy^2;
                T(i,j) = (term1+term2)/K;
            end
        end
        error = max(max(abs(T-T_old)));
        T_old = T;
        no_iteration = no_iteration+1;
    end

%PLOTTINNG

[X,Y] = meshgrid(x,y);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('yaxis')
title_text_1 = sprintf('STEADY STATE - JACOBI METHOD');
title_text_2 = sprintf('ITERATION NUMBER=%d',no_iteration);
title_text_3 = sprintf(' Time Step dt =%d', dt);
title({title_text_1,title_text_2});


end


%Gauss seidel method

if  method ==2
    while(error>tol)
        for i = 2:n-1
            for j = 2:n-1
                term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;
                term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;
                T(i,j) = (term1+term2)/K;
            end
        end
        error = max(max(abs(T_old-T)));
        T_old = T;
        no_iteration = no_iteration +1;
    end

%PLOTTINNG

[X,Y] = meshgrid(x,y);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('STEADY STATE - GAUSS SEIDEL METHOD');
title_text_2 = sprintf('ITERATION NUMBER=%d',no_iteration);
%title_text_3 = sprintf('CFL  =%d',cfl);
title_text_4 = sprintf(' Time Step dt =%d', dt);
title({title_text_1,title_text_2,title_text_4});

end


% SOR method

if method == 3
    while(error>tol)
        for i = 2:n-1
            for j = 2:n-1
                term1 =(T_old(i+1,j)+T(i-1,j))/dx^2;
                term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;
                T(i,j)=(1-Omega)*T_old(i,j) +(Omega*((term1+term2)/K));
            end
        end
        error = max(max(abs(T-T_old)));
        T_old = T;
        no_iteration = no_iteration+1;
    end

 %PLOTTINNG

[X,Y] = meshgrid(x,y);
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('yaxis')
title_text_1 = sprintf('STEADY STATE- SOR');
title_text_2 = sprintf('ITERATION NUMBER=%d',no_iteration);
%title_text_3 = sprintf('CFL  =%d',cfl);
title_text_4 = sprintf(' Time Step dt =%d', dt);
title({title_text_1,title_text_2,title_text_4});

end

Steady-State 2D Heat conduction

First inputs are defined into variables and the temperature array is created.
Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged.
Constants like alpha, omega, tolerance are defined to run the iteration.
Input must be provided by the user, to choose the iterative method.
The if-else loop is used to solve the equation using the selected(input from the user) iterative technique.
while loop is used to run the iterative technique till the required accuracy is achieved.
Two for loop is used to calculate the temperature values at every node at each iterative step.
Finally, the temperature distribution graph is obtained using 'conturf' command.

RESULT:

The program ran successfully.
The temperature distribution plot obtained for both the steady and unsteady conditions is the same.
Because in Transient analysis also, the solution approaches to steady-state.
The top, bottom, left, and right side of the graphs representing their respective boundary conditions.

The number of iterations taken to converge to the solution is compared among the 3 iterative techniques.
In both steady and transient analysis, the number of iteration by the Successive over-relaxation method is less,
The number of iteration taken by Successive over relaxaton<Number of iteration taken by gauss-siedel<Number of iteration taken by Jacobi.

Steady-State analysis