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

Aim- to solve the steady and unsteady state 2d heat conduction problem using iterative solvers- Jaobi, Gauss seidel method, Successive over-relaxation(SOR).

Objective-

• the boundary conditions given-

Top Boundary- 600k

Bottom boundary- 900k

Right Boundary- 800k

left Boundary - 400k

• absolute error consider as - 1e-4

Theory-

Conduction is a mode of heat transfer, that 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 the 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).

A Steady-state doesn't consider a time derivative. so, by rewrite the 2d steady state heat conduction problem is given by-

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

Steady-state 2d heat conduction equation

Transient Heat conduction equation in implicit form

Transient Heat conduction equation in explicit form

Now using the iterative solvers we will run the code and solve the above equation.

Steady State Analysis of a 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.

Jacobi iterations method-

`T(i,j)= 1/k[(T(i-1,j)+T(i+1,j))/(trianglex^2)+(T(i,j-1)+T(i,j+1))/(triangley^2)`

Gauss-Seidel Iterations method

`T(i,j)= 1/k[(T(i-1,j)+Told(i+1,j))/(trianglex^2)+(T(i,j-1)+Told(i,j+1))/(triangley^2)]`

SOR iterations method-

T(i,j)=((1-w)Told+omega)1/k[(T(i-1,j)+Told(i+1,j))/(trianglex^2)+(T(i,j-1)+Told(i,j+1))/(triangley^2)]

conditions for omega(w)= w>1 = over relaxation

w<1 = under relaxation

w= 1.1,1.2 = perfect solution

Matlab code for Steady-state analysis-

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;
no_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;
        no_iteration = no_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;
        no_iteration = no_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;
        no_iteration = no_iteration+1;
        
    end
end

% to visualize the result 

[X,Y] = meshgrid(x,y);
contourf(X,Y,T)
colorbar
set(gca,'ydir','reverse')

Unsteady 2d heat conduction (both implicit and explicit method)

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

Matlab 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
    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
    else
        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
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

% to visualize the result 

[X,Y] = meshgrid(x,y);
figure(1)
contourf(X,Y,T)
colorbar
set(gca,'ydir','reverse')

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.