Simulation of a 1D Super-sonic nozzle flow using Macormack Method

Aim:

1. To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle and derive both the conservation and non-conservation forms of the governing equations and solve them using MacCormack's technique. Also, to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
2. To obtain the following plots and draw suitable conclusions

• Steady-state distribution of primitive variables inside the nozzle
• Time-wise variation of the primitive variables
• Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
• Comparison of Normalized mass flow rate distributions of both forms

    3. Perform a grid independence test on the solution and explain the results  

 

Theory:

A Convergent-divergent nozzle consists of three main parts, namely, convergent section, throat, and divergent section. In the case of supersonic C-D nozzle in our problem, the flow entering through the convergent section will be subsonic in nature and this will change when it reaches the throat section, and subsequently the divergent section, ie, the flow changes to sonic and supersonic respectively. This problem can be solved by using the McCormick method using conservative or non-conservative equations. The continuity, momentum, and energy equations in the conservative and non-conservative form are given below. These equations are non-dimensionalized for easier coding of the equations.

 

Solving using the non-conservative method:

Terms used in the non-dimensionalizing process:

a) `x'=x/L`

b) `P'=P/P_0`

c) `V'=V/a_0`

d) `t'=t/(L/a_0)`

e) `A'=A/A^**`

Where, L = length of the nozzle

P0 = Inlet pressure

a0 = velocity of sound

`A^**` = Throat cross-sectional area

Non-conservative form of equations:

1) Continuity equation

`(delrho')/(delt')=-rho'(delV')/(delx')-rho'V'(del(lnA'))/(delx')-V'(delrho')/(delx')`

2) Momentum equation

`(delV')/(delt')=-V'(delV')/(delx')-1/(gamma)((delT')/(delx')+(T')/(rho')(delrho')/(delx'))`

3) Energy equation

`(delT')/(delt')=-V'(delT')/(delx')-(gamma-1)T'[(delV')/(delx')+V'(dellnA')/(delx')]`

To solve the above equations, the control volume is divided into n number of discrete grid points and the correct values of flow parameter is obtained using the steps given below

Predictor step:

Continuity equation:

Momentum equation:

Energy equation:

The flow parameters are predicted using the following formulas:

corrector step:

Continuity equation:

Momentum equation:

Energy equation:

Average time derivative is given by:

Hence corrected flow field variables is given by:

Boundary conditions:

Inlet boundary conditions:

`V_1=2V_2-V_3`

Density and temperature are independent of time

Outlet boundary conditions:

`V_N=2V_(N-1)-V_(N-2)`

`rho_N=2rho_(N-1)-rho_(N-2)`

`T_N=2T_(N-1)-T_(N-2)`

Initial conditions:

`rho=1-0.3146x`

`T=1-0.2314x`

`V=(0.1+1.09x)T^0.5`

 

Solving using the conservative method:

conservative form of equations:

1) Continuity equation

`(del(rho'A'))/(delt')+(del(rho'A'V'))/(delx')`

2) Momentum equation

`(del(rho'A'V'))/(delt')+(del[rho'A'V'^2+(1/gamma)p'A'])/(delx')=(1/gamma)p'(delA')/(delx')`

3) Energy equation

`(del[rho'((e')/(gamma-1)+(gamma/2)V'^2)A'])/(delt')+(del[rho'((e')/(gamma-1)+(gamma/2)V'^2)V'A'+p'A'V'])/(delx')=0`

Defining elements of the solution vector U, flux vector F, and source term J

Hence,

Non-dimensionalizing the flow parameters:

Writing F1, F2, F3, J2 in terms of U1, U2, U3:

`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_2U_3)/U_1-(gamma(gamma-1))/2U_2^3/U_1^2`

`J_2=(gamma-1)/gamma(U_3-gamma/2U_2^2/U_1)(dellnA')/(delx')`

Boundary conditions:

Initial conditions:

 

Calculating the time step:

Where C refers to the Courant number, for the solution to be stable, C<=1

The value of dt at ith node is given by:

The value dt at (i+1)th node is given by:

The least value of dt is chosen to be the final value of dt to satisfy C<=1

Nozzle shape:

`A=1+2.2(x-1.5)^2, 0<=x<=3`

 

Octave code:

a) Main code

clear all
close all
clc

% inputs variables
n = 31;
L = 3;
x = linspace(0,L,n);

% Courant number
c = 0.5;

% node number present at the throat
throat = ((n-1)/2)+1;

% Time steps
nt = 1400;

tic
% Function for non conservative form
[v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = non_conservative(n,L,c,nt,throat);

% Calculating the simulation time
Simulation_time_using_non_conservative_method = toc

tic
% Function for conservative form
[v2,rho2,t2,m2,mach2,pr2,v_t2,rho_t2,t_t2,m_t2,mach_t2,pr_t2] = conservative(n,L,c,nt,throat);
 
% Calculating the simulation time 
Simulation_time_using_conservative_method = toc
 
% Varaiation of flow parameters at the throat
figure(3)
a1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_t); 
hold on
plot(linspace(1,nt,nt),rho_t2,'color','r');
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.634,0.634,600),'color','k');
grid on;
ylabel("rho'",'Fontsize',12,'fontweight','bold');
h = legend(a1,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);
title('Timewise variation of flow parameters at throat','Fontsize',12,'color','r','fontweight','bold');

a2 = subplot(4,1,2);
plot(linspace(1,nt,nt),t_t); 
hold on
plot(linspace(1,nt,nt),t_t2,'color','r'); 
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.833,0.833,600),'color','k');
grid on;
ylabel("T'",'Fontsize',12,'fontweight','bold');
h = legend(a2,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);

a3 = subplot(4,1,3);
plot(linspace(1,nt,nt),pr_t); 
hold on
plot(linspace(1,nt,nt),pr_t2,'color','r'); 
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.528,0.528,600),'color','k');
grid on;
ylabel("Pr'",'Fontsize',12,'fontweight','bold');
h = legend(a3,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);

a4 = subplot(4,1,4);
plot(linspace(1,nt,nt),mach_t); 
hold on
plot(linspace(1,nt,nt),mach_t2,'color','r'); 
hold on
plot(linspace(nt-400,nt+200,600),linspace(1.0,1.0,600),'color','k');
grid on;
ylabel("Ma'",'Fontsize',12,'fontweight','bold');
xlabel('No of time steps','fontweight','bold');
h = legend(a4,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);

% Variation in flow parameters along x axis
figure(4)
p1 = subplot(4,1,1);
plot(x,mach); 
hold on
plot(x,mach2,'color','r','*', 'markersize', 3);
grid on;
ylabel("Ma'",'Fontsize',12,'fontweight','bold');
h = legend(p1,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')
title('Variation of Flow parameters along x axis','Fontsize',12,'color','r','fontweight','bold');

p2 = subplot(4,1,2);
plot(x,pr);
hold on
plot(x,pr2,'color','r','*', 'markersize', 3); 
grid on;
ylabel("Pr'",'Fontsize',12,'fontweight','bold');
h = legend(p2,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')

p3 = subplot(4,1,3);
plot(x,rho); 
hold on
plot(x,rho2,'color','r','*', 'markersize', 3);
grid on;
ylabel("rho'",'Fontsize',12,'fontweight','bold');
h = legend(p3,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')

p4 = subplot(4,1,4);
plot(x,t); 
hold on
plot(x,t2,'color','r','*', 'markersize', 3);
grid on;
ylabel("T'",'Fontsize',12,'fontweight','bold');
xlabel('Non-dimensional length','fontweight','bold');
h = legend(p4,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')

% mass flow rate comparison between two forms of equation
figure(5)
plot(x,m,'r', 'linewidth', 1.5)
hold on
plot(x,m2,'b', 'linewidth', 1.5)
hold on
plot(x,linspace(0.579,0.579,n),'color','k','--', 'markersize', 3)
ylabel('Non-dimensional mass flow rate','fontweight','bold','fontsize',12);
xlabel('Non-dimensional distance through nozzle(x)','Fontsize',12,'fontweight','bold');
legend({'Non-Conservative Form', 'Conservative Form','Exact analytical solution'});
title('Non-dimensional mass flow rate comparison','fontweight','bold','fontsize',12,'color','r'); 
axis([0 3 0.575 0.600]);
grid on

 

b) Code to solve non-conservative equations

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = non_conservative(n,L,c,nt,throat)

% inputs variables
x = linspace(0,L,n);
dx = x(2)-x(1);

% Initial conditions
rho = 1-0.3146*x;
t = 1-0.2314*x;
v = (0.1+1.09*x).*t.^0.5;
a = 1+2.2*(x-1.5).^2;
gamma = 1.4;

% calculating time step
for i = 1:n
  del_t(i) = c*(dx/(t(i)^0.5+v(i)));  
end
  dt = min(del_t);

% Time loop
for k = 1:nt

  rho_old = rho;
  v_old = v;
  t_old = t;

  % Predictor method
  for j = 2:n-1
     
    dvdx = (v(j+1)-v(j))/dx; 
    drhodx = (rho(j+1)-rho(j))/dx;
    dtdx = (t(j+1)-t(j))/dx;
    dlogadx = (log(a(j+1))-log(a(j)))/dx;
    
    % continuity equation
    drhodt_p(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;
  
    % momentum equation
    dvdt_p(j) = -v(j)*dvdx -(1/gamma)*(dtdx +(t(j)/rho(j))*drhodx);
    
    % energy equation
    dtdt_p(j) = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
        
    % solution update
    rho(j)=rho(j)+drhodt_p(j)*dt;
    v(j)=v(j)+dvdt_p(j)*dt;
    t(j)=t(j)+dtdt_p(j)*dt;
  end 

% corrector method
  for j = 2:n-1
     
    dvdx = (v(j)-v(j-1))/dx; 
    drhodx = (rho(j)-rho(j-1))/dx;
    dtdx = (t(j)-t(j-1))/dx;
    dlogadx = (log(a(j))-log(a(j-1)))/dx;
    
    % continuity equation
    drhodt_c(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;
  
    % momentum equation
    dvdt_c(j) = -v(j)*dvdx -(1/gamma)*(dtdx +(t(j)/rho(j))*drhodx);
    
    % energy equation
    dtdt_c(j) = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
  end 

% average derivative
drhodt = 0.5*(drhodt_p+drhodt_c);
dvdt = 0.5*(dvdt_p+dvdt_c);
dtdt = 0.5*(dtdt_p+dtdt_c);

% final solution update
for p = 2:n-1
    rho(p) = rho_old(p) + drhodt(p)*dt;
    v(p) = v_old(p) + dvdt(p)*dt;
    t(p) = t_old(p) + dtdt(p)*dt;
end
     
% boundary conditions

  % inlet BC
  v(1)=2*v(2)-v(3);
  
  % outlet BC
  v(n)=2*v(n-1)-v(n-2);
  rho(n)=2*rho(n-1)-rho(n-2);
  t(n)=2*t(n-1)-t(n-2);
  
  % calculating solutions 
   m = rho.*v.*a;
   mach = v./((t).^0.5);
   pr = rho.*t;
  
  % Throat values 
   rho_t(k) = rho(throat);
   t_t(k) = t(throat);
   pr_t(k) = pr(throat);
   mach_t(k) = mach(throat);
   v_t(k) = v(throat);
   m_t(k) = m(throat);
   
   % plotting mass flow rate at different time steps
     figure(1);
   if k  == 50
      plot(x, m, 'b', 'linewidth', 1.5)
    hold on
     elseif k == 100
      plot(x, m, 'r', 'linewidth', 1.5)
    hold on
     elseif k == 150
      plot(x, m, 'c', 'linewidth', 1.5)
    hold on
     elseif k == 200
      plot(x, m, 'g', 'linewidth', 1.5)
    hold on
     elseif k == 550
      plot(x, m, 'k', 'linewidth', 1.5) 
    hold on 
     elseif k == 700
      plot(x, m, 'y', 'linewidth', 1.5)
    hold on
      elseif k == 1400
      plot(x, m,'m', 'linewidth', 1.5)
      grid on
      ylabel('Non-dimensional mass flow rate','fontweight','bold','fontsize',12);
      xlabel('Non-dimensional distance through nozzle(x)','Fontsize',12,'fontweight','bold');
      title('Non-dimensional mass flow rate(Non-conservative form) along x, at different time steps','Fontsize',12,'fontweight','bold','color','r');
      legend({'50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps'});
   end
  end
end

 

c) Code to solve conservative equations

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = conservative(n,L,c,nt,throat)

% input variables
x = linspace(0,L,n);
dx = x(2)-x(1);

% Initial conditions
a = 1+2.2*(x - 1.5).^2;

for i = 1:n
    if (x(i)>=0 && x(i)<=0.5)
        rho(i) = 1;
        t(i) = 1;
    elseif (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);
    elseif (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

% Assuming steady state flow for the initial condition
  v = 0.59./(rho.*a);
  gamma = 1.4;

% Solution vector parameters
  pr = rho.*t;
  u1 = rho.*a;
  u2 = rho.*a.*v;
  u3 = u1.*((t./(gamma-1))+((gamma/2).*(v.^2)));

% calculating time step
     dt = min (c.*(dx./((t.^0.5)+v))); 
 
 % Time loop
for k = 1:nt
   
    % Copy of solution vectors   
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
        
    % flux terms
    f1 = u2;
    term = (((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
    f2 = ((u2.^2)./u1)+ term;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);

    % Predictor Method
    for i = 2:n-1
       
        df1dx(i) = (f1(i+1)-f1(i))/dx;
        df2dx(i) = (f2(i+1)-f2(i))/dx;
        df3dx(i) = (f3(i+1)-f3(i))/dx;
        dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
        
        % j_term(i) = term(i)*dlogadx(i); 
        dadx(i) = (a(i+1)-a(i))/dx;
        j2(i) = (1/gamma)*(rho(i)*t(i))*dadx(i); 
        
        % continuity equation
        du1dt_p(i) = -df1dx(i);
        
        % momentum equation
        du2dt_p(i) = -df2dx(i) + j2(i);
        
        % energy equation
        du3dt_p(i) = -df3dx(i);
        
        %Solution update
        u1(i) = u1(i) + du1dt_p(i)*dt;
        u2(i) = u2(i) + du2dt_p(i)*dt;
        u3(i) = u3(i) + du3dt_p(i)*dt;
    end
    
    % Decoding primitive terms to initialize the corrector step
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
    
    % flux terms for the corrector step
    f1 = u2;
    term1 = (((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
    f2 = ((u2.^2)./u1)+ term1;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
    % Corrector Method
    for i = 2:n-1

        df1dx(i) = (f1(i)-f1(i-1))/dx;
        df2dx(i) = (f2(i)-f2(i-1))/dx;
        df3dx(i) = (f3(i)-f3(i-1))/dx;
        dlogadx(i) = (log(a(i))-log(a(i-1)))/dx;
        
        % j_term1(i) = term1(i)*dlogadx(i); 
        dadx1(i) = (a(i)-a(i-1))/dx;
        j3(i) = (1/gamma)*(rho(i)*t(i))*dadx1(i); 
 
        % Continuity equation
        du1dt_c(i) = -df1dx(i);
        
        % Momentum equation
        du2dt_c(i) = -df2dx(i) + j3(i);
        
        %Energy equation
        du3dt_c(i) = -df3dx(i);
    end
    
    % Average derivative
    du1dt = 0.5*(du1dt_p + du1dt_c);
    du2dt = 0.5*(du2dt_p + du2dt_c);
    du3dt = 0.5*(du3dt_p + du3dt_c);
    
    % Final solution update
 for j = 2:n-1
    u1(j) = u1_old(j) + du1dt(j)*dt;
    u2(j) = u2_old(j) + du2dt(j)*dt;
    u3(j) = u3_old(j) + du3dt(j)*dt;
 end
    
% boundary conditions
    
  % inlet BC
    % u1 value remains same, because rho is not a floating variable and also 'a' doesnt change with time
    u1(1) = rho(1)*a(1);
    % velocity is the floating primitive term, which can be calculated by extrapolating u2
    u2(1) = 2*u2(2) - u2(3);
    % Floating variable v, is calculated for every time step
    v(1) = u2(1)/u1(1);
    % updating u3
    u3(1) = u1(1)*((t(1)/(gamma-1))+((gamma/2)*v(1)^2));
    
  % outlet BC
   u1(n) = 2*u1(n-1) - u1(n-2);
   u2(n) = 2*u2(n-1) - u2(n-2);
   u3(n) = 2*u3(n-1) - u3(n-2);
    
   % Decoding primitives terms 
   rho = u1./a; 
   v = u2./u1;
   t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
   
  % calculating solutions 
   m = rho.*v.*a;
   mach = v./((t).^0.5);
   pr = rho.*t;
    
   % Throat values 
   rho_t(k) = rho(throat);
   t_t(k) = t(throat); 
   pr_t(k) = pr(throat);
   mach_t(k) = mach(throat);
   v_t(k) = v(throat);
   m_t(k) = m(throat);
   
   % plotting mass flow rate at different time steps
    figure(2);
    if k  == 50
      plot(x, m, 'b', 'linewidth', 1.5)
      hold on
    elseif k  == 100
      plot(x, m, 'r', 'linewidth', 1.5)
      hold on
     elseif k  == 150
      plot(x, m, 'c', 'linewidth', 1.5)
      hold on
     elseif k  == 200
      plot(x, m, 'g', 'linewidth', 1.5)
      hold on
     elseif k == 550
      plot(x, m, 'y', 'linewidth', 1.5) 
      hold on  
     elseif k  == 700
      plot(x, m, 'k', 'linewidth', 1.5)
      hold on
     elseif k  == 1400
      plot(x, m,'m', 'linewidth', 1.5)
      grid on
      ylabel('Non-dimensional mass flow rate','fontweight','bold','fontsize',12);
      xlabel('Non-dimensional distance through nozzle(x)','Fontsize',12,'fontweight','bold');
      title('Non-dimensional mass flow rate(Conservative form) along x, at different time steps','Fontsize',12,'fontweight','bold','color','r');
      legend({'50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps'});
    end
  end
end

 

Outputs:

The code was run for courant number, c=0.5, and the number of grid points taken initially was 31(n). This Courant number was chosen to get a stable solution since we know that if c>1, the solution becomes unstable, and hence c is usually taken as less than or equal to one. The number of grid points is taken as an odd number, to get a node exactly at the throat of the nozzle

a) Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process (non-conservative)

Here, by observing the above plot, we can conclude that the solution is moving towards stability, ie, the solution is getting converged with the time steps performed, after the initial fluctuation. The plots of time step 550 and above are overlapping,  from which we can conclude that the solution gets converged at around 550-time steps.

b) Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process (conservative)

Here, by observing the above plot, we can conclude that the solution is moving towards stability, ie, the solution is getting converged with the time steps performed, after the initial fluctuation. The plots of time step 550 and above are overlapping,  from which we can conclude that the solution gets converged at around 550 time-steps. Also, we can see that the fluctuations in the solution solved using the conservative solution is much less when it is solved using non-conservative equations. 

c) Steady-state distribution of primitive variables inside the nozzle

The value of flow parameters solved using conservative and non-conservative equations gave similar trends.

Except for Ma', whose magnitude increases along the length of the nozzle, all other flow parameters simulated (Pr', Rho', T'), the magnitude decreases.

d) Time-wise variation of the primitive variables

The above plot gives the variation of flow parameters with the number of time steps. We can see that the solution tends to stabilize with the time steps performed, the solution is converged in case of both the equations. The conservative form of equations, even though they are more difficult to code, they tend to be more stable compared to the non-conservative form of equations.

e) Comparison of Normalized mass flow rate distributions of both forms

The conservative form of the equation is more stable compared to the non-conservative form of the equation. The conservative form of the equation gives almost constant mass flow rate throughout its length, whereas, in the case of the non-conservative form of the equation, the mass flow rate is fluctuating throughout the length.

f) Time taken to perform simulations

 

Grid independence test:

Whenever the cell sizes are reduced by increasing the number of grid points, the solution tends to be more accurate when they converge. If the difference between the original parameter and the updated parameter is less than a certain tolerance, the initial solution is said to be grid-independent. 

Even though smaller cell sizes tend to give more accurate results, the time steps need to be increased considerably for the solution to converge. This leads to an increase in simulation time, which is not desirable. Hence by performing the grid independence test, we are checking whether the difference between the initial and new parameters is within a certain tolerance limit or not. If they are within a certain tolerance limit, the solution is said to grid-independent, and the simulations can be run using the same number of grid points. If the difference is greater than a certain tolerance limit, the number of grid points can be increased and further simulations carried out.

In this problem, the number of grid points taken initially is 31, to perform the grid-independence test, we need to increase the number of grid points, say 61. Now that we have doubled the number of grid points, the number of time steps also must be increased to get a converged solution, 2800 time steps are taken. 

 

Outputs:

a) Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process (non-conservative)

b) Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process (conservative)

 

c) Steady-state distribution of primitive variables inside the nozzle

 

d) Time-wise variation of the primitive variables

 

e) Comparison of Normalized mass flow rate distributions of both forms

 

f) Time taken to perform simulations

 

Comparison:

Flow Parameters (Non-dimensionalized)	Non-conservative form		Conservative form		Exact solution for the flow parameter
	n = 31	n = 61	n = 31	n = 61
Density	0.639	0.638	0.652	0.639	0.634
Temperature	0.836	0.835	0.841	0.836	0.833
Pressure	0.534	0.532	0.548	0.533	0.528
Mach number	0.999	1.0	0.978	0.998	1.0
Simulation time(s)	24.721	75.650	24.096	82.195	-

 

 

 

 

 

 

After the code was re-run for 61 grid points, to perform a grid independence test, the obtained results are compared with the results obtained with 31 grid points. Since there is only a marginal variation in the results obtained, the solution obtained by solving using 31 grid points is said to be grid-independent.

 

Results and conclusions:

1. In a subsonic-supersonic C-D nozzle, all non-dimensionalized flow parameters except Mach number decreases along the length of the nozzle. The Mach number increases from subsonic to supersonic speed when it passes from inlet to outlet.
2. The results obtained using the conservative equations tend to be more stable compared to that of the non-conservative equations.
3. Since there is only a marginal difference in the solutions obtained using 31 and 61 grid points respectively, the solution obtained using 31 grid points is said to be grid-independent
4. The conservative method takes marginally more time to solve compared to the non-conservative method because of the presence of complex equations