Using equilibrate fct with virtual optimized kinetic schemes

[Alex Elias]

Hi, 

I'm working with virtual optimized kinetic schemes, composed by virtual reactions and virtual chemical species, which are built through optimization of both kinetic rate parameters and virtual species thermo-chemical properties in order to capture reference target flame quantities (for more details: Mélody Cailler, Nasser Darabiha, Benoît Fiorina - Development of a virtual optimized chemistry method. Application to hydrocarbon/air combustion - http://doi.org/10.1016/j.combustflame.2019.09.013).

I've been working with this methodology in CANTERA. One of the first tests to verify that the first optimization step (evaluate virtual species thermo-chemical properties) yields good results is to calculate the adiabatic flame temperature of a fuel/air mixture (say a stoichiometric methane/air at 300 K and 1 atm).

When I build a mechanism using two products (Opt2P.cti), I can use cantera's equilibrate function normally. 

import cantera as ct
gas = ct.Solution('Opt2P.cti')
gas.TP = 300, ct.one_atm
gas.set_equivalence_ratio(1.0, 'F', 'O2:1.0 N2:3.76')   # F here is equal to CH4 
gas.equilibrate('HP', solver='gibbs')

However, when I extend the formulation to four products (Opt4P.cti), the equilibrate function returns:

"CanteraError thrown by MultiPhase::equilibrate_MultiPhaseEquil:
 No convergence for T"

In both cases, when I implement an enthalpy balance, I get adiabatic temperature values consistent with those obtained with the reference mechanism (gri30.cti) using the equilibrate function.

I've tried other solver options, however the results are not consistent. Only the 'gibbs' solver gives me comparable values to 'gri30' model.

I have attached the optimized mechanisms, in case anyone could take a look. 

NASA coefficients and virtual species masses are optimized to recover averaged properties of the multi-component real ('gri30') mixture.

Does anyone have any idea why the equilibrate function works well with 2 virtual products and not with 4? Am I missing something obvious? 

I appreciate any suggestions!

Regards,

Alex

[Ray Speth]

Hi Alex,

Are the chemical potentials of the virtual species considered in the property optimization? One thing I noticed as a difference between your two mechanisms is that in the `opt4p` mechanism, the species have standard chemical potentials (accessible in Python as `gas.standard_gibbs_RT`) that are several orders of magnitude larger than any species in the GRI 3.0 mechanism, while the `opt2p` mechanism has values that are in the same range as the reference mechanism.

Since the goal of the equilibrium solver is find the composition that minimizes the Gibbs function for the mixture, these values are pretty critical.

Regards,
Ray

[Alex Elias]

Hi Ray,

No, chemical potentials of the virtual species are not considered in the optimization. 

This methodology was originally proposed using only enthalpic balances for 0D and 1D calculations (using an in-house chemkin-based code, if I am not mistaken).  Therefore, the 7th-NASA coefficient was not considered in the studies. Since I am applying this methodology together with CANTERA, it was necessary to extend the optimization process to include the 7th-NASA coefficient. This way, flame speed calculations, for example, could be done directly via FreeFlame.solve (I believe that without the entropy information these calculations are not possible, correct?)

Thank you for taking the time to review my question. I believe the problem may be related to your observations. For some reason that is not yet clear to me, the optimization using two products produces consistent results for the chemical potentials. However, when I use 3 or 4 products, as you pointed out, the chemical potentials skyrocket. From what I could see, the discrepancy comes mainly from the entropy values. So I think I need to add a constraint in the objective function for the values of the chemical potentials of the virtual species. Specifically, on the entropy values.

The first thing that comes to mind is to use some sort of limiter. Since the coefficients are optimized in order, when I optimize the 7th-NASA, I already know the result for all the previous ones. Then I can evaluate the chemical potential of the new species, and if Its potential is greater than a certain value (which could be based on the reference mechanism, such as [-500,500] for gri3.0), I penalize the objective function. I believe that this way the optimizer can converge to a solution with an 'acceptable' chemical potential.

Well, once again, thank you very much for your answer. If you have any ideas related to the constraint for the entropy values of virtual species, I'm all ears.

Regards,

Alex

[Ray Speth]

Hi Alex,

If you’re only using irreversible reactions, it should be possible to use the flame solver without having meaningful entropy (or Gibbs function) values. You would just need to skip the built-in function for generating the initial guess of the flame profile, as that uses an equilibrium calculation to determine the burned gas state. In Python, you can do this by creating a class that inherits from FreeFlame and overrides the set_initial_guess method. The implementation in the base class should provide a guide for how to generate a reasonable initial guess. You could just use the composition corresponding to complete combustion instead of equilibrium as the burned gas.

Regards,
Ray

​

[Alex Elias]

Hi Ray,

I appreciate your guidance. I will definitely work on this suggestion. It seems to me a quick way to solve the problem and continue my research. Later, I can look more carefully at the question of the chemical potentials of the virtual species.

Thank you very much for the feedback. I'll let you know the results here.

Regards,

Alex

[Alex Elias]

Hi Ray, 

I did what you've suggested. I was able to change the flame solver to not use the equilibrate function in evaluating the initial guess of the flame profile. However, I am still having some slight discrepancies in my optimization results for equilibrium temperature and flame speed, which were not observed in the original work I've cited above. In order to identify sources of error in my work, I was looking at CANTERA's source code and I wondered if it is completely ignoring the entropy values, which were absent in my modeling. When calculating the rate of progress, GasKinetics::updateROP() > update_rates_T() > updateKc(), getStandardChemPotentials is used. Isn't CANTERA  looking for the gibbs values of the virtual mechanism? Or is entropy definitely not used by CANTERA in the calculation of the flame speed?

Regards,

Alex