Solving Steady and Unsteady 2-D Heat Conduction Problem for a Square Plate

Aim:    To solve steady and unsteady 2D heat conduction problem

Objective: 

  1.  To use iterative solvers - Jacobi Method, Gauss-Seidel Method and Successive Over Relaxation(SOR) Method, to solve steady and unsteady 2D conduction equation
  2. To use Explicit Method for unsteady state and see the effect of time step on the solution and analyse CFL condition
  3. To use higher value of time step to check its stability on unsteady implicit solutions solved by the iterative solvers. 

 

Steady State Condition

The staedy state equation for 2D heat conduction is by:

 `alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))+f=0`

In the above equation `alpha` is thermal diffusivity, T is temperature and f is the source term.

For the convenience source term is neglected so the final equation is:

`alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

Finite difference equation using FTCS(Forward time Central Scheme) is given by:

`(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(/_\x)^2+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(/_\y)^2=0`

Unsteady State Condition

Unsteady state equation for 2D heat conduction is :

`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))+f`

t here is time

Again here f is neglected so equation is :

`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`

Finite difference equation is:

`(T_(i,j)^(n+1)-T_(i,j)^n)/(/_\t)=(T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(/_\x)^2+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(/_\y)^2`

let `k1=(alpha/_\t)/(/_\x)^2 and k2=(alpha/_\t)/(/_\y)^2`

`T_(i,j)^(n+1)-T_(i,j)^n=k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))`

For Explicit Approach the above equation is :

`T_(i,j)^(n+1)-T_(i,j)^n=[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^n`

`T_(i,j)^(n+1)=T_(i,j)^n+[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^n`

where n is the time index and i,j are the space indices in x-direction and y directions respectively 

For Implicit Approach the equation is:

`T_(i,j)^(n+1)-T_(i,j)^n=[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^(n+1)`

`T_(i,j)^(n+1)=T_(i,j)^n+[k1(T_(i+1,j)-2T_(i,j)+T_(i-1,j))+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))]^(n+1)`

On rearranging it becomes:

`T_(i,j)^(n+1)+2k1T_(i,j)^(n+1)+2k2T_(i,j)^(n+1)=T_(i,j)^n+[k1(T_(i+1,j)+T_(i-1,j))+k2(T_(i,j+1)+T_(i,j-1))]^(n+1)`

`T_(i,j)^(n+1)=(T_(i,j)^n+[k1(T_(i+1,j)+T_(i-1,j))+k2(T_(i,j+1)+T_(i,j-1))]^(n+1))/(1+2k1+2k2)`

n is the time index here and i,j are space indices

 

Main Program for Steady State Solutions

 

%% To solve Steady 2D Conduction problem
% d^2T/dx^2+d^2T/dy^2=0

clc
close all
clear

nx=10;             % Number of grid points in X-direction
ny=10;             % Number of grid points in Y-direction
T=300*ones(nx,ny); % Room Temp = 300K Taken
% Creating grids along X and Y directions
x=linspace(0,1,nx);
y=linspace(0,1,ny);
% Reversing the direction of y vector
y_rev=flip(y);

% Creating Mesh
[xn,yn]=meshgrid(x,y_rev);

% Assigning BC
T_lft=400;
T_top=600;
T_rgt=800;
T_bot=900;
T(2:end-1,1)=T_lft;
T(2:end-1,end)=T_rgt;
T(1,2:end-1)=T_top;
T(end,2:end-1)=T_bot;
T(1,1)=(T_lft+T_top)/2;
T(end,1)=(T_lft+T_bot)/2;
T(1,end)=(T_rgt+T_top)/2;
T(end,end)=(T_rgt+T_bot)/2;

% Initializing values
tol=1e-4;   % Convergence Criterion
T_old=T;
T_initial=T;

for iterative_solver=1 :3
    
  if iterative_solver==1
       err_jac=1;
       iter_jac=0;
     while(err_jac>tol)
        
        for j=2:ny-1
            for i=2:nx-1
             T(i,j)=(T_old(i-1,j)+T_old(i+1,j)+T_old(i,j-1)+T_old(i,j+1))/4;   % since dx=dy
            end                      
        end
            err_jac=max(max(abs(T-T_old)));
            iter_jac=iter_jac+1;
            figure(1)
            [c,d]=contourf(xn,yn,T);
            colormap(jet);
            colorbar;
            txt=sprintf('Jacobi Method \nIteration Number= %d',iter_jac);
            title(txt);
            xlabel('X-direction along square plate\rightarrow');
            ylabel('Y-direction along square plate\rightarrow');
            clabel(c,d);
            pause(0.003);
            
            T_old=T; % Updating old Temperature values
      end
  end
 
  if iterative_solver==2
      err_gs=1;
      iter_gs=0;
    while(err_gs>tol)
        
        for j=2:ny-1
            for i=2:nx-1
                
                T(i,j)=(T(i-1,j)+T_old(i+1,j)+T(i,j-1)+T_old(i,j+1))/4;
            end
        end
        
            err_gs=max(max(abs(T-T_old)));
            iter_gs=iter_gs+1;
            figure(2)
            [c,d]=contourf(xn,yn,T);
            colormap(jet);
            colorbar;
            txt=sprintf('Gauss-Seidel Method \nIteration Number= %d',iter_gs);
            title(txt);
            xlabel('X-direction along square plate\rightarrow');
            ylabel('Y-direction along square plate\rightarrow');
            clabel(c,d);
            pause(0.003);
            
            T_old=T; % Updating old Temperature Values
    end
  end

 if iterative_solver==3
    % Scaling Factor Omega for SOR Method
    omega=0.5:0.5:2;
    
    for k=1:length(omega)
        err_sor=1;
        iter_sor=zeros(1,4);
        T=T_initial;
        T_old=T_initial;
        while(err_sor>tol)
        
            for j=2:ny-1
               for i=2:nx-1
                
                T_gs=(T(i-1,j)+T_old(i+1,j)+T(i,j-1)+T_old(i,j+1))/4;
                T(i,j)=(1-omega(k))*T(i,j)+omega(k)*(T_gs);
               end
            end
             
              err_sor=max(max(abs(T-T_old)));
              iter_sor(k)=iter_sor(k)+1;
              figure(k+2)
              [c,d]=contourf(xn,yn,T);
              colormap(jet);
              colorbar;
              txt=sprintf('Succesive Over Relaxation(SOR) Method \n Iteration Number= %d\n Omega value=%0.1f',iter_sor(k),omega(k));
              title(txt);
              xlabel('X-direction along square plate\rightarrow');
              ylabel('Y-direction along square plate\rightarrow');
              clabel(c,d);
              pause(0.0003);
              
              T_old=T; % Updating Old Temperature values
            
       end
    end
 end
 
 % Initializing T values;
      T_old=T_initial;
      T=T_initial;
end

In the program number of grid points taken is 10 in both X and Y directions along the square plate domain. For SOR method omega (scaling factor) taken as array of (0.5,1,1.5 and 2).

In steady state conditions we are using iterative solvers to converge the solution such that the error falls below the tolerance level.

Plots for Steady State Conditions

1. From Jacobi Method

2. From Gauss-Seidel Method

3. From Successive Over Relaxation(SOR) Method

a) For Omega value equal to 0.5

b) For Omega value equal to 1

c) For Omega value equal to 1.5

d) For Omega value equal to 2

Plots Analysis of Steady State condition

  • All the above plots show same results except for the result from SOR Method (Omega Value=2).
  • The number of iterations in case of Jacobi method is 207,for Gauss-Seidel it is 111 which is less than the Jacobi method.So we can say that Gauss-Seidel method is more efficient than Jacobi method.
  • In case of SOR method it is observed that for Omega value< 1(Under-Relaxation), the solution took more number of iterations to converge than at Omega value =1. Moreover the SOR method with omega value =1 is same as Gauss-Seidel method and so their number of interations are also same i.e 111.
  • SOR method further proves to be more efficient than Gauss-Seidel method at Omega value>1(Over-Relaxation). For Omega value=1.5 it gives optimum solution with least number of iterations to reach the final solution. However at Omega value>1.5 the number of iterations again increases and at Omega value=2 the solutions blows up.

 

Main Program for Unsteady state Solutions

 

%% To solve Unsteady 2D conduction problem
% dT/dt=alpha(d^2T/dx^2+d^2T/dy^2)

clc
close all
clear
% Initializing Parameters
nx=10;
ny=10;
T=300*ones(nx,ny); % Room Temp = 300K Taken
time=0.3;
alpha=1;
x=linspace(0,1,nx);
y=linspace(0,1,ny);

dx=x(2)-x(1);
dy=y(2)-y(1);
% Time step size 1e-2 and 1e-3
dt=[1e-3 1e-2];
%% CFL calculation

CFL=alpha*(dt/dx^2+dt/dy^2);

% Number of time steps and constants k1 and k2 for dt=1e-3 and 1e-2
nt=1+(time./dt);
k1=dt/dx^2;
k2=dt/dx^2;
% Reversing the direction of y vector
y_rev=flip(y);

% Creating Mesh
[xn,yn]=meshgrid(x,y_rev);

% Assigning BC
T_lft=400;
T_top=600;
T_rgt=800;
T_bot=900;
T(2:end-1,1)=T_lft;
T(2:end-1,end)=T_rgt;
T(1,2:end-1)=T_top;
T(end,2:end-1)=T_bot;
T(1,1)=(T_lft+T_top)/2;
T(end,1)=(T_lft+T_bot)/2;
T(1,end)=(T_rgt+T_top)/2;
T(end,end)=(T_rgt+T_bot)/2;

% Initializing Variables

tol=1e-4; %  Convergence Criteria
T_old=T;
T_initial=T;
T_prev_dt=T;

     

     %% ######## Explicit Direct Method #######  
     
     
   % Loop for each time step 1e-2 and 1e-3     
  for m=1:length(dt)  
      
    iter_explicit=0;
    % Time Loop
    for t=1:nt(m)
    % Space Loop             
        for j=2:ny-1    
            for i=2:nx-1
              term1=1-2*k1(m)-2*k2(m);
              H=k1(m)*(T_old(i-1,j)+T_old(i+1,j));
              V=k2(m)*(T_old(i,j-1)+T_old(i,j+1));
              T(i,j)=T_old(i,j)*term1+H+V;
            end              
        end
            iter_explicit=iter_explicit+1;
        
            figure(m)
            [c,d]=contourf(xn,yn,T);
            colormap(jet);
            colorbar;
            txt=sprintf('Explicit Method \nIteration Number= %d\n\\Deltat=%0.4f',iter_explicit,dt(m));
            title(txt);
            xlabel('X-direction along square plate\rightarrow');
            ylabel('Y-direction along square plate\rightarrow');
            clabel(c,d);
            pause(0.003);
            % Updating Temperature values at the end of Time loop 
        T_old=T;  
    end
  end
  
  
%% ############  Implicit Method   ###########


% Specific dt values taken=1e-2 and corresponding nt to check whether solution blows or not 
 dt_imp= 1e-2;
 nt_imp=1+time/dt_imp;
 
 % Corresponding k1=alpha*dt/dx^2  and k2=alpha*dt/dy^2
 k1=(alpha*dt_imp)/dx^2;    
 k2=(alpha*dt_imp)/dy^2;
   
for iterative_solver=1:3
  
    % Initialising Temp at the beginning
    T=T_initial;
    T_old=T_initial;
    T_prev_dt=T_initial;
    
    if iterative_solver==1

   %% Jacobi Method ##########
   
     iter_jac=0;
   % Time Loop 
       for t=1:nt_imp
         err_jac=1;
        % Convergence Loop  
          while(err_jac>tol)
            % Space Loop
             for j=2:ny-1    
               for i=2:nx-1
                 term1=1+2*k1+2*k2;
                 H=k1*(T_old(i-1,j)+T_old(i+1,j));
                 V=k2*(T_old(i,j-1)+T_old(i,j+1));
                 T(i,j)=(T_prev_dt(i,j)+H+V)/term1;
               end              
             end
             iter_jac=iter_jac+1;
             err_jac=max(max(abs(T-T_old)));
        
             figure(3)
             [c,d]=contourf(xn,yn,T);
             colormap(jet);
             colorbar;
             txt=sprintf('Jacobi Method \nIteration Number= %d\n\\Deltat=%0.3f',iter_jac,dt_imp);
             title(txt);
             xlabel('X-direction along square plate\rightarrow');
             ylabel('Y-direction along square plate\rightarrow');
             clabel(c,d);
             pause(0.0003);
             % Updating Temperature values at the end of Convergence loop 
             T_old=T;  
          end
          % Updating previous time step Temperature values 
          T_prev_dt=T;
     
       end
   end

    
   if iterative_solver==2
    %% Gauss-Seidel Method ###   
       
     iter_gs=0;
      % Time Loop 
      for t=1:nt_imp
         err_gs=1;
         % Convergence Loop
          while(err_gs>tol)
        
             for j=2:ny-1    
                for i=2:nx-1
                  term1=1+2*k1+2*k2;
                  H=k1*(T(i-1,j)+T_old(i+1,j));
                  V=k2*(T(i,j-1)+T_old(i,j+1));
                  T(i,j)=(T_prev_dt(i,j)+H+V)/term1;
                end              
             end
            
          iter_gs=iter_gs+1;
          err_gs=max(max(abs(T-T_old)));
        
          figure(4)
          [c,d]=contourf(xn,yn,T);
          colormap(jet);
          colorbar;
          txt=sprintf('Gauss-Seidel Method \nIteration Number= %d\n\\Deltat=%0.3f',iter_gs,dt_imp);
          title(txt);
          xlabel('X-direction along square plate\rightarrow');
          ylabel('Y-direction along square plate\rightarrow');
          clabel(c,d);
          pause(0.0003);
            % Updating Temperature values at the end of Convergence loop 
          T_old=T;  
         end
     % Updating previous time step Temperature values 
         T_prev_dt=T;
      end
   end

  if iterative_solver==3
  %% SOR Method ####
  
  % Scaling Factor Omega for SOR Method
    omega=0.8:0.1:1.2;
    
    for k=1:length(omega)
       
        % Initializing for each omega value
        iter_sor=zeros(1,5);
        T=T_initial;
        T_old=T_initial;
        T_prev_dt=T_initial;
      % Time Loop  
        for t=1:nt_imp
            
          err_sor=1;
          % Convergence Loop
          while(err_sor>tol)
            % Space Loop
              for j=2:ny-1
                 for i=2:nx-1
                   term1=1+2*k1+2*k2;
                   H=k1*(T(i-1,j)+T_old(i+1,j));
                   V=k2*(T(i,j-1)+T_old(i,j+1));
                   T_gs=(T_prev_dt(i,j)+H+V)/term1;
                   T(i,j)=(1-omega(k))*T_prev_dt(i,j)+omega(k)*T_gs;
                 end
              end
            
              err_sor=max(max(abs(T-T_old)));
              
              iter_sor(k)=iter_sor(k)+1;
              figure(k+4)
              [c,d]=contourf(xn,yn,T);
              colormap(jet);
              colorbar;
              txt=sprintf('Successive Over Relaxation(SOR) Method \nIteration Number= %d\n\\Deltat=%0.3f\n Omega value=%0.1f',iter_sor(k),dt_imp,omega(k));
              
              title(txt);
              xlabel('X-direction along square plate\rightarrow');
              ylabel('Y-direction along square plate\rightarrow');
              clabel(c,d);
              
              pause(0.0003);
              T_old=T; % Updating Temperature values at the end of Convergence Loop
           end
         T_prev_dt=T; % Updating previous time step Temperature
        end
    end
  end 
end

In the program number of grid points taken is 10 in both X and Y directions along the square plate domain. For Explicit approach time step size used is 1e-2 and 1e-3 to check how the time step affects the solutions.

For Implicit approach time step size used is 1e-2 to check the stability of Implicit method at higher time step and also iterative solvers are used to converge the solutions such that the error falls below the set tolerance criteria.

For Implicit SOR method, Omega (scaling factor) is taken as array of (0.8 to 1.2).

Plots for Un-Steady State Conditions

1) From Explicit Method at time-step size=1e-3

2) From Explicit Method at time-step size=1e-2

Plots Analysis for Explicit Scheme

  • The Explicit method does not use iterative solvers to converge and reach the solution infact it is just one equation which is solved directly. There is no tolerance criteria to converge to hence compared to Implicit Method it is faster and requires less computational time. 
  • At time step size =1e-3, the solution converges to the same steady-state solution.
  • At time step size =1e-2, the solution blows up. This is due to the CFL condition which is given by- `alpha((/_\t)/(/_\x)^2+(/_\t)/(/_\y)^2)<=1/2`
  • Here `/_\x=/_\y`. So CFL condition is actually `alpha(/_\t)/(/_\x)^2<=1/4`
  • The CFL value is calculated for corresponding time step size using matlab and they are:

CFL value for 1e-2 is 1.62 which is greater than 0.25 and so the solution blows up but for 1e-3 =0.162 which is less than 0.25 so solution is stable.

 

 

 

3) From Implicit Scheme

a) Jacobi Method

b) Gauss-Seidel Method

c) SOR at omega value=0.8

d) SOR at omega value=0.9

e) SOR at omega value=1

f) SOR at omega value=1.1

g) SOR at omega value=1.2

 

Plots Analysis for Implicit Schemes

  • Compared to steady state or explicit methods, implicit shemes shows more number of iterations because the coupled equations formed in implicit shemes are very large and so it requires more computational time to show converging results.
  • All implicit scheme methods worked at `/_\t=`0.01 for which explicit scheme crashed. This proves that even at higher value of time step size implicit method is unconditionally stable.
  • Gauss-Seidel Implicit scheme is more efficient than Jacobi or SOR implicit methods with number of iterations=587
  • SOR method under relax condition (Omega value<1) although shows less iterations but the solutions are not same as steady state or explicit condition. The solutions from them are not correct.

SOR method with omega value>1 shows more number of iterations than SOR method with omega value=1. So SOR method does not seem to work well for implicit schemes. Although for some condition it works well when time step size is very small.


Projects by Md Imran Ansari

Aim  :    To simulate flow through a pipe upto the fully developed region Objective :  Take a wedge about the axis of the pipe and use axisymmetry condition to get the desired output Take length of the pipe to be entry length and simulate flow a Read more

In Finite Volume Method(FVM) the values of the variable are stored at the computational nodes (central position) of the control volume. The approximations of surface and volume integrals require values of the variables at the centre of the faces of the control volume. S Read more

Aim:  To simulate the flow through a backward facing step using icoFoam solver   Objective: To create create the geometry as given and then develop mesh for that Use different scale factor to refine mesh near the walls Solve the flow through the domain Read more

Aim: To solve quasi one-dimensional isentropic flow through convergent-divergent nozzle using Mac-Cormack\'s technique Objective: To use conservation and non-conservation forms of partial differential equations to          construct fin Read more

Fourth Order Approximation of Second Derivative of a Function    1. Central Difference Scheme In Central difference scheme the number of nodes required to get fourth order approximation for a second derivative of a function is =5 ( A+O-1) where A is order o Read more


Loading...

The End