4
$\begingroup$

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?
$\endgroup$
11
  • 1
    $\begingroup$ @uhoh NASA's CEA stands for Chemical Equilibrium with Applications. $\endgroup$
    – ecfedele
    Feb 5, 2022 at 1:03
  • $\begingroup$ @OrganicMarble I didn't readily see a provision for $c$ within Barrére et al.s' methods, so I didn't account for it. Temperature, respectively, was determined as a function $T(M)$ of the local Mach number at each axial position. $\endgroup$
    – ecfedele
    Feb 5, 2022 at 1:56
  • 1
    $\begingroup$ @rockets_go_boom I did at least attempt to correct for that - I did a simple quadratic regression of $\gamma$ at the three stations (C, T, E) and used that as a continuous form $\endgroup$
    – ecfedele
    Feb 6, 2022 at 11:24
  • 1
    $\begingroup$ Exciting to see CEA in use. That was a bullwark of the Apollo program. I did some work on that code with the people from NASA Lewis. I think all the authors are gone. Bonnie McBride was the last as I recall. Used to take a box of cards, a room full of computer to run it, and an 8 hour wait to get a set of answers. Now, a case is finished nearly as quickly as the key is released on my pc. $\endgroup$
    – tckosvic
    Feb 6, 2022 at 21:13
  • 1
    $\begingroup$ @ecfedele -- It is important to also understand how the thermo properties are developed and stored in the program. Your results hinge on those. In the Apollo days those properties received as much attention as the main code as it was all new. Theromodata used to be on a large mag tape. There was a separate code on how to sort out properties to get stored properly in the thermodata file. Perhaps now the properties have been firmly settled on. Did you compile from source code? $\endgroup$
    – tckosvic
    Feb 9, 2022 at 13:09

1 Answer 1

2
$\begingroup$
  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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.