Steady and Unsteady 2D heat conduction Problem

Aim: To solve the steady and unsteady/Transient 2D Heat conduction problem using iterative methods(Jacobi,Gauss seidal, Successive Over Relaxation).

1. Assume that the domain is a unit square

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

3. Boundary conditions for steady and transient case

          (i) Left : 400 K

          (ii) Top : 600 K

          (iii) Right : 800 K

           (iv) Bottom : 900 K

Introduction:

Heat conduction is a mode of heat transfer, which transfers the heat from one substance to other substance when they are in physical contact or one end of the substance to other end of the same substance. 

In solids,

      Conduction is accomplished by lattice vibrations and drift of electrons. 

In liquids, 

     Conduction is accomplished by collision of molecules and diffusion.

 Convection and Radiation are the other two modes of heat transfer. 

Generalized Heat conduction through 3 dimensional coordinates is 

`α({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} + {∂^2 T}/{∂z^2} ) +{dot q}/k = {∂T}/{∂t}`

By eliminating the internal heat generation, finally we get

`α({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} + {∂^2 T}/{∂z^2} ) = {∂T}/{∂t}`

From the above 3D governing equation it is clear that the conduction is a time dependent, the transient nature arises from 3 dimensions and eventually move towards steady state(the temperature doesnot vary with time). 

TRANSIENT/UNSTEADY SOLUTION:

Derivation of numerical scheme:

From the generalized Energy equation in 3 dimensions we make the following assumptions

          (i) convection and Radiation effects doesnot occurs

         (ii) No internal heat generation or dissipation

         (iii) constant properties of system i.e. `α != f(t,x,y)`

  Hence for 2 dimensional  condition we have:

`α({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} ) = {∂T}/{∂t}`

The second order of spatial terms is discretized by central differencing

`{∂^2 T}/{∂x^2} = (T_(i-1,j) -2*(T)_(i,j)+T_(i+1,j) )/{Δx^2}`

`{∂^2 T}/{∂y^2} = (T_(i,j-1) -2*(T)_(i,j)+T_(i,j+1) )/(Δy^2)`

The temporal terms is discretized by forward differencing

`{∂T}/{∂t} = ((T^(n+1))_(i,j) - (T^n)_(i,j))/{Δt}`

Hence we have:

`α((T_(i-1,j) - 2*(T)_(i,j)+T_(i+1,j) )/{Δx^2} + (T_(i,j-1) - 2*(T)_(i,j)+T_(i,j+1) )/{Δy^2}) = {(T^(n+1))_(i,j) - (T^n)_(i,j)}/{Δt}`

we define two constants K1 and K2 such that:

`K_1 =α {∆t}/{(Δx)^2}` and  `K_2 =α{∆t}/{(Δy)^2}`

 Here  Δx = Δy

Hence the equation:

`(T^(n+1))_(i,j) = (T^n)_(i,j) + K_1((T^(n+1))_(i-1,j) - 2(T^(n+1))_(i,j) + (T^(n+1))_(i+1,j)) + K_2((T^(n+1))_(i,j-1) - 2(T^(n+1))_(i,j) + (T^(n+1))_(i,j+1))`

`(T^(n+1))_(i,j) - K_1(T^(n+1))_(i-1,j) + 2 K_1(T^(n+1))_(i,j) - K_1 (T^(n+1)) _(i+1,j) - K_2(T^(n+1))_(i,j-1) + 2 K_2(T^(n+1))_(i,j) - K_2(T^(n+1))_(i,j+1) = (T^n)_(i,j)`

`(T^(n+1))_(i,j)(1+ 2 K_1 + 2 K_2) - K_1(T^(n+1))_(i-1,j) - K_1(T^(n+1))_(i+1,j) - K_2(T^(n+1)) _(i,j-1) - K_2(T^(n+1))_(i,j+1) = (T^n)_(i,j)`

`(T^(n+1))_(i,j) = ((T^n)_(i,j) + K_1(T^(n+1))_(i-1,j )+ K_1 (T^(n+1))_(i+1,j) + K_2(T^(n+1))_(i,j-1) + K_2(T^(n+1))_(i,j+1))/((1+ 2 K_1 + 2 K_2))`

 Here n denotes time index and (i, j) denotes space index

 The other parameters are assumed as: 

              Thermal diffusivity,`α = 1xx10^-1` 

              Time step, dt = 0.01 seconds

               number of grid points, nx =ny = 41

               size of the grid, dx=dy= 0.0250

               tolerance value = `10^-4`

               omega = 1.1

             ambient temperature = `30^oC`

SOLVING THE UNSTEADY STATE PROBLEM EXPLICITLY:

MATLAB CODE:

% Matlab code to solve 2D transient state Heat conduction equation explicitly

close all
clear all
clc

% Inputs
            % number of grid points in x direction, nx         
                 nx = 41
            % number of grid points in y direction, ny
                 ny = nx
            % thermal diffusivity, alpha
                alpha = 1e-4
            % Initializing time to start in seconds 
                t_start = 0
            % Initializing time to end in seconds
                t_end = 10000
            % time step in seconds
               dt = 0.01
% calculation
      % x array along the length of the domain in x direction
               x = linspace(0,1,nx)
      % y array along the length of the domain in y direction
               y = linspace(0,1,ny)
      % Grid spacing in x direction, dx 
               dx = 1/(nx-1)
       % Grid spacing in y direction, dx = 1/(nx-1)
               dy = 1/(ny-1)
               
% Initializing the temperatures in x and y direction
                  T = 303*ones(nx,ny)
% Assigning Boundary conditions
                   T(1,:) = 600; % Top temperature
                   T(end,:) = 900; % Bottom temperature
                   T(:,1) = 400; % Left temperature
                   T(:,end) = 800;  % Right temperature
% Edge Refining 
                   T(1,1) = (600+400)/2 ;
                   T(1,end) = (600+800)/2 ;
                   T(end,1) = (900+400)/2 ;
                   T(end,end) = (900+800)/2;

% calculation of temperature distribution explicitly using iterative solvers
               k1 = (alpha*dt)/(dx^2);
               k2 = (alpha*dt)/(dy^2);
               [xx,yy] = meshgrid(x,y);
               Told = T;
               T_prev_dt = Told;
tic;
counter = 1
error = 1e9
tol = 1e-4
for t = 1:10000
    for i = 2:nx-1
        for j = 2:ny-1
            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);
xlabel('xx');
ylabel('yy');
title(sprintf('2D Unsteady state Heat transfer (Explicit Method) n Number of iterations:%d  (Time: %0.2fs)', counter-1,time_counter));

RESULTS & CONCLUSION FOR UNSTEADY STATE 2D HEAT CONDUCTION(EXPLICITLY):

Here the timetaken for the iterations is very small because of the temperature at the time n+1 depends explicitly on the temperature at time n(i.e., Told). Here this is relatively simple and fast to compute the temperature but main drawback is stable solutions are obtained only if `0 < k{Delta t}/{Delta x^2} + k{Delta t}/{Delta y^2} le 0.5` . From the above plot it is clear that the temperature is a function of time. 

SOLVING THE  UNSTEADY STATE PROBLEM IMPLICITLY:

MATLAB CODE:

%% Matlab program to solve 2D Unsteady state Heat conduction equation implicitly

close all
clear all
clc

% Inputs
            % number of grid points in x direction, nx         
                 nx = 41
            % number of grid points in y direction, ny
                 ny = nx
            % thermal diffusivity, alpha
                alpha = 1e-4
            % Initializing time to start in seconds 
                t_start = 0
            % Initializing time to end in seconds
                t_end = 10000
            % time step in seconds
               dt = 0.01
               
% calculation
          % x array along the length of the domain in x direction
               x = linspace(0,1,nx)
          % y array along the length of the domain in y direction
               y = linspace(0,1,ny)
          % Grid spacing in x direction, dx 
               dx = 1/(nx-1)
          % Grid spacing in y direction, dx = 1/(nx-1)
               dy = 1/(ny-1)
          % No. of time steps
               nt = [t_start:dt:t_end]
% Initializing the temperatures in x and y direction
                  T = 303*ones(nx,ny)
% Assigning Boundary conditions
                   T(1,:) = 600; % Top temperature
                   T(end,:) = 900; % Bottom temperature
                   T(:,1) = 400; % Left temperature
                   T(:,end) = 800;  % Right temperature
% Edge Refining 
                   T(1,1) = (600+400)/2 ;
                   T(1,end) = (600+800)/2 ;
                   T(end,1) = (900+400)/2 ;
                   T(end,end) = (900+800)/2;
                   
% optimal value of over relaxation parameter (omega);
                   omega = 1.1
% calculation of temperature distribution
                   k1 =  (alpha*dt)/(dx^2);
                   k2 =  (alpha*dt)/(dy^2);
                   [xx, yy] = meshgrid(x,y);
                   Told = T;
                   T_prev_dt = Told;
                   
counter = 1;
iterative_solver = input('Enter the number for the dedicated solver type: n 1 for jacobi n 2 for gauss seidel n 3 for SOR n')
error = 1e9
tol = 1e-4
tic
   for k = 1:10000
        % Jacobi Iteration scheme
   if iterative_solver == 1
        term1 = (1+(2*k1)+(2*k2))^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9
      while (error > tol)    
            for i = 2:nx-1
              for j = 2:ny-1    
                  H=(Told(i-1,j)+Told(i+1,j));
                  V=(Told(i,j-1)+Told(i,j+1));
                  T(i,j)=((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
              end
            end
          error= max(max(abs(Told-T)))
          Told = T
          counter = counter + 1   
      end
   end
   
      % Gauss seidal scheme
   if iterative_solver == 2
        term1 = (1+(2*k1)+(2*k2))^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9
      while (error > tol)    
            for i = 2:nx-1
              for j = 2:ny-1    
                  H=(T(i-1,j)+Told(i+1,j));
                  V=(T(i,j-1)+Told(i,j+1));
                  T(i,j)=((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
              end
            end
          error= max(max(abs(Told-T)))
          Told = T
          counter = counter + 1
      end
   end
          
    % Successive over-relaxation scheme
   if iterative_solver == 3
        term1 = (1+(2*k1)+(2*k2))^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9
      while (error > tol)    
            for i = 2:nx-1
              for j = 2:ny-1    
                  H=(T(i-1,j)+Told(i+1,j));
                  V=(T(i,j-1)+Told(i,j+1));
                  T(i,j)=(omega)*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)) + (1-omega)*(T_prev_dt(i,j))
              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 = ["Jacobi","Gauss Seidel","Successful Over Relaxation"]; 
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(iterative_solver), k*dt, counter, execution_time));

             


RESULTS & CONCLUSION FOR UNSTEADY STATE 2D HEAT CONDUCTION(IMPLICITLY):

1. Jacobi Method

               

2. Gauss seidel Method

3. Successive over Relaxation

Here the timetaken for the iterations is very large because of the temperature at the time n+1 depends implicitly on the temperature at new time step n+1 at different nodes. Here this is relatively complicated and slow to compute the temperature From the above plots it is clear that the temperature is a function of time. 

STEADY SOLUTION:

Derivation of numerical scheme:

From the generalized Energy equation in 3 dimensions we make the following assumptions

          (i) convection and Radiation effects doesnot occurs

         (ii) No internal heat generation or dissipation

         (iii) For steady state, the `Tne f(time)` then 

                                   `{∂T}/{∂t} =0`

          Hence for 2D the equation becomes

                          `({ ∂^2 T}/{∂x^2} + {∂^2 T}/{∂y^2} ) = 0`

         The second order of spatial terms is discretized by central differencing

                  `{∂^2 T}/{∂x^2} = (T_(i-1,j) -2*(T)_(i,j)+T_(i+1,j) )/{Δx^2}`

                   `{∂^2 T}/{∂y^2} = (T_(i,j-1) -2*(T)_(i,j)+T_(i,j+1) )/(Δy^2)`

         Hence we have

`((T_(i-1,j) - 2*(T)_(i,j)+T_(i+1,j) )/{Δx^2} + (T_(i,j-1) - 2*(T)_(i,j)+T_(i,j+1) )/{Δy^2}) = 0`

Also Δx = Δy

 The other parameters are assumed as:

              Time step, dt = 0.01 seconds

               number of grid points, nx =ny = 41

               size of the grid, dx=dy= 0.0250

               tolerance value = `10^-4`

SOLVING THE STEADY STATE PROBLEM IMPLICITLY:

MATLAB CODE:

%% Matlab program to solve 2D steady state Heat conduction equation implicitly

close all
clear all
clc

% Inputs
            % number of grid points in x direction, nx         
                 nx = 41
            % number of grid points in y direction, ny
                 ny = nx
            % thermal diffusivity, alpha
                alpha = 1e-4
            % Initializing time to start in seconds 
                t_start = 0
            % Initializing time to end in seconds
                t_end = 10000
            % time step in seconds
               dt = 0.01
               
% calculation
          % x array along the length of the domain in x direction
               x = linspace(0,1,nx)
          % y array along the length of the domain in y direction
               y = linspace(0,1,ny)
          % Grid spacing in x direction, dx 
               dx = 1/(nx-1)
          % Grid spacing in y direction, dx = 1/(nx-1)
               dy = 1/(ny-1)
          % No. of time steps
               nt = [t_start:dt:t_end]
% Initializing the temperatures in x and y direction
                  T = 303*ones(nx,ny)
% Assigning Boundary conditions
                   T(1,:) = 600; % Top temperature
                   T(end,:) = 900; % Bottom temperature
                   T(:,1) = 400; % Left temperature
                   T(:,end) = 800;  % Right temperature
% Edge Refining 
                   T(1,1) = (600+400)/2 ;
                   T(1,end) = (600+800)/2 ;
                   T(end,1) = (900+400)/2 ;
                   T(end,end) = (900+800)/2;
                   
% optimal value of over relaxation parameter (omega);
                   omega = 1.1
% calculation of temperature distribution
                   k = 2*(1/(dx^2)+ 1/(dx^2));
                   [xx, yy] = meshgrid(x,y);
                   Told = T;
                   T_prev_dt = Told;
                   

                   iterative_solver = input('Enter the number for the dedicated solver type: n 1 for jacobi n 2 for gauss seidel n 3 for SOR n')
error = 1e9
tol = 1e-4
tic
        % Jacobi Iteration scheme
   if iterative_solver == 1
        jacobi_iter = 1;
      while (error > tol)    
            for i = 2:nx-1
              for j = 2:ny-1    
                  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
          jacobi_iter = jacobi_iter+1;    
       end
   end
    
      % Gauss seidal scheme
   if iterative_solver == 2
          gs_iter = 1;
      while (error > tol)    
            for i = 2:nx-1
              for j = 2:ny-1    
                 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)= term1+term2;
              end
            end
          error= max(max(abs(Told-T)))
          Told = T
          gs_iter = gs_iter+1; 
      end
   end
          
    % Successive Over Relaxation scheme
   if iterative_solver == 3
       sor_iter = 1
      while (error > tol)    
            for i = 2:nx-1
              for j = 2:ny-1    
                 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)= omega*(term1+term2)-(omega-1)*Told(i,j);
              end
            end
          error= max(max(abs(Told-T)))
          Told = T
          sor_iter = sor_iter+1;
      end
   end
 time_required = toc;

figure(1)
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if iterative_solver ==1
    title(sprintf('Temperature profile for steady state Heat conduction using Jacobi Method n No. of Iterations:%d n computation time:%f s',jacobi_iter, time_required));
else
    if iterative_solver == 2
    title(sprintf('Temperature profile for steady state Heat conduction using Gauss siedel Method n No. of Iterations:%d n computation time:%f s',gs_iter, time_required));
    else
       if iterative_solver == 3
       title(sprintf('Temperature profile for steady state Heat conduction using Successive over Relaxation Method n No. of Iterations:%d n computation time:%f s',sor_iter, time_required));
        end
    end
end
xlabel('xx');
ylabel('yy');

RESULTS AND CONCLUSION FOR STEADY STATE 2D HEAT CONDUCTION:

1. Jacobi Method  

2. Gauss seidel Method

       

3. Successive Over Relaxation Method

       

Jacobi method converged after 3187 iterations and Gauss-seidel method took 1711 iterations to converge where as successive over relaxation runs faster and took 1427 iterations to converge. The optimum value of overrelaxation constant omega is taken as 1.1 for Successive over Relaxation method, which gives quick convergence.

The contour plot show  the isothermal lines, If the isotherms are more close and more heat is transferred, if the isotherms are less close then less heat is transferred.

 


Projects by Shravankumar Nagapuri

Objectives: To perform a parametric study on the gate valve simulation. To obtain the mass flow rates at the outlet for 5 design points. To show the cut section view for different design points, and also to show the gate disc lift and fluid volume extraction. To s Read more

Objectives: To Write a Matlab program that can generate the computational blockMesh file  automatically for a wedge angle of 3 degrees and compare the results obtained for Symmetry BC and Wedge BC. To Write a Matlab program that takes an angle as input and gene Read more

1Q. Briefly explain about the possible types of combustion simulations in FLUENT. A:  Combustion or burning is a high temperature exothermic chemical reaction between a fuel and an oxidant accompanied by the production of heat, light and unburnt gases in the fo Read more

Cyclone Separator: Principle and Working: A high speed rotating (air)flow is established within a cylindrical or conical container called a cyclone. Air flows in a spiral pattern, beginning at the top (wide end) of the cyclone and ending at the bottom (narrow) end bef Read more

Objective: To analyse the flow pattern of the fluid inside the gearbox, for two different clearances of the same geometry and each geometry flow is analyzed by two different fluids (Water and Oil) as lubricants. Gearbox sloshing effect:  Slosh refers Read more

Objectives: To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and Grading scheme. To calculate, length of the pipe using the entry length formula for laminar flow through a pipe To show that entry length is suf Read more

Objectives: To find out the maximum temperature attained by the processor. To Prove that the simulation has achieved convergence with appropriate images and plots. To Find out the heat transfer coefficient at appropriate areas of the model. To Identify potential h Read more

Finite volume Method: The finite volume method (FVM) is a method for representing and evaluating partial differential equations in the form of algebraic equations. Similar to the finite difference method (FDM) or finite element method (FEM), values are calculated at Read more

Q1. What are some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves? In your own words, explain how these mathematical models have been adapted for CFD calculations  Kelvin–Helmholtz instability: The Kelvin Read more

Objectives: To Generate the BlockMesh file for the given Geometry using OpenFOAM and use the icoFoam solver to simulate the flow through a backward-facing step. To create multiple meshes and will be comparing the results obtained from each mesh. To show the velocit Read more


Loading...

The End