SIMULATION OF 2-D HEAT CONDUCTION EQUATION SOLUTIONS IN TRANSIENT AND STEADY STATES

SIMULATION OF 2-D HEAT CONDUCTION EQUATION SOLUTIONSIN TRANSIENT                                                            AND STEADY STATES

Basic Introduction About Conduction

• Thermal conduction is the transfer of heat by molecular collisions of particles and movement of electrons with high kinetic energy within a body.
• It is an important mode of heat transfer that occurs between two bodies in contact. Conduction requires an temperature gradient between two bodies to take place, the other important modes of heat transfer are convection and radiation.
• Conduction is essentially a time-dependant phenomena, the transient nature arises from the governing energy equation in 3 dimensions and transient conduction eventually moves towards steady state. (i.e. when the temperature gradients that have been established within a body do not change with time)

NOTE: MS Word Equations in form of pictures have been used instead of Math Equations as the dialog box was not opening when clicking Insert Math Option. Same has been notified to Skill Lync via Email.

DERRIVATION OF TRANSIENT STATE NUMERICAL SCHEME

From the generalized Energy equation in 3 Dimensions we make the following assumptions:

• No convection occurs
• No internal heat formation or heat loss
• Thermal diffusibility factor \'alpha\' is independent of space and time and therefore, system is isolated and hence, system properties are constant

In 2 Dimensions, the equation becomes:

Descretization of space terms (x,y) by Central Differencing Scheme. (O(h^2))

Discretization of time derrivative by forward differencing scheme. (O(h))

Therefore the resultant equation becomes:

Defining 2 constants K1 and K2:

Here Delta(x) = Delta(y)

Hence the equatuion finally becomes: 

DERRIVATION OF STEADY STATE NUMERICAL SCHEME

MATLAB INPUT CODE(s) WITH EXPLANATIONS

Here n = no. of grid points

acc = absolute error value as specified in the problem

%Code for TRANSIENT 2D heat conduction (SOLVED IMPLICITLY)

clear all
close all
clc

%Initialization and grid formation
nx=51;
ny=nx;
x=linspace(0,1,nx); y=linspace(0,1,ny);        %Equally spacing a square domain into 51 horizontal and 51 vertical strips
dx=1/(nx-1); dy=1/(ny-1);                      %Specifying mesh element size in x and y directions
alpha = 1e-4;                                  %Specifying Value of Thermal Diffusibility
t_start = 0;                                   %Specifying time at the start of the observtion/study
t_end = 1/alpha;                               %Specifying time at the start of the observtion/study
dt = 0.01;                                     %Specifying differential time or the time step
error = 1e9;                                   %Initializing error value
acc = 1e-4;                                    %Stop the iteration when achieved 4 decimal accurate convergence
CFL_Number = alpha*dt/(dx^2);                  %Courant Fredreich Lewy Condition for convergence
Nt = [t_start:dt*50:t_end];                    %Time spacing


%SPECIFICATION OF BOUNDARY CONDITIONS 
T=298*ones(nx,ny);              %Setting standard temperature at all nodes initially
T(1,:) = 900;                   % Bottom Temperature
T(end,:) = 600;                 % Top Temperature
T(:,1) = 400;                   % Left Temperature
T(:,end) = 800;                 % Right Temperature

% Refinement of the created mesh at the edges
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%Temperature Caliberation using Explicit Approach by iterative method
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
counter = 1;
solver = input(\'Enter Number:    \')
tic;

for k = 1:20000
   %GAUSS JACOBI SOLVER
    if solver == 1
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > acc)
            for i = 2:nx-1
                for j = 2:ny-1
                    %Discretized Terms of the Equation
                    H = (Told(i-1,j)+Told(i+1,j));
                    V = (Told(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    
    %GAUSS SIEDEL SOLVER
    if solver == 2
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > acc)
            for i = 2:nx-1
                for j = 2:ny-1
                    %Discretized Terms of the Equation
                    H = (T(i-1,j)+Told(i+1,j));
                    V = (T(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    
    %SOR SOLVER
        a = 2/(1 + sin(pi*dx));                                            %Relaxation Parameter
        if solver == 3
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > acc)
            for i = 2:nx-1
                for j = 2:ny-1
                    %Discretized Terms of the Equation
                    H=(a*T(i-1,j)+(1-a)*Told(i-1,j)+Told(i+1,j));
                    V=(a*T(i,j-1)+(1-a)*Told(i,j-1)+Told(i,j+1));
                    T(i,j)=((T_prev_dt(i,j))*term1)+(term2*H)+(term3*V);
                end
            end 
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    T_prev_dt = T;
end
execution_time = toc;

figure(1);
str_array = [\"Gauss-Jacobi\",\"Gauss-Siedel\",\"SOR\"];
formatSpec = \"Unsteady state 2D Heat Transfer by implicit: (%s Method) \\n Time: %0.3fs \\n Iteration Counter: %d \\n Computation time: %0.3fs\";

[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel(\'xx\');
ylabel(\'yy\');
title(sprintf(formatSpec,str_array(solver),k*dt,counter,execution_time));

%Code for TRANSIENT 2D heat conduction (SOLVED EXPLICITLY)

clear all
close all
clc

%Initialization and grid formation
nx=51;
ny=nx;
T=298*ones(nx,ny);                             %Setting standard temperature at all nodes initially
x=linspace(0,1,nx);                            %Equally spacing a square domain into 51 horizontal strips
y=linspace(0,1,ny);                            %Equally spacing a square domain into 51 vertical strips
dx=1/(nx-1);                                   %Specifying mesh element size in x directions
dy=1/(ny-1);                                   %Specifying mesh element size in y directions
alpha = 1e-4;                                  %Specifying Value of Thermal Diffusibility
t_start = 0;                                   %Specifying time at the start of the observtion/study
t_end = 1/alpha;                               %Specifying time at the start of the observtion/study
dt = 0.01;                                     %Specifying differential time or the time step
error = 1e8;                                   %Initializing error value
acc = 1e-4;                                    %Stop the iteration when achieved 4 decimal accurate convergence
CFL_Number = alpha*dt/(dx^2);                  %Courant Fredreich Lewy Condition for convergence


%SPECIFICATION OF BOUNDARY CONDITIONS 
T(1,:) = 900;                   % Bottom Temperature
T(end,:) = 600;                 % Top Temperature
T(:,1) = 400;                   % Left Temperature
T(:,end) = 800;                 % Right Temperature

% Refinement of the created mesh at the edges
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%Temperature Caliberation using Explicit Approach by iterative method
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
tic;
counter = 1;

for z = 1:20000
    for i = 2:nx-1
        for j = 2:ny-1
            %Discretized Terms of the Equation
            term1 = Told(i,j);
            term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
            term3 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
            T(i,j) = term1+term2+term3;
        end
    end
    Told = T;
    counter = counter+1;
    
end
time_counter = toc;

figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf(\'2D Unsteady state heat conduction (Explicit Method) \\n Number of iterations: %d (Time: %0.2fs)  \\n Computation Time: %0.4f \\n CFL Number: %0.4f\',counter-1,z*dt,time_counter,CFL_Number));
xlabel(\'xx\');
ylabel(\'yy\');

%Code for steady state 2D heat equation (SOLVING IMPLICITLY)

clear all
close all
clc

%Initialization and grid formation
nx=51;
ny=nx;
x=linspace(0,1,nx); y=linspace(0,1,ny);        %Equally spacing a square domain into 51 horizontal and 51 vertical strips
dx=1/(nx-1); dy=1/(ny-1);                      %Specifying mesh element size in x and y directions
k = 2*(1/(dx^2)+1/(dy^2));                     %Coefficient of temperature of the current node under observation
acc = 1e-4;                                    %Stop the iteration when achieved 4 decimal accurate convergence1
error = 1e8;                                   %Initializing error value


%SPECIFICATION OF BOUNDARY CONDITIONS 
T=298*ones(nx,ny);              % Setting standard temperature at all nodes initially
T(1,:) = 900;                   % Bottom Temperature
T(end,:) = 600;                 % Top Temperature
T(:,1) = 400;                   % Left Temperature
T(:,end) = 800;                 % Right Temperature

% Refinement of the created mesh at the edges
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%Design of the iterative solver
x=y;
[xx,yy] = meshgrid(x,y);
Told = T;

SOLVER = input(\'Enter Number:    \')

tic;

%GAUSS JACOBI SOLVER
if SOLVER == 1
    iteration_count_jacobi = 1;
    while (error > acc)
        for i = 2:nx-1
            for j = 2:ny-1
                %Discretized Terms of the Equation
                term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
                term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
                T(i,j) = term1+term2;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iteration_count_jacobi = iteration_count_jacobi + 1;
    end
end

%GAUSS SIEDEL SOLVER
if SOLVER == 2
    iteration_count_siedel = 1;
    while (error > acc)
        for i = 2:nx-1
            for j = 2:ny-1
                %Discretized Terms of the Equation
                term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
                term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
                T(i,j) = term1+term2;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iteration_count_siedel = iteration_count_siedel + 1;
    end
end


%SUCCESSIVE OVER RELAXATION SOLVER
if SOLVER == 3
    a = 2/(1+sin(pi*dx));                                                  %Relaxation Parameter
    iteration_count_SOR = 1;
    while (error > acc)
        for i = 2:nx-1
            for j = 2:ny-1
                %Discretized Terms of the Equation
                term1 = (T(i-1,j)+Told(i+1,j))/(k*(dx^2));
                term2 = (T(i,j-1)+Told(i,j+1))/(k*(dy^2));
                 T(i,j) = a*(term1+term2)-(a-1)*Told(i,j);
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iteration_count_SOR = iteration_count_SOR + 1;
    end
end
time_elapsed = toc;

%sprintf(formatSpec,A1,...,An) formats the data in arrays A1,...,An according to formatSpec in column order, and returns the results to str.

figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if (SOLVER == 1)
    title(sprintf(\'Temperature profile of steady state heat conduction using GAUSS JACOBI METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.3f s.\',iteration_count_jacobi,time_elapsed));
end
if (SOLVER == 2)
        title(sprintf(\'Temperature profile of steady state heat conduction using GAUSS SIEDEL METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.3f s.\',iteration_count_siedel,time_elapsed));
end
if (SOLVER == 3)
    title(sprintf(\'Temperature profile of steady state heat conduction using SOR method. \\n No. of Iterations: %d. \\n Computation Time: %0.3f s. \\n Relaxation Parameter: %f\',iteration_count_SOR,time_elapsed,real(a)));
end
xlabel(\'xx\');
ylabel(\'yy\');

ANALYSIS RESULTS

    1. TRANSIENT STATE ANALYSYS

      2. STEADY STATE