Constant pressure reactor - with liquids

[Chris Neal]

I'm looking into using the Cantera constant pressure homogeneous reactor to model a problem.

I'm trying to assess if it would be feasible to use the constant pressure homogenous reactor to simulate gas-phase reactions from a mixture that starts off as two liquids. The liquids are included in the mechanism as new species, and the conversion from the liquid to the vapor is modeled as an Arrhenius reaction that converts the liquid species to a corresponding vapor species at a particular rate. 

I'd like to understand how the mixture would be handled for a situation like the one described above. Is something like Amagat's law used?

 

Is the ideal gas eos the only supported equation of state for the species? Is there any way to provide an implementation for an alternative eos for a species?

[Ray Speth]

Hi Chris,

Cantera’s reactor models work with any of the equations of state implemented as versions of the ThermoPhase class. However, each reactor can only contain a single, homogeneous phase. If you specify an ideal gas model for the phase, then all species in the phase will be treated as part of an ideal gas.

A better way to handle this might be to try out a feature which was recently merged into the main branch, which is the ability to augment reactor governing equations in your own Python code. There is an Python example at https://github.com/Cantera/cantera/blob/main/interfaces/cython/cantera/examples/reactors/custom2.py, and documentation in the development version docs at https://cantera.org/documentation/dev/sphinx/html/cython/zerodim.html#extensiblereactor.

Regards,
Ray

​

[Chris N]

Thanks for the clarification Ray. To make sure I understand, let's say I go into my Cantera input .yaml file and specify thermo: Peng-Robinson under the "phases" key. Activating this for the phase will require that I have Peng-Robinson parameters for every species in my input file?   There's no way to have some species tied to the peng-robinson with the rest defaulting to the ideal gas?

[Ray Speth]

Hi Chris,

Yes, that’s essentially correct. If you want to use the Peng-Robinson model, you need to provide either the Peng-Robinson parameters (a, b, and acentric factor) or the critical properties (critical temperature, critical pressure, and acentric factor) for each species.

However, you can get individual species to behave as though they were ideal gasses by specifying an arbitrarily low critical temperature and an arbitrarily high critical pressure. This can even be done programmatically by modifying the species data after it’s been loaded (a new feature for Cantera 2.6 that’s in the development version, which you are presumably already using if you’re working with the Peng-Robinson model).

import cantera as ct
species = ct.Species.listFromFile('gri30.yaml')
for S in species:
    if 'critical-parameters' in S.input_data or 'equation-of-state' in S.input_data:
        continue
    S.update_user_data(
        {'critical-parameters':
          {'critical-temperature': 0.001,
           'critical-pressure': 1e20,
           'acentric-factor': 0.0
          }
        }
    )

pr = ct.Solution(thermo='Peng-Robinson', species=species)

In this example, I added a check to skip species where there is already either a critical-parameters field or an equation-of-state field (under the presumption that it contains Peng-Robinson parameters, though you could make that more specific).

One thing I haven’t given thought to so far is whether the binary interaction parameters work out correctly in this limit, once you add some species with “real” Peng-Robinson coefficients.

Regards,
Ray

​

[Chris N]

Thanks Ray,

I came by a similar approach(I think) from examining the Peng-Robinson class' file parsing method.  In the Cantera input yaml file under the `phases` key I changed the `thermo` key to Peng-Robinson, and for all the species that I wanted to be treated as ideal gases I added the following key-value pairs to the individual species entries.

equation-of-state:
      model: Peng-Robinson
      a: 0
      b: 0
      acentric-factor: 0

By setting those constants to zero, the Peng-Robinson equation of state devolves to the ideal-gas equation of state. For the species that I wanted to be treated with the Peng-Robinson equation of state I did not add this section, and instead added an entry for those species into the critical-proprieties.yaml file.

[Ingmar Schoegl]

Interesting. The question whether "the binary interaction parameters work out correctly in this limit, once you add some species with “real” Peng-Robinson coefficients." is an interesting one. In case there aren't major concerns from a technical viewpoint, I wonder whether this may be something where the Peng-Robinson phase parser could just make assumptions if `equation-of-state` is missing.

-ingmar-

[Chris N]

I might be barking up the wrong tree here. I'm primarily interested in a problem that can be summarized as follows.

1.) Take a constant pressure reactor (no heat transfer, no inlets, no outlets)

2.) Set up an initial mixture of two species(call them materialA_liquid and materialB_liquid).

The two liquid materials vaporize at the reactor's set pressure(call the vapor species materialA and materialB). This vaporization process is modeled by adding two reactions to the mechanism file that are forward reactions that take the materialA_liquid and produce materialA with an Arrhenius rate that has parameters that are tuned to match empirical vaporization rates.

Two things are in my mind currently:
1.) At the start I have the two liquids together, and  the Peng-Robinson class appears to take the vapor root for the mixture(densities should be in the hundreds, but I am seeing the density on the order of 1 kg/m^3 in the output). 

2.) On top of that, as the liquid species are transformed to the vapor species, I have a mixture of liquids and gases. But the Peng-Robinson coding appears to compute a mixture-level set of roots and picks either the vapor or liquid root, which wouldn't make any sense for a case where the constituents are both vapor and liquids. I would need some other form of combining the densities of the individual species to get the mixture density.

Would you agree that this might be a problem that I should not use Cantera for?

[Ray Speth]

I'm not sure what useful default assumptions could be made that are better than those currently in place. If the `equation-of-state` field with Peng-Robinson parameters is not present, Cantera next looks for the `critical-parameters` field to read the critical temperature and pressure. If that is not present, it then checks the `critical-properties.yaml` database for a species with the same name. For any of these cases, the Peng-Robinson model gives default formulas for the binary species parameters, which can be overridden by entries in the `equation-of-state` field. I think Chris's use case is a bit of a special one, and it would not make sense generally to assume a = b = omega = 0 for species where no parameters were provided.

Regards,

Ray

[Ray Speth]

Hi Chris,

I think it might help to see a concrete example of what you're implementing, and how you're trying to set the state of the system.

The Peng-Robinson model is definitely able to represent states in the two-phase region. However, once you have a state with multiple phases, the properties that the PengRobinson class will provide are averages for the phase, and don't provide an indication of how the two phases differ (for example, the mole fractions of each species will be different in different phases). As such, I don't think you can calculate reaction rates in such a phase using the GasKinetics model.

From what you've said, I'd reiterate my earlier suggestion of trying the new ExtensibleReactor class, which would give you the flexibility to say that your reactor is basically a gas phase reactor with a couple of additional equations for your liquids and the rates at which they are evaporating.

Regards,

Ray

[Chris N]

Sure thing Ray. I've attached an example C++ file.  It is based on one of the kinetics examples on the Cantera website.

The example was compiled with the following command
g++ -I/home/neal/software/cantera/cantera2.5.1/cantera_build/include/  liquid_reactor.cpp -L/home/neal/software/cantera/cantera2.5.1/cantera_build/lib -lstdc++ -lcantera -pthread -o kinetics

It outputs a .csv file that contains time, temperature, pressure, density, and species mole fractions. I have it set up here to only be pure oxygen, and the density that is coming out in the output file is around 5kg/m^3.  At 80 Kelvin and 101325 Pascals, the NIST chemistry webbook has the density being around 1200 kg/m^3. That's the first issue that I'm running up against.

[Ray Speth]

Hi Chris,

Thanks for providing the example. It seems like Cantera's implementation of the Peng-Robinson model does have some interesting behavior in the two-phase region. I think this comes from one of the early use cases for the closely-related Redlich-Kwong model in Cantera being an interest in the metastable states beyond the saturation curve. So what it does is it tends to stick to the solution of the equation of state that's closest to the current molar volume. If you want to force it onto the liquid branch, you would need to set the temperature much lower first, and then work your way back up. I think this is something that should probably be revisited for Cantera's Peng-Robinson and Redlich-Kwong models.

But beyond that, I still think you need to look at a modification of the reactor governing equations. The existing Reactor classes are for homogeneous mixtures, and you are specifically looking at an inhomogeneous mixture. You hadn't mentioned that you were working in C++, but in this case the route to implementing your own equations is fairly straightforward (insomuch as anything in C++ is straightforward) -- you just need to implement a class derived from `Reactor` and implement the modified governing equations that way.

Another option might be to represent your liquids as being separate reactors, connected to the gas phase by mass flow controllers. You would then use the evaporation rate to set the mass flow from the liquid to vapor phase. In this case, you could actually use a different (condensed phase) equation of state for the liquid phase, but with the same species names for the gas phase.

Regards,

Ray

[Chris N]

Thanks for the suggestions Ray. I'll look into them. I'm definitely committing modeling crimes by my treatment of the liquid phases as just being another reactive species in a homogeneous mixture. For the time being, I'm ok with that. I'm looking to track a mixture fraction of the system as it evolves over time in the homogeneous reactor.


The text below contains some of my observations I had while trying to understand the logical flow through the Phase classes.

I put some print statements into the PengRobinson class and saw that the densityCalc method was always being called with the phaseRequested being the gas phase FLUID_GAS(=-1). So I tracked the calls to the densityCalc method to be coming from the class that the PengRobinson class derived from,  MixtureFugcaityTP. The setPressure method of the MixtureFugcacityTP class appears to take a pressure and the current temperature and calls the densityCalc method. The initialization of the MixtureFugacityTP class sets the forcedState_ variable to be FLUID_UNDEFINED and so the code always passes through the first section of the if-else block in that method. The iState_ variable is initialized to be FLUID_GAS, and so when the densityCalc method is called, the iState_ variable is passed which requests that the gas root from the cubic EOS be returned, and that density is set as the new density.

Moving to the phaseState method in the same class

I placed a print statement in the setPhase method and this is the output that I see during execution. Note, the "t, rho, rhoMid" line is the one in the setPhase method. The other outputs are in the setPressure method. I don't know if it matters, but the initial density is 0.001, which appears to be the value that the density variable is initialized to when an object is instantiated. I don't have a grasp on the process that the objects go through from instantiation to finally being used, but it appears that the setPhase result being called with such a low initial density would guarantee that the iState_ variable would be set as a FLUID_GAS. I'm probably wrong, but it felt like there was a bit of circularity in some of the logic if a state started out as a gas, it seems like subsequent calls to these methods would be self-fulfilling and have a hard time swapping to the liquid. Which might just be a re-wording of what you were saying about it having a gas bias.

Constant-pressure ignition of a Liquid  mixture
beginning at T = 300 K and P = 1 atm.
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 298.15, 0.001, 699.886, 823.843, 30.2258, -0.00116778, 567, 265.441
We're in the setPressure method
We're in the FLUID_UNDEFINED branch
T, P, rhoNow, iState_, forcedState_: 298.15, 101325, 0.001, -1, -3
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 298.15, 0.001, 699.886, 823.843, 30.2258, -0.00116778, 567, 265.441
Num Solutions: 3
Phase Requested: -1
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 298.15, 1.94429, 699.886, 823.843, 30.2258, -4137.76, 567, 265.441
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 300, 1.94429, 695.561, 822.843, 30.2258, -4166.27, 567, 265.441
We're in the setPressure method
We're in the FLUID_UNDEFINED branch
T, P, rhoNow, iState_, forcedState_: 300, 101325, 1.94429, -1, -3
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 300, 1.94429, 695.561, 822.843, 30.2258, -4166.27, 567, 265.441
Num Solutions: 3
Phase Requested: -1

I also tested the coding with O2 at 80K and 1atm (crit temperature is 154K), and the results are similar to my other case. The main difference is that the algorithm seems to think that the O2 is a supercritical fluid at the state that is set.

Constant-pressure ignition of a Liquid O2 mixture
beginning at T = 80 K and P = 1 atm.
We're in the setPressure method
We're in the FLUID_UNDEFINED branch
T, P, rhoNow, iState_, forcedState_: 300, 101325, 0.001, -2, -3
Num Solutions: 1
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 80, 0.653131, 665.693, 1506.76, 0.00133693, -272.429, 154.58, 408.431
We're in the setPressure method
We're in the FLUID_UNDEFINED branch
T, P, rhoNow, iState_, forcedState_: 80, 101325, 0.653131, -1, -3
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 80, 0.653131, 665.693, 1506.76, 0.00133693, -272.429, 154.58, 408.431
Num Solutions: 3
Phase Requested: -1
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 80, 5.22694, 665.693, 1506.76, 0.00133693, -15362.6, 154.58, 408.431
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 80, 5.22694, 665.693, 1506.76, 0.00133693, -15362.6, 154.58, 408.431
t, rho, rhoMid, densLiqTmid, densGasTmid, dpdv, tcrit, rhocrit: 80, 5.22694, 665.693, 1506.76, 0.00133693, -15362.6, 154.58, 408.431

[Chris N]

I've had a little more time to examine the classes at play here.  What I'm trying to do is to take a mixture of non-ideal gasses, and assume that they make up an ideal mixture. So the cross-interactions between the species is not considered. I believe not specifying the binary interaction parameters will default to this type of case.

I threw together a simple test-case for just the PengRobinson class to set a TP and make a density query for a species of pure O2. I based it on the c++ unit tests for the PengRobinson class. Attached is the code and files for running the test case.

I got the following results:

T=300K, P=101325 Pa,   rho_preos = 0.653 kg/m^3.......NIST chemistry webook:  rho = 1.3 kg/m^3
T=1000K, P=101325 Pa, rho_preos = 0.195464 kg/m^3........NIST chemistry webbook: rho = 0.38986 kg/m^3.....ideal gas: rho=0.704 kg/m^3

I thought that there would be more agreement between the two results. Am I not properly specifying my problem?

[Chris N]

Disregard that last message. I had the mole fractions vector mixed up and was querying CH4 instead of O2. When I switched to the O2 the densities matched up for the vapor phase.

[Chris N]

Ok. I've spent more time looking into this. I just want get confirmation that the Peng Robinson eos as implemented can not return any information regarding states that are "under the bell".  Is that a correct assessment?

[Chris N]

I want to step back and try and re-visit my original question with some additional context.

The attached document details the modeling approach that I'm interested in trying to replicate using Cantera. The elements that I believe to be most relevant are shown below.

• Section 2.1.3 on page 33 of the document explains the governing equation for the homogeneous reactor, which seems to match the forms used in Cantera for the constant pressure and constant volume reactors.

• The section titled "Homogeneous Liquid-Gas Model using Amagat's Law" on page 34 details how the mixture properties are determined for the homogeneous reactor. It also explains how the liquids are handled with respect to the equation of state for them.
• Section 2.3 on page 41 of the document explains the liquid vaporization Arrhenius rate approach to handling vaporization in the homogeneous reactor.

So the governing reactor equations look to be a match to Cantera's. 

I crafted the Arrhenius rate vaporization into the yaml mechanism file, and the it seems to be handled by Cantera with the liquid phase being converted to the vapor phase at a particular rate. It's mainly the density of the mixture that is off.

In my attempt to set up the test-case problem I ran into the issue of the initial density of the liquid phase being reported as being 2 kg/m^3 instead of 880 kg/m^3 when using the Peng-Robinson equation of state. I believe that this might be related to the algorithm in the MixtureFucaityTP class that determines what phase(liquid/gas) the system is currently in.  I did see that if say at 1 atm, and 293K, the density is 880 kg/m^3 that the PREOS in Cantera would return liquid values once the temperature was set to 130K. But marching back up to 293K using the "setState_TP" resulted in the gas density being reported once the temperature got over 132K.  

[Ray Speth]

Hi Chis,

I think your assessment regarding the Peng-Robinson implementation is correct -- the way the algorithm is currently structured, I think it is able to calculate the homogeneous metastable states under the dome, but not the states where the mixture has actually separated into two phases.

While the overall governing equations for conservation of species / enthalpy are the same, if you want to use a different implementation of the equation of state based on Amagat's law, then you're going to have to implement something yourself to do that. I don't think Cantera's implementation of Peng-Robinson is what you're looking for.

Regards,

Ray

[Chris N]

Thanks Ray. I think we've moved away from the Amagat's law approach.

I have been looking at the MixtureFugacityTP class with the goal of maybe adjusting the logic in there that picks states to maybe have it temporarily return the liquid roots consistently.  Two general coding questions that I had regarding the class structure are:

1.)  Is there a conscious design choice to have methods that utilize changing the value of arguments passed to them over methods that simply return results via the standard "return" feature? 

2.) Is there a reason for why there is coding that only updates states if they change? Seems like this accounting for the current state and having to update if changed is a bit of a pain.

[Ray Speth]

Hi Chris,

Regarding your first question, I think all the cases where the MixtureFugacityTP class uses mutable arguments are for functions that need to calculate several values. C++ doesn't really have a nice syntax for returning multiple values (until you start using C++17 features, but Cantera doesn't require a C++17-compatible compiler, so we can't use those features yet).

As to your second question, you're right that keeping track of what values have and haven't been updated can be a bit tricky, but it's pretty common to have to deal with in Cantera. The reason is that there are many cases where the same intermediate properties need to be queried for a given value of the state. For example, in a reactor calculation, at each evaluation of the reactor governing equations you would generally need to evaluate the pressure, the internal energy, and also the species chemical potentials (to calculate reverse rate constants), and it's more efficient if you don't recalculate all of these values in such cases.

Regards,

Ray