Solving the steady and unsteady 2D heat conduction problem

Observation :-

1. Assume that the domain is a unit square.

2. Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

3. Boundary conditions for steady and transient case

1. Left:400K

2. Top:600K

3. Right:800K

4. Bottom:900K

4. Steady state analysis & Transient State Analysis

Explanation of Steady state :-

we have the assumptions as there is n covenction . no internal heat generation and no change of temperature gradient with respect to tiome therefore the generalised energy 2d equation is :

Solving this equation does not involve in time marching therefore the temperature gradient do not change with respect to time .

Steady state Analysis :-

Programe for steady state

clear all 
close all
clc

% steady state condition
n = 1;
% length of domain in unit square 
l = 1;
nx = 15;
ny = nx;

%2dimensional array
x = linspace(0,l,nx);
y = x;
T = 273*ones(nx,ny);
% initial temperature condition
T_left = 400;
T_right = 800;
T_top = 600;
T_bottom = 900;
% boundary condition for left 
T(2:ny-1,1) = T_left;
% boundary condition for right
T(2:ny-1,nx) = T_right;
% boundary condition for top
T(1,2:nx-1) = T_top;
% boundary condition for bottom
T(nx,2:ny-1) = T_bottom;

% calculating average temperature at the matrix corner
T(1,1) = (T_top + T_left)/2;
T(nx,ny) = (T_right + T_bottom)/2;
T(1,ny) = (T_top + T_right)/2;
T(nx,1) = (T_left + T_bottom)/2;

%inputs for jacobi,gauss & succesive iteration 
T_old = T;
T1 = T;
T2 = T;
T3 = T;

% factor to be considered
[xx,yy] = meshgrid(x,y);
%three iteration inputs
jacobi = 1; gauss = 2; successive = 3;
dx = 1/(nx-1);
dy = dx;
k1 = (dy^2)/(2*(dx^2+dy^2));
k2 = (dx^2)/(2*(dx^2+dy^2));
for Method = 1:3
    if Method == 1
    %jacobi iteration 
    error1 = 9e9;
    jacobi_iteration = 0;
    while error1 > 1e-4
        jacobi_iteration = jacobi_iteration + 1;
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                            T1(i,j) = (k1*(T_old(i-1,j)+T_old(i+1,j)))+(k2*(T_old(i,j+1)+T_old(i,j-1)));
            end
        end
        error1 = max(max(abs(T1-T_old)));
        T_old = T1;%updating t_old values
        subplot(3,1,1);
        surf(x,y,T1);
        title(['jacobi num of iteration =' num2str(jacobi_iteration)]);
        xlabel('x-axis');
        ylabel('y-axis');
        zlabel('Temperature');
        colorbar;
        subplot(3,1,2);
        pcolor(x,y,T);
        shading interp;
        xlabel('x-axis');
        ylabel('y-axis');
        zlabel('Temperature');
        subplot(3,1,3);
        contour(x,y,T1);
        xlabel('x-axis');
        ylabel('y-axis');
        zlabel('Temperature');
        pause(0.01)
    end
    end   

    %Gauss iteration
 
        if Method == 2
            error2 = 9e9;
            gauss_iteration = 0;
            while error2 > 1e-4
                gauss_iteration = gauss_iteration + 1;
                for i = 2:(nx-1)
                    for j = 2:(ny-1)
                        T2(i,j) = (k1*(T2(i-1,j)+T_old(i+1,j)))+(k2*(T2(i,j+1)+T_old(i,j-1)));
                    end
                end
                error2 = max(max(abs(T2-T_old)));
                T_old = T2;
                figure(2)
                subplot(3,1,1);
                surf(x,y,T2);
                title(['gauss num of itration =' num2str(gauss_iteration)]);
                xlabel('x-axis');
                ylabel('y-axis');
                zlabel('Temperature');
                colorbar
                subplot(3,1,2);
                pcolor(x,y,T);
                xlabel('x-axis');
                ylabel('y-axis');
                zlabel('Temperature');
                subplot(3,1,3);
                contour(x,y,T2);
                xlabel('x-axis');
                ylabel('y-axis');
                zlabel('Temperature');
                pause(0.01)
            end
        end
 
            
            % Successive over relaxation
            
        if Method == 3
            c = cos(pi/nx)+cos(pi/ny);
            % relaxation parameter
            w = 4/(2+sqrt(4-c^2));
            successive_iteration = 0;
            error3 = 9e9;
            while error3 > 1e-4
                successive_iteration = successive_iteration + 1;
                for i = 2:(nx-1)
                    for j = 2:(ny-1)
                        T3(i,j) = (1-w)*T3(i,j)+w*(T3(i-1,j)+T_old(i+1,j)+T3(i,j-1)+T_old(i,j+1))/4;
                    end
                end
                error3 = max(max(abs(T3-T_old)));
                T_old = T3;
                figure(3);
                subplot(3,1,1);
                surf(x,y,T3);
                title([' successive num of iteration =' num2str(successive_iteration)]);
                xlabel('x-axis');
                ylabel('y-axis');
                zlabel('Temperature');
                colorbar
                subplot(3,1,2);
                pcolor(x,y,T);
                xlabel('x-axis');
                ylabel('y-axis');
                zlabel('Temperature');
                subplot(3,1,3);
                contour(x,y,T3);
                xlabel('x-axis');
                ylabel('y-axis');
                zlabel('Temperature');
                pause(0.01);
            end
        end
end

Explanation of Steady - state :-

Steady state is when heat being generated in a system is in equilibrium with the heat being dissipated or taken away from the system.

Consider a steady state system where heat is being generated, and is being conducted and / radiated away or, in addition, also being radiated away. Steady state condition for such system is when its temperature at each and every point in the system attains a unchanging fixed level. Heat being generated at any instant is the heat being dissipated.

Result for steady state Analysis :-

1. Jacobi

2. Gauss-seidel :-

3. Successive over-relaxation :-

Output :-

1.Jacobi :-

Num of iteration : 474

Error1 = 9.7281e-05

2.Gauss-seidel :-

Num of iteration : 334

Error2 = 9.8304e-05

3. Successive over-relaxation :-

Num of iteration : 49

Error3 = 9.8760e-05

Explantion for programe and plots :-

• Steady state form we have initiated a time step method for jacobi, gauss, successive .
• Here I given the tolerence value 1x10-4 have the maximum error .
• After that initiated the temperature condition for top, Bottom , left , right.
• By using meshgrid command for 2D array.
• Inputs have given for T1,T2,T3 are T & T_initial also given as T. then error value has 9x10^9
• Then using FOR loop iteration for jacobi , gauss , successive . we could calculated time step and error value by using the formulas.
• After that plot has been open in surf command , subplot command , contou plot command in that we have seen step iteration one by one.

Workspace for steady state :-

Transient State Analysis :-

Therefore

In explicit method the temperature at next level is calculated by means of surrounding temperature values and increasing the time steps i.e dt shall resut in oscillation of result .

Explanation for unsteady state :-

• Transient means unsteady.
• When there is a change in temperature of the body with time then it is said to transient or unsteady state.
• Whereas in case of steady state, no change in the temperature of the body with time.
• Transient heat conduction is more complex than steady state conduction.
• transient conduction is a mode of heat transfer, where it's properties will change with respect to time. It need not to be temperature alone. It can be any property.if it varies with respect to time,we can say that as in transient condition.

Programe for Impicit :-

clear all 
close all
clc

% variables
l = 1;
nx = 15;
ny = nx;
% total time
t = 18000;
% num of time step 
nt = 1800;
% time step
dt = t/nt;

%2dimensional array
x = linspace(0,l,nx);
y = x;
T = 273*ones(nx,ny,nt);

%Thermal diffusivity of gold
alpha = 127e-6;
dx = 1/(nx-1);
dy = dx;

% initial temperature condition
for k = 1:nt+1
T(1,:,k) = 900;
T(end,:,k) = 600;
T(:,1,k) = 400;
T(:,end,k) = 800;
T(1,1,k) = (900+400)/2;
T(1,end,k) = (900+800)/2;
T(end,1,k) = (400+600)/2;
T(end,end,k) = (800+600)/2; 
end

% factor to be considered
[xx,yy] = meshgrid(x,y);
k1 = (alpha*dt)/dx^2;
k2 = (alpha*dt)/dy^2;
g = (1+(2*k1)+(2*k2))^-1;
T_initial = T(:,:,1);
T1 = T; T2 = T;T3 = T;
T1_old = T(:,:,1);T2_old = T(:,:,1);T3_old = T(:,:,1);
error = 9e9;
jacobi_iteration = 0;
%Loop for Iterative solvers
for Method = 1:3
   %Jacobi Iterative Method
   if Method == 1
   %Time Loop
   for k = 1:nt
      %Convergence Loop
     if error > 1e-4 
        jacobi_iteration = jacobi_iteration+1;
    %Nodal Loops
    for i = 2:nx-1
        for j = 2:ny-1
           T1(i,j,k+1) = (T1_old(i,j)+(k1*(T1(i+1,j,k)+T1(i-1,j,k)))+(k2*(T1(i,j+1,k)+T1(i,j-1,k))))*g;
        end
    end 
   %Plotting Jacobi Method
    figure(1)
    subplot(2,1,1); 
    surf(xx,yy,T_initial);
    xlabel('X-axis')
    ylabel('Y-axis')
    zlabel('Temperature')
    title('Initial Temperature')
    subplot(2,1,2)
    surf(xx,yy,T1(:,:,k))
    print1 = sprintf('Jacobi number of Iteration = %d',jacobi_iteration);
    title({'Final Temperature';print1})
    xlabel('X-axis')
    ylabel('Y-axis')
    zlabel('Temperature')
    error =max(max(abs(T1(:,:,k)-T1(:,:,k+1))));
    T1_old = T1(:,:,k+1);
    
     end
   end
   end

    %Gauss Iterative Method
    if Method == 2
       error = 9e9;
       gauss_iteration = 0;
    for k = 1:nt
    %Convergence Loop
    if error > 1e-4
        gauss_iteration = gauss_iteration + 1;
    %Nodal Loop    
    for i = 2:nx-1
        for j = 2:ny-1
            T2(i,j,k+1) = (T2_old(i,j)+(k1*(T2(i+1,j,k)+T2(i-1,j,k+1)))+(k2*(T2(i,j+1,k)+T2(i,j-1,k+1))))*g;
        end
    end 
     %Plotting Gauss Method
     figure(2)
     subplot(2,1,1)
     surf(xx,yy,T_initial)
     xlabel('X-axis')
     ylabel('Y-axis')
     zlabel('Temperature')
     title('Initial Temperature')
     subplot(2,1,2)
     surf(xx,yy,T2(:,:,k))
     xlabel('X-axis')
     ylabel('Y-axis')
     zlabel('Temperature')
     print2 = sprintf('Gauss number of iteration = %d',gauss_iteration);
     title({'Final Temperature';print2})
     xlabel('X-axis')
     ylabel('Y-axis')
     zlabel('Temperature')
     error =max(max(abs(T2(:,:,k)-T2(:,:,k+1))));
     T2_old = T2(:,:,k+1);

    end
    end
    end

   
      %SOR Iterative Method
      if Method == 3
      z = cos(pi/nx)+cos(pi/ny);
      w = 4/(2+sqrt(4-z^2));
      error = 9e9;
      successive_iteration = 0;
      %Time Loop
for k = 1:nt
    %Convergence Loop
    if error > 1e-4
        successive_iteraton = successive_iteration + 1;
    %Nodal Loops    
    for i = 2:nx-1
        for j = 2:ny-1
            T3(i,j,k+1)=(1-w)*T3_old(i,j)+w*(T3_old(i,j)/(1+4*k1)+(k2/(1+4*k2))*(T3(i,j-1,k)+T3(i,j+1,k)+T3(i-1,j,k)+T3(i+1,j,k)));
        end
    end 
   %Plotting SOR Method
   figure(3)
   subplot(2,1,1)
   surf(xx,yy,T_initial)
   xlabel('X-axis')
   ylabel('Y-axis')
   zlabel('Temperature')
   title('Initial Temperature')
   subplot(2,1,2)
   surf(xx,yy,T3(:,:,k))
   xlabel('X-axis')
   ylabel('Y-axis')
   zlabel('Temperature')
   print3 = sprintf('successive number of iteration = %d',successive_iteration);
   title({'Final Temperature';print3})
   xlabel('X-axis')
   ylabel('Y-axis')
   zlabel('Temperature')
   error =max(max(abs(T3(:,:,k)-T3(:,:,k+1))));
   T3_old = T3(:,:,k+1);
    end
end
     end
end

Explantion for programe and plots :-

• Implicit form we have initiated a time step method for jacobi, gauss, successive .
• Here I given the thermal diffusivity for gold alpha = 127x10^-6
• After that initiated the temperature condition for top, Bottom , left , right.
• By using meshgrid command for 2D array.
• Inputs have given for T1,T2,T3 are T & T_initial also given as T. then error value has 9x10^9
• Then using FOR loop iteration for jacobi , gauss , successive . we could calculated time step and error value by using the formulas.
• After that plot has been open in surf command in that we have seen time step iteration one by one.

Results for Implicit :-

1.Jacobi :-

2.Gauss-seidal :-

3. Successive over-relaxation :-

Output :-

1.Jacobi :-

Num of iteration : 899

Error = 9.7970e-05

2.Gauss-seidel :-

Num of iteration : 696

Error = 9.7970e-05

3. Successive over-relaxation :-

Num of iteration : 0

Error = 9.7970e-05

Workspace for Implicit :-

clear all
close all
clc
%Defining variables
nx = 15;
ny = nx;
%Total time
t = 1800;
%num of timesteps
nt = 180;
%timestep
dt = t/nt;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,y);      
T = 273*ones(nx,ny,nt);
%Thermal diffusivity of gold
alpha = 127e-6;
dx = 1/(nx-1);
dy = dx;
c = (alpha*dt)/(dx^2);
%Assigning boundary conditions
for k = 1:nt+1
T(1,:,k) = 900;
T(end,:,k) = 600;
T(:,1,k) = 400;
T(:,end,k) = 800;
T(1,1,k) = (900+400)/2;
T(1,end,k) = (900+800)/2;
T(end,1,k) = (400+600)/2;
T(end,end,k) = (800+600)/2; 
end
k1 = (alpha*dt)/dx^2;
k2 = (alpha*dt)/dy^2;
%Initial temperature as per boundary condition
T_initial = T(:,:,1);
%Time Loop
for k = 1:nt
    %Nodal Loops
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j,k+1) = (T(i,j,k)+(k1*(T(i+1,j,k)+T(i-1,j,k)-(2*T(i,j,k))))+(k2*(T(i,j+1,k)+T(i,j-1,k)-(2*T(i,j,k)))));
        end
    end 
end
T_final = T(:,:,nt);
%Plotting the temperature change
figure(1)
subplot(2,2,1)
surf(xx,yy,T_initial)
title('Initial Temperature along the surface of the plate')
colorbar;
subplot(2,2,2)
surf(xx,yy,T_final)
title('Final Temperature along the surface of the plate')
colorbar;
subplot(2,1,2)
ct = 1;
for k = 1:nt
    Time = sprintf('Time = %d%s',k*20,'Sec');
    T_time = T(:,:,k);
    surf(xx,yy,T_time)
    set(gcf,'units','normalized','outerposition',[0 0 1 1])
    shading interp
    title({'Temperature change with respect to time';Time})
    colorbar;
    M(ct)=getframe(gcf);
    pause(1.0)
    ct=ct+1;
end

Explanation for explicit :-

• Explicitt form we have initiated a initial and final temperature has calculated t seconds.
• Here I given the thermal diffusivity for gold alpha = 127x10^-6
• After that initiated the temperature condition for top, Bottom , left , right.
• By using meshgrid command for 2D array.
• Inputs have given for T1,T2,T3 are T & T_initial also given as T. then error value has 9x10^9
• Then using FOR loop iteration. we could calculated time loop and error value by using the formulas.
• After that plot has been open in surf command in that we have seen time step iteration one by one.
• In this plot temperature will change with respect to time .

Result for Explicit :-

Workspace for explicit -