Fix HEOS fugacity for mixtures in the two-phase region

[friederikeboehm]

Description of the Change

I noticed that when calculating the fugacity for mixtures in Peng-Robinson EoS the results did not quite made sense. Since the fugacity of a component in the mixture must be the same in liquid and gas, I checked the saturation states which both had the same fugacity but different from the overall state. As I am not quite sure the two-phase region is currently appropriately handled in this case, I suggest this PR which returns the sat-state fugacity. It now returns the weighted sum of both sat-states, which is of course equal to the single sat-states (should I just decide on liquid or gas?).

Benefits

Get better fugacity values for two-phase region

Possible Drawbacks

/

Verification Process

Tested before and after change fugacity values for AbstractStates: EoS "HEOS" and "PR":

• components: "Methane&Ethane&Propane"
• fractions: (0.25, 0.25, 0.5)
• PQ inputs: p = 500000, q = 0.5
• get molar fractions for sat states -> update State appropriately
• compare differences in fugacities for overall, gas and liquid state

after the changes, overall and gas/liquid returned the same values

[ibell]

I don't think this is a good idea; mixture fugacity doesn't linearly scale between saturated liquid and vapor values and I thus I think we should close this PR.

[friederikeboehm]

Umm, correct me if I'm wrong, but there is nothing to scale as the values have to be the same.
But I am fine with abandoning the current idea, but then there should be an error message instead of a value that makes no sense at all. I also think in this case it would be nice to have a function similar to AbstractState_keyed_output_satState in CoolPropLib to be able to get the value without an additional flash calculation. Happy to implement it :)
Should I close this PR and make a new one with the proposed changes?

[ibell]

Yes, you are right that the fugacity coefficients (and likewise the chemical potentials and fugacities) should be the same in both phases. They will not be exactly the same due to the tolerance in the phase equlibrium solver. If two phase, I think I would leave your solution, certainly that would be better than yielding the homogeneous result inside the two phase dome. Overall I think the proposed solution of adding the routines to expose the values for the coexisting phases in c++, python, and julia is the best one.