Challenge

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

DDhilipan Sivasamy1 Aug 2022

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

`(del^2 T)/(del x^2) + (del^2T)/(del y^2) + (del^2 T)/(del z^2) + q/k = 1/alpha (delT)/(del tau)`

Where,`T = text(Temperature)`

`q = text(Heat Generation)`

`alpha = text(Thermal diffusivity)`

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).

`((del^2 T)/(del x^2) + (del^2 T)/(del y^2)) alpha = (del T)/(del tau)`

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

`((del^2 T)/(del x^2) + (del^2 T)/(del y^2)) = 0`

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

Steady State 2D heat conduction

`T_p = 1/K [((T_R + T_L)/(Delta x^2)) + ((T_T+T_B)/(Delta y^2))]`

`K = 2((Delta x^2+Delta y^2)/(Delta x^2Delta y^2))`

Transient 2D heat conduction (explicit)

`T_p^(n+1) = T_p^n + K_1(T_R^n - 2T_p^n + T_L^n) + K_2(T_T^n -2T_p^n + T_B^n)`

`K_1 = alpha((Delta t)/(Deltax^2))`

`K_2 =alpha((Delta t)/(Delta y^2))`

Transient 2D heat conduction (implicit)

`T_p^(n+1) = 1/((1 + 2K_1 + 2K_2))[T_p^n+K_1(T_R^(n+1) + T_L^(n+1)) + K_2(T_T^(n+1) + T_B^(n+1))]`

`K_1 = alpha((Delta t)/(Deltax^2))`

`K_2 =alpha((Delta t)/(Delta y^2))`

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

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.

Code :

clear all
close all 
clc

% Domain and dicretization length
number_node = 10;
domain_length = 1;

%value of x and y at nodes in a discretization plane
x = linspace(0,domain_length,number_node);
y = linspace(0,domain_length,number_node);
dx = domain_length/(number_node-1);
dy = domain_length/(number_node-1);


% Let define a temperature value at all nodes with initial guess 300k
T = 300*ones(number_node,number_node);

% Initia boundary conditions at left and right
T(1:end,1) = 400;
T(1:end,end) = 800;

% Initia boundary conditions at bottom and top
T(end,1:end) = 900;
T(1,1:end) = 600;

%average values of T at corner boundaries
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;

%copy of Temperature to compute to compute
T_old = T;

%constant
K = (2*(dx^2+dy^2))/(dx^2*dy^2);
omegha = 1.2;

% to enter while loop initialize error value randomly
error = 10;
tol = 1e-4;
iteration = 0;

method = input('Enter(1) for jacobi,(2) for gauss siediel,(3) for sos :');
if method == 1
    while(error>tol)
        for i = 2:number_node-1
            for j = 2:number_node-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_old-T)));
        
        % to upadte the previous velcity with new velocity
        %got from the iteration
        T_old = T;
        iteration = iteration+1;
        
    end
elseif method == 2
    while(error>tol)
        for i = 2:number_node-1
            for j = 2:number_node-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)));
        
        % to upadte the previous velcity with new velocity
        %got from the iteration
        T_old = T;
        iteration = iteration+1;
    end
    
else
    while(error>tol)
        for i = 2:number_node-1
            for j = 2:number_node-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-omegha)*T_old(i,j))+(omegha*((term1+term2)/K));
            end
        end
        error = max(max(abs(T_old-T)));
        
        % to upadte the previous velcity with new velocity
        %got from the iteration
        T_old = T;
        iteration = iteration+1;
        
    end
end

% to visualize the result 
 for i = 1:method
[X,Y] = meshgrid(x,y);
contourf(X,Y,T);
[C,h]=contourf(X,Y,T);
colorbar;
colormap('jet');
clabel(C,h);
xlabel('x-axis');
ylabel('y-axis');
ii = {'Jacobi Process','Gauss Siediel Process','SOS Process'};
str = strcat({'Method '}, num2str(i),' is : ',ii(i));
title_text = sprintf('iterantion number =%d',iteration);
title(str,title_text,'fontweight','bold');
set(gca,'ydir','reverse')
 end

Result :

1. jacobi Method :

2. Gauss Siedel Method :

3. Successive over relaxation Method :

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.

Code :

clear all
close all 
clc

% Domain and dicretization length
number_node = 10;
domain_length = 1;

%value of x and y at nodes in a discretization plane
x = linspace(0,domain_length,number_node);
y = linspace(0,domain_length,number_node);
dx = (domain_length/number_node-1);
dy = (domain_length/number_node-1);
dt = 1/10;


% Let define a temperature value at all nodes with initial guess 300k
T = 300*ones(number_node,number_node);

% Initia boundary conditions at left and right
T(1:end,1) = 400;
T(1:end,end) = 800;
% Initia boundary conditions at bottom and top
T(end,1:end) = 900;
T(1,1:end) = 600;
%average values of T at corner boundaries
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;

%copy of Temperature to compute for previous space terms
T_old = T;

% copy of T for previous time_step
T_prev = T;

%constants 
alpha = 1.5;
omegha = 1.15;
K1 = (alpha*dt)/(dx^2);
K2 = (alpha*dt)/(dy^2);
iteration = 0;


approach = input('Enter 1 for implicit or Enter 2 for explicit : ')
if approach == 1;
else
     for nt = 1:1500
        for i = 2:number_node-1
            for j = 2:number_node-1
                term1 = K1*(T_prev(i+1,j)-2*T_prev(i,j)+T_prev(i-1,j));
                term2 = K2*(T_prev(i,j+1)-2*T_prev(i,j)+T_prev(i,j-1));
                T(i,j) = T_prev(i,j)+term1+term2;
            end
        end
        % to upadte the previous timestep value with new timestep value
        %got from the iteration
        T_prev = T;
        iteration = iteration+1;
    end
end
method = input('Enter 1 for jacobi or 2 for gauss siedel or 3 for sos : ')
if method == 1;
    for nt = 1:1500
        error = 10;
        tol = 1e-4;
        while(error>tol)
            for i = 2:number_node-1
                for j = 2:number_node-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
            iteration = iteration+1;
            error = max(max((abs(T_old-T))));
            T_old = T;
        end
        T_prev = T;
    end
elseif method == 2;
    for nt = 1:1500
        error = 10;
        tol = 1e-4;
        while(error>tol)
            for i = 2:number_node-1
                for j = 2:number_node-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
            iteration = iteration+1;
            error = max(max(abs(T_old-T)));
            T_old = T;
            
        end
        T_prev = T;
    end
elseif method == 3;
    for nt = 1:1500
        error = 10;
        tol = 1e-4;
        while(error>tol)
            for i = 2:number_node-1
                for j = 2:number_node-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-omegha)*T_old(i,j)+(omegha*term4);
                end
            end
            iteration = iteration+1;
            error = max(max(abs(T_old-T)));
            T_old = T;
        end
        T_prev = T;
    end
end


% to visualize the result 
for i = 1:approach
    for j = 1:method
[X,Y] = meshgrid(x,y);
figure(1)
contourf(X,Y,T);
[C,h]=contourf(X,Y,T);
colorbar;
colormap('jet');
clabel(C,h);
xlabel('x-axis');
ylabel('y-axis');
ii = {'implicit','explicit'};
iii = {'Jacobi Process','Gauss Siediel Process','SOS Process'};
str = strcat({'Approach '}, num2str(i),{' is :- '},ii(i),{' ; Methode is : '},iii(j));
title_text = sprintf('iterantion number = %d',iteration);
title(str,title_text,'fontweight','bold');
set(gca,'ydir','reverse')
    end
end

Result :

1. Implicit Approach :

1. Jacobi Method :