Code for quasi 1D nozzle flows

Question

The quasi-one-dimensional model describing the flow of compressible gas in rocket nozzles is very common. The corresponding equations have a divergent non dimensional form

$ \frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}=\mathbf{J}$

where

$ U_1=\rho A, U_2=\rho A V, U_3=\rho A (\frac{e}{\gamma -1}+\frac{\gamma}{2}V^2)$

$ F_1=\rho A V, F_2=\rho A V^2 +\frac{1}{\gamma}p A, F_3=\rho V A (\frac{e}{\gamma -1}+\frac{\gamma}{2}V^2)$

$ J_1=J_3=0, J_2=\frac{1}{\gamma}p\frac{\partial A}{\partial x}, p=\rho T.$

Here $e=T$ is an internal energy, $T$ is a temperature, $\rho$ is a density, $V$ is velocity, and $A$ is an area of nozzle cross section.

There are 7 repository results on GitHub for quasi 1D nozzle flows with 3 codes in MATLAB, 2 in Python, 1 in C++, and 1 in Java:

To simulate quasi-one-dimensional compressible Euler equations C++

CFD solution of quasi-1D nozzle flows using finite difference approach (MATLAB)

MATLAB code for computing quasi-1D nozzle flow and sensitivities using CSA

Numerical Solutions of quasi-1d nozzle flows (Python)

Quasi-1D Subsonic-Supersonic Isentropic Convergent-Divergent Nozzle Flow Simulation with MacCormack's Technique (Python)

MATLAB Codes to simulate inference for Transient (Conservative) and Non-Transient (Non-conservative) conditions and plot different graphs of mass flow rates, temperatures, pressures, and densities respectively

Numerical solution of quasi-one-dimensional nozzle flows given in the book "Computational Fluid Dynamics: The Basics with Application" by John D. Anderson, Jr (Java)

But there is nothing for Mathematica. In this post, we demonstrate codes translated for Mathematica and testing using data from the NASA portal. First code translated from here is given by

(*# Global variables*)
gamma = 1.4(* # Specific heat ratio*);
c = 0.5 (*# Courant number*);

(*# Grid generation*)
dx = 0.2 (*# Grid size*);
xstart = 0.0 (* # Starting x-position*);
xstop = 10.0 (* # Stopping x-position*);
x = Range[xstart, xstop, dx];
n = Length[x](* # Number of grid points*);

(*# Nozzle area distribution*)
 A0[x_] := 1 + 2.2*((x - 1.5)^2)(* # Anderson Appendix B*);
half = Round[n/2];
A = Table[If[i <= half, 1.75 - 0.75*Cos[(0.2*x[[i]] - 1.0) Pi],
   1.25 - 0.25*Cos[(0.2*x[[i]] - 1.0)*Pi]], {i, n}]; ListPlot[A]         
(*# Flow field variables
# Note that these are dimensionless*)
rho = 0.9608492 - 
   0.07526785000260212`*
    x (*# Density normalized by reservoir density,rho0 in case of \
Pout/Pin=.16; *);
T = 1 - 0.02314*
    x  (*# Temperature normalized by reservoir temperature,T0*);
V = (0.1 + 0.1*x)*
   Sqrt[T](* # Velocity normalized by stagnation acoustic speed,a0*);
ListPlot[{rho, T, V}]
(*# Time history data*)maxt = 600 (*# Maximum number of time-steps*);
hd = Table[0., {n}, {maxt}];
rho1 = hd; V1 = hd; T1 = hd; p1 = hd; M1 = hd;


zeros = Table[0., {n}];
ddtrho = zeros;
ddtV = zeros;
ddtT = zeros;
ddtrhobar = zeros;
ddtVbar = zeros;
ddtTbar = zeros;

ddtrhoavg = zeros;
ddtVavg = zeros;
ddtTavg = zeros;
dtn = c*dx/(Sqrt[T] + V);
dt = Min[dtn](*# Time-step size determination*);
Do[
  
  (*# Prediction step*)
  
  Do[ddtrho[[i]] = -rho[[i]]*((V[[i + 1]] - V[[i]])/dx) - 
     rho[[i]]*V[[i]]*(Log[A[[i + 1]]] - Log[A[[i]]])/dx - 
     V[[i]]*(rho[[i + 1]] - rho[[i]])/dx;
   ddtV[[
     i]] = -V[[i]]*(V[[i + 1]] - V[[i]])/
       dx - (1/gamma)*(((T[[i + 1]] - T[[i]])/
          dx) + (T[[i]]/rho[[i]])*((rho[[i + 1]] - rho[[i]])/dx));
   ddtT[[
     i]] = (-V[[i]]*(T[[i + 1]] - T[[i]])/dx) - ((gamma - 1)*
       T[[i]]*(((V[[i + 1]] - V[[i]])/
           dx) + (V[[i]]*(Log[A[[i + 1]]] - Log[A[[i]]])/dx)));, {i, 
    2, n - 1}];
  
  (*# Predicted values at next time-step*)
  
  rhobar = rho + (ddtrho*dt);
  Vbar = V + (ddtV*dt);
  Tbar = T + (ddtT*dt);
  
  (*# Correction step*)
  
  
  
  Do[ddtrhobar[[i]] = -rhobar[[i]]*((Vbar[[i]] - Vbar[[i - 1]])/dx) - 
     rhobar[[i]]*Vbar[[i]]*(Log[A[[i]]] - Log[A[[i - 1]]])/dx - 
     Vbar[[i]]*(rhobar[[i]] - rhobar[[i - 1]])/dx;
   ddtVbar[[
     i]] = (-Vbar[[i]]*(Vbar[[i]] - Vbar[[i - 1]])/dx) - (1/
        gamma)*(((Tbar[[i]] - Tbar[[i - 1]])/
          dx) + (Tbar[[i]]/
           rhobar[[i]])*((rhobar[[i]] - rhobar[[i - 1]])/dx));
   ddtTbar[[
     i]] = (-Vbar[[i]]*(Tbar[[i]] - Tbar[[i - 1]])/dx) - ((gamma - 1)*
       Tbar[[i]]*(((Vbar[[i]] - Vbar[[i - 1]])/
           dx) + (Vbar[[i]]*(Log[A[[i]]] - Log[A[[i - 1]]])/
            dx)));, {i, 2, n - 1}];
  
  ddtrhoavg = 0.5*(ddtrho + ddtrhobar);
  ddtVavg = 0.5*(ddtV + ddtVbar);
  ddtTavg = 0.5*(ddtT + ddtTbar);
  
  (*# Update flow field variables*)
  
  rho = rho + ddtrhoavg*dt;
  V = V + ddtVavg*dt;
  T = T + ddtTavg*dt;
  
  (*# Compute pressure from the equation of state*)
  p = rho*T;
  
  (*# Compute boundary points*)
  (*# Subsonic inflow*)
  
  V[[1]] = (2*V[[2]]) - V[[3]];
  (*# Supersonic outflow*)
  rho[[n]] = (2*rho[[n - 1]]) - rho[[n - 2]];
  V[[n]] = (2*V[[n - 1]]) - V[[n - 2]];
  T[[n]] = (2*T[[n - 1]]) - T[[n - 2]];
  
  (*# Compute Mach number*)
  M = V/Sqrt[T];
  
  (*# Data collection*)
  
  rho1[[All, j]] = rho;
  V1[[All, j]] = V;
  T1[[All, j]] = T;
  p1[[All, j]] = p;
  M1[[All, j]] = M;, {j, 1, 600}];

To test this code we use data

data = Drop[
   Import["https://www.grc.nasa.gov/www/wind/valid/cdv/axial.Pex.p16.\
gen", "List"], 7];

p016=Interpreter[DelimitedSequence["Number"]][#] & /@ data;
data1 = Drop[
   Import["https://www.grc.nasa.gov/www/wind/valid/cdv/axial.Mex.p16.\
gen", "List"], 7];

M016=Interpreter[DelimitedSequence["Number"]][#] & /@ data1;

Finally, we show computation results (solid lines) together with p016, M016 (red points) for the Mach number (left) and pressure (right)

{Show[ListLinePlot[Transpose[{x, M}], PlotRange -> {All, {0, 2}}, 
   Frame -> True, ImageSize -> 300, PlotLabel -> "Pin/Pout=0.16", 
   PlotStyle -> Blue, FrameLabel -> {"x", "M"}], 
  ListPlot[M016, PlotStyle -> Red]], 
 Show[ListLinePlot[Transpose[{x, p}], PlotRange -> {All, {0, 1}}, 
   Frame -> True, ImageSize -> 300, PlotLabel -> "Pin/Pout=0.16", 
   PlotStyle -> Blue, FrameLabel -> {"x", "p"}], 
  ListPlot[p016, PlotStyle -> Red]]}

As we can see this is a good correlation for this simple case of supersonic outflow. The problem comes when we test flow with a shock wave in the diffuser. The corresponding data are given by

data = Import[
   "https://www.grc.nasa.gov/www/wind/valid/cdv/axial.Mex.p75.gen", 
   "List"];
M075=Interpreter[DelimitedSequence["Number"]][#] & /@ Drop[data, 7];

data1 = Import[
   "https://www.grc.nasa.gov/www/wind/valid/cdv/axial.Pex.p75.gen", 
   "List"];
p075=Interpreter[DelimitedSequence["Number"]][#] & /@ Drop[data1, 7];

For this problem, we need code with shock capturing method implementation like here. Therefore we translated this code as well for Mathematica users:

(*Note: All parameters are nondimensionalized.*)
dx = 0.2;(*grid size*)
c = 0.5;(*courant number*)
gamma = 1.4;(*Ratio of specific heats*)

(*Discretization of nondimensional distance along the nozzle*)
x = Range[0, 10., dx];
k = Length[x]; nth = Round[k/2];
jmax = 200;

Cx = .05;(*% artificial viscosity arbitrary parameter*)

(*% Parabolic area distribution along the nozzle*)
n = Length[x](* # Number of grid points*);

(*# Nozzle area distribution*)
 half = Round[n/2];
A = Table[If[i <= half, 1.75 - 0.75*Cos[(0.2*x[[i]] - 1.0) Pi],
   1.25 - 0.25*Cos[(0.2*x[[i]] - 1.0)*Pi]], {i, k}]; ListPlot[A] 
(*% Initial Condition*)
P = p075[[All, 2]];(*pressure,p/p0*)
rho = P^(1/gamma);
T = P/rho;
V = M075[[All, 2]]*Sqrt[T];(*Velocity,V/a0*)

mf0 = rho*A*V;(*mass flow rate@t=0*)

ListPlot[{rho, P, T, V, mf0}]
gamma1 = 1/gamma;
gamma2 = gamma - 1;

(* Initial condition of solution vectors*)
U1 = rho*A;
U2 = rho*A*V;
uu3 = T/gamma2 + 0.5 gamma V^2;
U3 = rho*A;
U3 = U3*uu3;

(*% Initial condition of flux vectors*)
F1 = U2;
F2 = (U2^2/U1) + (1 - gamma1)*(U3 - 0.5*gamma*U2^2/U1);
F3 = (gamma*U2*U3/U1) - 0.5*gamma*gamma2*U2^3/U1^2;


{ListPlot[{U1, U2, U3}, PlotRange -> All], ListPlot[{F1, F2, F3}]}

dta = (c*dx)/(V + Abs[T]^0.5);
dt = Min[dta];
resmax = 10^-3;(*% maximum error*)
res = 1;
t = 0;
nstep = 0; M = Table[0, {k}]; rho1 = 
 P1 = Table[0, {k}, {jmax}]; M1 = P1; V1 = P1; T1 = P1;

Do[If[res <= resmax, Break[]];
  
  
  nstep = nstep + 1;
  
  (*% Predictor Step*)
  
  Do[
   J2p[i] = gamma1*rho[[i]]*T[[i]]*((A[[i + 1]] - A[[i]])/dx);
   dU1dtp[i] = -(F1[[i + 1]] - F1[[i]])/dx;
   dU2dtp[i] = -(F2[[i + 1]] - F2[[i]])/dx + J2p[i];
   dU3dtp[i] = -(F3[[i + 1]] - F3[[i]])/dx;, {i, 2, k - 1}];
  Do[num = Abs[P[[i + 1]] - 2*P[[i]] + P[[i - 1]]];
   den = P[[i + 1]] + 2*P[[i]] + P[[i - 1]];
   S1p[i] = (Cx*num*(U1[[i + 1]] - 2*U1[[i]] + U1[[i - 1]]))/den;
   S2p[i] = (Cx*num*(U2[[i + 1]] - 2*U2[[i]] + U2[[i - 1]]))/den;
   S3p[i] = (Cx*num*(U3[[i + 1]] - 2*U3[[i]] + U3[[i - 1]]))/den;, {i,
     2, k - 1}];
  
  
  Do[
   U1p[i] = U1[[i]] + (dU1dtp[i]*dt) + S1p[i];
   U2p[i] = U2[[i]] + (dU2dtp[i]*dt) + S2p[i];
   U3p[i] = U3[[i]] + (dU3dtp[i]*dt) + S3p[i];, {i, 2, k - 1}];
  
  U1p[1] = U1[[1]];
  U2p[1] = U2[[1]];
  U3p[1] = U3[[1]];
  U1p[k] = U1[[k]];
  U2p[k] = U2[[k]];
  U3p[k] = U3[[k]];
  
  
  Do[
   rhop[i] = U1p[i]/A[[i]];
   Vp[i] = U2p[i]/U1p[i];
   Tp[i] = gamma2*(U3p[i]/U1p[i]) - 0.5*gamma*Vp[i]^2;
   Pp[i] = rhop[i]*Tp[i];
   F1p[i] = U2p[i];
   F2p[i] = (U2p[i]^2/
       U1p[i]) + (1 - gamma1)*(U3p[i] - 0.5*gamma*U2p[i]^2/U1p[i]);
   F3p[i] = (gamma*U2p[i]*U3p[i]/U1p[i]) - (0.5*gamma*gamma2*
       U2p[i]^3/U1p[i]^2);, {i, 1, k}];
  
  
  
  (*% corrector Step*)
  
  Do[
   J2c = gamma1*rhop[i]*Tp[i]*((A[[i]] - A[[i - 1]])/dx);
   dU1dtc[i] = -(F1p[i] - F1p[i - 1])/dx;
   dU2dtc[i] = -(F2p[i] - F2p[i - 1])/dx + J2c;
   dU3dtc[i] = -(F3p[i] - F3p[i - 1])/dx;, {i, 1, k - 1}];
  
  
  Do[
   num = Abs[(Pp[i + 1] - 2*Pp[i] + Pp[i - 1])];
   den = Pp[i + 1] + 2*Pp[i] + Pp[i - 1];
   S1c[i] = (Cx*num*(U1p[i + 1] - 2*U1p[i] + U1p[i - 1]))/den;
   S2c[i] = (Cx*num*(U2p[i + 1] - 2*U2p[i] + U2p[i - 1]))/den;
   S3c[i] = (Cx*num*(U3p[i + 1] - 2*U3p[i] + U3p[i - 1]))/den;, {i, 2,
     k - 1}];
  
  
  (*% Average time derivatives*)
  
  Do[
   dU1dtav[i] = 0.5*(dU1dtp[i] + dU1dtc[i]);
   dU2dtav[i] = 0.5*(dU2dtp[i] + dU2dtc[i]);
   dU3dtav[i] = 0.5*(dU3dtp[i] + dU3dtc[i]);, {i, 2, k - 1}];
  
  Do[
   U1[[i]] = U1[[i]] + (dU1dtav[i]*dt) + S1c[i];
   U2[[i]] = U2[[i]] + (dU2dtav[i]*dt) + S2c[i];
   U3[[i]] = U3[[i]] + (dU3dtav[i]*dt) + S3c[i];, {i, 2, k - 1}];
  
  
  (*% Boundary condition@first node*)
  
  (*U1[[2]]=A[[2]];
  U2[[2]]=2*U2[[3]]-U2[[4]];
  ve=U2[[2]]/U1[[2]];
  U3[[2]]=U1[[2]]*((1/gamma2)+0.5*gamma*ve^2);
  
  (*% Boundary condition@last node*)
  
  U1[[k-1]]=2*U1[[k-2]]-U1[[k-3]];
  U2[[k-1]]=2*U2[[k-2]]-U2[[k-3]];
  U3[[k-1]]=2*U3[[k-2]]-U3[[k-3]];(*(0.6784*A[[k]]/gamma2)+(0.5*gamma*
  U2[[k]]^2/U1[[k]]);*)*)
  
  
  Do[
   F1[[i]] = U2[[i]];
   F2[[i]] = (U2[[i]]^2/
       U1[[i]]) + (1 - gamma1)*(U3[[i]] - 0.5*gamma*U2[[i]]^2/U1[[i]]);
   F3[[i]] = (gamma*U2[[i]]*U3[[i]]/U1[[i]]) - 
     0.5*gamma*gamma2*U2[[i]]^3/U1[[i]]^2;, {i, 1, k}];
  
  
  
  (*% corrected values of primitive variables*)
  
  Do[
   rho[[i]] = U1[[i]]/A[[i]];
   V[[i]] = U2[[i]]/U1[[i]];
   T[[i]] = gamma2*((U3[[i]]/U1[[i]]) - 0.5*gamma*V[[i]]^2);
   M[[i]] = V[[i]]*T[[i]]^-0.5;
   P[[i]] = rho[[i]]*T[[i]];, {i, 1, k}];
  
  (*% mass flow rate after t+dt*)
  
  mf = rho*A*V;
  
  res1 = Abs[dU1dtav[nth]];
  res2 = Abs[dU2dtav[nth]];
  res = If[ res1 > res2,
    res1,
    res2];
  
  
  
  Do[
   rho1[[i, j]] = rho[[i]];
   T1[[i, j]] = T[[i]];
   P1[[i, j]] = P[[i]];
   M1[[i, j]] = M[[i]];
   V1[[i, j]] = V[[i]];, {i, 1, k}];
  , {j, 1, jmax}];

Finally, we check this numerical solution using the exact solution M075,p075

{Show[ListLinePlot[Transpose[{x, M}], PlotRange -> {All, {0, 2}}, 
   Frame -> True, ImageSize -> 300, PlotLabel -> "Pin/Pout=0.75", 
   PlotStyle -> Blue], ListPlot[M075, PlotStyle -> Red]], 
 Show[ListLinePlot[Transpose[{x, P}], PlotRange -> {All, {0, 1.2}}, 
   Frame -> True, ImageSize -> 300, PlotLabel -> "Pin/Pout=0.75", 
   PlotStyle -> Blue], ListPlot[p075, PlotStyle -> Red]]}  

As we can see the correlation is not so good even if we start from the exact solution. The question is how we can solve this problem using maybe code written for Mathematica from scratch and with NDSolve?

Update 1 To make research with variable dx in the last code we use initial conditions in a form

pi = Interpolation[p075, InterpolationOrder -> 1];

mi = Interpolation[M075, InterpolationOrder -> 1];

(*% Initial Condition*)
P = pi[x];(*pressure,p/p0*)
rho = P[[1]]^(1/gamma) - (P[[1]]^(1/gamma) - P[[k]]^(1/gamma))/10*x ;
T = P[[1]]/rho[[1]] - (P[[1]]/rho[[1]] - P[[k]]/rho[[k]])/10*x; M = 
 mi[x];
V = M[[1]]*
    Sqrt[T[[1]]] - (M[[1]]*Sqrt[T[[1]]] - M[[k]]*Sqrt[T[[k]]])/10*
    x;(*Velocity,V/a0*)

mf0 = rho*A*V;(*mass flow rate@t=0*) 

Then for dx=0.1, Cx=0.0005, jmax=1200 we have pictures below

Update 2. In a case of subsonic flow we retrieve data from

data = Drop[
   Import["https://www.grc.nasa.gov/www/wind/valid/cdv/axial.Pex.p89.\
gen", "List"], 7];

p089=Interpreter[DelimitedSequence["Number"]][#] & /@ data;
data1 = Drop[
   Import["https://www.grc.nasa.gov/www/wind/valid/cdv/axial.Mex.p89.\
gen", "List"], 7];

M089=Interpreter[DelimitedSequence["Number"]][#] & /@ data1;

Input data for the last code

pi = Interpolation[p089, InterpolationOrder -> 1];

mi = Interpolation[M089, InterpolationOrder -> 1];

(*% Initial Condition*)
P = pi[x];(*pressure,p/p0*)
rho = P[[1]]^(1/gamma) - (P[[1]]^(1/gamma) - P[[k]]^(1/gamma))/10*x ;
T = P[[1]]/rho[[1]] - (P[[1]]/rho[[1]] - P[[k]]/rho[[k]])/10*x; M = 
 mi[x];
V = M[[1]]*
    Sqrt[T[[1]]] - (M[[1]]*Sqrt[T[[1]]] - M[[k]]*Sqrt[T[[k]]])/10*
    x;(*Velocity,V/a0*)

mf0 = rho*A*V;(*mass flow rate@t=0*)

Final results after 690 iteration

Taken into account that we solve quasi 1D problem using McCormack's method this correlation is not bad. But we need something better, since for supersonic outflow we have nice correlation.

Update 3. Third code we translated from here. It is good working for supersonic outflow, and for subsonic outflow as well

(*#setting up initial conditions*)
gamma = 1.4; dx = 0.2;
X = Range[0, 10, 
  dx]; M016 = {{0.0000000, 0.2395428}, {0.2000000, 
   0.2401525}, {0.4000000, 0.2419912}, {0.6000000, 
   0.2450882}, {0.8000000, 0.2494931}, {1.000000, 
   0.2552769}, {1.200000, 0.2625332}, {1.400000, 
   0.2713805}, {1.600000, 0.2819642}, {1.800000, 
   0.2944594}, {2.000000, 0.3090735}, {2.200000, 
   0.3260485}, {2.400000, 0.3456640}, {2.600000, 
   0.3682380}, {2.800000, 0.3941266}, {3.000000, 
   0.4237210}, {3.200000, 0.4574398}, {3.400000, 
   0.4957166}, {3.600000, 0.5389794}, {3.800000, 
   0.5876219}, {4.000000, 0.6419665}, {4.200000, 
   0.7022194}, {4.400000, 0.7684249}, {4.600000, 
   0.8404243}, {4.800000, 0.9178230}, {5.000000, 
   0.9997698}, {5.200000, 1.049278}, {5.400000, 1.099540}, {5.600000, 
   1.150450}, {5.800000, 1.201643}, {6.000000, 1.252754}, {6.200000, 
   1.303419}, {6.400000, 1.353288}, {6.600000, 1.402030}, {6.800000, 
   1.449334}, {7.000000, 1.494914}, {7.200000, 1.538515}, {7.400000, 
   1.579906}, {7.600000, 1.618885}, {7.800000, 1.655275}, {8.000000, 
   1.688924}, {8.200000, 1.719701}, {8.400000, 1.747494}, {8.600000, 
   1.772210}, {8.800000, 1.793771}, {9.000000, 1.812114}, {9.200000, 
   1.827187}, {9.400000, 1.838950}, {9.600000, 1.847373}, {9.800000, 
   1.852435}, {10.00000, 1.854124}}; p016 = {{0.`, 0.9608492`}, {0.2`,
    0.9606546`}, {0.4`, 0.9600656`}, {0.6`, 0.9590641`}, {0.8`, 
   0.9576204`}, {1.`, 0.95569`}, {1.2`, 0.9532128`}, {1.4`, 
   0.9501103`}, {1.6`, 0.946282`}, {1.8`, 0.9416005`}, {2.`, 
   0.9359082`}, {2.2`, 0.9290103`}, {2.4`, 0.920669`}, {2.6`, 
   0.9105966`}, {2.8`, 0.8984537`}, {3.`, 0.8838449`}, {3.2`, 
   0.8663263`}, {3.4`, 0.8454205`}, {3.6`, 0.8206479`}, {3.8`, 
   0.7915765`}, {4.`, 0.7578946`}, {4.2`, 0.7194992`}, {4.4`, 
   0.6765875`}, {4.6`, 0.6297246`}, {4.8`, 0.5798628`}, {5.`, 
   0.5286468`}, {5.2`, 0.4983042`}, {5.4`, 0.4686222`}, {5.6`, 
   0.4395775`}, {5.8`, 0.4114938`}, {6.`, 0.3846425`}, {6.2`, 
   0.3592383`}, {6.4`, 0.3354384`}, {6.6`, 0.3133435`}, {6.8`, 
   0.2930009`}, {7.`, 0.2744158`}, {7.2`, 0.2575561`}, {7.4`, 
   0.2423654`}, {7.6`, 0.2287681`}, {7.8`, 0.2166777`}, {8.`, 
   0.2060031`}, {8.2`, 0.1966533`}, {8.4`, 0.1885405`}, {8.6`, 
   0.1815824`}, {8.8`, 0.1757052`}, {9.`, 0.1708438`}, {9.2`, 
   0.1669425`}, {9.4`, 0.1639555`}, {9.6`, 0.1618473`}, {9.8`, 
   0.1605926`}, {10.`, 0.1601759`}};
k = Length[p016]; gamma = 1.4;
nx = Length[X];
(*% Parabolic area distrubition along the nozzle*)
(*# Nozzle area distribution*)
 half = Round[nx/2];
A = Table[If[i <= half, 1.75 - 0.75*Cos[(0.2*X[[i]] - 1.0) Pi],
    1.25 - 0.25*Cos[(0.2*X[[i]] - 1.0)*Pi]], {i, nx}];

pi = Interpolation[p016, InterpolationOrder -> 1];
mi = Interpolation[M016, InterpolationOrder -> 1];
(*% Initial Condition*)
P = pi[X];(*pressure,p/p0*)k = nx;
rho = P[[1]]^(1/gamma) - (P[[1]]^(1/gamma) - P[[k]]^(1/gamma))/10*X ;
T = P[[1]]/rho[[1]] - (P[[1]]/rho[[1]] - P[[k]]/rho[[k]])/10*X; M = 
 mi[X];
V = M[[1]]*
    Sqrt[T[[1]]] - (M[[1]]*Sqrt[T[[1]]] - M[[k]]*Sqrt[T[[k]]])/10*
    X;(*Velocity,V/a0*)

mf0 = rho*A*V;(*mass flow rate@t=0*)

dt = Min[0.5*dx/(T^0.5 + V)];
nt = 1200;

(*#Variable preDef*)
ini0 = Table[0., {nx}];        (*#predicted values*)
drhodtp = rhop = ini0;

Vp = ini0;
dVdtp = ini0;

Tp = ini0;
dTdtp = ini0;

drhodtc = ini0;     (*#Corrected Values*)
dVdtc = ini0;
dTdtc = ini0;

TinTime = XTP = Table[0., {nx}, {nt}] ;        (*#Used for plotting *)
rhoinTime = XTP;
VinTime = XTP;
PinTime = XTP;
MinTime = XTP;

err = 1;
iterNo = 0;

Do[If[ err <= 10^-4, Break[]];
      iterNo = iterNo + 1;
      rhon = rho;
      Vn = V;
      Tn = T;
      Mn = M;
      
  
      Do[
   drhodtp[[i - 1]] = -rho[[i - 1]]*(V[[i]] - V[[i - 1]])/dx - 
     rho[[i - 1]]*V[[i - 1]]*(Log[A[[i]]] - Log[A[[i - 1]]])/dx - 
     V[[i - 1]]*(rho[[i]] - rho[[i - 1]])/dx;    
   dVdtp[[i - 
      1]] = -V[[i - 1]]*(V[[i]] - V[[i - 1]])/
       dx - (1/gamma)*((T[[i]] - T[[i - 1]])/
         dx + (T[[i - 1]]/rho[[i - 1]])*(rho[[i]] - rho[[i - 1]])/
          dx);  dTdtp[[
     i - 1]] = -V[[i - 1]]*(T[[i]] - T[[i - 1]])/dx - (gamma - 1)*
      T[[i - 1]]*( (V[[i]] - V[[i - 1]])/dx + 
        V[[i - 1]]*(Log[A[[i]]] - Log[A[[i - 1]]])/dx);, {i, 3, nx}];
      rhop = rho + drhodtp*dt;
          Vp = V + dVdtp*dt; 
             Tp = T + dTdtp*dt;
        Do[
   drhodtc[[i]] = -rhop[[i]]*(Vp[[i]] - Vp[[i - 1]])/dx - 
     rhop[[i]]*Vp[[i]]*(Log[A[[i]]] - Log[A[[i - 1]]])/dx - 
     Vp[[i]]*(rhop[[i]] - rhop[[i - 1]])/dx;    
   dVdtc[[i]] = -Vp[[i]]*(Vp[[i]] - Vp[[i - 1]])/
       dx - (1/gamma)*((Tp[[i]] - Tp[[i - 1]])/
         dx + (Tp[[i]]/rhop[[i]])*(rhop[[i]] - rhop[[i - 1]])/dx);  
   dTdtp[[i]] = -Vp[[i]]*(Tp[[i]] - Tp[[i - 1]])/dx - (gamma - 1)*
      Tp[[i]]*( (Vp[[i]] - Vp[[i - 1]])/dx + 
        Vp[[i]]*(Log[A[[i]]] - Log[A[[i - 1]]])/dx);, {i, 2, 
    nx - 1}];
      rho = rho + .5 (drhodtp + drhodtc)*dt;
          V = V + .5 (dVdtp + dVdtc)*dt; 
             T = T + .5 (dTdtp + dTdtc)*dt; 
     
      (*rho[[1]]=1;    T[[1]]=1;*)
  
     
      (*V[[1]]=2*V[[2]]-V[[3]];
      V[[-1]]=2*V[[-2]]-V[[-3]];*)
  
      (*T[[-1]]=2*T[[-2]]-T[[-3]]; rho[[-1]]=2*rho[[-2]]-rho[[-3]];*)
  
  
      
      P = rho*T;
      PinTime[[All, j]] = P;
      
      M = V/Abs[T]^0.5;
      MinTime[[All, j]] = M;
  rhoinTime[[All, j]] = rho; VinTime[[All, j]] = V; 
  TinTime[[All, j]] = T;
      err = Total[Abs[M - Mn]];, {j, nt}];
    
{err, iterNo}

{Show[ListLinePlot[Transpose[{X, M}], PlotRange -> {All, {0, 2}}, 
   Frame -> True, ImageSize -> 300, PlotLabel -> "Pin/Pout=0.16", 
   PlotStyle -> Blue, FrameLabel -> {"x", "M"}], 
  ListPlot[M016, PlotStyle -> Red]], 
 Show[ListLinePlot[Transpose[{X, P}], PlotRange -> {All, {0., 1.2}}, 
   Frame -> True, ImageSize -> 300, PlotLabel -> "Pin/Pout=0.16", 
   PlotStyle -> Blue, FrameLabel -> {"x", "p"}], 
  ListPlot[p016, PlotStyle -> Red]]}

If we put in the last code dt = Min[0.75*dx/(T^0.5 + V)] (it looks like 0.75 is special Courant number for subsonic flow) and use p089, M089 from above, then we have nice fit with exact solution with error of $10^{-4}$ - see picture below

Answer 1

After much experimentation, I decided that I needed to build the code from the scratch using LinearSolve. This code based on implicit step in time and 3-4 order DifferentiationMatrix in space. As test example we take famous Sod shock tube problem. Let consider the Euler equations in a divergent form

As initial state we have $\rho=p=1, 0\le x \le 2 $, and $\rho=0.125, p=0.1, 2\le x\le 4$. This problem has exact solution; therefore, we can compare numerical solution with exact one. As well known, the exact solution defined in 5 regions separated by shock and a centered rarefaction wave - see for example here and here. This code generates numerical solution (run it first!)

Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];

sols[p0_, rho0_, p1_, rho1_] := 
  Module[{L = 4, dt = 0.002, nt = 300, nx = 401, \[Gamma] = 1.4, 
    nu = .0025}, {rvar0, mvar0, evar0, pvar0} = 
    Table[ConstantArray[0, {nx, nt}], {4}];
   {rvar1, mvar1, evar1, pvar1} = MakeVariables[{rr, mm, ee, pp}, nx];
    xgrid = Range[0, L, L/(nx - 1)];
   rvar0[[All, 1]] = rho0; rvar0[[Round[nx/2] ;; nx, 1]] = rho1; 
   pvar0[[All, 1]] = p0; pvar0[[Round[nx/2] ;; nx, 1]] = p1; 
   evar0[[All, 1]] = p0/(\[Gamma] - 1); 
   evar0[[Round[nx/2] ;; nx, 1]] = p1/(\[Gamma] - 1);
   dx2 = 
    NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, 
      DifferenceOrder -> 4]@"DifferentiationMatrix"; 
   dx1 = NDSolve`FiniteDifferenceDerivative[Derivative[1], xgrid, 
      DifferenceOrder -> 3]@"DifferentiationMatrix";
   Do[{rvar, mvar, evar} = {rvar0[[All, i]], mvar0[[All, i]], 
      evar0[[All, i]]}; cdtr = (rvar1 - rvar)/dt; 
    cdtm = (mvar1 - mvar)/dt;
    cdte = (evar1 - evar)/dt;
    eq1 = cdtr + dx1 . mvar1;
    eq2 = 
     cdtm + dx1 . (mvar1 mvar/rvar + pvar1) - nu (dx2 . (mvar1/rvar));
     eq3 = cdte + dx1 . ( mvar/rvar (evar1 + pvar1)) - 
      0 nu dx2 . evar1;
    eq4 = pvar1 - (evar1 - 1/2 mvar1 mvar/rvar ) (\[Gamma] - 1);
    
    eq1[[1]] = rvar1[[1]] - rho0; eq1[[-1]] = rvar1[[-1]] - rho1; 
    eq2[[1]] = mvar1[[1]]; eq2[[-1]] = mvar1[[-1]]; 
    eq4[[1]] = pvar1[[1]] - p0; eq4[[-1]] = pvar1[[-1]] - p1;
    eqs = Join[eq1, eq2, eq3, eq4]; 
    var = Join[rvar1, mvar1, evar1, pvar1]; {vec, mat} = 
     CoefficientArrays[eqs, var]; sol = LinearSolve[mat, -vec];
    rvar0[[All, i + 1]] = Take[sol, nx]; 
    mvar0[[All, i + 1]] = Take[sol, nx + 1 ;; 2 nx]; 
    evar0[[All, i + 1]] = Take[sol, 2 nx + 1 ;; 3 nx]; 
    pvar0[[All, i + 1]] = Take[sol, 3 nx + 1 ;; 4 nx];, {i, nt - 1}];
   ];

sols[1, 1, .1, .125] // AbsoluteTiming

Visualization

dt = .002; mmax = 300; tm = 
 dt mmax; {ListLinePlot[Transpose[{xgrid, rvar0[[All, -1]]}], 
  AxesLabel -> {"x", "\[Rho]"}, PlotLabel -> Row[{"t = " , mmax dt}]],
  ListLinePlot[Transpose[{xgrid, mvar0[[All, -1]]/rvar0[[All, -1]]}], 
  PlotRange -> All, AxesLabel -> {"x", "u"}], 
 ListLinePlot[Transpose[{xgrid, pvar0[[All, -1]]}], PlotRange -> All, 
  AxesLabel -> {"x", "p"}], 
 ListLinePlot[
  Transpose[{xgrid, pvar0[[All, -1]]/rvar0[[All, -1]]/(1.4 - 1)}], 
  PlotRange -> All, AxesLabel -> {"x", "\[Epsilon]"}]}

time = Range[0, mmax dt, dt]; rho = 
 Table[{xgrid[[i]], time[[j]], rvar0[[i, j]]}, {i, Length[xgrid]}, {j,
    Length[time] - 1}];
ListContourPlot[Flatten[rho, 1], Contours -> 20, 
 ColorFunction -> "Rainbow", ContourStyle -> None, 
 FrameLabel -> {"x", "t"}, PlotLabel -> "\[Rho]", 
 PlotLegends -> Automatic]

This code can be used to reproduce exact solution

\[Gamma] = 1.4; \[Rho]1 = 1; p1 = 1; \[Rho]5 = .125; p5 = .1; u1 = 0; \
u5 = 0;

cs1 = Sqrt[\[Gamma] p1/\[Rho]1]; cs5 = 
 Sqrt[\[Gamma] p5/\[Rho]5]; a = (\[Gamma] - 1) (cs5/cs1) (P - 1);
b = Sqrt[2*\[Gamma]*(2*\[Gamma] + (\[Gamma] + 1)*(P - 1))];
P = P /. 
  FindRoot[ 
   P - p1/p5*(1 - a/b)^(2.*\[Gamma]/(\[Gamma] - 1.)) == 0, {P, .5}]

cshock = 
 u5 + cs5*Sqrt[((\[Gamma] - 1 + P*(\[Gamma] + 1))/(2*\[Gamma]))]

(*region 2*)

u2[x_, t_] := 
 2/(\[Gamma] + 1) (cs1 + (x - x0)/t); \[Rho]2[x_, 
  t_] := \[Rho]1 (1 - (\[Gamma] - 1)/2 u2[x, t]/cs1)^(2/(\[Gamma] - 
       1)); p2[x_, t_] := 
 p1 (1 - (\[Gamma] - 1)/2 u2[x, t]/cs1)^(2 \[Gamma]/(\[Gamma] - 1));

(*region 3*)

gg = \[Gamma]; alpha = (gg + 1)/(gg - 1);
ccontact = u1 + 2*cs1/(gg - 1)*(1 - (P*p5/p1)^((gg - 1.)/2/gg));
\[Rho]3 = \[Rho]1*(P*p5/p1)^(1/gg);
u3 = ccontact;
p3 = P*p5;

(*region 4*)

p4 = p3; u4 = u3; \[Rho]4 = (1 + alpha*P)/(alpha + 
     P)*\[Rho]5; cfanright = ccontact - Sqrt[gg*p3/\[Rho]3]    

To compare exact solution with numerical solution shown above we define separation lines as follows

x0 = N@xgrid[[Round[Length[xgrid]/2]]]; x1[t_] := x0 - cs1 t; 
x2[t_] := x0 + cfanright t; x4[t_] := x0 + cshock t; 
x3[t_] := x0 + u3 t;

Finally, we define numerical solution and compare with numerical one (red points)

u[x_, t_] := 
  Piecewise[{{u1, x <= x1[t]}, {u2[x, t], x1[t] < x <= x2[t]}, {u3, 
     x2[t] < x <= x3[t]}, {u4, x3[t] < x <= x4[t]}, {u5, x > x4[t]}}];

p[x_, t_] := 
  Piecewise[{{p1, x <= x1[t]}, {p2[x, t], x1[t] < x <= x2[t]}, {p3, 
     x2[t] < x <= x3[t]}, {p4, x3[t] < x <= x4[t]}, {p5, x > x4[t]}}];

\[Rho][x_, t_] := 
  Piecewise[{{\[Rho]1, x <= x1[t]}, {\[Rho]2[x, t], 
     x1[t] < x <= x2[t]}, {\[Rho]3, x2[t] < x <= x3[t]}, {\[Rho]4, 
     x3[t] < x <= x4[t]}, {\[Rho]5, x > x4[t]}}];

{Show[Plot[u[x, tm], {x, 0, 4}, PlotRange -> {0, 1.}, 
   Exclusions -> None, AxesLabel -> {"x", "u"}], 
  ListPlot[Transpose[{xgrid, mvar0[[All, -1]]/rvar0[[All, -1]]}], 
   PlotRange -> All, PlotStyle -> Red]], 
 Show[Plot[p[x, tm], {x, 0, 4}, PlotRange -> {0, 1.}, 
   Exclusions -> None, AxesLabel -> {"x", "p"}], 
  ListPlot[Transpose[{xgrid, pvar0[[All, -1]]}], PlotRange -> All, 
   PlotStyle -> Red]], 
 Show[Plot[\[Rho][x, tm], {x, 0, 4}, PlotRange -> {0, 1.}, 
   Exclusions -> None, AxesLabel -> {"x", "\[Rho]"}], 
  ListPlot[Transpose[{xgrid, rvar0[[All, -1]]}], PlotRange -> All, 
   PlotStyle -> Red]],
 Show[Plot[p[x, tm]/\[Rho][x, tm]/(gg - 1), {x, 0, 4}, 
   PlotRange -> {1, 3}, Exclusions -> None, 
   AxesLabel -> {"x", "\[Epsilon]"}], 
  ListPlot[
   Transpose[{xgrid, pvar0[[All, -1]]/rvar0[[All, -1]]/(gg - 1)}], 
   PlotRange -> {0, 3}, PlotStyle -> Red]]}