Talk to a counsellor
CoursesIIT ProgramsProjectsPlacementsFor BusinessBlogLog inTalk to a counsellor
Challenge

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

SSubash Karan25 Sept 2023
MATLAB

AIM:

  • To address the 2D heat conduction problem in both steady and unsteady states through the application of iterative techniques such as Jacobi, Gauss-Seidel, and Successive Over-Relaxation

  • 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 Conduction, a way heat moves, mainly occurs in solids when temperature differences result from electrons and atoms vibrating. We can describe it with a general conduction equation.
  • Thermal diffusivity is a property of materials that tells us how quickly they can conduct heat. It's calculated by dividing the material's thermal conductivity by its heat capacity. In this 2D scenario without any heat generation, we can ignore those terms. Consequently, the equation describes unsteady 2D heat conduction, which has a parabolic nature
  • For steady-state, time derivative is not considered, their fore below equation represents the 2D steady-state Heat conduction (elliptic).
  •  We derive difference equations to calculate the numerical solutions for both the steady-state and unsteady-state 2D conduction equations 
  • Difference equations for the steady-state and unsteady-state 2D conduction equation are derived to obtain the numerical solution.

 

 

 


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

Steady-State analysis

  • Jacobi

  • Gauss-Siedel

  • Successive over relaxation

Transient analysis

  • Out of a number of iterations,1500 iterations is due to time step
  • Jacobi

  • Gauss-Siedel

  • Successive over-relaxation

  • Explicit
  • 1500 iterations are due to time step
  • According to CFL stability criteria, K1+K2<0>

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.