-
-
Notifications
You must be signed in to change notification settings - Fork 348
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Update / Add High pressure transport #1704
base: main
Are you sure you want to change the base?
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1704 +/- ##
==========================================
+ Coverage 73.10% 73.67% +0.56%
==========================================
Files 383 383
Lines 54632 54900 +268
Branches 9106 9122 +16
==========================================
+ Hits 39940 40448 +508
+ Misses 11692 11431 -261
- Partials 3000 3021 +21 ☔ View full report in Codecov by Sentry. |
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. |
819f302
to
720e601
Compare
With regards to diffusion coefficients at high pressure, the Funazukuri, 1992 paper method that Poling mentions in the book might be worth checking out. |
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. |
@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. |
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. |
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? |
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. |
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. |
@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). |
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. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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 runscons 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.
//These are the parameters that are needed to calculate the viscosity using the Lucas method. | ||
struct LucasMixtureParameters | ||
{ | ||
double FQ_mix_o; | ||
double FP_mix_o; | ||
double Tr_mix; | ||
double Pr_mix; | ||
double Pc_mix; | ||
double Tc_mix; | ||
double MW_mix; | ||
double P_vap_mix; | ||
}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is there a way to group variables together in the doxygen system?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is (although we don't use it anywhere yet). Not sure whether we should start.
@param x,y,z These variables do something
src/transport/TransportFactory.cpp
Outdated
@@ -46,6 +46,7 @@ TransportFactory::TransportFactory() | |||
addDeprecatedAlias("water", "Water"); | |||
reg("high-pressure", []() { return new HighPressureGasTransport(); }); | |||
addDeprecatedAlias("high-pressure", "HighP"); | |||
reg("high-pressure-chung", []() { return new ChungHighPressureGasTransport(); }); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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".
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"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.
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. |
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. |
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. |
I have completed the documentation/description of the @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? |
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. |
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 |
As a minor comment, this PR would finally (!) close issue #267 ... 🎉 |
834d08d
to
da50a0c
Compare
Some thoughts:
|
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.
I agree that this is worth resolving. I think the place to do it is in an override of
|
675c8d0
to
1b334f8
Compare
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. |
…vate class variables
…ies.yaml, if needed
…st that would never work
…the method as well
… high pressure transport model
…common species to crit database, and updated keywords file to prevent CI failure for new example
2079bf4
to
d0e2f75
Compare
Yes, that's likely what happened, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi @wandadars, thanks for these updates. I think this PR is very nearly done. Most of my comments below are related to a few minor tweaks to the input data files and the new example.
In addition, there are a couple of documentation updates that need to be made:
- Please add an entry for the
high-pressure-Chung
model todoc/sphinx/reference/transport/index.md
and revise the entry forhigh-pressure
to explain the differences between the two models. - Please add an entry for the new model in the YAML documentation,
doc/sphinx/yaml/phases.rst
, under thetransport
section.
Update: The one other request I'd make is for a few additional tests to improve coverage of the newly introduced methods. Right now, there are no tests of any of the newly introduced methods except incidentally through the tests that check that the input data is being handled as expected. While you've shown that it's hard to get a good "ground truth" reference value to compare to due to the limitations of the thermo models, we could at least introduce some regression tests, with the comparison values marked clearly as such and perhaps reference values from the textbook or NIST with a short explanation of why we don't expect to be able to match them exactly. If you can do this for both the new class as well as the HighPressureGasTransport
class, I think we can finally mark #267 as completed.
I'm going to get rid of the |
I added two new tests to |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the updates, @wandadars. Most of my comments here are about the documentation. I'd suggest looking through the generated HTML files to check these issues as you fix them.
I think the new tests are good -- a regression test that calculates values similar to what we're showing in the new example is a reasonable choice, and that example shows the degree of discrepancy between these models and the NIST data.
The one significant issue seems to be the one with the mixture-averaged diffusion coefficients. Right now, these appear to be identical to the ideal gas model because of the way the calculations are structured.
And one minor issue: please check the line lengths in HighPressureGasTransport.cpp
to keep them to our 88 character limit -- there are some sections where this is being broken by a lot.
* Where @f$ \eta_i^* @f$ is the referred to as the dilute gas viscosity and is the | ||
* component that is temperature dependent described in @cite ely-hanley1981 | ||
* (see elyHanleyDilutePureSpeciesViscosity() ), | ||
* @f$ MW_i @f$ is the molecular weight of species i, @f$ f_int @f$ is a constant |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
* @f$ MW_i @f$ is the molecular weight of species i, @f$ f_int @f$ is a constant | |
* @f$ MW_i @f$ is the molecular weight of species i, @f$ f_{int} @f$ is a constant |
* with @f$ h_{ij} @f$ defined by Equation 11 in @cite ely-hanley1983 : | ||
* | ||
* @f[ | ||
* h_{ij} = \frac{1}{8} (h_i^{1/3} + h_j^{1/3})^3 (1 - l_{ij}) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
* h_{ij} = \frac{1}{8} (h_i^{1/3} + h_j^{1/3})^3 (1 - l_{ij}) | |
* h_{ij} = \frac{1}{8} \left( h_i^{1/3} + h_j^{1/3} \right)^3 (1 - l_{ij}) |
* Where @f$ MW_0 @f$ is the molecular weight of the reference fluid, @f$ MW_m @f$ | ||
* is the molecular weight of the mixture, @f$ f_{m,0} @f$ and @f$ h_{m,0} @f$ are |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it would be good to be very clear that
* @note The correlations for the reference fluid viscosity return a value of | ||
* micro-grams/cm/s and the parts of the thermal conductivity that utilize | ||
* the correlation fit parameters give values in units of mW/m/K. Appropriate | ||
* conversions were made in the relevant functions. | ||
* @see elyHanleyDiluteReferenceViscosity() | ||
* @see elyHanleyReferenceThermalConductivity() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you want to refer to these other two methods within the @note
, you shouldn't use the special @see
directive -- this causes them to be formatted in the HTML output as just general "see also" methods associated with the whole method, not specifically this note.
* viscosity at high pressure is computed using the pure-fluid high pressure | ||
* relation of Lucas with the difference being that mixture values of the | ||
* model parameters are used. These mixture values are computed using the mixing | ||
* rules described in @see computeMixtureParameters() . |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please check the resulting HTML documentation -- I don't think you want to use the @see
directive here.
phase entry to `high-pressure`. Implemented by class {ct}`HighPressureGasTransport`. | ||
|
||
Chung high-pressure Gas | ||
: A model for high-pressure gas transport properties that usesthe method of Chung |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
: A model for high-pressure gas transport properties that usesthe method of Chung | |
: A model for high-pressure gas transport properties that uses the method of Chung |
# Create the gas object using a subset of the gri30 mechanism for the methane-co2 | ||
# system |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
# Create the gas object using a subset of the gri30 mechanism for the methane-co2 | |
# system | |
# %% | |
# Create the gas object for the methane-CO2 system using reference state thermo from | |
# the GRI 3.0 mechanism, species critical properties from `critical-properties.yaml`, | |
# and the Peng-Robinson equation of state. |
non-physical values. These non-physical values should trigger an error | ||
in the high-pressure transport model. | ||
""" | ||
gas = ct.Solution(test_data_path / 'methane_co2_noCritProp.yaml') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The test data path is automatically added to the Cantera search path, so you can just write:
gas = ct.Solution(test_data_path / 'methane_co2_noCritProp.yaml') | |
gas = ct.Solution('methane_co2_noCritProp.yaml') |
and remove the fixture from the method.
The case where you need test_data_path
is for accessing files that aren't read by Cantera's AnyMap/YAML functions.
thermal_conductivities.append(gas.thermal_conductivity) | ||
diffusion_coefficients.append(gas.mix_diff_coeffs) | ||
|
||
data = np.empty((len(viscosities), 3+len(diffusion_coefficients[0]))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
data = np.empty((len(viscosities), 3+len(diffusion_coefficients[0]))) | |
data = np.empty((len(viscosities), 3+gas.n_species)) |
} | ||
|
||
|
||
void HighPressureGasTransport::getBinaryDiffCoeffs(const size_t ld, double* const d) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I noticed that according to Codecov, this method doesn't get called, despite the fact that you've added a test that calculates the mixture-averaged diffusion coefficients which rely on the binary diffusion coefficients. I think the problem is that GasTransport::getMixDiffCoeffs
just operates on the m_bdiff
array directly, which here only contains the low-pressure binary diffusion coefficients.
Probably the most direct solution to this would be to override getMixDiffCoeffs
, getMixDiffCoeffsMass
, and getMixDiffCoeffsMole
in HighPressureGasTransport
so they use getBinaryDiffCoeffs
.
This PR is a WIP
Changes proposed in this pull request
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
scons build
&scons test
) and unit tests address code coverage