Example 7 from RP-1311

Note

The python script for this example is available in the source/bind/python/cea/samples directory of the CEA repository.

Here we describe how to run example 7 from RP-1311 [1] using the Python API. This is a shock problem example using a mixture of H₂, O₂, and Ar.

First import the required libraries:

import numpy as np
import cea

Use cea.units for unit conversions.

Declare the reactant names, and the moles of each:

reac_names = [b"H2", b"O2", b"Ar"]
moles = np.array([0.05, 0.05, 0.9])

Define the unshocked temperature and pressure in SI units, as well as the inital shock velocity in m/s.

p0 = cea.units.mmhg_to_bar(10.0)  # Unshocked pressure (bar)
T0 = 300.0   # Unshocked temperature (K)
u1 = np.array([1100, 1200, 1250, 1300, 1350, 1400])  # Initial velocity (m/s)

Now instantiate the reactant and product Mixture objects. To create the product Mixture, we pass the list of reactant names along with the flag products_from_reactants=True, which will return the full set of possible product species.

reac = cea.Mixture(reac_names)
prod = cea.Mixture(reac_names, products_from_reactants=True)

Initialize the ShockSolver and ShockSolution objects. We want to compute properties across the incident and reflected shock waves, so we set reflected=True when we initilaize the ShockSolution object. This needs to be done at this point to properly size the variable arrays, but note that this option will need to be set again later when we call solve().

solver = cea.ShockSolver(prod, reactants=reac)
solution = cea.ShockSolution(solver, reflected=True)

Next we convert the reactant moles to weights:

weights = reac.moles_to_weights(moles)

Initialize arrays for storing the solution variables:

n = len(u1)
mach1 = np.zeros(n)
u1_out = np.zeros(n)
P1_out = np.zeros(n)
T1_out = np.zeros(n)
rho1 = np.zeros(n)
H1 = np.zeros(n)
U1 = np.zeros(n)
G1 = np.zeros(n)
S1 = np.zeros(n)
M1 = np.zeros(n)
cp1 = np.zeros(n)
gamma1 = np.zeros(n)
v_sonic1 = np.zeros(n)

u2 = np.zeros(n)
P2 = np.zeros(n)
T2 = np.zeros(n)
rho2 = np.zeros(n)
H2 = np.zeros(n)
U2 = np.zeros(n)
G2 = np.zeros(n)
S2 = np.zeros(n)
M2 = np.zeros(n)
cp2 = np.zeros(n)
gamma2 = np.zeros(n)
v_sonic2 = np.zeros(n)
P2_P1 = np.zeros(n)
T2_T1 = np.zeros(n)
M2_M1 = np.zeros(n)
rho2_rho1 = np.zeros(n)
v2 = np.zeros(n)

u5 = np.zeros(n)
P5 = np.zeros(n)
T5 = np.zeros(n)
rho5 = np.zeros(n)
H5 = np.zeros(n)
U5 = np.zeros(n)
G5 = np.zeros(n)
S5 = np.zeros(n)
M5 = np.zeros(n)
cp5 = np.zeros(n)
gamma5 = np.zeros(n)
v_sonic5 = np.zeros(n)
P5_P2 = np.zeros(n)
T5_T2 = np.zeros(n)
M5_M2 = np.zeros(n)
rho5_rho2 = np.zeros(n)
u5_p_v2 = np.zeros(n)
mole_fractions_incd = {}
mole_fractions_refl = {}

Loop over the initial shock velocities and solve the shock problem, noting again that we set reflected=True to compute the equilibrium reflected shock solution.

i = 0
for u in u1:
    # Solve the shock problem
    solver.solve(solution, weights, T0, p0, u1=u, reflected=True)

    # Store the output
    u1_out[i] = u
    P1_out[i] = p0
    T1_out[i] = T0
    if solution.converged:
        mach1[i] = solution.Mach[0]
        rho1[i] = solution.density[0]
        H1[i] = solution.enthalpy[0]
        U1[i] = solution.energy[0]
        G1[i] = solution.gibbs_energy[0]
        S1[i] = solution.entropy[0]
        M1[i] = solution.M[0]
        cp1[i] = solution.cp_eq[0]
        gamma1[i] = solution.gamma_s[0]
        v_sonic1[i] = solution.sonic_velocity[0]

        u2[i] = solution.velocity[1]
        P2[i] = solution.P[1]
        T2[i] = solution.T[1]
        rho2[i] = solution.density[1]
        H2[i] = solution.enthalpy[1]
        U2[i] = solution.energy[1]
        G2[i] = solution.gibbs_energy[1]
        S2[i] = solution.entropy[1]
        M2[i] = solution.M[1]
        cp2[i] = solution.cp_eq[1]
        gamma2[i] = solution.gamma_s[1]
        v_sonic2[i] = solution.sonic_velocity[1]
        P2_P1[i] = solution.P21
        T2_T1[i] = solution.T21
        M2_M1[i] = solution.M21
        rho2_rho1[i] = 1.0/solution.rho12
        v2[i] = solution.v2

        u5[i] = solution.velocity[2]
        P5[i] = solution.P[2]
        T5[i] = solution.T[2]
        rho5[i] = solution.density[2]
        H5[i] = solution.enthalpy[2]
        U5[i] = solution.energy[2]
        G5[i] = solution.gibbs_energy[2]
        S5[i] = solution.entropy[2]
        M5[i] = solution.M[2]
        cp5[i] = solution.cp_eq[2]
        gamma5[i] = solution.gamma_s[2]
        v_sonic5[i] = solution.sonic_velocity[2]
        P5_P2[i] = solution.P52
        T5_T2[i] = solution.T52
        M5_M2[i] = solution.M52
        rho5_rho2[i] = solution.rho52
        u5_p_v2[i] = solution.u5_p_v2

        if i == 0:
            for prod in solution.mole_fractions:
                    mole_fractions_incd[prod] = np.array([solution.mole_fractions[prod][1]])
                    mole_fractions_refl[prod] = np.array([solution.mole_fractions[prod][2]])
        else:
            for prod in mole_fractions_incd:
                mole_fractions_incd[prod] = np.append(mole_fractions_incd[prod], solution.mole_fractions[prod][1])
            for prod in mole_fractions_refl:
                mole_fractions_refl[prod] = np.append(mole_fractions_refl[prod], solution.mole_fractions[prod][2])

    i += 1

Print out the results:

print()
print("INITIAL GAS (1)")
print()

print("Mach Number1    ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(mach1[i]), end=" ")
    else:
        print("{0:10.3f}".format(mach1[i]))

print("u1, m/s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(u1_out[i]), end=" ")
    else:
        print("{0:10.3f}".format(u1_out[i]))

print("P, bar          ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(P1_out[i]), end=" ")
    else:
        print("{0:10.4f}".format(P1_out[i]))

print("T, K            ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(T1_out[i]), end=" ")
    else:
        print("{0:10.3f}".format(T1_out[i]))

print("Density, kg/m^3 ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4e}".format(rho1[i]), end=" ")
    else:
        print("{0:10.4e}".format(rho1[i]))

print("H, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(H1[i]), end=" ")
    else:
        print("{0:10.3f}".format(H1[i]))

print("U, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(U1[i]), end=" ")
    else:
        print("{0:10.3f}".format(U1[i]))

print("G, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(G1[i]), end=" ")
    else:
        print("{0:10.3f}".format(G1[i]))

print("S, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(S1[i]), end=" ")
    else:
        print("{0:10.3f}".format(S1[i]))

print("M, (1/n)        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(M1[i]), end=" ")
    else:
        print("{0:10.3f}".format(M1[i]))

print("Cp, kJ/kg-K     ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(cp1[i]), end=" ")
    else:
        print("{0:10.3f}".format(cp1[i]))

print("Gamma_s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(gamma1[i]), end=" ")
    else:
        print("{0:10.3f}".format(gamma1[i]))

print("Son. Vel., m/s  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(v_sonic1[i]), end=" ")
    else:
        print("{0:10.3f}".format(v_sonic1[i]))

print()
print("SHOCKED GAS (2) - Incident, Equilibrium")
print()
print("u2, m/s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(u2[i]), end=" ")
    else:
        print("{0:10.3f}".format(u2[i]))

print("P, bar          ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(P2[i]), end=" ")
    else:
        print("{0:10.4f}".format(P2[i]))

print("T, K            ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(T2[i]), end=" ")
    else:
        print("{0:10.3f}".format(T2[i]))

print("Density, kg/m^3 ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4e}".format(rho2[i]), end=" ")
    else:
        print("{0:10.4e}".format(rho2[i]))

print("H, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(H2[i]), end=" ")
    else:
        print("{0:10.3f}".format(H2[i]))

print("U, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(U2[i]), end=" ")
    else:
        print("{0:10.3f}".format(U2[i]))

print("G, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(G2[i]), end=" ")
    else:
        print("{0:10.3f}".format(G2[i]))

print("S, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(S2[i]), end=" ")
    else:
        print("{0:10.3f}".format(S2[i]))

print("M, (1/n)        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(M2[i]), end=" ")
    else:
        print("{0:10.3f}".format(M2[i]))

print("Cp, kJ/kg-K     ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(cp2[i]), end=" ")
    else:
        print("{0:10.3f}".format(cp2[i]))

print("Gamma_s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(gamma2[i]), end=" ")
    else:
        print("{0:10.3f}".format(gamma2[i]))

print("Son. Vel., m/s  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(v_sonic2[i]), end=" ")
    else:
        print("{0:10.3f}".format(v_sonic2[i]))

print()
print("P2/P1           ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(P2_P1[i]), end=" ")
    else:
        print("{0:10.3f}".format(P2_P1[i]))

print("T2/T1           ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(T2_T1[i]), end=" ")
    else:
        print("{0:10.3f}".format(T2_T1[i]))

print("M2/M1           ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(M2_M1[i]), end=" ")
    else:
        print("{0:10.3f}".format(M2_M1[i]))

print("rho2/rho1       ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(rho2_rho1[i]), end=" ")
    else:
        print("{0:10.3f}".format(rho2_rho1[i]))

print("v2, m/s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(v2[i]), end=" ")
    else:
        print("{0:10.3f}".format(v2[i]))

print()
print("MOLE FRACTIONS")
print()
trace_species = []
for prod in mole_fractions_incd:
    if np.any(mole_fractions_incd[prod] > 5e-6):
        print("{0:15s}".format(prod), end=" ")
        for j in range(n):
            if j < n-1:
                print("{0:10.5g}".format(mole_fractions_incd[prod][j]), end=" ")
            else:
                print("{0:10.5g}".format(mole_fractions_incd[prod][j]))
    else:
        trace_species.append(prod)

if len(trace_species) > 0:
    print()
    print("TRACE SPECIES:")
    max_cols = 8
    nrows = (len(trace_species) + max_cols - 1) // max_cols
    for i in range(nrows):
        print(" ".join("{0:15s}".format(trace_species[j]) for j in range(i * max_cols, min((i + 1) * max_cols, len(trace_species)))))

print()
print("SHOCKED GAS (5) - Reflected, Equilibrium")
print()
print("u5, m/s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(u5[i]), end=" ")
    else:
        print("{0:10.3f}".format(u5[i]))

print("P, bar          ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(P5[i]), end=" ")
    else:
        print("{0:10.4f}".format(P5[i]))

print("T, K            ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(T5[i]), end=" ")
    else:
        print("{0:10.3f}".format(T5[i]))

print("Density, kg/m^3 ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4e}".format(rho5[i]), end=" ")
    else:
        print("{0:10.4e}".format(rho5[i]))

print("H, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(H5[i]), end=" ")
    else:
        print("{0:10.3f}".format(H5[i]))

print("U, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(U5[i]), end=" ")
    else:
        print("{0:10.3f}".format(U5[i]))

print("G, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(G5[i]), end=" ")
    else:
        print("{0:10.3f}".format(G5[i]))

print("S, kJ/kg        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(S5[i]), end=" ")
    else:
        print("{0:10.3f}".format(S5[i]))

print("M, (1/n)        ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(M5[i]), end=" ")
    else:
        print("{0:10.3f}".format(M5[i]))

print("Cp, kJ/kg-K     ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(cp5[i]), end=" ")
    else:
        print("{0:10.3f}".format(cp5[i]))

print("Gamma_s         ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(gamma5[i]), end=" ")
    else:
        print("{0:10.3f}".format(gamma5[i]))

print("Son. Vel., m/s  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(v_sonic5[i]), end=" ")
    else:
        print("{0:10.3f}".format(v_sonic5[i]))

print()
print("P5/P2           ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(P5_P2[i]), end=" ")
    else:
        print("{0:10.3f}".format(P5_P2[i]))

print("T5/T2           ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(T5_T2[i]), end=" ")
    else:
        print("{0:10.3f}".format(T5_T2[i]))

print("M5/M2           ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(M5_M2[i]), end=" ")
    else:
        print("{0:10.3f}".format(M5_M2[i]))

print("rho5/rho2       ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(rho5_rho2[i]), end=" ")
    else:
        print("{0:10.3f}".format(rho5_rho2[i]))

print("u5+v2, m/s      ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3f}".format(u5_p_v2[i]), end=" ")
    else:
        print("{0:10.3f}".format(u5_p_v2[i]))

print()
print("MOLE FRACTIONS")
print()
trace_species = []
for prod in mole_fractions_refl:
    if np.any(mole_fractions_refl[prod] > 5e-6):
        print("{0:15s}".format(prod), end=" ")
        for j in range(n):
            if j < n-1:
                print("{0:10.5g}".format(mole_fractions_refl[prod][j]), end=" ")
            else:
                print("{0:10.5g}".format(mole_fractions_refl[prod][j]))
    else:
        trace_species.append(prod)

print()
print("TRACE SPECIES:")
max_cols = 8
nrows = (len(trace_species) + max_cols - 1) // max_cols
for i in range(nrows):
    print(" ".join("{0:15s}".format(trace_species[j]) for j in range(i * max_cols, min((i + 1) * max_cols, len(trace_species)))))

This results in the following output to the terminal:

INITIAL GAS (1)

Mach Number1         3.353      3.658      3.810      3.962      4.115      4.267
u1, m/s           1100.000   1200.000   1250.000   1300.000   1350.000   1400.000
P, bar              0.0133     0.0133     0.0133     0.0133     0.0133     0.0133
T, K               300.000    300.000    300.000    300.000    300.000    300.000
Density, kg/m^3 2.0126e-02 2.0126e-02 2.0126e-02 2.0126e-02 2.0126e-02 2.0126e-02
H, kJ/kg             1.062      1.062      1.062      1.062      1.062      1.062
U, kJ/kg           -65.182    -65.182    -65.182    -65.182    -65.182    -65.182
G, kJ/kg         -1556.265  -1556.265  -1556.265  -1556.265  -1556.265  -1556.265
S, kJ/kg             5.191      5.191      5.191      5.191      5.191      5.191
M, (1/n)            37.654     37.654     37.654     37.654     37.654     37.654
Cp, kJ/kg-K          0.574      0.574      0.574      0.574      0.574      0.574
Gamma_s              1.625      1.625      1.625      1.625      1.625      1.625
Son. Vel., m/s     328.087    328.087    328.087    328.087    328.087    328.087

SHOCKED GAS (2) - Incident, Equilibrium

u2, m/s            666.380    575.193    558.801    546.820    536.956    527.905
P, bar              0.1093     0.1642     0.1872     0.2104     0.2343     0.2591
T, K              1528.516   1816.425   1930.785   2040.819   2147.535   2249.341
Density, kg/m^3 3.3222e-02 4.1988e-02 4.5020e-02 4.7847e-02 5.0600e-02 5.3374e-02
H, kJ/kg           384.032    555.644    626.184    696.556    768.446    841.973
U, kJ/kg            54.945    164.500    210.327    256.838    305.453    356.577
G, kJ/kg         -8121.784  -9579.928 -10165.744 -10731.423 -11281.059 -11805.155
S, kJ/kg             5.565      5.580      5.589      5.600      5.611      5.623
M, (1/n)            38.618     38.612     38.603     38.589     38.566     38.530
Cp, kJ/kg-K          0.588      0.610      0.629      0.659      0.703      0.765
Gamma_s              1.578      1.550      1.529      1.499      1.463      1.422
Son. Vel., m/s     720.739    778.703    797.279    811.953    822.963    830.836

P2/P1                8.200     12.318     14.043     15.781     17.572     19.432
T2/T1                5.095      6.055      6.436      6.803      7.158      7.498
M2/M1                1.026      1.025      1.025      1.025      1.024      1.023
rho2/rho1            1.651      2.086      2.237      2.377      2.514      2.652
v2, m/s            433.620    624.807    691.199    753.180    813.044    872.095

MOLE FRACTIONS

Ar                 0.92305    0.92289     0.9227    0.92236    0.92179    0.92093
H               1.2167e-07  7.121e-06 2.5615e-05 7.6407e-05 0.00019702 0.00044351
H2              2.6593e-06 4.9279e-05 0.00012308 0.00026816 0.00052464  0.0009274
H2O               0.051235    0.05093   0.050603   0.050071   0.049255   0.048091
O               2.8709e-06 5.5898e-05 0.00014199 0.00031455 0.00062675  0.0011323
O2                0.025619    0.02549   0.025366   0.025185   0.024937   0.024628
OH               8.642e-05 0.00057729  0.0010447  0.0017282  0.0026641  0.0038441

TRACE SPECIES:
H2O2            HO2             O3              H2O(L)          H2O(cr)

SHOCKED GAS (5) - Reflected, Equilibrium

u5, m/s            252.707    303.716    315.779    325.157    333.059    339.920
P, bar              0.2590     0.4832     0.5832     0.6880     0.8009     0.9248
T, K              2108.375   2576.766   2712.530   2829.429   2936.955   3038.970
Density, kg/m^3 5.7006e-02 8.6377e-02 9.8543e-02 1.1083e-01 1.2352e-01 1.3693e-01
H, kJ/kg           740.707   1120.094   1267.269   1411.248   1557.667   1707.473
U, kJ/kg           286.301    560.643    675.433    790.470    909.312   1032.110
G, kJ/kg        -11015.930 -13316.978 -13970.349 -14527.742 -15038.011 -15520.148
S, kJ/kg             5.576      5.603      5.617      5.633      5.651      5.669
M, (1/n)            38.578     38.296     38.107     37.896     37.664     37.413
Cp, kJ/kg-K          0.680      1.078      1.277      1.462      1.628      1.768
Gamma_s              1.481      1.310      1.274      1.253      1.239      1.232
Son. Vel., m/s     820.448    855.925    868.334    881.790    896.418    912.129

P5/P2                2.369      2.942      3.115      3.270      3.418      3.570
T5/T2                1.379      1.419      1.405      1.386      1.368      1.351
M5/M2                0.999      0.992      0.987      0.982      0.977      0.971
rho5/rho2            1.716      2.057      2.189      2.316      2.441      2.566
u5+v2, m/s         180.914    321.091    375.421    428.022    479.986    532.175

MOLE FRACTIONS

Ar                 0.92209    0.91534    0.91084     0.9058    0.90023    0.89425
H               0.00012694  0.0029053  0.0055041  0.0088064   0.012802   0.017411
H2              0.00038641  0.0033392  0.0049636  0.0064716   0.007792  0.0088491
H2O               0.049672   0.041577   0.037066     0.0325   0.027947   0.023548
O               0.00045802  0.0045891  0.0074745   0.010813   0.014638   0.018923
O2                0.025057   0.023285    0.02251   0.021717   0.020846   0.019865
OH               0.0022103  0.0089647   0.011639   0.013892   0.015742   0.017153

TRACE SPECIES:
H2O2            HO2             O3              H2O(L)          H2O(cr)

[1]

McBride, B.J., Gordon, S., "Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications II. Users Manual and Program Description: Users Manual and Program Description - 2", NASA RP-1311, 1996. [NTRS](https://ntrs.nasa.gov/citations/19960044559)