Add isothermal_compressibility and thermal_expansion_coefficient for Redlich-Kwong and Peng-Robinson EOS

[corykinney]

Abstract

Implementing the required functions to calculate the isothermal compressibility and thermal expansion coefficient for the Peng-Robinson EOS added in Cantera/cantera#1047.

Motivation

I'm currently developing a Python package that solvers for shock conditions (in a shock tube) using a real gas EOS. Currently, the implementation of Peng-Robinson does not define the isothermal_compressibility and thermal_expansion_coefficient properties that the algorithm requires.

Description

I found this paper that includes the derivation of the derivative of molar mass with respect to temperature at constant pressure, which can be used for the thermal expansion coefficient. I couldn't find a similar resource relevant to the isothermal compressibility factor. I have attached my attempt at a derivation of the isothermal compressibility factor in case there are any errors in my understanding of the thermodynamic principles or the derivatives.

I have not yet started on any code, as I have limited C++ experience and expect to have some trouble implementing it. I will document my progress here, and if anyone is willing to help, that would be great!

[corykinney]

This process is intended to be used for relatively major features, e.g. adding a new thermodynamic model or other rewriting one of the language interfaces, that would require significant development time. Small feature additions like introducing a single function or adding a new example are more simply handled by creating an issue or pull request in the Cantera/cantera repository directly.

Would it be better to move this issue to the main repository as it is smaller?

[bryanwweber]

@corykinney Nope, this is a great place for this issue. I'd be happy to have a Zoom/Webex/etc. to help you get started on the contribution, if it'd be helpful. Feel free to reach out to the email in my profile!

[bryanwweber]

Also ping @gkogekar @decaluwe for their thoughts on the thermo.

[decaluwe]

thanks, @bryanwweber - I had noticed this, yesterday. This would definitely be a great extension of that model. I'll try and remember to give it a look once I've finished grading for the semester and recovered that portion of my brain. @corykinney please feel free to bug me if I forget, if you want feedback (but I also trust you to proceed, without any input from me).

[corykinney]

@bryanwweber I will definitely be reaching out to you and asking some questions, I appreciate the offer!

@decaluwe I'll go ahead and start trying to implement it next week. I'll likely proceed fairly slowly as I figure out the C++ side of things, but I would definitely appreciate if you had time to double check the equations sometime along the way!

[corykinney]

Starting to work on it at corykinney:Peng-Robinson

[ischoegl]

@corykinney … are you planning to post a PR for this?

[corykinney]

@ischoegl ... Yes, I plan on posting a PR eventually. It's just taken a back seat to some other things, so it might be a month or two before I get around to it.

[corykinney]

@corykinney … are you planning to post a PR for this?

Unfortunately I won't be able to continue working on this, but I'll leave this open in case anyone wants to finish it from the derivations I had included in the initial post.

[corykinney]

@ischoegl I've been working on real gas stuff more recently and had time to finish the implementations over on corykinney:real-gas-derivatives. I just need to add test cases to make sure that I didn't make any mistakes. For Peng Robinson I will probably compare to CoolProp. Any ideas for how to test the Redlich Kwong since CoolProp doesn't have that implemented?

[ischoegl]

@corykinney … sounds great! However, I myself don’t have much experience with non-ideal gas models. @decaluwe should be able to help though as his group has been driving some of these efforts lately.

[decaluwe]

Hi @corykinney. Here are my recommendations:

1. For most tests, my preference is to hard-code the EoS results into the test. I.e. input the a and b and any other required coefficients directly into the test, and calculate the EoS prediction manually to compare against Cantera's result. I would probably prefer this even over an above comparing to external software. External software calcs obviously provide a good "sanity check" on our results, but since we can't guarantee that other software packages are free of error, there is some very tiny risk there. Not a huge preference, and consulting external software can certainly be easier but all things equal...
2. A lot of times you can use fundamental identities to check internal consistency. I.e. if the compressibility is equal to a derivative of other properties, you can use finite differencing to estimate that derivative and check your result against it, with some small tolerance for error.
3. If neither of those is a realistic or presents an undue burden for some reason, the bare minimum would be a regression test - I.e. print out some current values returned by the software, enter them manually into the test and check that the results never deviate from these values. This is definitely sub-optimal (how do we know that these values are right, in the first place?), but at least checks against unintended consequences from other changes in the code base.

Thanks!
Steven

[corykinney]

1. For most tests, my preference is to hard-code the EoS results into the test. I.e. input the a and b and any other required coefficients directly into the test, and calculate the EoS prediction manually to compare against Cantera's result. I would probably prefer this even over an above comparing to external software. External software calcs obviously provide a good "sanity check" on our results, but since we can't guarantee that other software packages are free of error, there is some very tiny risk there. Not a huge preference, and consulting external software can certainly be easier but all things equal...

@decaluwe If I calculate the result manually and hard code the value, wouldn't that be open to an error in my derivation of the equation? I suppose if I verify the manual calculation matches the external software it would be a good option.

2. A lot of times you can use fundamental identities to check internal consistency. I.e. if the compressibility is equal to a derivative of other properties, you can use finite differencing to estimate that derivative and check your result against it, with some small tolerance for error.

This seems like it would be a bit more straightforward, as I could implement a central finite difference for a small change in either temperature or pressure and the corresponding change in molar volume. I suppose I would just need to select the right tolerance for error.

[decaluwe]

If I calculate the result manually and hard code the value, wouldn't that be open to an error in my derivation of the equation? I suppose if I verify the manual calculation matches the external software it would be a good option.

At some point somewhere, one has to assume that somebody has done something correctly. 😁 I don't see any way around this, in a test. But yes, a one-time verification of your derivation against external software would be a good way to be double-sure. It would probably also be a good standard procedure to always get at least two sets of eyes on any derivation that isn't already an established part of the literature. I'll try and make some time for this, between now and when this PR is ready.

could implement a central finite difference for a small change in either temperature or pressure and the corresponding change in molar volume.

My only caution here is: what would differentiate what P-R returns for any given identity, vs. ideal gas, vs. R-K, etc? I guess it would check against gross/large errors in the derivation or implementation, but still wouldn't give a guarantee that the theory was done correctly. Again, a workable solution, and yes more straightforward, but it would be good if tests like these were supplemented by others. Also, thermodynamic identities, such as g = h - Ts can be useful here, in that they are exact and can catch errors in one property, so long as those errors aren't also propagated to the calculation of the other properties in the identity.

[ischoegl]

A lot of times you can use fundamental identities to check internal consistency. I.e. if the compressibility is equal to a derivative of other properties, you can use finite differencing to estimate that derivative and check your result against it, with some small tolerance for error.

As a comment on thermodynamic consistency, @speth introduced comprehensive checks in Cantera/cantera#1299, which implements #114. While presumably not available for all properties, the framework is in place.

[corykinney]

At some point somewhere, one has to assume that somebody has done something correctly.

@decaluwe That's a fair point, perhaps I was going a little overboard 😁

My only caution here is: what would differentiate what P-R returns for any given identity, vs. ideal gas, vs. R-K, etc? I guess it would check against gross/large errors in the derivation or implementation, but still wouldn't give a guarantee that the theory was done correctly.

I implemented a quick test for the work on the SRK EOS #159 to check the implementations for the same functions, and it seems like a central finite difference can give at least a decent check on accuracy. Here is some sample output for TPX of 1000, 300e5, "CO2: 1" with the prototype SRK implementation:

SRK [Cantera]
Density = 145.14896320725558
Isothermal compressibility = 3.036143602111974e-08
Finite difference approx = 3.036143608806698e-08
Thermal expansion coeff = 0.0009779386240779259
Finite difference approx = 0.0009779386255063225

Ideal [Cantera]
Density = 158.7919821922599
Isothermal compressibility = 3.333333333333334e-08
Finite difference approx = 3.3333333275688315e-08
Thermal expansion coeff = 0.001
Finite difference approx = 0.0009999999999998508

the check is pretty simple

def isothermal_compressibiliy_central_diff(gas, dP):
    T, P = gas.TP
    gas.TP = T, P - dP
    rho1 = gas.density_mass
    gas.TP = T, P + dP
    rho2 = gas.density_mass
    gas.TP = T, P
    return - gas.density_mass * (1 / rho2 - 1 / rho1) / (2 * dP)

and something similar for thermal_expansion_coeff, where dP and dT were $1\,Pa$ and $1\,K$, respectively. The difference might not be so apparent at lower pressures though. Of course having other checks would be good as well.

As a comment on thermodynamic consistency, @speth introduced comprehensive checks in Cantera/cantera#1299, which implements #114. While presumably not available for all properties, the framework is in place.

@ischoegl I see that isothermalCompressibility and thermalExpansionCoeff were on the initial list of checks in the enhancement proposal, were those two specifically implemented? If so, how do I run this test suite and is there a need for additional tests to be written or is that suite intended to be sufficient?

[speth]

Thanks for sharing your derivation of the Peng-Robinson isothermal compressibility. I think there is actually a simpler approach to calculating all of these quantities, by applying a couple of partial derivative identities.

Peng-Robinson

Writing the EoS as:
$$p = \frac{RT}{v-b} - \frac{a\alpha(T)}{v^2 + 2 b v - b^2}$$

Isothermal compressibility, $1/v\cdot(\partial v/\partial p)_T$

We can compute the isothermal compressibility by noting that
$$\left(\frac{\partial v}{\partial p}\right)_T = 1 \bigg/ \left(\dfrac{\partial p}{\partial v}\right)_T$$
where the derivative on the right hand side can be computed directly from the EoS:
$$\left(\dfrac{\partial p}{\partial v}\right)_T = -\frac{RT}{(v-b)^2} + \frac{2a(v+b)}{(v^2+2bv-b^2)^2}$$

I tested this with your PR, and the results are identical using this more compact formula.

Thermal expansion coefficient, $1/v\cdot (\partial v/\partial T)_p$

Using the triple product rule:
$$\left(\frac{\partial v}{\partial T}\right)_p \left(\frac{\partial p}{\partial v}\right)_T \left(\frac{\partial T}{\partial p}\right)_v = -1$$
The first term is what we're looking for. The second is what we already calculated to get the isothermal compressibility, and the inverse of the third can be computed directly from the EoS:
$$\left(\frac{\partial p}{\partial T}\right)_v = \frac{R}{v-b} - \frac{d(a\alpha)}{dT} \frac{1}{v^2 + 2bv-b^2}$$

This again gives identical results to what's implemented in your PR, but I think this formula is a bit simpler.

Redlich-Kwong

Isothermal compressibility

Again, differentiating the EoS directly gives:
$$\left(\dfrac{\partial p}{\partial v}\right)_T = - \frac{RT}{(v-b)^2} + \frac{a}{\sqrt{T}}\frac{2v+b}{(v^2+vb)^2}$$

Thermal expansion coefficient

The same process applies, and we just need the derivative:
$$\left(\frac{\partial p}{\partial T}\right)_v = \frac{R}{v-b} - \frac{1}{v(v+b)}\left(\frac{da}{dT} \frac{1}{\sqrt{T}} - \frac{a}{2 T^{3/2}}\right)$$

I compared results using both of these to the formulas with the implementation in your PR (is there an accompanying derivation somewhere?) and found that they give identical results.

[corykinney]

@speth Using the partial derivative identities is a great idea! I can update the code with the simplified equations and implement the equations into the comments.