Discrepancy in local Mach behavior between CEA and other sources

Question

I'm exploring some designs using the NASA CEA (Chemical Equilibrium with Applications) code and I ran into a discrepancy between CEA predicted values for local Mach number at the exit plane and values given by Mach flows through a duct. I was wondering if anyone here has ideas as to why.

First, I tried running a CEA simulation of a RD-170, as there is a great deal of data about both its nozzle and operating point, but also its combustion chamber:

              THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM

            COMPOSITION DURING EXPANSION FROM FINITE AREA COMBUSTOR

Pin =  3553.4 PSIA
Ac/At =  2.6000      Pinj/Pinf =  1.030413
CASE = RD170__________

            REACTANT                    WT FRACTION      ENERGY      TEMP
                                        (SEE NOTE)     KJ/KG-MOL      K  
FUEL        RP-1                         1.0000000    -24717.700    298.150
OXIDANT     O2(L)                        1.0000000    -12979.000     90.170

O/F=    2.63000  %FUEL= 27.548209  R,EQ.RATIO= 1.294930  PHI,EQ.RATIO= 1.294930

                INJECTOR  COMB END  THROAT     EXIT
Pinj/P            1.0000   1.0637   1.7875   392.12
P, BAR            245.00   230.33   137.06  0.62480
T, K             3859.00  3842.28  3646.34  1810.35
RHO, KG/CU M    1.8325 1 1.7309 1 1.0984 1 1.0662-1
H, KJ/KG         -781.08  -823.45 -1492.37 -6335.03
U, KJ/KG        -2118.04 -2154.16 -2740.17 -6921.06
G, KJ/KG        -42411.4 -42313.3 -40866.4 -25883.6
S, KJ/(KG)(K)    10.7878  10.7982  10.7982  10.7982

M, (1/n)          23.999   24.007   24.297   25.685
(dLV/dLP)t      -1.03490 -1.03479 -1.02985 -1.00005
(dLV/dLT)p        1.5747   1.5754   1.5221   1.0016
Cp, KJ/(KG)(K)    5.3586   5.3746   5.1909   1.9225
GAMMAs            1.1434   1.1430   1.1401   1.2032
SON VEL,M/SEC     1236.4   1233.3   1192.7    839.7
MACH NUMBER        0.000    0.236    1.000    3.969

TRANSPORT PROPERTIES (GASES ONLY)
CONDUCTIVITY IN UNITS OF MILLIWATTS/(CM)(K)

VISC,MILLIPOISE   1.1735   1.1699   1.1295  0.69246

WITH EQUILIBRIUM REACTIONS

Cp, KJ/(KG)(K)    5.3586   5.3746   5.1909   1.9225
CONDUCTIVITY     11.9520  11.9587  11.0581   1.9898
PRANDTL NUMBER    0.5261   0.5258   0.5302   0.6690

WITH FROZEN REACTIONS

Cp, KJ/(KG)(K)    2.0402   2.0394   2.0310   1.8523
CONDUCTIVITY      3.6656   3.6523   3.4912   1.8594
PRANDTL NUMBER    0.6531   0.6533   0.6571   0.6898

PERFORMANCE PARAMETERS

Ae/At                      2.6000   1.0000   36.870
CSTAR, M/SEC               1814.8   1814.8   1814.8
CF                         0.1604   0.6572   1.8365
Ivac, M/SEC                4862.1   2238.9   3508.7
Isp, M/SEC                  291.1   1192.7   3332.9

I estimated the Mach number using a Newton-Raphson approach and the Mach relations for flow through an expansion nozzle. This equation was given by Barrére et al. (1960). Rocket Propulsion, Sec. 2.2.8, and is as follows:

$$ \varepsilon = \frac{1}{M}\Big[\frac{2(1+\frac{\gamma-1}{2}M^2}{\gamma + 1}\Big]^\frac{\gamma+1}{2(\gamma-1)}$$

This model agrees extremely well with CEA regarding values in the chamber - at the combustion end, there was a discrepancy of only $\Delta M = 2\cdot10^{-4}$ or $0.085\%$. The error at the exit plane, however was much higher - $\Delta M = 0.231$ or $5.825\%$.

I reran CEA with different propellants and different parameters in order to see how this error would change - using liquid methane/LOX at $p_c=100\,\text{bar}$, $\varphi = 3.6$, and $\varepsilon=60$, the chamber error was about $0.190\%$ but the exit plane error had drifted to $9.088\%$.

That said, my question is threefold:

1. What is the primary mechanism for this error?
2. Is there any way to sensibly adjust the above Barrére approach to better make it match the CEA estimates (or the other way around)?
3. If it is not suitable to adjust these values, which values would be considered more reliable from a design context?

Answer 1

1. The primary mechanism for your error is likely the assumption of constant $\gamma$ (specific heat ratio) in the isentropic flow equation you listed. In an equilibrium calculation, the composition of the flow at various axial locations will vary, since the variations in temperature and pressure along the nozzle change the equilibrium constants for various chemical reactions. The changes in composition will slightly affect your $\gamma$. More importantly, the 'GAMMAs' value listed in the CEA output is not the exact same as $\gamma$ in your equation. I don't have the derivation/source available to me, but from previous lecture materials, I can tell you that the GAMMAs listed in the output accounts for variations in the speed of sound along the nozzle. To convert: $$\gamma = \frac{-\gamma_s}{(dLV/dLP)t}$$

Accounting for that, the gamma is roughly 1.18 to 1.20 throughout the nozzle. Solving the isentropic flow equation with an average 1.19 value, you get an exit M of 4.128, which is 4% off of the listed CEA output. I believe the remaining discrepancy is because of a constant $\gamma$ value used in the equation and the equilibrium flow condition.

2 and 3. The isentropic flow equation (Barrére) is based around the assumption of a constant $\gamma$ value, which is generally reasonable. The isentropic flow equations are purely thermodynamics/fluids equations, they do not consider any sort of chemical reactions. CEA will be the more accurate of the two. To get an even more accurate computation, you would use a code/application that includes chemical kinetics (factoring in reaction rates and timescales, since the velocity of the flow means that not all reactions reach equilibrium). Essentially, the most accurate representation of reactions in the flow falls somewhere between CEA's 'frozen' (no reactions in nozzle) and 'equilibrium' (all reactions reach equilibrium) cases. CEA is a thermochemistry code only, and does not include kinetics.

To answer question #2, if you wanted to use the isentropic flow equations only, you could improve your accuracy by using an averaged $\gamma$ value, as described above. Also, the equation you listed is a simplification of the full Mach number equation for the case where the first area is A_throat, and the corresponding Mach number is 1 (at the throat). You could divide the nozzle into smaller axial increments, and calculate the exit Mach number at each of those using the $\gamma$ value for that increment, instead of a constant for the whole nozzle. In this case, you would transform the $\epsilon$ in your equation into the ratio of the initial and exit area of that increment. You can do this with the relation $$\frac{A_2}{A_t}*\frac{A_t}{A_1} = \frac{A_2}{A_1} $$. Iterate across small increments with more accurate $\gamma$ values in those increments, and you can probably get a more accurate result. I don't think there's any value in this, though, since CEA is doing that, and more, for you.