Fix downstream shock Mach numbers (#122)

* fixed CO derivatives by updated element convergence check

* corrected shock mach values which were incorrectly based on mw ratios

Author: markleader

CHANGELOG.md

@@ -7,6 +7,11 @@ All notable user-visible changes to this project are documented here.
 ### Changed
 
 ### Fixed
+- Corrected downstream shock Mach numbers to use finalized gas speeds and
+  reported sound speeds, including frozen and retained last-valid states.
+  Documented shock reference frames and the legacy reflected-frozen sound-speed
+  basis; `M21`/`M52` remain molecular-weight ratios, with no thermochemical or
+  legacy output changes.
 - Fixed premature convergence of minor elements such as carbon in database Air
   by checking element-relative balances and gas updates. This can add Newton
   iterations and slightly change equilibrium results.

docs/source/examples/shock/example7.rst

@@ -154,7 +154,7 @@ is still not a converged shock point.
             v_sonic2[i] = solution.sonic_velocity[1]
             P2_P1[i] = solution.P21
             T2_T1[i] = solution.T21
-            M2_M1[i] = solution.M21
+            M2_M1[i] = solution.M21  # Legacy molecular-weight ratio, not a Mach ratio
             rho2_rho1[i] = 1.0/solution.rho12
             v2[i] = solution.v2
 
@@ -172,7 +172,7 @@ is still not a converged shock point.
             v_sonic5[i] = solution.sonic_velocity[2]
             P5_P2[i] = solution.P52
             T5_T2[i] = solution.T52
-            M5_M2[i] = solution.M52
+            M5_M2[i] = solution.M52  # Legacy molecular-weight ratio, not a Mach ratio
             rho5_rho2[i] = solution.rho52
             u5_p_v2[i] = solution.u5_p_v2
 

docs/source/examples/shock_examples.rst

@@ -3,7 +3,45 @@ Shock Examples
 
 Here we describe the example shock problems.
 
+Velocity and Mach conventions
+-----------------------------
+
+The Fortran, C, Python, MATLAB, and Excel shock interfaces report stations
+``[initial, incident, reflected]``, numbered ``[1, 2, 5]`` in legacy CEA.
+Speeds are positive magnitudes. The unshocked gas and reflecting wall are
+stationary in the laboratory frame.
+
+* ``velocity[0] = u1`` is the upstream gas speed in the incident-shock frame,
+  also the incident wave speed in the laboratory frame.
+* ``velocity[1] = u2`` is the downstream gas speed in the incident-shock frame.
+  The laboratory-frame gas speed is ``v2 = u1 - u2``.
+* ``velocity[2] = u5`` is the downstream gas speed in the reflected-shock frame.
+  That gas is stationary at the wall, so ``u5`` also equals the reflected
+  wave-speed magnitude in the laboratory frame (the wave travels back toward
+  the unshocked gas). The upstream reflected-frame gas speed is
+  ``u5_p_v2 = u5 + v2``. It is not ``u2``.
+
+Thus the reflected mass balance is ``rho2*(v2 + u5) = rho5*u5`` and
+``u5 = v2/(rho5/rho2 - 1)``. ``Mach`` is ``velocity/sonic_velocity`` in these
+respective shock frames, not a vector of laboratory-frame Mach numbers.
+The upstream reflected-shock Mach number would use ``(v2 + u5)/a2`` with
+the appropriate upstream sound speed, not the reported downstream entry.
+
+Sound speed is reported as ``sqrt(R*T*gamma_s/M)``. Initial and incident-frozen
+states use frozen heat capacities; equilibrium states use the equilibrium
+isentropic exponent. Reflected-frozen output retains the legacy incident-state
+frozen heat-capacity basis and molecular weight, evaluated with the reflected
+temperature in the sound-speed expression. For temperature-dependent heat
+capacities this can differ from a locally evaluated frozen sound speed;
+``Mach`` uses the reported convention. Retained last-valid states use the same
+definition but remain unconverged and need not satisfy shock conservation.
+
+``M21`` and ``M52`` are **molecular-weight ratios**, not Mach-number ratios.
+The legacy ``M2/M1`` and ``M5/M2`` output labels and example variables retain
+this meaning. Mach ratios, if needed, must be formed from ``Mach`` explicitly,
+with due regard to the different reference frames.
+
 .. toctree::
    :maxdepth: 1
 
-   shock/example7
+   shock/example7

source/bind/c/cea_enum.h

@@ -111,6 +111,15 @@
     CEA_ROCKET_FROZEN_PRANDTL           = 25, \
     CEA_ROCKET_EQUILIBRIUM_PRANDTL      = 26
 
+/* Shock arrays are ordered [initial, incident, reflected] (stations 1, 2, 5).
+ * VELOCITY gives gas speeds in the incident frame (1, 2) and reflected frame (5).
+ * Station 5 is at rest at the wall; its speed also equals the reflected wave-speed
+ * magnitude in the wall frame. V2 = u1-u2 is wall-frame gas speed at station 2;
+ * U5_P_V2 is upstream gas speed relative to the reflected wave, not the wall.
+ * MACH = VELOCITY/SONIC_VELOCITY. Reflected-frozen sound speed retains the legacy
+ * incident frozen-Cp basis at reflected T. M21/M52 are molecular-weight ratios,
+ * not Mach ratios. See docs/source/examples/shock_examples.rst for conventions.
+ */
 #define CEA_SHOCK_PROPERTY_TYPE_ENUM \
     CEA_SHOCK_TEMPERATURE              = 0, \
     CEA_SHOCK_PRESSURE                 = 1, \

source/bind/python/CEA.pyx

@@ -3391,7 +3391,8 @@ cdef class ShockSolver:
         p0 : float
             Unshocked pressure in bar
         u1 : float, optional
-            Shock velocity in m/s (mutually exclusive with Mach1)
+            Incident wave speed relative to the stationary unshocked gas in m/s
+            (mutually exclusive with Mach1)
         Mach1 : float, optional
             Shock Mach number (mutually exclusive with u1)
         reflected : bool, default True
@@ -3453,7 +3454,12 @@ cdef class ShockSolution:
     Solution object containing shock wave calculation results.
 
     Stores shock properties at initial, incident, and optionally
-    reflected shock states.
+    reflected shock states (legacy stations 1, 2, and 5).
+
+    Gas speeds and Mach numbers use the incident-shock frame at stations 1
+    and 2 and the reflected-shock frame at station 5; they do not share a
+    single laboratory frame. The unshocked gas and reflecting wall are at
+    rest in the laboratory frame.
 
     Parameters
     ----------
@@ -3509,7 +3515,13 @@ cdef class ShockSolution:
 
     property velocity:
         """
-        Velocity at each shock state in m/s.
+        Gas speed in m/s, in each station's shock frame.
+
+        Entries 0 and 1 are upstream/downstream speeds relative to the incident
+        shock. Entry 2 is downstream speed relative to the reflected shock.
+        Reflected downstream gas is stationary at the wall, so entry 2 also
+        equals the magnitude of the reflected wave speed in the wall frame.
+        The upstream reflected-frame gas speed is ``u5_p_v2``, not entry 1.
 
         Returns
         -------
@@ -3521,7 +3533,12 @@ cdef class ShockSolution:
 
     property Mach:
         """
-        Mach number at each shock state.
+        Gas Mach number: ``velocity / sonic_velocity`` at each shock state.
+
+        Uses the incident-shock frame for entries 0 and 1 and the reflected-shock
+        frame for entry 2. Entry 2 is neither a wall-frame gas Mach number (zero)
+        nor the upstream reflected-shock Mach number. See ``sonic_velocity``
+        for the retained reflected-frozen sound-speed convention.
 
         Returns
         -------
@@ -3533,7 +3550,14 @@ cdef class ShockSolution:
 
     property sonic_velocity:
         """
-        Sonic velocity at each shock state in m/s.
+        Reported sound speed in m/s: ``sqrt(R * T * gamma_s / M)``.
+
+        Initial and incident-frozen states use frozen heat capacities;
+        equilibrium states use the equilibrium isentropic exponent.
+        For compatibility, reflected-frozen output uses the incident state's
+        frozen heat-capacity basis and molecular weight at the reflected
+        temperature. This legacy convention can differ from the local frozen
+        sound speed for a gas with temperature-dependent heat capacities.
 
         Returns
         -------
@@ -3545,7 +3569,7 @@ cdef class ShockSolution:
 
     property rho12:
         """
-        Density ratio across the incident shock (ρ2/ρ1).
+        Density ratio across the incident shock (ρ1/ρ2).
 
         Returns
         -------
@@ -3617,43 +3641,48 @@ cdef class ShockSolution:
 
     property M21:
         """
-        Mach number ratio across the incident shock (M2/M1).
+        Molecular-weight ratio across the incident shock (M2/M1), not a Mach ratio.
 
         Returns
         -------
         float
-            Mach number ratio across the incident shock
+            Legacy CEA2 molecular-weight ratio across the incident shock
         """
         def __get__(self):
             return self._get_scalar_property(SHOCK_M21)
 
     property M52:
         """
-        Mach number ratio across the reflected shock (M5/M2).
+        Molecular-weight ratio across the reflected shock (M5/M2), not a Mach ratio.
 
         Returns
         -------
         float
-            Mach number ratio across the reflected shock
+            Legacy CEA2 molecular-weight ratio across the reflected shock
         """
         def __get__(self):
             return self._get_scalar_property(SHOCK_M52)
 
     property v2:
         """
-        Velocity at the incident shock state in m/s.
+        Station-2 gas speed in the wall/unshocked-gas frame in m/s.
+
+        Equals ``velocity[0] - velocity[1]``; this is not a shock wave speed.
 
         Returns
         -------
         float
-            Velocity at the incident shock state
+            Gas speed behind the incident shock relative to the wall
         """
         def __get__(self):
             return self._get_scalar_property(SHOCK_V2)
 
     property u5_p_v2:
         """
-        Velocity of the reflected shock wave in m/s.
+        Upstream gas speed relative to the reflected shock in m/s.
+
+        Equals ``velocity[2] + v2``. This is also the reflected wave speed
+        relative to station-2 gas, not the wave speed relative to the wall.
 
         Returns
         -------

source/bind/python/cea/lib/libcea.pyi

@@ -689,18 +689,25 @@ class ShockSolver:
 
 
 class ShockSolution:
+    """Stations [1, 2, 5]; speeds/Mach use incident, incident, reflected shock frames."""
     last_error: int
     def __init__(self, solver: ShockSolver, reflected: bool = True) -> None: ...
     @property
     def T(self) -> ScalarOrArray: ...
     @property
     def P(self) -> ScalarOrArray: ...
     @property
-    def velocity(self) -> ScalarOrArray: ...
+    def velocity(self) -> ScalarOrArray:
+        """Gas speeds [u1, u2, u5] in the shock frames; u5 is also the wall-frame wave speed magnitude."""
+        ...
     @property
-    def Mach(self) -> ScalarOrArray: ...
+    def Mach(self) -> ScalarOrArray:
+        """velocity / sonic_velocity, including the legacy reflected-frozen sound-speed basis."""
+        ...
     @property
-    def sonic_velocity(self) -> ScalarOrArray: ...
+    def sonic_velocity(self) -> ScalarOrArray:
+        """sqrt(R*T*gamma_s/M); reflected-frozen uses incident frozen Cp at reflected T."""
+        ...
     @property
     def rho12(self) -> ScalarOrArray: ...
     @property
@@ -714,13 +721,21 @@ class ShockSolution:
     @property
     def T52(self) -> ScalarOrArray: ...
     @property
-    def M21(self) -> ScalarOrArray: ...
+    def M21(self) -> ScalarOrArray:
+        """Incident molecular-weight ratio M2/M1, not a Mach ratio."""
+        ...
     @property
-    def M52(self) -> ScalarOrArray: ...
+    def M52(self) -> ScalarOrArray:
+        """Reflected molecular-weight ratio M5/M2, not a Mach ratio."""
+        ...
     @property
-    def v2(self) -> ScalarOrArray: ...
+    def v2(self) -> ScalarOrArray:
+        """Station-2 gas speed in the wall frame, u1-u2 [m/s]."""
+        ...
     @property
-    def u5_p_v2(self) -> ScalarOrArray: ...
+    def u5_p_v2(self) -> ScalarOrArray:
+        """Upstream gas speed in the reflected-shock frame [m/s]."""
+        ...
     @property
     def volume(self) -> ScalarOrArray: ...
     @property

source/bind/python/cea/samples/rp1311/example7.py

@@ -116,7 +116,7 @@
         v_sonic2[i] = solution.sonic_velocity[1]
         P2_P1[i] = solution.P21
         T2_T1[i] = solution.T21
-        M2_M1[i] = solution.M21
+        M2_M1[i] = solution.M21  # Legacy molecular-weight ratio, not a Mach ratio
         rho2_rho1[i] = 1.0/solution.rho12
         v2[i] = solution.v2
 
@@ -134,7 +134,7 @@
         v_sonic5[i] = solution.sonic_velocity[2]
         P5_P2[i] = solution.P52
         T5_T2[i] = solution.T52
-        M5_M2[i] = solution.M52
+        M5_M2[i] = solution.M52  # Legacy molecular-weight ratio, not a Mach ratio
         rho5_rho2[i] = solution.rho52
         u5_p_v2[i] = solution.u5_p_v2