STUDY ON BEHAVIOR OF ITERATIVE TECHNIQUES FOR TWO DIMENSIONAL CONDUCTION DISTRIBUTION

DESCRIPTION:

In this code, we are going to solve the two dimesnional convection equation defining a temperture distribution. The equation can be solved by 2 methods.

1. Implicit method

2. Explicit method

In explicit method, we will be solving the equation sequencially for the successive term by using only the nth terms  in the RHS. i.e.., by using the known values only

In Implicit method, we will be solving the equation sequencially for the successive terms by using the values of surrounding terms where some will be known and some will be unknown values.

Hence, solving is done by assuming some values initially and solving it iteratively to get the best values from it.

There are three iterative solvers. They are,

1. Jacobi solver

2. Guass seidel solver

3. successive over / under relxation solver

 

MAIN CODE:

clear all
close all
clc

%input data
nx = 10;
ny = 10;
x = linspace(0,1,nx);
dx = x(2) - x(1);
y = x;
dy = dx;

%Initialising the temperature array
temp = ones(nx,ny);

%Boundary conditions
temp(1,:) = 400;
temp(nx,:) = 900;
temp(:,1) = 400;
temp(:,ny) = 800;

%Define state and solver

%state =1 --- STEADY STATE
%state =2 --- UNSTEADY STATE

%solver =1--- IMPLICIT-JACOBI METHOD
%solver =2--- IMPLICIT-GUASS SEIDEL METHOD
%solver =3--- IMPLICIT-SOR METHOD
%solver =4--- EXPLICIT METHOD

state = 2;
solver = 4;

sor = 1.1;      %Relaxation factor
alpha = 1.4;    %Thermal diffusivity

%make a copy
tempold = temp;

%error &tol
error = 90;
tol = 1e-4;
iter = 1;
k =(2/(dx^2))+(2/(dy^2));


%steady state
if state == 1
 %jacobi
      if solver == 1
          disp('Steady state -- Implicit -- Jacobi method')
          
          while error > tol
              for i = 2:nx-1
                  for j = 2:ny-1
                      %central differencing method
                      temp(i,j) = ((tempold(i-1,j) + tempold(i+1,j) + tempold(i,j-1) + tempold(i,j+1))/k)/(dx^2);
                  end
              end
              
              error = max(max(abs(tempold-temp)));
              tempold = temp;
              contourf(x,y,temp,'ShowText','on')
              colorbar
              title(sprintf('STEADY STATE - JACOBI METHOD  n ITERATION : %d',iter-1))
              iter = iter + 1;
          end
      end


      %gauss seidel method
      if solver == 2
         disp('Steady state -- Implicit -- Guass seidel method')
          
          while error > tol
              for i = 2:nx-1
                  for j = 2:ny-1
                      %central differencing method
                      temp(i,j) = ((temp(i-1,j) + tempold(i+1,j) + temp(i,j-1) + tempold(i,j+1))/k)/(dx^2);
                  end
              end
              
              error = max(max(abs(tempold-temp)));
              tempold = temp;
              contourf(x,y,temp,'ShowText','on')
              colorbar
              title(sprintf('STEADY STATE - GUASS SEIDEL METHOD  n ITERATION : %d',iter-1))
              iter = iter + 1;
          end
      end 
      
  %successive over iteration method (sor >=1)
      
     if solver == 3
         disp('Steady state -- Implicit -- SOR method')
         
          while error > tol
              for i = 2:nx-1
                  for j = 2:ny-1
                      %central differencing method
                      temp(i,j) = (tempold(i,j)*(1-sor)) + (((temp(i-1,j) + tempold(i+1,j) + temp(i,j-1) + tempold(i,j+1))/k)/(dx^2))*(sor);  
                  end
              end
             
              error = max(max(abs(tempold-temp)));
              tempold = temp;
              iter = iter + 1;
              
          end
          
           contourf(x,y,temp,'ShowText','on') %contour plots lines ; contourf plots shaded areas
           colorbar
           title(sprintf('STEADY STATE - SOR METHOD  n ITERATION : %d',iter-1))
     end
    
     
end
     

%transient state
if state == 2
        
        %transient state inputs
        time = 0.6;
        dt = 2.1e-3;
        nt = time/dt;
        k1 = (alpha*dt)/(dx^2);
        k2 = (alpha*dt)/(dy^2);
        kf = 1/(1+ 2*k1 + 2*k2);
        temp_dt = tempold;
        
        cfl = alpha*dt*((1/(dx*dx))+(1/(dy*dy)));

        %implicit jacobi method
        if solver == 1
        disp('Transient state -- Implicit -- Jacobi method')
        
        %timeloop
        for timestep = 1:nt
            error = 1e9;
            
            %convergence loop
            while error > tol            
                
                for i = 2 : nx-1
                    for j = 2 : ny-1

                        temp(i,j) = kf*(temp_dt(i,j) + (k1*(tempold(i-1,j) + tempold(i+1,j))) + (k2*(tempold(i,j-1) + tempold(i,j+1))));
                    end
                end
                error = max(max(abs(tempold-temp)));
                tempold = temp;
                iter = iter + 1;                        
            end
        temp_dt = tempold;
        end

    contourf(x,y,temp,'ShowText','on')
    title(sprintf('TRANSIENT STATE - JACOBI METHOD  n ITERATION : %d',iter-1))
    colorbar  
    
        end
        
    %implicit - gauss seidel method 
      if solver == 2
          disp('Transient state -- Implicit -- Guass seidel method')
          
          for timestep = 1:nt
              error = 1e9;
             
              while error > tol
                  for i = 2:nx-1
                      for j = 2:ny-1
                          %central differencing method
                          temp(i,j) = kf*(temp_dt(i,j) + (k1*(temp(i-1,j) + tempold(i+1,j))) + (k2*(temp(i,j-1) + tempold(i,j+1))));
                      end
                  end

                  error = max(max(abs(tempold-temp)));
                  tempold = temp;

                  iter = iter+1;
              end
          temp_dt = tempold;
          end
          
      contourf(x,y,temp,'ShowText','on')
      colorbar
      title(sprintf('TRANSIENT STATE - GAUSS SEIDEL METHOD  n ITERATION : %d',iter-1))
      
      end
      
      if solver == 3
          disp('Transient state -- Implicit -- SOR method')
          
          for timestep = 1:nt
              error = 1e9;
              
              while error > tol
                  for i = 2:nx-1
                      for j = 2:ny-1
                          %central differencing method
                          temp(i,j) = ((1-sor)*tempold(i,j)) + (sor*(kf*(temp_dt(i,j) + (k1*(temp(i-1,j) + tempold(i+1,j))) + (k2*(temp(i,j-1) + tempold(i,j+1))))));
                      end
                  end

                  error = max(max(abs(tempold-temp)));
                  tempold = temp;
                  iter = iter+1;

              end

              temp_dt = tempold;
          end
          
          contourf(x,y,temp,'ShowText','on')
          colorbar
          title(sprintf('TRANSIENT STATE - SOR METHOD  n ITERATION : %d',iter-1))
          
      end


      if solver == 4
          disp('Transient state -- Explicit method')
          for timestep = 1:nt
              error = 1e9;
              
             
                  for i = 2:nx-1
                      for j = 2:ny-1
                          %central differencing method
                          temp(i,j) = tempold(i,j) + (k1*(temp(i-1,j) + temp(i+1,j) + temp(i,j-1) + temp(i,j+1) - (4*tempold(i,j))) );
                      end
                  end
                  error = max(max(abs(tempold-temp)));
                  tempold = temp;
                  iter = iter+1;
              end                     
               contourf(x,y,temp,'ShowText','on')
               title(sprintf('TRANSIENT STATE - EXPLICIT METHOD n ITERATION : %d',iter-1))
               colorbar 
      end
           
end

PLOT:

 

1. STEADY STATE:

 

1.1 JACOBI SOLVER:

 

1.2 GUASS SEIDEL SOLVER:

 

1.3 SOR SOLVER:

 

 

2. TRANSIENT STATE:

 

2.1 JACOBI SOLVER:

 

2.2 GUASS SEIDEL SOLVER:

 

2.3 SOR SOLVER:

 

2.4 EXPLICIT METHOD:

 

 

TIME-STUDY:

Based on the number of iterations on each method, we can compare them for computational time as the iterations are propotional to the computational time.

 

STEADY STATE:

 

In Steady state, we can observe that the Computational time for SOR method is the least and Jacobi was the greatest of all the three methods.

 

JACOBI METHOD: 214

GUASS-SEIDEL METHOD: 114

SOR METHOD: 95

 

 

TRANSIENT STATE:

 

In Transient state, we can observe that the Computational time for Explicit method is the least and Jacobi was the greatest of all the three methods. Also, SOR method is the next least one.

 

JACOBI METHOD: 3455

GUASS-SEIDEL METHOD: 2855

SOR METHOD: 2452

EXPLICIT METHOD: 2440

 

 

STABILITY ANALYSIS:

 

The Courant–Friedrichs–Lewy condition is a necessary condition for convergence while solving certain partial differential equations numerically. To determine the stability of the simulation, we can use courant's number.

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

 

1. EXPLICIT METHOD:

 

When Courant number = 0.4763, we can see that the solution was stable.

 

 

When Courant number = 0.5171, we can observe that the solution has been blown up.

 

So, this states the importance of Courant number in determining the stability of the simulation.

 

IMPLICIT METHOD:

 

When Courant number = 0.5171, We can see the implicit solution doesnt blow up. Also, for any of the Courant number the implicit solution is highly stable. Hence, we can state that when we are simulating for observing the effect of timestep, then we can go with explicit method and when we are going to run for the steady state, then we can go for Implicit method.


Projects by Yokesh R

Demo
Yokesh R · 2020-03-04 11:49:40

Hello Read more

Mixing tee project
Yokesh R · 2019-08-22 06:06:22

Case-1:       CASE-2:           CASE-3:         Read more

MAIN CODE: clear all close all clc A=[5 1 2 ; -3 9 4; 1 2 -7]; B=[10 ;-14; 33]; u = [0 1 2; 0 0 4; 0 0 0]; l = [0 0 0; -3 0 0; 1 2 0]; di = [5 0 0; 0 9 0; 0 0 -7]; mag = [0.25,0.5,0.75,1,1.5,2]; for i = 1 : length(mag) d = di*mag(i); Read more

DESCRIPTION: Using the central and skewed scheme methods initially, we have to develop the Taylor table, which should be converted into a matrix and solved for the constants. These constants help in solving the second order differential equation. Finally, the error res Read more

DESCRIPTION: In this code, we are going to compare the results of the 1D convection distribution equation describing the velocity for a range of time-step values. As a result, we will be getting the plots on the final velocity profile on time marching and computational Read more

DESCRIPTION:      In this code, we are going to discretize the function for the first-order derivative using first, second and fourth order approximation techniques and studying the error bounced by every method by a bar chart.   CODE: clear all Read more

DESCRIPTION:      Here in this code, we are using a range of dx terms for discretizing the first order derivative and studying the results of the dx vs error plot  CODE: clear all close all clc x = pi/3; dx = linspace(pi/4,pi/4000,30); %y Read more

STOICHIOMETRIC COMBUSTION   DESCRIPTION:   Stoichiometric combustion defines the ideal combustion that needs to take place for a given chemical equation. This is actually the general form of equation for alkanes. The general form might differ for e Read more

DESCRIPTION:         In this report, we will be studying the effects of the grid points in the first order convection equation describing one dimensional velocity distribution.   `frac{delu}{delt} = c*frac{delu}{delx}`   CODE: cl Read more


Loading...

The End