Cantera-CEA verification

[Aleksei Ermashkevich]

Hi all! I seek your expert advice on a challenging issue I have encountered in my research on the thermodynamic calculation of the O2(L) + RP-1 fuel mixture, specifically at an O/F ratio of 0.5. For my analysis, I utilized the NASA-9 polynomial coefficients obtained from NASA Thermobuild and the Cantera 3.0.0 module for the calculations.

Upon comparing my results with those generated by the Chemical Equilibrium with Applications (CEA) program, I observed significant discrepancies in the combustion product composition, particularly in the quantity of solid carbon (graphite, C(gr)), and a substantial difference in the specific heat capacity (Cp).

For instance, at a pressure of 10 MPa and an O/F ratio of 0.5, the results are as follows:

Mass fraction of C(gr) = 0.10838 (using Cantera)
Mass fraction of C(gr) = 0.29633 (using CEA)
Cp = 4571.4176 J/(kg*K) (using Cantera)
Cp = 9150.9 J/(kg*K) (using CEA)
These differences raise concerns about the accuracy of my calculations or the potential misunderstanding of the methodologies used. I am keen to understand the root cause of these discrepancies and seek guidance on how to address this issue.

My source code is:

```python
def calculate_thermo(press, temp_ox, temp_fu, k_m):
    import cantera as ct
    # Initialize pure oxidizer solution object in Cantera.
    oxidizer_pure = ct.Solution('reactants.yaml', 'O2(L)')
    oxidizer_pure.TP = temp_ox, press  # Set temperature and pressure for pure oxidizer
   
    # Initialize pure fuel solution object in Cantera.
    fuel_pure = ct.Solution('reactants.yaml', 'RP-1')
    fuel_pure.TP = temp_fu, press  # Set temperature and pressure for pure fuel
   
    molar_ratio = k_m * fuel_pure.mean_molecular_weight / oxidizer_pure.mean_molecular_weight
    moles_oxidizer = molar_ratio / (1 + molar_ratio)
    moles_fuel = 1 - moles_oxidizer
   
    # Initialize gas state for two pure propellants
    gas = ct.Solution('nasa9.yaml')
   
    # Create a mixture of the liquid phases with the gas-phase model, with the number of moles for fuel and oxidizer
    mix = ct.Mixture([(oxidizer_pure, moles_oxidizer), (fuel_pure, moles_fuel), (gas, 0)])
    # Solve for the equilibrium state, at constant enthalpy and pressure
    mix.equilibrate('HP')

    temp = gas.T
    enthalpy = gas.enthalpy_mass
    entropy = gas.entropy_mass
    composition = {species: {'mass_fraction': gas.Y[i], 'mol_fraction': gas.X[i]} for i, species in
                            enumerate(gas.species_names) if gas.Y[i] and gas.X[i] > 1e-09}
    return temp, enthalpy, entropy, composition
```
I also attached 2 files reactants.yaml and nasa9.yaml.

Thank you for your time and consideration. I look forward to your expert advice.

[Ray Speth]

Hi Alex,

The problem with the C(gr) composition is due to including it as a gas phase species, which is non-physical. You need to define a separate solid phase for the graphite and include it in your Mixture definition. For example,

def calculate_thermo(press, temp_ox, temp_fu, k_m): import cantera as ct # Initialize pure oxidizer solution object in Cantera. oxidizer_pure = ct.Solution('reactants.yaml', 'O2(L)') oxidizer_pure.TP = temp_ox, press # Set temperature and pressure for pure oxidizer # Initialize pure fuel solution object in Cantera. fuel_pure = ct.Solution('reactants.yaml', 'RP-1') fuel_pure.TP = temp_fu, press # Set temperature and pressure for pure fuel molar_ratio = k_m * fuel_pure.mean_molecular_weight / oxidizer_pure.mean_molecular_weight moles_oxidizer = molar_ratio / (1 + molar_ratio) moles_fuel = 1 - moles_oxidizer # Initialize gas state for two pure propellants gas = ct.Solution('nasa9-twophase.yaml') # excludes C(gr) from the gas phase species list carbon = ct.Solution('graphite.yaml') # Phase containing C(gr), built into Cantera # Create a mixture of the liquid phases with the gas-phase model, with the number of moles for fuel and oxidizer mix = ct.Mixture([(oxidizer_pure, moles_oxidizer), (fuel_pure, moles_fuel), (gas, 0), (carbon, 0)]) # Solve for the equilibrium state, at constant enthalpy and pressure mix.equilibrate('HP') temp = gas.T enthalpy = gas.enthalpy_mass entropy = gas.entropy_mass phase_masses = [mix.phase_moles(i) * mix.phase(i).mean_molecular_weight for i in range(mix.n_phases)] total_mass = sum(phase_masses) Y_carbon = phase_masses[3] / total_mass return Y_carbon calculate_thermo(10e6, 60, 298, 0.5)

For me, this returns 0.29598, which is in pretty good agreement with CEA.

Regarding the difference in specific heats, this is a question of definitions. Cantera’s definition of cp is the partial derivative of enthalpy with respect to temperature at constant pressure and composition. The value returned by CEA is the derivative while holding the system at equilibrium. You can compute a value like what’s returned by CEA by computing dh/dT using a finite difference method where the mixture is equilibrated at constant T and P at each point before calculating the enthalpy. You can see a related example here: https://cantera.org/examples/python/thermo/sound_speed.py.html.

Regards,
Ray

​