@corykinney Awesome, thanks for taking this on! Please feel free to reach out directly to me if you want to set up a meeting to go over the code structure or how to make additions to C++ 😄

@bryanwweber Thanks! I will most likely take you up on that once I have a chance to experiment a bit more and figure out the right questions to ask.

Regarding the input file structure, I see on the example yaml files for the PR and RK equations of state that the a and b values are specified by the user with the model hard-coded as either Redlich-Kwong or Peng-Robinson. Is it expected that additional models would follow this approach?

I would think it would be most flexible if the critical temperature, pressure (and optionally acentric factor) were included in the yaml file and the user could set the equation of state upon loading the mechanism so that the files could be reused as much as possible, the a and b values required by each model would then be calculated from the critical properties. Otherwise, if I wanted to compare predictions from multiple models, I'd have to make multiple input files originating from the same critical property data to begin with.

Hi @corykinney Glad you are tacklling this! The code is set up to calculate a and b, as needed, from critical properties. You can see here where the input file is processes and looks for either an a and b value first, then looks for critical-parameters here. We also created a database of known or estimated critical parameters (crit-properties.yaml) for a wide range of species, so if no inputs are given we try to look them up as a last-ditch effort.

@decaluwe Thanks for the explanation and the links to the code segments. I was aware of the crit-properties.yaml database lookup, but didn't realize that instead of specifiying a and b for the model that you could also specify the critical-parameters in the file instead based on looking at the example files.

Do you have any thoughts on reusing input files and specifying equation of state when the Solution is initialized? I think it would be very useful to enable comparisons across EOS. For example, performing IDT simulations with a mechanism and doing something like:

gas_ideal = ct.Solution("input.yaml", eos="ideal")
gas_PR = ct.Solution("input.yaml", eos="PR")

I certainly see the value there. It would be a slight departure from current Cantera norms. Currently, to do this, you'd have a phase declaration for each phase, and just load that phase. Which is not a huge amount of extra work, if you just have the repeat the phase declarations in the input file, but can re-use the same species definitions. I actually don't know if this is possible, currently - can the ideal gas model use a species defined for a real gas phase? I'll try and look into it, when I've got a minute.

Likewise thanks for taking this on! I am linking #20 and #107 as there was some discussion that is pertinent. YAML does support multiple thermo models in a single file (the implemented version is option 4 in #20)

YAML does support multiple thermo models in a single file (the implemented version is option 4 in #20)

@ischoegl So if the currenty YAML implementation supports multiple thermo models, would the current approach require copies of the file with different values for the thermo key at the top of the file? So if I wanted to compare equations of state, I'd either have to keep multiple copies of the file or do some file IO to automatically modify it before initializing the Solution objects?

Also, as a side note, I'll take an approach similar to @gkogekar did in documenting her approach with PengRobinson, which I will try to keep updated at corykinney/Soave-Redlich-Kwong. Due to the similarities between the SRK and PR equations, it looks like I should be able to reuse a lot of the code that was written there.

Hi @corykinney -- as in @gkogekar's input file, there, I think the prefered approach is to just have one input file. In this file, each species is defined just once, and each species definition can contain all the necessary info for multiple EoS options.

Then in the phases section, you can have a phase declaration for each type of phase. In your model file, then, you can load each phase separately:

gas_RK = ct.Solution('input.yaml', 'RK_phase')
gas_SRK = ct.Solution('input.yaml', 'SRK_phase')

where RK_phase and SRK_phase, respectively, are the names given in the yaml file for the Redlich-Kwong and Soave-Redlich-Kwong phase declarations, for example.

Hopefully I'm understanding your question correctly - let me know, if not!

To provide a concrete example, the following is how you can currently define a species that can be used in either a Redlich-Kwong or Peng-Robinson thermo model:

- name: H2O
  composition: {H: 2, O: 1}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 1000.0, 3500.0]
    data:
    - [4.19864056, -0.0020364341, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12,
      -30293.7267, -0.849032208]
    - [3.03399249, 0.00217691804, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14,
      -30004.2971, 4.9667701]
  critical-parameters:
    critical-temperature: 647.096
    critical-pressure: 22.064 MPa
    acentric-factor: 0.344
  equation-of-state:
  - model: Peng-Robinson
    binary-a:
      H2: 4 bar*cm^6/mol^2
      CO2: 7.897e7 bar*cm^6/mol^2
  - model: Redlich-Kwong
    a: [1.14972e+14, 0]
    b: 166.7180412

This definition would also work if either the critical-parameters or equation-of-state field was removed (but not both). The model-specific coefficients are used preferentially, with the critical state parameters being used as a fallback. This is the format that is described in the YAML API Docs. You should be able to achieve the same for SRK if you follow structure of PengRobinson::initThermo.

For phase definitions, to enable the use like in @decaluwe's example, the phases section of the input file would look something like:

phases:
- name: RK_phase
  thermo: Redlich-Kwong
  kinetics: gas
  state:
    T: 300.0
    P: 1.01325e+05
- name: SRK_phase
  thermo: Soave-Redlich-Kwong
  kinetics: gas
  state:
    T: 300.0
    P: 1.01325e+05

Given the relatively small differences between the R-K and SRK equations of state, I hope this implementation can be done in a way that keeps as much functionality as possible in a common base class and avoids copy-pasting large chunks of the existing R-K implementation.

@decaluwe @speth Thanks for the information, I didn't know it could be done that way! That's more or less what I was asking if was possible (aside maybe the explicit definition of compatible phases, but that's not inconvenient).

Given the relatively small differences between the R-K and SRK equations of state, I hope this implementation can be done in a way that keeps as much functionality as possible in a common base class and avoids copy-pasting large chunks of the existing R-K implementation.

Right now I'm thinking about getting it working through large amounts of copying with the idea of writing test cases and then working back toward maybe a CubicBase class that could remove the overlap between existing PengRobinson and RedlichKwongMFTP classes (on a related note, I think dropping the MFTP and just having RedlichKwong would be nice, but I'm sure there was probably a reason it wasn't named that way to start with), but that would be a ways away. The SRK is really closer to the PR than the original RK in terms of implementation since they share the very similar temperature dependent parameter $\alpha (T)$.

As you may have already noticed, there is a common base class used by the PengRobinson and RedlichKwongMFTP classes, which is MixtureFugacityTP. As confusing as this name is (at least to me), I think it's probably a fine place for any common functionality that's relevant for all of these models. For instance, this is already where the solver for the cubic EOS is implemented.

I'm not opposed to dropping the incomprehensible suffix from RedlichKwong (it just describes the class inheritance, i.e. mixture fugacity thermo phase), but that does require some care to provide a smooth transition, and I doubt it's worth the effort.

@speth so it can be assumed that the MixtureFugacityTP will only be used by cubic equations of state and thus it would be safe to put all common code there?

At present, I think that's probably a safe assumption.