Solving A Quasi-1 Dimensional Supersonic Incompressible Isentropic Nozzle Flow Using Maccormack Explict Method In MATLAB

Abstract:

This project presents the Modelling of the  behaviour of  fluid and its properties when it passes through a convergent-divergent nozzle using a Finite Difference Method. However, to simplify the understanding of the dynamics of the flow in the nozzle, only 1-dimention flow variation has been considered where the properties vary only in 1 (x-axis) direction making it "Quasi 1-D subsonic-supersonic isentropic nozzle flow". Here, MATLAB codes are used to understand how the governing equation  of the flow is transformed into a numeric form (Finite Difference Form).

Solver used : Maccormack predictor-corrector explicit finite difference method.

Flow Type   : Incompressible-Unsteady-inviscid(Ideal) and Isentropic Flow.

 

Introduction:

We consider the steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by p0 , and To, respectively. The cross-sectional area of the reservoir is large ( theoretically, A -> infinity), and hence the velocity is very small ( V -> O). Thus, Po and T0 are the stagnation values or total pressure and total temperature, respectively. The flow expands 1sentrop1cally to supersonic speeds at the nozzle exit, where the exit pressure, temperature, velocity and Mach number are denoted by Pe, Te, V and M respectlvely. The flow is locally subsonic in the convergent section if the nozzle, some at the throat (minimum area), and supersonic at the divergent section. The sonic flow (M = 1) at the throat means that the local velocity at this location is equal the local speed of sound. Using an asterisk to denote sonic flow values, we have at the throat V = V* = a*. Similarly, the sonic flow values of pressure and temperature are denoted by p* and T*, respectively. The area of the sonic throat is denoted by A*. We assume that at a given section, where the cross-sectional area is A, the flow properties are uniform across that section. Hence, although the area of the nozzle changes as a function of distance along the nozzle, x, and therefore in reality the flow field is two dimensional (the flow varies in.the two-dimensional .ry space), we make the assumption that the flow properties-vary only with x; this is tantamount to
assuming uniform flow properties across any given cross section. Such flow is
defined as quasi-one-dimensional flow.

 

 

 

The governing continuity, momentum, and energy equations for this quasione-
dimensional, steady, isentropic flow can be expressed, respectively, as

Continuity equation:

 

Momentum equation:

Energy equation:

here 1 and 2 denotes the quantities at different location of the nozzle.

The perfect gas equation stated as:

where

ρ = Density.

R = Universal gas constant.

T = Temperature.

h = Enthalpy.

A = Area of cross-section.

V = Velocity.

(1) Non-Conservation form of the governing equation:

(1.1)The initial conditions and the nozzle shape :-

The nozzle shape is kept constant over time. The parabolic area distribution is given as:-

At x = 1.5, is the throat of the nozzle. The convergent region is at the region in between [0 , 1.5] whereas the divergent section occupies the region in between [1.5 , 3]. The nozzle shape is shown below as:

The initials conditions of the fluid properties are as follows:

ρ = 1 - (0.3146*x)

T = 1 - (0.2314*x)

V = (0.9 + 1.09*x)*T^0.5

 

(1.2) The boundary conditions:-

The subsonic (inlet) flow condition:- Here we choose velocity(V) to float because on a physical basis we know the mass flow through the nozzle must be allowed to adjust to the proper steady state, and allowing V to float makes the most sense as part of this adjustment.

The independent(fixed) parameters are:-

  ρ1 = 1 bar.

  T1 = 1 k

The supersonic (outlet) flow conditions:-Here we choose to keepe all the parameters float. Linear extrapolation based on flow field values is used at inter points.

V(n) = 2V(n-1) - V(n-2)

ρ(n) = 2ρ(n-1) - ρ(n-2)

T(n) = 2T(n-1) - T(n-2)

 

(1.3)The non-dimentionalised non-conservative form of the above governing equations are:-

Continuity equation:

Momentum equation:

 

Equation equation:

  • The prime variables indicates the non-domentionalised form.
  • The above equation is the the non-conservative form of the goerning equation.
  • We perform Maccormack predictor-corrector explicit method to determine the dependent variables.
  • In predictor, forward difference scheme of 1st order is used to determine the variable for the next time step and update it.
  • In corrector, backward difference scheme of 1st order is used to determine the variable and finnaly their average is taken which is 2nd order finite difference accurate. Similarly, for next time step the same process will be implemented.   

Discretization of the PDE using Maccormack explicit method:

Predictor Step (1st order forward difference scheme):

The predicted value of ρ, V ,T  for the next time step is calculated as:

 Corrector step(1st order backward difference scheme):

 

The average time derivative for all the primitive is calculated as:-

Therefore, the corrected values of flow field for [t+△t] is given as:

Hence, from the abeove equation the values of primitive variables can be calculated.

 

(1.4)The non-dimentionalised Conservation form of governing equation:-

Initial Conditions:

 

 

Continuity equation

Momentum Equation

Energy equation

The above equations are the nondimensional conservation form of the continuity, momentum, and energy equations for quasi-one-dimensional flow, respectively.

The equations for quasi-onedimensional flow can be expressed in a generic form. Let us define the elements of the solutions vector "U", the flux vector "F", and the source term "J" as follows:-

Therefore, the above governing equation can be written as:-

Before directly discretizing the above generaic form, we must express the primitive values in terms of the solution vectors and the solution vections has to be converted in terms of the flux terms.

 

MATLAB CODES:-

 

MAIN SCRIPT:

% main script for conservative 1-D quasi supersonic
% nozzle flow.

clear all
close all
clc

% 1-D grid defination.
n     = 31;
x     = linspace(0,3,n);
dx    = x(2)-x(1); % grid mesh.
A     = 1 + 2.2*(x - 1.5).^2;
gamma = 1.4;
c     = 0.5;       % CFL Number.
n_throat = find(A-1<0.001);
%%--------------------------------------------------------------Conservative---------------------------------------------------------------------%%
[rho,V,T,P,mass_flow_c,Mach,rho_throat,T_throat,Mach_throat,P_throat,iter] = conservation_1_D_quasisupersonic_nozzle_flow(n,x,dx,A,gamma,n_throat,c);

fprintf('mim. number of iteration to convergencefor conservation form: %fn',iter);

% Steady-properties of the nozzle.
figure(3)
subplot(4,1,1);
plot(x,rho,'g','linewidth',1);
title('Steady state properties at the nozzle (conservation Form)');
xlabel('x/L');
ylabel('rho/rho_{0}');
legend('Density ratio');


subplot(4,1,2);
plot(x,T,'b','linewidth',1);
xlabel('x/L');
ylabel('Temperature Ratio');
legend('Temperature ratio');


subplot(4,1,3);
plot(x,P,'y','linewidth',1);
xlabel('x/L');
ylabel('Pressure Ratio');
legend('Pressure ratio');


subplot(4,1,4);
plot(x,Mach,'r','linewidth',1);
xlabel('x/L');
ylabel('Mech Number');
legend('Mech');



subplot(4,1,4);
plot(x,V,'r','linewidth',1);
xlabel('x/L');
ylabel('Mech Number');
legend('Mech');

% Transient properties at the throat of nozzle

figure(4);
subplot(4,1,1);
plot(rho_throat);
title('Transient properties variation at the nozzle throat (Conservation Form)');
xlabel('Time Steps');
ylabel('Density');
%ylabel('rho/rho_{0}');
grid on;

% 
subplot(4,1,2);
plot(P_throat);
xlabel('Time Steps');
ylabel('Pressure');
%ylabel('P/P_{0}');
grid on;


subplot(4,1,3);
plot(T_throat);
xlabel('Time Steps');
%ylabel('T/T_{0}');
ylabel('Temperature');
grid on;


subplot(4,1,4);
plot(Mach_throat);
xlabel('Time Steps');
%ylabel('Mach Number');
ylabel('Mach Number');
grid on;


%%-------------------------------------------------------------Non-Conservative---------------------------------------------------------------%%
% main script for non- conservative 1-D quasi supersonic
% nozzle flow.


[rho,V,T,A,P,mass_flow,Mach,rho_throat,T_throat,P_throat,Mach_throat,iter,k1] = Main_Script_Non_conservation(n,x,dx,A,gamma,n_throat,c);

fprintf('mim. number of iteration to convergence for Non-conservation form: %fn',iter);

% Steady-properties of the nozzle.
figure(5)
subplot(4,1,1);
plot(x,rho,'g','linewidth',1);
title('Steady state properties at the nozzle (Non-conservation Form)');
xlabel('x/L');
ylabel('rho/rho_{0}');
legend('Density ratio');
grid minor;


subplot(4,1,2);
plot(x,T,'b','linewidth',1);
xlabel('x/L');
ylabel('Temperature Ratio');
legend('Temperature ratio');
grid minor;


subplot(4,1,3);
plot(x,P,'y','linewidth',1);
xlabel('x/L');
ylabel('Pressure Ratio');
legend('Pressure ratio');
grid minor;


subplot(4,1,4);
plot(x,Mach,'r','linewidth',1);
xlabel('x/L');
ylabel('Mech Number');
legend('Mech');
grid minor;


% Transient properties at the throat of nozzle

figure(6);
subplot(4,1,1);
plot(rho_throat);
title('Transient properties variation at the nozzle throat (Non-Conservation Form)');
xlabel('Time Steps');
%ylabel('rho/rho_{0}');
ylabel('Density');
grid on;

% 
subplot(4,1,2);
plot(P_throat);
xlabel('Time Steps');
%ylabel('P/P_{0}');
ylabel('Pressure');
grid on;


subplot(4,1,3);
plot(T_throat);
xlabel('Time steps');
%ylabel('T/T_{0}');
ylabel('Temperature');
grid on;


subplot(4,1,4);
plot(Mach_throat);
xlabel('Time steps');
%ylabel('Mach Number');
ylabel('Mach Number');
grid on;



figure(7)
plot(x,mass_flow_c)
hold on
plot(x,mass_flow)
xlabel('Length Of Nozzle');
ylabel('Mass Flow Rate');
legend('Conservation Mass Flow Rate','Non-Conservation Mass Flow Rate');

 

USER DEFINED FUNCTIONS:

Conservation Form Script:

% 1-D quasi-Isentropic Supersonic Nozzle flow.
% Conservative form approach.

function [rho,V,T,P,mass_flow_c,Mach,rho_throat,T_throat,Mach_throat,P_throat,iter] = conservation_1_D_quasisupersonic_nozzle_flow(n,x,dx,A,gamma,n_throat,c);



% Initial Conditions.
for i = 1:length(x);
    
    if x(i)>= 0 && x(i)<= 0.5
        
        rho(i) = 1;
        T(i)   = 1;
        
    end
    
    if x(i)>= 0.5 && x(i)<= 1.5
        
        rho(i) = 1 - 0.366*(x(i) - 0.5);
        T(i)   = 1 - 0.167*(x(i) - 0.5);
        
    end
    
    if x(i)>= 1.5 && x(i)<= 3.5
        
        rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
        T(i)   = 0.833 - 0.3507*(x(i) - 1.5);
        
    end
    
end

V = 0.59./(rho.*A);  % Velocity. 
P = rho.*T;          % Pressure.

% Initial Conditions for the solution vectors.
U_1 = rho.*A;
U_2 =(rho.*V).*A;
U_3 = rho.*A.*(((T./(gamma-1)) + (gamma/2)*(V.^2)));

% Initial Conditions for Flux Vectors.
F_1 = U_2;
F_2 = ((U_2.^2)./U_1) + (((gamma-1)/gamma)*(U_3 - ((gamma/2)*(U_2.^2./U_1))));
F_3 = (gamma*(U_2.*U_3)./U_1) - ((gamma*(gamma-1)/2).*((U_2.^3)./U_1.^2));

% % Storing the initial values.
% U_1_old = U_1;
% U_2_old = U_2;
% U_3_old = U_3;

% F_1_old = F_1;
% F_2_old = F_2;
% F_3_old = F_3;


% Initializing Numerical discretization using  Maccormack explicit method.
% Initialzing time loop.
% Initializing predictor-Corrector step.
nt = 5000;
iter = 0;
for k1 =1:nt
   
    iter = iter+1;
    
    U_1_old = U_1;
    U_2_old = U_2;
    U_3_old = U_3;
    
    
    F_1 = U_2;
    F_2 = ((U_2.^2)./U_1) + (((gamma-1)/gamma)*(U_3 - ((gamma*U_2.^2)./(2*U_1))));
    F_3 = (gamma*(U_2.*U_3)./U_1) - ((gamma*(gamma-1)/2)*((U_2.^3)./(U_1.^2)));
   
    % CFL Time loop.
    a = sqrt(T);
    timestep = c*dx./(a + V);
    dt = min(timestep);
    
    % Initializing predictor Step.
    
    for j = 2:n-1
        
        % Source vector.
        J_2(j)      = (1/gamma)*rho(j)*T(j)*((A(j+1) - A(j))/dx);
        
        % Continuity Equation.
        dU1_dt_p(j) = -((F_1(j+1) - F_1(j))/dx);
        
        % Momentum Equation.
        dU2_dt_p(j) = -((F_2(j+1) - F_2(j))/dx) + J_2(j);
        
        % Energy Equation.
        dU3_dt_p(j) = -((F_3(j+1) - F_3(j))/dx);
        
        % Updating the Solution Vactors.
        
        U_1(j) = U_1(j) +  (dU1_dt_p(j)*dt);
        U_2(j) = U_2(j) +  (dU2_dt_p(j)*dt);
        U_3(j) = U_3(j) +  (dU3_dt_p(j)*dt);
        
    end
    
    % Updating premitive variables.
    
    rho = U_1./A;
    V   = U_2./U_1;
    T   = (gamma-1)*((U_3./U_1) - ((gamma/2)*(V.^2)));
    P   = rho.*T;
    
    % Updating Flux vectors.
    F_1 = U_2;
    F_2 = ((U_2.^2)./U_1) + (((gamma-1)/gamma)*(U_3 - ((gamma*U_2.^2)./(2*U_1))));
    F_3 = (gamma*((U_2.*U_3)./U_1)) - ((gamma*(gamma-1)/2)*((U_2.^3)./(U_1.^2)));
    
    
    
    
    
   
    % Initializing corrector step.
    
    for k = 2:n-1
        
        % Source Vector.
        J_2(k)      = (1/gamma)*rho(k)*T(k)*((A(k) - A(k-1))/dx);
        
        % Continuity Equation.
        dU1_dt_c(k) = -((F_1(k) - F_1(k-1))/dx);
        
        % Momentum Equation.
        dU2_dt_c(k) = -((F_2(k) - F_2(k-1))/dx) + J_2(k);
        
        % Energy Equation.
        dU3_dt_c(k) = -((F_3(k) - F_3(k-1))/dx);
        
    end
    
    
    % Computing the Average Time Step.
    
    dU1_dt_avg = 0.5*(dU1_dt_p + dU1_dt_c);
    
    dU2_dt_avg = 0.5*(dU2_dt_p + dU2_dt_c);
    
    dU3_dt_avg = 0.5*(dU3_dt_p + dU3_dt_c);
    
    
    % Computing final solution vectors.
    
   for k = 2:n-1
       
       U_1(k) = U_1_old(k) + (dU1_dt_avg(k))*dt;
    
       U_2(k) = U_2_old(k) + (dU2_dt_avg(k))*dt;
       
       U_3(k) = U_3_old(k) + (dU3_dt_avg(k))*dt;
       
   end
   
   
   % Applying Boundary conditions
   % Inlet conditions.
   
   % Primitive variables.
   rho(1) = 1;   % Fixed Boundary condition.
   T(1)   = 1;   % Fixed Boundary condition.
   V(1)   = 2*V(2) - V(3);
   
   % Solution Vectors.
   U_1(1) = rho(1)*A(1);
   U_2(1) = 2*U_2(2) - U_2(3);
   U_3(1) = U_1(1)*((T(1)/(gamma -1)) + (gamma/2)*V(1)^2);
   
   
   % Outlet conditions.
   U_1(n) = 2*U_1(n-1) - U_1(n-2);
   U_2(n) = 2*U_2(n-1) - U_2(n-2);
   U_3(n) = 2*U_3(n-1) - U_3(n-2);
   
   
   % Computing final primitive variables.
   
   rho       = U_1./A;                                        % Non-dimentional density.
   V         = U_2./U_1;                                      % Non-dimentional velocity.
   T         = (gamma-1)*((U_3./U_1) - ((gamma/2)*(V.^2)));   % Non-dimentional Temperature.
   P         = rho.*T;                                        % Non-dimentional pressure.
   Mach      = V./sqrt(T);                                    % Mach number.
   mass_flow_c = (rho.*V).*A;                                 % Non-dimentional mass flow rate.
   
   
   % Transient throat values.
   Mach_throat(k1) = Mach(n_throat);
   rho_throat(k1)  = rho(n_throat);
   P_throat(k1)    = P(n_throat);
   T_throat(k1)    = T(n_throat);
   
   % Transient mass flow rate.
   
   Transientmassflow_c = figure(1);
   hold on;
   
   if iter == 200;
       
       plot(x,mass_flow_c,'r');
       
   else if iter ==500;
           
           plot(x,mass_flow_c,'b','linewidth',1);
           
       else if iter ==700;
               
               plot(x,mass_flow_c,'g','linewidth',2);
               
           end
           
       end
       
   end
   
end
figure(1)
title('Transient quasi 1-D nozzle mass flow (conservative form)');
xlabel('x/L');
ylabel('Massflowrate');
legend('200Deltat','500deltat','700deltat');
grid on  
   
     
end

 

Non-Conservation Form Script:

function [rho,V,T,A,P,mass_flow,Mach,rho_throat,T_throat,P_throat,Mach_throat,iter,k1] = Main_Script_Non_conservation(n,x,dx,A,gamma,n_throat,c);
  
% Initializing Numerical Discretization Method Using Maccormack Principle.
% Initializing time marching loop.
% Initial Conditions
rho   = 1 - 0.3146*x;                % Density.
T     = 1 - 0.2314*x;                % Temperature.                
V     = ((0.1 + 1.09*x).*(T)).^0.5;  % Velocity.
nt = 1000;
iter = 0;
for k1 = 1:nt
    
    iter = iter + 1;
    
    V_old   = V;
    T_old   = T;
    rho_old = rho;
    
    % CFL Time loop.
    a = sqrt(T);
    timestep = c*dx./(a + V);
    dt = min(timestep);
    
    % Initializing Predictor Step.
   
   for j = 2:n-1
        
       % Continuity equation.
       dV_dx            = (V(j+1) - V(j))/dx;
       
       dT_dx            = (T(j+1) - T(j))/dx;
       
       dlogA_dx         = (log(A(j+1)) - log(A(j)))/dx;
       
       drho_dx          = (rho(j+1) - rho(j))/dx;
       
       
       drho_dt_P(j)     = -rho(j)*dV_dx - rho(j)*V(j)*dlogA_dx - V(j)*drho_dx;
       
       % Momentum equation
        
       dV_dt_P(j)       = -V(j)*dV_dx  -  ((gamma^-1)*(dT_dx + ((T(j)/rho(j))*drho_dx)));
       
       % Energy equation.
       
       
       dT_dt_P(j)       = -V(j)*dT_dx - ((gamma-1)*T(j)*(dV_dx + (V(j)*dlogA_dx)));
       
       rho(j)           = rho(j) + drho_dt_P(j)*dt;
       V(j)             = V(j)   + dV_dt_P(j)*dt;
       T(j)             = T(j)   + dT_dt_P(j)*dt;
       
      
   end
   
   
   
   % Initializing Corrector Step
      
   for j = 2:n-1
          
      % Continuity equation.
       dV_dx            = (V(j) - V(j-1))/dx;
       
       dT_dx            = (T(j) - T(j-1))/dx;
       
       dlogA_dx         = (log(A(j)) - log(A(j-1)))/dx;
       
       drho_dx          = (rho(j) - rho(j-1))/dx;
       
       
       drho_dt_C(j)     = -rho(j)*dV_dx - rho(j)*V(j)*dlogA_dx - V(j)*drho_dx;
       
       % Momentum equation
        
       dV_dt_C(j)       = -V(j)*dV_dx  -  ((gamma^-1)*(dT_dx + ((T(j)/rho(j))*drho_dx)));
       
       % Energy equation.
       
       
       dT_dt_C(j)       = -V(j)*dT_dx - ((gamma-1)*T(j)*(dV_dx + (V(j)*dlogA_dx)));
       
       
     
   end
          
      % Time Average Value.
      
      drho_dt           = 0.5*(drho_dt_P + drho_dt_C);
      dV_dt             = 0.5*(dV_dt_P   + dV_dt_C);
      dT_dt             = 0.5*(dT_dt_P   + dT_dt_C);
      
      % Final output solution.
      
   for k = 2:n-1
          
      V(k)              = V_old(k)   + dV_dt(k)*dt;
      T(k)              = T_old(k)   + dT_dt(k)*dt;
      rho(k)            = rho_old(k) + drho_dt(k)*dt;
        
   end
      
   P         = rho./T;
   Mach      = V./sqrt(T);                                    % Mach number.
   mass_flow = (rho.*V).*A;                                   % Non-dimentional mass flow rate.
   
   % Boundary Condition.
   % Inlet Conditions.
   V(1)                 = 2*V(2) - V(3);
   % Outlet Conditions.
   V(n)                 = 2*V(n-1) - V(n-2);
   T(n)                 = 2*T(n-1) - T(n-2);
   rho(n)               = 2*rho(n-1) - rho(n-2);
   
   % Transient throat values.
   Mach_throat(k1) = Mach(n_throat);
   rho_throat(k1)  = rho(n_throat);
   P_throat(k1)    = P(n_throat);
   T_throat(k1)    = T(n_throat);
   
   % Transient mass flow rate.
   
   Transientmassflow = figure(2);
   hold on;
   
   if iter == 100;
       
       plot(x,mass_flow,'r');
       
   else if iter ==300;
           
           plot(x,mass_flow,'b','linewidth',1);
           
       else if iter ==500;
               
               plot(x,mass_flow,'g','linewidth',2);
               
           else if iter ==700;
                   
                   plot(x,mass_flow,'--r','linewidth',1);
                   
               end
               
           end
           
       end
       
   end
   
end
figure(2)
title('transient quasi 1-D nozzle mass flow (Non-conservative form)');
xlabel('x/L');
ylabel('Massflowrate');
legend('100deltat','300deltat','500deltat','700deltat');
grid on  
      
      
end

Results:

(A) Mass Flow Rate at different time iteration for both conservative and non-conservative form of governing equation. @ 31 grid points

 

(B) Flow field variation at end of time step @ 31 grid points. (steady-state).

 

 

(C) Fluid properties variations at different time at nozzle throat @ 31 grid points. (transient-state).

(D) Mass Flow Rate of Conservative form vs non-conservative form @31 grids.

 

(E) Mass Flow Rate at different time iteration for both conservative and non-conservative form of governing equation. @ 61 grid points

 

(F) Flow field variation at end of time step @ 61 grid points. (steady-state).

(G) Fluid properties variations at different time at nozzle throat @ 61 grid points. (transient-state).

(H) Mass Flow Rate of Conservative form vs non-conservative form @61 grids.

(I)Residual terms of time derivative properties at Throat @ 61 grid points.

(j)Residual terms of time derivative properties at Throat @ 31 grid points.

 

Observation:

  • The first plot gives the insite of non-dimetionalised mass-flow rate vs time . 3/4 different curves are shown, each for a different time during the course of the time-marching procedure. The dashed curve is the variation of p, V,A which pertains to the initial conditions, and hence it is labeled 0△t. As tthe time step increases, the curve becomes more staight essentially converges to a steady state value.
  • The non-dimentional mass flow rate  has 2 plots. one is conservative and another is non-conservative. both the plots shows the curves of the steady-state mass flow rate using different forms of governing equations.
  • The next plots shows the flow-field variations at the end of the time steps.
  • Plot "E" and "H" shows the variation of mass flow rate over nozzle length for both conservative and non-conservative form. The comparison  illustrates a general advantage of the conservation form of the equations. The conservation form does a better job of preserving mass throughout the flow field, mainly because the mass flow itself is one of the dependent variables in the equations~the mass flow is a primary result from these equations. In contrast, the dependent variables in the nonconservation form of the equations are the primitive variables, and the mass flow is obtained only
    as a secondary result. Because the conservation form of the equations does a better job of conserving mass throughout the flow.
  • The nonconservation form leads to smaller residuals. The amount by which the residuals decay is often used as an index of "quality" of the numerical algorithm.In this sense, the nonconservation form does a better job.
  • There is no clear superiority of either form in terms of accuracy of the results.
    The nonconservation form seems to produce slightly more accurate results for
    the primitive variables, and the conservation form seems to produce slightly
    more accurate results for the flux variables. The results in either case are certainly satisfactory.
  • It is interesting to examine the variation of the time derivatives as a function of time itself from (plot "H" and "I"), or equivalently as a function of the number of time steps. Once again focusing on the nozzle throat ( at grid point i = 16) gives the variation of the time derivatives of nondimensional density and velocity as a function of the number of time steps.
  • From the plot H and I we can observe that at early times, the time derivatives are large, and they oscillate in value. These oscillations are associated with various unsteady compression and expansion waves which propagate through the nozzle during the transient process.
  • Higher the number of grid points, the more there will be residuals of thetime derivatives and more time it will take to become steady-state which seems to be the characteristics of Maccormack's technique.
  • The steady is achieved at around 400 time steps for 31 grid points and around 700 time steps for 61 grid points.

 

 

Grid Dependency

                                                                          Conditions At Nozzle Throat

                                                                        ρ/ρ*        T/T*         P/P*        

Case 1: 31 grid points                                      0.6397     0.8365      0.7636       

Case 2:  61 grid points                                     0.6383     0.8357      0.7638

Case 3:  91 grid points                                     0.6290     0.8307      0.7572

From the above grid dependency table we can conclude that higher the number of of grid points, the lesser will be the truncation error and the numerical values will be closer to the actual results but however it will take more simulation time to achieve steady state conditions.

But if we increase the grid points to a number where the solutions are no longer sensitive to the grid points, then at that point it becomes "grid independent".

 


Projects by Ankit Chakrabarty

challange 4
Ankit Chakrabarty · 2019-11-15 12:00:52

Challanges 4 open edges: Open edges or free edges error are the one that is joined to only one faces. 2 joined faces shares 2 vertices for a common edges. It is shown by colour of green.   Read more

 Turbocharger: Literature Review: A turbocharger is mechanical rotating device which forcebly induces air at high pressure into the internal combustion engine, in turn increases the power and effeciency (thermal and volumetric effeciency) of the engine. A turb Read more

  (A) Problem defination:- Solve the following Linear equation using Gauss-Jacobi, Gauss-Seidal and SOR iterative method and understand the following parameters:- Spectral radius and the impact of spectral radius on the convergence of the solution. Effect of d Read more

CASE STUDY:-  A Thermodynamic Analysis Of A Tractor Engine For Use In India With Given Engine Performance Targets. A Brief Literature Review: A tractor is an engineered vehicle specifically designed to deliver Power and traction/tractive effort to haul heavy equi Read more

Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.   Explicit Scheme : We say that a scheme is explit when the information at time level "n+1" depends on previous time Read more

  CASE STUDY: Effect of Burn Rate on engine performance Of A 6-Cyl-Turbocharged-DI Engine Using A Non-predictive Combustion model (DI WIEBE) Abstract: Analytical functions approximating the burn rate in internal combustion engines are useful and cost-effective to Read more

Solving 1-D Linear Convection Using First-Order Backward Difference And Forward Difference method Using MATLAB LITERATURE REVIEW: WHAT IS CONVECTION? Convection is the sum of bulk transport of the  fluid  and brownian/ osmotic dispersion of fluid constituen Read more

Deriving the 4th order approximation for 2nd order derivative of a function using Taylor table for centered, skewed right-sided stencil and skewed left_sided stencil and compare the error results for the same using MATLAB.   Given function, f(x)    =& Read more

Objective:- (1) Determine the engine parameters running at 1800 rpm and record the same. (2) Further increase the power output obtained at 3600 rpm by 10%.   Model:-  The GT model Of the 1-cylinder engine is shown below:- Observations (@1800rpm):-  Read more


Loading...

The End