"Multicomponent" critical point calculation immediately hits NaN

[ImagineBaggins]

Hi there,

First of all, I'm very impressed with FeOs so far! I've found it to be very functional and convenient to use in what I've tried with it so far. However, I am having an issue with calculating a multicomponent mixture critical point.

The mixture technically has four components (N₂, C₁, nC₄, nC₁₄), but the feed composition is actually a binary ([0.934, 0, 0, 0.066]).

I am calling critical_point as follow:

cp = State.critical_point(
    pcsaft,
    moles=z_SI,
    initial_temperature=T_crit_est
)

where pcsaft is a typically configured dft.HelmholtzEnergyFunctional (used successfully with many other systems), z_SI = [0.934, 0, 0, 0.066] mol, and T_crit_est = 163.6088 K (which is admittedly not a great estimate in this case, but it doesn't seem to matter for this issue).

When I attempt to run this, I get a RuntimeError stating the critical point calculation failed to converge. Running it with iteration-level verbosity produces the following:

 iter |    residual    |   temperature   |       density        
----------------------------------------------------------------
    0 |                |  163.60880000 K |   7.34725956 kmol/m³ 
    1 |            NaN | NaN K |   2.44908652  mol/m³
    2 |            NaN | NaN K |   2.44908652  mol/m³
    3 |            NaN | NaN K |   2.44908652  mol/m³
    4 |            NaN | NaN K |   2.44908652  mol/m³
    5 |            NaN | NaN K |   2.44908652  mol/m³
    .
    .
    .

^{(some spacing missing in the temperature column when T = NaN K, btw)}

I admit I may be making some sort of mistake here since this is my first time working with a multicomponent mixture in FeOs, but in any case, State.critical_point_binary seems to find critical composition at the given T/P, while State.critical_point seems to find the Tc/Pc at the given composition (which is what I want in this case).

(On this topic, is there an easy way to calculate the critical composition (à la critical_point_binary) for a multicomponent mixture at a fixed T? I haven't found built-in functionality like this yet.)

Anyway, when running the exact same critical point calculation with pcsaft configured as a true binary (with the zero-mole fraction components removed from the system) and z_SI = [0.934, 0.066] mol, the calculation is successful, with the following iteration history:

 iter |    residual    |   temperature   |       density        
----------------------------------------------------------------
    0 |                |  163.60880000 K |   7.34725956 kmol/m³ 
    1 |   7.69913528e0 |  204.51100000 K |   6.89349730 kmol/m³ 
    2 |   4.97543261e0 |  224.29238325 K |   6.15877134 kmol/m³ 
    3 |   3.83972289e0 |  220.57958195 K |   5.42404539 kmol/m³ 
    4 |   3.61214874e0 |  233.47841300 K |   6.15877134 kmol/m³ 
    5 |   3.54237113e0 |  227.26476762 K |   5.42404539 kmol/m³ 
    .
    .
    .
   31 |  1.58171667e-1 |  431.45279611 K |  17.17966068 kmol/m³ 
   32 |  4.11533652e-2 |  442.85147626 K |  17.26818896 kmol/m³ 
   33 |  1.96759453e-3 |  443.70148728 K |  17.25522211 kmol/m³ 
   34 |  2.66398192e-6 |  443.70278176 K |  17.25519451 kmol/m³ 
   35 | 2.29771757e-12 |  443.70278176 K |  17.25519451 kmol/m³

The initial estimate of Tc does not seem to matter for the failure of the multicomponent call; even copying the converged binary Tc and initializing the multicomponent with that yields the same behavior (the density values are even the same).

Any help/guidance here would be greatly appreciated! :)

---

A second unrelated question while I'm here:

Is there any sort of build of FeOs that is capable of using DGT to compute surface tension for a planar interface? I know that is not a trivial problem to solve reliably, but I've seen at least one paper by Prof. Gross's group (doi: 10.1016/j.fluid.2017.11.032) that compares DFT with DGT for many mixtures including multicomponent ones. However, I can't seem to find any such functionality in FeOs (or any other publicly-available thermo package, besides one, which unfortunately had a lot of other issues). At one point I tried PlanarInterface.from_pdgt, but if I recall correctly, that that did not support multicomponent mixtures.

[ImagineBaggins]

Actually, I found a workaround to my main issue. Adding a small non-zero value to the zeros causes critical_point to solve as normal, as below:

ε = 1e-8
z_SI = [0.934-ε, ε, ε, 0.066-ε]

I suppose something is being divided by zero when a mole fraction is exactly zero.

[prehner]

Hi, thank you for using FeOs and the positive feedback. Often we get to hear what doesn't work, so it's nice to hear that the framework has been useful to you.

Regarding the problem with the critical points: I think you are on point with the issue arising from the 0 in the composition. I thought, I had tried it out for binary mixtures over the entire composition range, but I'll need to recheck that. If it converges with a small value instead of 0 that means it should be somewhat easy to fix (hopefully).

Multicomponent mixtures at a fixed temperature: that is actually not straightforward, you need to specify N-1 intensive degrees of freedom to uniquely determine the critical point of an N-component mixture. For a pure component there are no dofs to specify, for a binary mixture you could chose the composition, the temperature or the pressure, but for a multicomponent mixture the only generic specification you can make is the composition of all the components. Other cases like the temperature plus some ratios of certain components are valid but very specific to a certain application.

Finally, if you are familiar with DGT, you know that for pure components it reduces to a simple integration, whereas for mixtures we need to solve for the composition of the mixture along the interface. Because the kind of equations that is being solved is somewhat different to what we are doing in DFT, and none of the developers or contributors is currently actively working on DGT, we only included the simplest (pure component) case. It could be a nice standalone application built on top of feos, but I'm not aware of any efforts going in that direction.

[ImagineBaggins]

Hello again. I apologize for the lapse in communication.

You are absolutely right regarding multicomponent critical points. I had forgotten to specify the additional constraints I am imposing to deal with the additional DOF. I was interested in implementing the procedure of Shardt & Elliott (10.1016/j.ces.2020.116095), but their publication was very inconsistent with its description of their critical point algorithm, hence why I asked. After some time, however, I did find a functional description of their algorithm, so I have that working on my end now. They allow one component's mole fraction to vary freely, while the remaining fractions always take the same relative proportions as in the initial mixture, and iterate flashes until a maximum pressure is reached, which they use as a pseudo-critical point.

Regarding DGT, I am definitely aware of the difficulties in solving the general multicomponent interface problem. That is actually exactly why I came to ask here; because I am trying to avoid having to implement it myself in the interest of time. I totally understand if the FeOs team has no plans to implement such a method, however the main reason I asked here is because some of the publications from Prof. Gross's group - namely Mairhofer & Gross, 2018 (10.1016/j.fluid.2017.11.032) and your own dissertation - seem to show a working implementation of (P)DGT + PCP-SAFT for multicomponent systems. Specifically, Fig. 5.6 from your dissertation depicts exactly what I'd like to be able to do.

If you know of any publicly-available version of this code or are willing to share something, I'd love to talk about this more. Perhaps via email would be better than GitHub.

[g-bauer]

Regarding DGT/surface tensions: Are you interested in DGT in particular or in surface tensions for multi-component mixtures in general? If you "just" need surface tensions for mixtures - that's possible using DFT. There are several examples here.

[ImagineBaggins]

Thank you for the link. I am not necessarily limiting myself to just DGT, and I have run a lot of systems with FeOs's DFT now. It definitely performs very well in a lot of cases, but also not good enough in just as many other cases. Occasionally I run into calculation issues (trivial flashes, CP failures, DFT failures), but more often the model does "work", just with a poor match to the experimental data. So I believe that may be more a shortcoming of DFT itself rather than FeOs.

Whether or not that would change by using DGT is doubtful (particularly since the complication of untangling the choice of thermodynamic model, pure params., and BIPs from the performance of the surface tension model itself would only be exacerbated by introducing several more params. with DGT), but having more options is always good. I have investigated other models besides these two, but they seem to be among the best for semi-/fully-predictive models applicable at high pressures, esp. when the system is above the critical point of individual pure components.