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

SShivbrat Singh23 Nov 2023
MATLAB

SOLVING STEADY AND UNSTEADY 2D HEAT CONDUCTION EQUATION USING ITERATIVE METHOD

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

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