Simulation of a 1D Super-sonic nozzle flow simulation using MacCormack s Method

Objective: To write code to solve the 1D supersonic nozzle flow equations using the Macormack Method.

1. Implement both conservative and non-conservative forms of the governing equations.
2. Perform a grid dependence test.
3. Write separate functions for conservative and non-conservative forms.
4. Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions.
5. Which method is faster?
6. Minimum number of cycles for which the simulation should be run in order for convergence

Theory:

Maccormack Technique(Predictor & Corrector method):

MacCormack\'s technique is a variant of the Lax-Wendroff approach but is much simpler in its application. Maccormack method is also an explicit finite-difference technique which is second-order-accurate in both space and time. The results obtained by using MacCormack\'s method are perfectly satisfactory for many fluid flow applications.

we assume that the flow field at each grid point is known at time t, and we proceed to calculate the flow-field variables at the same grid points at time `t+Deltat`.
First, consider the density at grid point (i) at time `t+Deltat`.

where, `((delrho)/(delt))_(avg)`  is mean value of `((delrho)/(delt))` between time t and `Deltat`

similarly, other flow field variable can be written as:

Here we are going to consider 1D supersonic fluid flow through the nozzle as illustrated in fig.

We are assuming here is flow along nozzle is only in one direction i.e. X-axis that\'s why all flow-field variables calculated above only for one direction i.e X-axis.

SETUP:

we will set up three echelons of equations as follows:

1. The governing equation in terms of partial differential equations suitable for the time-marching solution of quasi-1D  flow.
2. The finite-difference expressions pertaining to MacCormack\'s technique as applied to this problem will be set up. 
3. Other details for the numerical solution (such as the calculation of the time step and the treatment of boundary conditions) will be formulated. 

1. Governing Equations:

Non-conservative form of the equation for flow field variable:

Continuity equation:

Momentum Equation:

Energy Equation:

`e= C_v* T` where Cv is Specific heat at constant volume and T is temperature

All the above flow field variables are non-dimensionalized i.e.

where,

variable with subscripts 0 is flow filed variables at a reservoir

Conservative form of the equation for flow field variable:

Continuity equation:

 `where, U1=rhoA, F1= rhoAv`

Momentum Equation:

`where, U2=rhoAv, F2= rhoAv^2+(1/(gamma))PA, J2=1/(gamma)*P((delA)/(delx))`

P(Pressure)= rho * T.

Energy Equation:

where,

`U3=[rhoA(e/(gamma-1)+gamma/2v^2)], F3=[rhovA(e/(gamma-1)+gamma/2v^2)+Pav]`

2. Finite Difference Expressions:

In Predictor method- All Spatial terms we setup as FORWARD DIFFERENCES
In Corrector method- All Spatial terms we setup as BACKWARD DIFFERENCES

In the Predictor step: 
we will get predicted values as

From the above equations flow-field variables known at time t.

In the Corrector step: 

all these values are calculated using BACKWARD differences for spatial nodes

Average time-derivative will be given by:

Finally, we have for the corrected values of the flow-field variables at time `t+Deltat`

3. Details for the numerical solution

Calculation of TIME-STEP:

The governing equations we are using is Hyperbolic with respect to time.

applicable to both Conservative and Non-conservative forms of equations.

Boundary Conditions:

For Non-conservative form

Subsonic Inflow BC:

The slope of the linear extrapolation line is determined from points 2 and 3 as

SLOPE=(V3-V2)/dx

V1 by linear extrapolation

V1=2*V2 - V3

Supersonic Outflow BC:

We again choose to use linear extrapolation based on the flow field values at all the internal points

V(N) = 2*V(N-1) - V(N-2)

P(N) = 2*P(N-1) - P(N-2)

T(N)= 2*T(N-1) - T(N-2)

For Conservation Form

Subsonic In-Flow

U1= rho* A

U2(N)=  2u2(N-1) - U2(N-2)

Supersonic Outflow 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)

where N denotes grid points at the outflow boundary

NOZZLE SHAPE AND INITIAL CONDITIONS:

For Non-conservation form

A= 1 + (2.2 * (x - 1.5)^2)

rho= 1-0.3146 * x;
t = 1-0.2314 * x; 
v= (0.1+1.09*x) * T^0.5;

For Conservation form

A= 1 + (2.2 * (x - 1.5)^2)

rho= 1 }
T = 1   }   for 0 <= x <= 0.5

rho = 1 - 0.366(x - 0.5)  }
T = 1 - 0.l67 (x\' - 0.5)     }  for 0.5<= x<= l.5

rho = 0.634 - 0.3879(x - 1.5) }
T = 0.833 - 0.3507(x\' _ 1.5)  }  for 1.5 <= x <= 3.5 

Main Program:

clear all
close all
clc

% initial conditions
n=31; % grid points
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma = 1.4;


C=0.5;
nt=1000;


% For Non-conservation form of governing equation
tic;
[rho,v,t,mass_flow_rate,p_ncn,M_ncn,mass_flow_rate_th,p_ncn_th,M_ncn_th,rho_th,t_th]=non_conservative(n,x,dx,gamma,C,nt)
t_ncn=toc;

sprintf(\'Computation time for Non-conservative form is %d\',t_ncn)

% For Conservation form of governing equation
tic;
[rho_c,v_c,t_c,mass_flow_rate_c,p_cn,M_cn,mass_flow_rate_c_th,p_cn_th,M_cn_th,rho_c_th,t_c_th]=conservative(n,x,dx,gamma,C,nt) 
t_cn=toc;

sprintf(\'Computation time for Conservative form is %d\',t_cn)


% comparing mass flow rates
figure(1)
plot(x,mass_flow_rate)
hold on
plot(x,mass_flow_rate_c)
legend (\'Non-conservative\',\'Conservative\')
xlabel(\'Length of Nozzle x(m)\')
ylabel(\'Mass flow rate across lengh od nozzle\')
title(\'Normalised mass flow rate comparision along length of nozzle\')

figure(2)
subplot(3,1,1)
plot(x,v,\'b\')
hold on
plot(x,v_c,\'k*\')
legend (\'Non-conservative form\',\'Conservative form\')
xlabel(\'Length of Nozzle x(m)\')
ylabel(\'Velocity variation\')
title(\'Normalised Velocity comparision along length of Nozzle\')

subplot(3,1,2)
plot(x,rho,\'b\')
hold on
plot(x,rho_c,\'k*\')
legend (\'Non-conservative form\',\'Conservative form\')
xlabel(\'Length of Nozzle x(m)\')
ylabel(\'Density variation\')
title(\'Normalised Density comparision along Length of Nozzle\')

subplot(3,1,3)
plot(x,t,\'b\')
hold on
plot(x,t_c,\'k*\')
legend (\'Non-conservative form\',\'Conservative form\')
xlabel(\'Length of Nozzle x(m)\')
ylabel(\'Temperature variation\')
title(\'Normalised Temperature comparision along Length of Nozzle\')

Non_conservative Part:

function [rho,v,t,mass_flow_rate,p_ncn,M_ncn,mass_flow_rate_th,p_ncn_th,M_ncn_th,rho_th,t_th]=non_conservative(n,x,dx,gamma,C,nt)


% calculate initial profiles

rho= 1-0.3146*x;
t = 1-0.2314*x;  % t = temperature
v=  (0.1+1.09*x).*t.^0.5;

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

throat_ncn= find(a==1);


[t_step_ncn]=time_step_ncn(x,v,dx,n,C,t)
dt= t_step_ncn;
 
 for k=1:nt
  
rho_old = rho;
  v_old = v;
  t_old = t;

    % Predictor method (FDS)

     for j=2:n-1
      dvdx=(v(j+1)-v(j))/dx;
      dlogadx=(log(a(j+1))-log(a(j)))/dx;
      drhodx= (rho(j+1)-rho(j))/dx;
      dtdx=(t(j+1)-t(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;

     endfor

  % Corrector method

     for j=2:n-1
      dvdx=(v(j)-v(j-1))/dx;
      dlogadx=(log(a(j))-log(a(j-1)))/dx;
      drhodx= (rho(j)-rho(j-1))/dx;
      dtdx=(t(j)-t(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);  


     endfor

     % To compute Average time derivative
   
   drhodt = 0.5.* (drhodt_p + drhodt_c);
   dvdt = 0.5.* (dvdt_p + dvdt_c);
   dtdt = 0.5.* (dtdt_p + dtdt_c);

  for i=2:n-1
     
  % Final solution update
    rho(i) = rho_old(i) + drhodt(i).*dt;
    v(i) = v_old(i) + dvdt(i).*dt;
    t(i) = t_old(i) + dtdt(i).*dt;

  endfor
  
  % Apply Boundary condition
   % Inlet
     v(1) = 2*v(2) - v(3);
   
   % Outlet
     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); 
     
   
% Nondimentional flow properties
% masss flow rate
mass_flow_rate= rho.*a.*v;

% Pressure
p_ncn= rho.*t;  

% Mach no.
M_ncn= v./(t.^0.5);


% flow properties @ throat
% masss flow rate
mass_flow_rate_th= mass_flow_rate(throat_ncn);

% Pressure
p_ncn_th= p_ncn(throat_ncn);  

% Mach no.
M_ncn_th= M_ncn(throat_ncn);

% Density
rho_th= rho(throat_ncn);

% Temperarue
t_th= t(throat_ncn);
    

endfor

endfunction

Conservative Form:

function [rho_c,v_c,t_c,mass_flow_rate_c,p_cn,M_cn,mass_flow_rate_c_th,p_cn_th,M_cn_th,rho_c_th,t_c_th]=conservative(n,x,dx,gamma,C,nt)
  

% calculate initial profiles
for i=1:n
if x(i)>=0 && x(i)<=0.5
    rho_c(i)=1;
    t_c(i)=1;
elseif x(i)>0.5 && x(i)<=1.5
rho_c(i)=1-(0.366*(x(i)-0.5));
t_c(i)=1-(0.167*(x(i)-0.5));

elseif x(i)>1.5 && x(i)<=3.5
rho_c(i)=0.634-(0.3879*(x(i)-1.5));
t_c(i)=0.833-(0.3507*(x(i)-1.5));

end
end

    
%area
a_c = 1+(2.2*(x-1.5).^2);

throat_cn = find(a_c==1);
%Velocity
v_c = 0.59./(rho_c.*a_c);

% Solution vector
u_1= rho_c.*a_c; 
u_2= rho_c.*a_c.*v_c;
u_3= u_1.*((t_c./(gamma-1))+(gamma/2).*v_c.^2);
p = rho_c.*t_c;


for k=1:nt

u1_old= u_1;
u2_old= u_2;
u3_old= u_3;

[t_step_cn] = time_step_cn(x,v_c,dx,n,C,t_c);
dt=t_step_cn;

% flux terms 
f_1 = u_2;
f_2 = ((u_2.^2)./u_1)+((gamma-1)/gamma).*(u_3-(0.5*gamma*((u_2.^2)./u_1)));
f_3 = (gamma.*u_2.*u_3./u_1)-0.5*gamma*(gamma-1).*((u_2.^3)./(u_1.^2));

% Predictor method (FDS) 

 for j=2:n-1

 %source term
 
     dF1_dx_p=(f_1(j+1)-f_1(j))/dx;
     dF2_dx_p=(f_2(j+1)-f_2(j))/dx;
     dF3_dx_p=(f_3(j+1)-f_3(j))/dx;
     dA_dx=((a_c(j+1))-(a_c(j)))/dx;
     dlogA_dx=(log(a_c(j+1))-log(a_c(j)))/dx;

%continuity equation

     du_1dt_p(j)= -dF1_dx_p;
%Momentum equation
     
     j_2(j) = dlogA_dx*((gamma-1)/gamma)*(u_3(j)-(0.5*gamma*u_2(j)^2)/(u_1(j)));
     du_2dt_p(j) = -dF2_dx_p + j_2(j);

% energy equation

     du_3dt_p(j)= -dF3_dx_p;
 
% solution update
    u_1(j)= u_1(j) + du_1dt_p(j)*dt;
    u_2(j)= u_2(j) + du_2dt_p(j)*dt;
    u_3(j)= u_3(j) + du_3dt_p(j)*dt;

 end
 
f_1 = u_2;
f_2 = ((u_2.^2)./u_1)+((gamma-1)/gamma).*(u_3-(0.5*gamma*((u_2.^2)./u_1)));
f_3 = (gamma.*u_2.*u_3./u_1)-0.5*gamma*(gamma-1).*((u_2.^3)./(u_1.^2));


rho_c = u_1./a_c;
v_c = u_2./u_1;
t_c = (gamma-1).*((u_3./u_1)-((0.5*gamma).*(u_2./u_1).^2));
 
 
% Corrector method (BDS)

for j=2:n-1
 
     dF1_dx_c=(f_1(j)-f_1(j-1))/dx;
     dF2_dx_c=(f_2(j)-f_2(j-1))/dx;
     dF3_dx_c=(f_3(j)-f_3(j-1))/dx;
     dA_dx=((a_c(j))-(a_c(j-1)))/dx;
     dlogA_dx=(log(a_c(j))-log(a_c(j-1)))/dx;
   
% source term
     j_2(j) = dlogA_dx*((gamma-1)/gamma)*(u_3(j)-(0.5*gamma*u_2(j)^2)/(u_1(j)));
 
% continuity equation

     du_1dt_c(j)= -dF1_dx_c;

% Momentum equation

     du_2dt_c(j)= -dF2_dx_c + j_2(j);

% Energy equation

     du_3dt_c(j)= -dF3_dx_c;

end

% Average values

du_1dt_avg= 0.5.*(du_1dt_p + du_1dt_c);
du_2dt_avg= 0.5.*(du_2dt_p + du_2dt_c);
du_3dt_avg= 0.5.*(du_3dt_p + du_3dt_c);

% Final solution update

 for i= 2:n-1
   
   u_1(i)= u1_old(i)+ du_1dt_avg(i).*dt;
   u_2(i)= u2_old(i)+ du_2dt_avg(i).*dt;
   u_3(i)= u3_old(i)+ du_3dt_avg(i).*dt;

end   

% Inlet BC
u_1(1) = a_c(1); 
u_2(1) = 2*u_2(2)-u_2(3);
u_3(1) = u_1(1)*((t_c(1)/(gamma-1)) + ((gamma/2)*(v_c(1)^2)));

% Outlet Bc
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);  


% updating values
% Density
rho_c= u_1./a_c;

% Velocity
v_c= u_2./u_1;

% temperature
t_c= (gamma-1)*((u_3./u_1)-((gamma.*v_c.^2)/2));


% Nondimentional flow properties
% masss flow rate
mass_flow_rate_c= rho_c.*a_c.*v_c;

% Pressure
p_cn= rho_c.*t_c;  

% Mach no.
M_cn= v_c./((t_c).^0.5);


% flow properties @ throat
% masss flow rate
mass_flow_rate_c_th= mass_flow_rate_c(throat_cn);

% Pressure
p_cn_th= p_cn(throat_cn);  

% Mach no.
M_cn_th= M_cn(throat_cn);

% Density
rho_c_th= rho_c(throat_cn);

% Temperature
t_c_th= t_c(throat_cn); 



end

Time_step Programs:

Conservative part

function [t_step_cn]=time_step_cn(x,v_c,dx,n,C,t_c)
  
  for i=1:n
    
    t_step(i) =  C*(dx / (((t_c(i)).^0.5)+v_c(i)));
    t_step_cn = min(t_step(i));
    
  endfor

  endfunction

Non_Conservative part

 function [t_step_ncn]=time_step(x,v,dx,n,C,t)
  
  for i=1:n
    
    t_step(i) =  C*(dx / (((t(i)).^0.5)+v(i)));
    t_step_ncn = min(t_step(i));
    
  endfor

  endfunction

Results are as below:

Here, we perform a grid independence test to find out whether our solution is a function of grid points or not. If we try to test our solution, increase nos. of grid points that mean we are decreasing the value of increment i.e. dx if the value of flow field variables we obtain from the second calculation is different than first calculation then our solution is grid dependence. But in this case, the Values of all Flow field variables are somewhat closer to each other for different grid points.

Graph of Normalised Mass flow rate comparison between Conservative and Non_conservative form. 

Conclusion:

As we can see conservation form does a better job of preserving the mass throughout the flow field, mainly because mass flow itself is one of the dependent variables in the conservation form of governing equation. 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.
Although conservation does a better job for conserving the mass throughout the length of the nozzle that doesn\'t define its superiority, let\'s take a look at the values of primitive variables.

1. The conservation form yields a better mass flow distribution. The conservation form simply does a better job of conserving mass.
2. 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.
3. There is no clear superiority of either form in terms of the 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.
4. Comparing the amount of calculational effort to achieve a solution, we note that the solution of the conservation form requires marginally more work. Most of this is due to the need to decode the primitive variables from the flux variables; such decoding is not necessary when you are solving the nonconservation form.