-
-
Notifications
You must be signed in to change notification settings - Fork 361
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
Add isothermal_compressibility and thermal_expansion_coefficient for Redlich-Kwong and Peng-Robinson EOS #1421
Conversation
Codecov Report
@@ Coverage Diff @@
## main #1421 +/- ##
==========================================
- Coverage 70.78% 69.51% -1.27%
==========================================
Files 373 377 +4
Lines 53872 57993 +4121
Branches 17579 20693 +3114
==========================================
+ Hits 38131 40313 +2182
- Misses 13339 14730 +1391
- Partials 2402 2950 +548
... and 301 files with indirect coverage changes 📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
2baf401
to
c494173
Compare
For reference, the following finite difference script was used to calculate the values in the test and set reasonable tolerances: import cantera as ct
import numpy as np
def isothermal_compressibility(state, dP):
T, P = state.TP
state.TP = T, P - dP
rho1 = state.density_mass
state.TP = T, P + dP
rho2 = state.density_mass
state.TP = T, P
return - state.density_mass * (1 / rho2 - 1 / rho1) / (2 * dP)
def thermal_expansion_coeff(state, dT):
T, P = state.TP
state.TP = T - dT, P
rho1 = state.density_mass
state.TP = T + dT, P
rho2 = state.density_mass
state.TP = T, P
return state.density_mass * (1 / rho2 - 1 / rho1) / (2 * dT)
gas = ct.Solution("thermo-models.yaml", "CO2-PR")
gas.X = "CO2: 1"
betaT_ct = []
betaT_fd = []
alphaV_ct = []
alphaV_fd = []
P = 5e5
for i in range(10):
T = 300 + 25 * i
gas.TP = T, P
betaT_ct.append(gas.isothermal_compressibility)
alphaV_ct.append(gas.thermal_expansion_coeff)
betaT_fd.append(isothermal_compressibility(gas, P * 1e-5))
alphaV_fd.append(thermal_expansion_coeff(gas, T * 1e-5))
betaT_error = np.array(betaT_fd) - np.array(betaT_ct)
alphaV_error = np.array(alphaV_fd) - np.array(alphaV_ct)
print(f"Isothermal compressibility\nCantera: {betaT_ct}")
print(f"Finite difference error: ≤ {max(abs(betaT_error)):.5e}\n")
print(f"Thermal expansion coefficient\nCantera: {alphaV_ct}")
print(f"Finite difference error: ≤ {max(abs(alphaV_error)):.5e}") which gives the following output for
Note: For the It seems that there might be a slight difference in critical parameters used for carbon dioxide, as the latest version of CoolProp gives slightly different values for density than were used in the original |
@ischoegl Sorry to tag you again, but the finite difference tests are now implemented as thermodynamic consistency tests as you suggested using the framework by @speth, which is great because now I don't have to do any extra work for #1423 😄. Everything passes |
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.
@corykinney ... thanks - while I'm not familiar with the theory (and don't have time to catch up at the moment), I looked over the code; the only thing I found was a LaTeX docstring typo. @decaluwe - the actual changes look like they are well aligned with recent work by your group - could you review and approve as appropriate?
Regarding your question about squashing commits: as you introduce tests that are subsequently removed it makes sense to squash corresponding commits.
eaeb52a
to
4c155a6
Compare
Sorry for the delay - taking a look at the derivation, just to get another set of eyes on it... Stay tuned :) |
@decaluwe I appreciate you taking the time to review it! The derivations in the original enhancement proposal were a bit outdated, as I didn't originally realize that |
@speth For the consistency tests with finite differences implemented here, I followed the tolerance pattern used for checking the final result of many of the other properties, in this case it looked like:
I would think |
|
@speth Thank you for the explanation. For some reason my impression was the opposite, which is that I tried to select some reasonable numbers for the two quantities with units
@speth Do you think these tolerances are reasonable, or are they perhaps too restrictive? |
Actually I realized the values were swapped, so now the tests pass and the only question is whether the |
@decaluwe Have you had a chance to look at the final version (also linked above) of the derivations and the implementation in the last couple weeks? |
@ischoegl Any chance this might get reviewed soon? |
@corykinney … I was under the impression that @decaluwe would review the theory. |
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.
Thank you for this PR, @corykinney. I think both the implementation of these properties for the P-R and R-K models and the additional consistency tests for these properties across all phase models are useful contributions.
I had some comments on the derivation of these formulas that I shared in Cantera/enhancements#122 that I think can be used to simplify the calculations here. Other than that, I have just a few minor comments shown below.
//! Returns the isothermal compressibility. Units: 1/Pa. | ||
virtual double isothermalCompressibility() const; | ||
|
||
//! Return the volumetric thermal expansion coefficient. Units: 1/K. | ||
virtual double thermalExpansionCoeff() const; |
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 this would be a good place to document the formulas for these properties. Likewise in RedlichKwongMFTP.h
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.
Added the formulas in a9b5945
Unfortunately I'm not certain how to build the documentation to verify, but it should be formatted correctly
test/thermo/consistency.cpp
Outdated
double P1 = phase->pressure(); | ||
double mv1 = phase->molarVolume(); | ||
|
||
double P2 = P1 * (1 + 1e-7); |
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 use the same relative perturbation of 1e-6 that is used on other finite difference tests, I think you can set much tighter absolute tolerances for these tests. At least locally, with this larger perturbation, I was able to pass the test suite with atol_c = 1e-14
and atol_e = 1e-18
.
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.
Adjust tolerances in a9b5945 as mentioned
test/thermo/consistency.cpp
Outdated
double mv_mid = 0.5 * (mv1 + mv2); | ||
double alphaV_fd = 1 / mv_mid * (mv2 - mv1) / (T2 - T1); | ||
|
||
EXPECT_NEAR(alphaV_fd, alphaV_mid, max({rtol_fd * alphaV_mid, rtol_fd * alphaV_fd, atol_e})); |
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 wrap long lines (over 88 characters).
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.
Wrapped both lines in a9b5945
@speth Currently working on testing the simplified equations mentiond in your other comment. Regarding the applicability of the inverse function theorem, the Jacobian just needs to be invertible at whichever point we're evaluating the derivatives at in order to use these simplifications, correct? |
@speth Switching to the derivative identity approach allowed for leveraging existing pressure derivative functions, so now the implementation is extremely simple. Let me know if everything looks good, and I can squash the changes. |
Add finite difference checks for `isothermalCompressibility` and `thermalExpansionCoeff` for thermo models.
1b28545
to
37409c0
Compare
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 update to this PR, @corykinney. The revised implementation looks good to me -- I wasn't even expecting it to be that simple!
Besides the minor doc update, please add yourself to the AUTHORS
file, and then I think this will be good to go.
//! Returns the isothermal compressibility (\f$\beta_T\f$). Units: 1/Pa. | ||
/*! | ||
* \f[ | ||
* \beta_T = -\frac{1}{v} \left(\frac{\partial v}{\partial p}\right)_T | ||
* \f] | ||
*/ | ||
virtual double isothermalCompressibility() const; |
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's no need to repeat documentation for overrides unless there's some specific differences in the documentation, which isn't the case here.
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 didn't realize the equations were already documented in ThermoPhase
! The unnecessary comments are now removed.
Add `isothermalCompressibility` and `thermalExpansionCoeff` definitions to the `RedlichKwongMFTP` and `PengRobinson` models
37409c0
to
931de24
Compare
@speth Let me know if there are any other changes :) |
Changes proposed in this pull request
Add
isothermalCompressibility
andthermalExpansionCoeff
definitions to theRedlichKwongMFTP
andPengRobinson
models.Add corresponding thermo consistency tests
betaT_eq_minus_dmv_dP_const_T_div_mv
andalphaV_eq_dmv_dT_const_P_div_mv
based on finite differences.Closes Cantera/enhancements#122
Checklist
scons build
&scons test
) and unit tests address code coverage