Critical point calculations, initializations, and hyperdual numbers in Python

[ImagineBaggins]

I've been able to generate some nice diagrams with FeOs (see below), and I want to express again how useful this project is (and your replies to my previous questions)! I've used it quite a bit in my work/dissertation. Though I am experiencing a few difficulties with critical points and near-critical phase equilibrium, and so I have a few questions.

(the above two diagrams calculated with FeOs, using Peng-Robinson since accuracy of critical properties are important; a volume-translated version of SAFT would've been interesting to try too)

As I'm sure you're aware (#236), critical point calculations can sometimes find the "wrong" critical point at given T or P. Specifically, I've had a lot of trouble generating the liquid-liquid critical line (blue) in particular. I've simply been picking a high P and stepping down to ~zero, and the calculation often fails (sometimes a trivial solution, and sometimes an error something like Unexpected T=NaN in validation) or finds the vapor-liquid critical line (red) instead. Using critical_point_binary with both initial_temperature and initial_molefracs from the previous successful point, I can usually extend the convergence range, but it still seems to be extremely sensitive to the initialization and doesn't always extend all the way down to the UCEP.

In creating the Txy diagram above, despite my best efforts, I could not get FeOs to find the LL critical point at P = 73 bar; it would always find one of the other three, no matter the T init (somewhat understandably though, since this is an exceptional case with four isobaric critical points!). I ended up crudely approximating it with a series of bubble-point flashes, making micro-steps in T (𝒪(10^-8) K), until even that would not converge anymore. The VLE envelopes also had lots of holes that needed manual recalculation to patch.

I also tried implementing my own naïve approach for the LL critical line which simply performs repeated bubble-point calculations, and averages L and V compositions until they stop changing, but that usually finds the azeotropic line, if it exists (although, it also starts failing near the VL critical azeotrope).

---

So, with that context, my questions are:

• Is there a better way to initialize CP calculations in FeOs to target LL solutions? Any way to impose bounds on the solution (e.g. maximum temperature)?
• Any method in FeOs for related points such as UCEP/LCEP, critical azeotrope?
• To implement a more rigorous solution to some of the above problems, high-order (4th+) composition and density derivatives are required, which I cannot seem to reverse-engineer from the properties currently exposed in the Python version of FeOs. I was wondering if there is any present way to calculate an arbitrary derivative of the Helmholtz or Gibbs energies in FeOs using hyperdual numbers, i.e., without having to resort to finite difference approximations. I saw the hyperdual section in the docs, but it only appears to show a way to view the dual numbers, not use them.

I totally understand if there is no current solution to any of this, since this is a hard topic with no perfect approach.

[prehner]

Thanks for the questions and for sharing the results you have. It is always good to see the specific cases where calculations don't work. Some things come to mind:

• I usually use PC-SAFT to test algorithms, so maybe some heuristics (e.g., for starting values) might be a bit biased
• More importantly: The critical point algorithm always uses 0.3*max_density to initialize the density. For SAFT models max_density usually relates to a packing fraction of 0.5, for cubics it's related to the value of b. That's a very good initial guess for many systems (including mixtures), but of course only for the critical point, not for UCEP/LCEP. Passing the initial density as an optional argument should be a straightforward addition to the solver.
• A critical azeotrope has the same conditions as a critical point, right? After all, they're all on the same locus. I haven't thought much about how a dedicated algorithm for that would look like to be honest. For the UCEP/LCEP better initial values for the density might already be enough?
• Internally FeOs calculates arbitrary derivatives using recusive (hyper-) dual numbers. However, the generics in Rust don't interact well with the Python interface. Therefore, only specific derivatives are part of the Python library. If you can narrow it down to a set of specific properties/derivatives, that you would need to develop custom algorithms in Python, we can add those to the interface.

[ImagineBaggins]

Thank you for the reply! I should say that the liquid-liquid critical point solution is my main interest here; my 2nd and 3rd questions are more curiosities than actual feature requests as they are rather niche problems. But to address your points:

• Normally I would use PC-SAFT as well, but for PT diagrams specifically, I opted for Peng-Robinson because, as I'm sure you know, the base PC-SAFT model does not exactly reproduce critical properties. Below is the same diagram using PC-SAFT (black stars = PC critical points, red line = calculated VL critical line):

(I don't know how I generated that LL critical line in the original diagram, because I cannot reproduce it now.)

I also just noticed that the figure title is incorrect; this is indeed carbon dioxide and n-butane, not methane.

• Since I haven't gone to the effort of implementing such solution algorithms myself yet, I can't really comment on the best way to solve for a CEP. However, the book "High-Pressure Fluid Phase Equilibria" by Deiters & Kraska writes the following:

where $G_{\left( \mathrm{\#} x \right)} = \left( \frac{\partial^{\mathrm{\#}} G_{\mathrm{m}}}{\partial x_1^{\mathrm{\#}}} \right)_{P,T}$.

This sounds to me like a bit of a headache to implement, so it's probably not worth it for such a niche use case.

• As far as I understand it, yes, a critical azeotrope is a single stable phase (like an ordinary critical point), but also a non-PC endpoint of the azeotropic line. Again, I won't pretend to know the best way to compute such a thing, but Deiters & Kraska have the following to say:

where: $A_{\left( \mathrm{\#} V \right)} = \left( \frac{\partial^{\mathrm{\#}} A_{\mathrm{m}}}{\partial V_{\mathrm{m}}^{\mathrm{\#}}} \right)_{V,T}$,

and: $A_{\left( Vx \right)} = \left( \frac{\partial^2 A_{\mathrm{m}}}{\partial V_{\mathrm{m}} \partial x_1} \right)_T$.

• I understand the limitation of the Rust - Python interface, so my ideal scenario is not possible without a lot of work then. The algorithm I was interested in, again, from Deiters & Kraska, is a variant of the Heidemann - Khalil method (Algorithm 5 in §5.8 of D&K).

The binary case requires $\frac{\partial^2 G}{\partial x_1^2}$ and $\frac{\partial^3 G}{\partial x_1^3}$ (@ constant T & P). Or, in Helmholtz form (which D&K recommend), $\frac{\partial^2 A}{\partial x_1^2}$, $\frac{\partial^3 A}{\partial x_1^3}$, $\frac{\partial^2 A}{\partial V^2}$, $\frac{\partial^3 A}{\partial V^3}$, $\frac{\partial^2 A}{\partial V \partial x_1}$, $\frac{\partial^3 A}{\partial V^2 \partial x_1}$, $\frac{\partial^3 A}{\partial V \partial x_1^2}$ (all @ constant T & V/ρ (the latter only if applicable)). $\frac{\partial^4 G}{\partial x_1^4}$ is also needed for a stability check. The multicomponent form requires similar terms, but now including cross derivatives w.r.t. different mole fractions. Given the large number of distinct derivatives needed, unless there's some way of implementing a Python interface for a generic hyperdual derivative of a generic property, I don't think it's worth the effort of adding all these specific cases.

The further down the phase diagram rabbit hole you go, though, the worse the derivatives get, so technically there are even higher derivatives of interest. For example, the tricritical point (occurring only at very specific combinations of EOS parameters) requires the 1st- through 6th-order mole fraction derivatives. Though that's not something I think most people are interested in calculating.

So, basically, you don't worry about all these situations. I was just curious if there happened to be any procedures in place for them. Besides the fact that the flash/critical point algorithms in FeOs can sometimes have convergence issues near these points, routines that solve for them are more "would be nice to have" features for phase diagram generation only.

[ImagineBaggins]

To come back to the issue of the algorithm for LL critical points, I think this CO₂ - n-Butane system might be a good test case for troubleshooting. I know that this is a type II system, so it should have a continuous VL critical line connecting the two PC critical points, and a separate LL critical line at lower temperatures). Because of the weird shape of this particular VL critical line, the following regions of critical point multiplicity exist (Peng-Robinson values, but qualitatively the same for PC-SAFT):

• $P_{c,\mathrm{VL,max}} < P$: LL only
• $P = P_{c,\mathrm{VL,max}} \approx 80.98 \mathrm{bar}$: LL and VL double critical point (DCP)
• $P_{c,\mathrm{CO_2}} < P < P_{c,\mathrm{VL,max}}$: LL and 2x VL
• $P = P_{c,\mathrm{CO_2}} \approx 73.83 \mathrm{bar}$: LL, pure-component, and 2x VL
• $P_{c,\mathrm{VL,min}} < P < P_{c,\mathrm{CO_2}}$: LL and 3x VL (and vapor pressure curve now starts to have solutions)
• $P = P_{c,\mathrm{VL,min}} \approx 71.67 \mathrm{bar}$: LL, VL ordinary CP, and VL DCP
• $P_{c,n\mathrm{C_4}} < P < P_{c,\mathrm{VL,min}}$: LL and VL
• $P = P_{c,n\mathrm{C_4}} \approx 37.96 \mathrm{bar}$: LL and pure-component
• $P_{\mathrm{UCEP}} < P < P_{c,n\mathrm{C_4}}$: LL only
• $P_{\mathrm{UCEP}} \approx 2.8 \mathrm{bar}$ (very approximate; simply the visual endpoint of the LL curve in my original diagram, assuming that is correct/accurate): LL CEP; subcritical vapor in equilibrium with LL critical point (in type II systems)
• $0 < P < P_{\mathrm{UCEP}}$: no critical points; three-phase equilibrium present

As a quick numerical example with FeOs, I ran a coarse scan down from 100 bar in increments of 1 bar. I am simply calling State.critical_point_binary(eos, Pi*BAR) in hopes of finding the LL solution. Above 80 bar, I get a RuntimeError: `Critical point` did not converge within the maximum number of iterations (as expected). At 80 bar, the first solution occurs, at T = 340.13 K; this is a VL solution. After that, it misses a few points, then finds the remainder of the VL solutions down to the CP of pure n-butane, and nothing thereafter. Trying various combinations of T and x inits., both fixed and based on previous valid solutions, all result in various combinations of failing and finding a VL solution.

I also have an alternate algorithm that simply performs repeated PhaseEquilibrium.bubble_point flashes, averages the "liquid" and "vapor" compositions for the next feed, and repeats, until the differences approach zero. This can find the upper parts of the LL critical line (or what appears to be the LL critical line, anyway), though it starts to fail around 36.5 bar and finds a solution approximately on the CO₂ VP curve (not sure why it's not on it exactly; may be tolerance/precision issue, since I don't think it is finding the three-phase line). I cannot remember what I did to get the "full" LL line in my original diagram.

I don't know if this is helpful for the FeOs algorithm, but D&K write the following:

Only critical curves with almost constant mole fractions might pose a problem, namely if their mole fractions fall between the mole fractions of the $x_1$ search grid, which can sometimes happen with l=l critical curves. It is therefore important to use a narrow $x_1$ spacing in regions of interest. Alternatively, one might do the $V_{\mathrm{m}}$ scan in the outer loop and the $x_1$ scan in the inner loop.

I also don't know the details of the FeOs binary critical point algorithm, so I don't know if the issue is T init., density init, step size restriction, or something else entirely, but I should point out that the particular LL critical line I am trying to solve for is at lower temperature than the VL line, though this need not necessarily be the case for different phase diagram types; in some cases, a VL and LL critical line may even cross in P,T-space, differing only in composition. This is probably only worth pursuing if you want FeOs to be completely robust for phase diagram generation, which is a herculean task. But if you are interested, I'll note some papers here that I've found useful in understanding the situation: 10.1016/j.cherd.2013.06.026, 10.1016/j.fluid.2008.08.017, 10.1016/j.fluid.2008.08.018, 10.1021/acs.iecr.1c04703, 10.1016/j.supflu.2006.03.011.