Update / Add High pressure transport

[wandadars]

This PR is a WIP

Changes proposed in this pull request

• Updates to the current version of the HighPressureGasTransport class including:
  • cleaner coding
  • switching to mixture average properties as the multicomponent was incomplete and lacking implementations of required base class methods.
  • New class that implements Chung's method from Poling's 2001 book.

I will add examples when the primary coding is complete. The original class had no tests or examples and was marked for deprecation, so I'm hoping that this effort will save it from deletion.

Checklist

• The pull request includes a clear description of this code change
• Commit messages have short titles and reference relevant issues
• Build passes (scons build & scons test) and unit tests address code coverage
• Style & formatting of contributed code follows contributing guidelines
• The pull request is ready for review

[codecov]

Codecov Report

Attention: Patch coverage is 8.77193% with 364 lines in your changes missing coverage. Please review.

Project coverage is 75.51%. Comparing base (026c5c7) to head (834d08d).
Report is 102 commits behind head on main.

Files	Patch %	Lines
src/transport/HighPressureGasTransport.cpp	6.95%	361 Missing ⚠️
...clude/cantera/transport/HighPressureGasTransport.h	40.00%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1704      +/-   ##
==========================================
- Coverage   75.69%   75.51%   -0.19%     
==========================================
  Files         443      445       +2     
  Lines       60971    61719     +748     
  Branches     9552     9662     +110     
==========================================
+ Hits        46154    46607     +453     
- Misses      11786    12054     +268     
- Partials     3031     3058      +27     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[wandadars]

One thing that I noticed is that, while we may utilize the critical-properties.yaml file in the cubic equation of state to obtain critical properties, the input YAML file needs to also have its own "critical-paramters' section with the acentric-factor defined. This is the current way that the transport properties class obtains the quantity.

Also, the difference between using the tabulated values directly from the YAML data versus making a query to the cubic eos asking it for the mixture critical properties Tc, Pc, Vc, does seem to have a strong effect. In the Poling book, example 9-1 for Sulfur Dioxide, the tabulated critical volume is 122 cm^3/mol, but the Peng-Robinson equation of state gives a critical volume of 97.6 cm^3/mol. So we have a 20% difference in the critical properties between tabulated values and EOS-derived values.

[wandadars]

With regards to diffusion coefficients at high pressure, the Funazukuri, 1992 paper method that Poling mentions in the book might be worth checking out.

[wandadars]

One thing that might need to be checked is whether a user has provided an acentric-factor, because if they don't have it included in the "critical-properties" section in the YAML file, a default value of 0 is used quietly I think. That would be undesirable. So, even if a user provides the "acentric-factor" in the "equation-of-state" section, it won't be used by the transport model and a value of 0 is used.

I'm also working on checking the case of what happens if a user only has a subset of species defined as having critical properties (for cases where critical properties of the other mixture gases may not be known). I want to make sure that the mixing rules don't cause any FPEs when the EOS is asked for the critical values of Pc, Tc, Vc for a species that is ideal. I know for the PREOS, you can provide a=0 and b=0 to recover ideal behavior. I just want to make sure it all works out in that situation.

[wandadars]

@speth Is there a general guiding philosphy established for what types of statements should be in loglevel print statements? I have code blocks in the high pressure transport that I'm using for evaluating the operation of the model, but they're just cout blocks and I have them commented out. I was thinking that having them accessible via a loglevel might be useful.

[speth]

It looks like most of these blocks are just printing out parameter values on input to a function, right? I'd be inclined to just delete them. I think once you have this working, they will have served their purpose and won't be needed.

[wandadars]

This implementation is equation-heavy with references to the Poling book. I am trying to find a balance between having all equations replicated in the documentation versus referring to the book. Any thoughts on how I should proceed here?

[speth]

That's a good question. I guess the ideal would be to include all the equations that are needed to explain the calculations. That would provide users with the ability to understand what Cantera's doing without needing to consult a somewhat difficult to obtain reference, and makes it possible to more easily verify that the implementation corresponds to the equations. On the other hand, that may be an excessive number of equations and effort, and anything at all would be an improvement over the current version of this class.

[wandadars]

One issue I have is with regards to a unit test. I want to try and compare the result to something from the Poling book, but the critical volume estimated from the cubic EOS is 20% off of the tabulated/known value, so that generates an error in the viscosity or thermal conductivity.

[wandadars]

@speth I think this is mostly in a working state now. I know you guys have tons of other things to review. I think this is ready if anyone wants to start the process of reviewing the changes.

The general approach that I have taken with these Lucas and Chung method is that often in the Poling book, the mixture-based calculation is identical to the pure-species equation with the only difference being that mixture-average values of the parameters in the equations are used. So I implemented "blind" methods that are the pure-species versions of the low pressure viscosity and high pressure viscosity methods and then I left the calculation of the mixture values and the passing of those mixture values to the "blind" methods up to the main method (viscosity()).

I put the mixture values into a struct so they could be cleanly computed in one place and passed around if needed.

We have an updated Lucas method (The HighPressureGas class) and a new Chung method (ChungHighPressureGas) class.

Both of these classes use the same Takahasi method, so I moved that into a separate standalone function that both classes could call. I also have a standalone Neufeld collision integral function. This is used in the Poling book for estimating the Omega parameter in the viscosity equations. I know Cantera has it's own calculation of that value, but I wasn't sure how to go about getting it properly, so I put the standalone function in for now.

I tried to keep the comments clear and inside the methods, make reference to the equations in the Poling book (pdfs can be found online of the book).

[wandadars]

When looking into the Takahashi result, we can not interpolate the coefficients from the paper as it produces wild results that look nothing like the plots. I do believe that the authors used some sort of curve fit the data points in the table in order to generate the plots from their paper. For us, using the results from the table presented in the Takahashi paper produces piecewise linear segments between the data points. That's the best that can be done I think, without us introducing our own smooth curve fit to the data in the table, which I've always been taught to not do.

[speth]

Thanks for sharing this work-in-progress, @wandadars. I'm glad to see this capability getting some attention, and I like where this is going in terms of providing a well-documented implementation of these methods. I've left a variety of comments, some of which I realize are addressing issues that you've inherited from the existing implementation. I think this would be a good opportunity to clean these up. Besides the in-line comments, I had a couple of general points:

• Please format all lines with a maximum of 88 characters.
• I'd strongly recommend getting Doxygen up and running on your local machine so you can quickly see all the warning messages it emits and view the generated HTML output. I usually just run python -m http.server 8080 -d .\build\doc\html and leave that running so I can view the generated docs at http://localhost:8080 and just use "Ctrl+Shift+R" each time I run scons doxygen.
• Have you given any thought to an example? Comparing properties calculated using these models with the ideal/low-pressure model might be interesting, and would give us a rare example focused on just transport properties.

[speth]

I'm a little confused by this struct. If these are all properties of the mixture, would it be simpler to just make them member variables of the HighPressureGasTransport class?

[wandadars]

Yes, that is true. I was just leaning towards encapsulating those variables into a group to make it clear what variables are being set by the method that computes them and returns the struct. I tend to dislike methods that change the state of an object across a series of variables. I know that essentially the same thing is done when the struct is updated.

[wandadars]

Is there a way to group variables together in the doxygen system?

[ischoegl]

There is (although we don't use it anywhere yet). Not sure whether we should start.

@param x,y,z  These variables do something

[speth]

Should we consider giving the high-pressure model a more specific name (and if so, what would it be)? From my reading of the Poling book, it's not clear to me what set of models you'd want as a "default".

[wandadars]

This is definitely open to suggestions. The high-pressure model is a selection of 3 different models lumped together. Lucas for the viscosity, Takahashi for the diffusion coefficient, and that Ely and Hanley method for the thermal conductivity. It's hard to slap a name on that. With the Chung class, at least it's 2/3rd of Chung's work, so the name makes a little more sense, but even that might be unclear because it's using Takahashi's method for diffusion coefficinets (there is no Chung method in the Poling book for diffusion coefficients sadly).

[wandadars]

Maybe high-pressure-LET and high-pressure-CCT, or high-pressure-LTE and high-pressure-CT Where we just use the first letter of the model name for the viscosity, thermal conductivity, and diffusion coefficient.

The transport class is essentially 3 smaller transport classes in a trench coat.

[speth]

"3 smaller transport classes in a trench coat" is a very apt description. In the long run (but definitely not right now), it might be worth thinking of a different code structure that was more amenable to mixing and matching the individual models. I'm OK with just using high-pressure and high-pressure-Chung (note the capitalization) as the names for now.

[wandadars]

The Ely and Hanley method used by the original high pressure transport class appears to be this one. So at least we know now where it came from. I will update that method now that we have a source for it.

[wandadars]

Thanks for sharing this work-in-progress, @wandadars. I'm glad to see this capability getting some attention, and I like where this is going in terms of providing a well-documented implementation of these methods. I've left a variety of comments, some of which I realize are addressing issues that you've inherited from the existing implementation. I think this would be a good opportunity to clean these up. Besides the in-line comments, I had a couple of general points:

• Please format all lines with a maximum of 88 characters.
• I'd strongly recommend getting Doxygen up and running on your local machine so you can quickly see all the warning messages it emits and view the generated HTML output. I usually just run python -m http.server 8080 -d .\build\doc\html and leave that running so I can view the generated docs at http://localhost:8080 and just use "Ctrl+Shift+R" each time I run scons doxygen.
• Have you given any thought to an example? Comparing properties calculated using these models with the ideal/low-pressure model might be interesting, and would give us a rare example focused on just transport properties.

This was really helpful. Just as the text is, I think it would be great to have it included in that "contributing to Cantera" documentation section. I can't find the section now, but it was that one that you linked to previously that had notes about how to structure an example and things like that.

[wandadars]

I've got most of the suggestions completed. I added a new example that plots Lucas, Chung, and Ideal viscosities (Lucas is just the standard HighPressureGasTransport class). I moved those repetitive calculations up into an init method for the critical properties as well as the other fluid properties that are used in the mixing rules.

[wandadars]

I have completed the documentation/description of the thermalConductivity() method that was implemented in the original HighPressureGasTransport class.

@speth With regards to an example, I was thinking of maybe taking Methane and running a case at 350K over a range of pressures from 101325 Pa up to 60MPa using the high-pressure and high-pressure-chung models. Then also having a data file with a set of NIST data for methane over that range and plotting them together. What's the proper way to handle a data file for an example? Put the file in the samples directory? Something else?

[wandadars]

It took some time, but I added what I think is a detailed description of the Ely and Hanley model to a doxygen header comment.

[speth]

With regards to an example, I was thinking of maybe taking Methane and running a case at 350K over a range of pressures from 101325 Pa up to 60MPa using the high-pressure and high-pressure-chung models. Then also having a data file with a set of NIST data for methane over that range and plotting them together. What's the proper way to handle a data file for an example? Put the file in the samples directory? Something else?

I think that would be a decent example. For the data, the best way to handle it depends a bit on what it is. We've started moving some of the YAML input data used for samples into a separate repository (https://github.com/Cantera/cantera-example-data), with it being installed so you can import it with something like ct.Solution("example_data/whatever.yaml"). For experimental reference data, if its not too large, you can also just embed it in the .py file directly, like we do in the continuous_reactor.py example. A separate file alongside the .py file doesn't work very well as there isn't a good way to make it downloadable with the example script.

[ischoegl]

As a minor comment, this PR would finally (!) close issue #267 ... 🎉

[wandadars]

I ran the chung-high-pressure method and compared it with NIST data for pure methane at the same conditions.

[wandadars]

Some thoughts:

• We can only use this model in conjunction with an equation of state that can provide critical properties.

• The unity lewis number assumption might be applicable for this transport model as well. I'm not totally sure about this, but the lewis number assumption seems like more of a modifier than a full transport model.

• There's that issue of the cubic eos' getting their data from the critical-properties.yaml file, but the transport models read their acentric factor from a section in the standard Cantera input yaml file, and if one isn't found, a default value is quietly used.

• Both high-pressure transport models use that Takahashi method. If we split the models up into separate files, that shared function between them needs to be located elsewhere. I don't want to create an entire header just for that one function.

[speth]

The unity lewis number assumption might be applicable for this transport model as well. I'm not totally sure about this, but the lewis number assumption seems like more of a modifier than a full transport model.

That's true, but I also feel like if you're using a Le=1 transport model, you can't be very concerned about accuracy, so why bother with the extra complexity of the high-pressure model for other properties.

There's that issue of the cubic eos' getting their data from the critical-properties.yaml file, but the transport models read their acentric factor from a section in the standard Cantera input yaml file, and if one isn't found, a default value is quietly used.

I agree that this is worth resolving. I think the place to do it is in an override of getTransportData() method, which can be made virtual. In the high pressure versions of this function, what I'd suggest is to extend the set of cases for places that it looks for a value of the acentric factor, in order:

• a non-zero value in the species GasTransportData object
• A value in the critical-parameters section of the species input
• A value from critical-properties.yaml. For this, I'd suggest a similar structure to how the critPropsDb/dbSpecies variables are handled in RedlichKwongMFTP::initThermo, where we avoid reading the database unless there's a species where we have no data.

[wandadars]

If I have 2 yaml input files that I want to use to test the critical properties parsing, where is the best place to put them? Into test/data?

[ischoegl]

If I have 2 yaml input files that I want to use to test the critical properties parsing, where is the best place to put them? Into test/data?

Yes, correct. That is the place we use to store input data for unit tests.

[wandadars]

Another comparison, looking at all the models.

I'm going to look again at the Chung model to make sure there is not error in it anywhere because I think it should be closer to the NIST data. But overall, they're tracking the nist data trend fairly well.

[wandadars]

I'm looking again at this handling of critical properties for ideal gases in the mixtures. One thing I tried was to pick the values such that the critical properties maintain their respective ratios as they were in the original implementations. The lower values of say Tcrit going from 1e6 to 1e-6 fixes some issues in the Lucas model of viscosity/thermal conductivity, but then introduces issues into the Chung models. Another issue I noticed is that in the original implementation, the speciesCritTemperature() method can return a critical temperature of 0 for a=0 case, which would then cause divide-by-zero FPEs in multiple places in the equation of state coding. I think it's just been safe because often when a=0, b is also zero, and the first branch of the if-statement gets executed returning a crit temperature of 1e6.

I feel like I've opened up a can or worms with this.

[speth]

I feel like I've opened up a can or worms with this.

Yeah, this does seem like a bit of one. I think you are on to something in thinking about how picking values for the critical properties as $a \to 0$ and/or $b \to 0$ relates to the critical compressibility $Z_c$. Specifically, for the Peng-Robinson EoS, the critical properties should always give a constant value for $Z_c$. One fairly clean way to ensure this is the case would be to always calculate $T_c$, $P_c$, and $V_c$ according to the values of $a$ and $b$, but first clip $a$ and $b$ so they are slightly non-negative. This would eliminate the weird order-of-operations dependence of the current checks for negative $a$ and $b$, and ensure that $Z_c$ always has the expected value. The only adjustments would be to how small the values of $a$ and $b$ can set without sending any of the calculations that depend on these values into strange directions.

A similar update should work for RedlichKwongMFTP::calcCriticalConditions as well.

[wandadars]

I'm starting to wonder if having these ideal gases in the models for the viscosity is a good idea with just how sensitive they are to pure species having non-physical critical properties. It's also very strange that it is so hard to locate any references for how to incorporate ideal gases into the real-gas mixture rules.

mixture_vis = X_real * mu_real + X_ideal*mu_ideal

where mu_real = real-gas model estimate of a mixture viscosity of a mixture of all the species that have critical properties. In the model, the mixture fraction would be the fractions of the real gases amongst themselves.

mu_ideal = typical gas dynamics estimates of the viscosity of the ideal mixture of species, where the values are weighted by say the mole fraction of the ideal species amongst themselves.

And X_real is the mole fraction of the real gas species versus all of the species, and X_ideal is the mole fraction of the ideal gases species versus all of the species.

[speth]

mixture_vis = X_real * mu_real + X_ideal*mu_ideal

This is kind of a strong assertion about how the different species behave, and I doubt that it's consistent with the viscosity formulas that are implemented here. Perhaps for now it's best to just stick with the status quo, where the real gas models require parameters for all species (either the critical state or model specific parameters).

[wandadars]

mixture_vis = X_real * mu_real + X_ideal*mu_ideal

This is kind of a strong assertion about how the different species behave, and I doubt that it's consistent with the viscosity formulas that are implemented here. Perhaps for now it's best to just stick with the status quo, where the real gas models require parameters for all species (either the critical state or model specific parameters).

I agree. I don't think these models were crafted with the possibility of a species that has non-physical critical properties. Or something like a Wilke rule would need to be used instead.