VT flash: water + a tiny amount of hydrogen?

[B-LIE]

I am puzzled by the following... assume a tank with pure water and the following data:

using Clapeyron
#
# content of a cathode separation tank
n_H2O_c = 0.648e4
V_c = 0.35

# models
mod_pr = PR(["water"])
mod_grg = GERG2008(["water"])
mod_phsft = pharmaPCSAFT(["water"])
mod_mf = MultiFluid(["water"])

# temperature
T = 70 + 273.15

If I do VT flash, I get:

julia> vt_flash(mod_pr,V_c,T,[n_H2O_c])
Flash result at T = 343.15, p = 28583.8 with 2 phases:
 (x = [1.0], β = 6477.92, v = 2.2099e-5)
 (x = [1.0], β = 2.07864, v = 0.0995097)

julia> vt_flash(mod_grg,V_c,T,[n_H2O_c])
Flash result at T = 343.15, p = 31207.5 with 2 phases:
 (x = [1.0], β = 6477.46, v = 1.84249e-5)
 (x = [1.0], β = 2.54407, v = 0.0906631)

and similar for the two other models. Fine, no problem.

Then I try to add a relatively small amount of hydrogen:

n_H2_c = 251
mod_pr = PR(["water", "hydrogen"])
mod_grg = GERG2008(["water", "hydrogen"])
mod_phsft = pharmaPCSAFT(["water", "hydrogen"])
mod_mf = MultiFluid(["water", "hydrogen"])

With this addition, there appears to be only a single phase only for most models, unless I reduce the H2 amount. With the PR model:

julia> vt_flash(mod_pr,V_c,T,[n_H2O_c, n_H2_c])
Flash result at T = 343.15, p = 2.7025e6 with 1 phase:
 (x = [0.96271, 0.0372902], β = 6731.0, v = 2.22071e-5)

julia> vt_flash(mod_pr,V_c,T,[n_H2O_c, 0.1n_H2_c])
Flash result at T = 343.15, p = 2.53105e6 with 1 phase:
 (x = [0.996141, 0.00385851], β = 6505.1, v = 2.20996e-5)

julia> vt_flash(mod_pr,V_c,T,[n_H2O_c, 0.01n_H2_c])
Flash result at T = 343.15, p = -2.0756e8 with 2 phases:
 (x = [0.999997, 3.12906e-6], β = 8.49378e-10, v = 2.33729e-5)
 (x = [0.999746, 0.000253753], β = 6482.51, v = 3.87907e-5)

julia> vt_flash(mod_pr,V_c,T,[n_H2O_c, 0.001n_H2_c])
Flash result at T = 343.15, p = 2.53067e6 with 2 phases:
 (x = [0.999961, 3.87331e-5], β = 6480.25, v = 2.20887e-5)
 (x = [0.908839, 0.0911615], β = 2.87781e-13, v = 2.24414e-5)

I'm puzzled. Adding relatively little H2 seems to (a) change the system from 2 phase water to 1 phase (gas) H2O+H2 (and the pressure goes up from ca. 0.29 bar for pure water, to 27 bar for the mixture) . If I (b) add only 10% of the H2 (10% refers to n_H2_c), still 1 phase but lower pressure (which makes sense). (c) Adding only 1% H2 gives a negative pressure and two phases -- doesn't really make sense?? (d) Adding still less, 0.1% H", gives a positive pressure and two phases.

Essentially, I'm puzzled by how little H2 I can add and still have 2 phases, + that I get a negative pressure in one case.

---

If I use the other models...

julia> vt_flash(mod_grg,V_c,T,[n_H2O_c,n_H2_c])
Flash result at T = 343.15, p = NaN with 1 phase:
 (x = [0.96271, 0.0372902], β = 6731.0, v = NaN)

with similar result for all 4 cases of added H2. The MultiFluid model behaves similar to the GERG2008.

julia> vt_flash(mod_phsft,V_c,T,[n_H2O_c, n_H2_c])
Flash result at T = 343.15, p = NaN with 1 phase:
 (x = [0.96271, 0.0372902], β = 6731.0, v = NaN)

julia> vt_flash(mod_phsft,V_c,T,[n_H2O_c, 0.1n_H2_c])
Flash result at T = 343.15, p = NaN with 1 phase:
 (x = [0.996141, 0.00385851], β = 6505.1, v = NaN)

julia> vt_flash(mod_phsft,V_c,T,[n_H2O_c, 0.01n_H2_c])
Flash result at T = 343.15, p = -2.28029e8 with 2 phases:
 (x = [0.999998, 1.56926e-6], β = 1.49156e-9, v = 1.93223e-5)
 (x = [0.99983, 0.000169777], β = 6482.51, v = 3.1487e-5)

julia> vt_flash(mod_phsft,V_c,T,[n_H2O_c, 0.001n_H2_c])
Flash result at T = 343.15, p = 2.32886e6 with 2 phases:
 (x = [0.999961, 3.87331e-5], β = 6480.25, v = 1.84188e-5)
 (x = [0.813639, 0.186361], β = 2.87781e-13, v = 1.96121e-5)

So the PharmaSAFT model gives sort of similar results as the PR (i.e., in the sense of 1 phase vs. 2 phase, but does not compute the pressure nor the molar volume in the one phase case -- which PR does).

---

Questions:

1. Does my results make sense? At the outset, I tried to compute how much H2O and H2 I could have in the tank and still have ca. 30 bar pressure and ca. 30% liquid in the tank. Perhaps I should do PT flash to find that?
2. What does the negative pressure mean (PR, PharmaSAFT, some cases of H2 content).
3. What does it mean that some models can not compute pressure or molar (gas) volume?

[B-LIE]

I did some more testing, and tried to specify the pressure, temperature, and amount using the tp_flash algorithm. Then I used the volume function to compute the volumes of the phases. Next, I used the computed total volume + temperature + amounts, and hoped I should get back the pressure and number of phases from tp_flash. But that didn't happen...

Here is what I did (continueing from the Julia code above):

# Defining a utility function tpv_flash
function tpv_flash(mod, p, T, n=SA[1.0])
    # Input arguments
    # - mod: model of system, including substances
    # - p: pressure in Pa 
    # - T: temperature in K 
    # - n: amount vector in mol, consistent with substances in mod
    # Return variables
    # - 1: mole fraction matrix, first row is liquid, second is gas
    # - 2: amount matrix, first row is liquid, second row is gas, in mol
    # - 3: vector of liquid volume and gas volume, in m3
    # - 4: vector of liquid molar volume and gas molar volume, m3/mol
    #
    tp_sol = tp_flash(mod, p, T, n)
    if size(tp_sol[1],1) == 1
        println("1-phase")
        V = volume(mod, p, T, tp_sol[2][1,:])
        return tp_sol[1], tp_sol[2], V, V/sum(tp_sol[2][1,:])
    else
        VL = volume(mod, p, T, tp_sol[2][1,:]; phase=:liquid)
        VG = volume(mod, p, T, tp_sol[2][2,:]; phase=:gas)
        return tp_sol[1], tp_sol[2], [VL, VG], [VL/sum(tp_sol[2][1,:]), VG/sum(tp_sol[2][2,:])]
    end
end

Then running some code:

julia> p = 3e6 # 30 bar
julia> tp_sol = tpv_flash(mod_pr,p_c,T,[n_H2O_c,n_H2_c])   # volume vector is 3rd return value
julia> vt_flash(mod_pr, sum(tp_sol[3]),T,[n_H2O_c,n_H2_c])
Flash result at T = 343.15, p = 2.7025e6 with 1 phase:
 (x = [0.96271, 0.0372902], β = 6731.0, v = 2.22071e-5)

I'm confused. The tp_flash routine gives a two phase solution at 30 bar and 70 C with the given amount. But when I back-calculate with the vt_flash routine, I get a 1 phase solution and a lower pressure.

---

I tried to use tp_flash on the GERG2008 and the MultiFluid models, but I get the same error message:

julia> tp_flash(mod_mf,p_c,T,[n_H2O_c,n_H2_c])
ArgumentError: Collection has multiple elements, must contain exactly 1 element
...

If I use the pharmaSAFT model, I get a result, it is strange -- finds the mole fractions, but not the mole numbers (which should be trivial), and not the Gibbs energy:

julia> tp_flash(mod_phsft,p_c,T,[n_H2O_c,n_H2_c])
([0.9810006741922844 0.01899932580771557; 1.955278169111263e-5 0.9999804472183089], [NaN NaN; NaN NaN], NaN)

[longemen3000]

Hello,

The anomalous results are definitely a bug. Negative pressures just mean that pressure(model,v,T,z) < 0, but that does mean that the volume point is not a valid one.

The main source of bugs in this case is the initial point generator. In the case of VT-flash, the initial point is calculated via estimating the pressure from temperature (via `Pproperty), and then performing K-value separation at specified p and T.

Pproperty works by calculating bubble and dew poi ts, and finding where is the dpecified property (in this case, volume). In the case of hydrogen + water, the bubble pressure fails (as there is no amount of pressure that would create a main phase with that amount of liquid hydrogen). Im working on improving the robustness there.

If you know your system a little bit more, you could pass an initial point. In this case, passing p0 = my_p) could help.

[B-LIE]

Thanks for suggestion on providing initial guess of pressure in the vt_flash algorithm. Here are some results:

julia> # without initial pressure guess
julia> vt_flash(mod_pr,V_c,T,[n_H2O_c,n_H2_c])
Flash result at T = 343.15, p = 2.7025e6 with 1 phase:
 (x = [0.96271, 0.0372902], β = 6731.0, v = 2.22071e-5)

julia> # Observe: for PR, 1 phase with pressure 27 bar
julia> vt_flash(mod_pr,V_c,T,[n_H2O_c,n_H2_c]; p0=2e6)
Flash result at T = 343.15, p = 3.5223e6 with 2 phases:
 (x = [0.999946, 5.37728e-5], β = 6478.04, v = 2.20846e-5)
 (x = [0.0091312, 0.990869], β = 252.962, v = 0.00081805)

Kind of weird... if I make an initial guess that is lower than the vt_flash algorithm found when incorrectly suggesting a single phase, the vt_flash algorithm now correctly finds a two-phase solution with a pressure of 35 bar.

Providing initial guess for the other (and better?) models also works, and the algorithm seems to be relatively robust wrt. initial guess. These other algorithms does suggest a lower pressure, though: ca. 32 bar.

---

Still, tp_flash fails for all methods except the PR model. More precisely, the pharmaSAFT method finds a mole fraction matrix, but fails for the rest:

julia> tp_sol = tp_flash(mod_phsft,p_c,T,[n_H2O_c,n_H2_c])
julia> tp_sol[1]
2×2 Matrix{Float64}:
 0.981001    0.0189993
 1.95528e-5  0.99998
julia> tp_sol[2]
2×2 Matrix{Float64}:
 NaN  NaN
 NaN  NaN
julia> tp_sol[3]
NaN

---

I should add: when I use an initial guess for p (specifying p0), the result is consistent for the case of the PR algorithm:

julia> p_c=3e6
julia> tp_sol = tpv_flash(mod_pr,p_c,T,[n_H2O_c,n_H2_c])
julia> # mole fractions
julia> tp_sol[1]
2×2 Matrix{Float64}:
 0.999954   4.58651e-5
 0.0105417  0.989458
julia> # amounts in the two phases
julia>sum(tp_sol[2][1,:]), sum(tp_sol[2][2,:])
(6477.626098760676, 253.37390123932533)
julia> # molar volumes
julia> tp_sol[4]
2-element Vector{Float64}:
 2.208675409551813e-5
 0.0009588718235683983

julia> V_c_tpv = sum(tp_sol[3])
julia> vt_flash(mod_pr,V_c_tpv,T,[n_H2O_c,n_H2_c]; p0=p_c)
Flash result at T = 343.15, p = 3.0e6 with 2 phases:
 (x = [0.999954, 4.58651e-5], β = 6477.63, v = 2.20868e-5)
 (x = [0.0105417, 0.989458], β = 253.374, v = 0.000958872)

[B-LIE]

Some results I get when I multiply $n_\mathrm{H_2}$ by a factor $\chi_\mathrm{H_2}$ varying from 0.1 to 6. Note logarithmic scale along the abscissa.

H2O liquid fraction:

Remarks:

• MultiFluid and Gerg2008 give virtually the same result for all properties I compute, but fail at high pressures.
• pharmaSAFT is closest to MultiFluid and Gerg2008 in many cases
• the quick routine is my own code, based on ideal gas law + Raoult's law for water and Henry's law for hydrogen. Surprisingly, the quick results appear to be indistinguishable from those of MultiFluid and Gerg2008, except for the liquid phase mole fractions and liquid total amount.

Next, H2O gas fraction:

Next, liquid molar volume:

Next, gas molar volume:

Next, total pressure (note: log-log scale).

Conclusions: (?)

• The Gerg2008 and MultiFluid results are perhaps most correct?? Unfortunately, they fail for high pressures in my tests.
• The quick routine seems to give decent performance, except for the liquid phase.

[longemen3000]

I'm gonna check what is happening with those multiparameter EoS. probably the initial point for those can be improved (as there are ancillaries for the saturation properties, and one can search those directly instead of performing saturation calculations).

Also, i forgot to mention it, but you can also provide, in addition to p0, a vector of K-values for the initial point:

help?> GeneralizedXYFlash
search: GeneralizedXYFlash GeneralizedSchur GeneralizedEigen GeneralizedSVD

  GeneralizedXYFlash{T}(;kwargs...)

  Method to solve non-reactive multicomponent, two-phase flash problem, using a generalized formulation.

  Only two phases are supported. if K0 is nothing, it will be calculated via fugacity coefficients at p,T conditions.

  Keyword Arguments:
  ––––––––––––––––––

    •  equilibrium (optional) = equilibrium type ":vle" for liquid vapor equilibria, ":lle" for liquid liquid equilibria, :unknown if not specified

    •  p0 (optional), initial guess pressure, ignored if pressure is one of the flash specifications.

    •  T0 (optional), initial guess temperature, ignored if temperature is one of the flash specifications.

    •  K0 (optional), initial guess for the constants K

    •  x0 (optional), initial guess for the composition of phase x

    •  y0 = optional, initial guess for the composition of phase y

    •  vol0 = optional, initial guesses for phase x and phase y volumes

    •  atol = absolute tolerance to stop the calculation

    •  rtol = relative tolerance to stop the calculation

    •  max_iters = maximum number of iterations

[longemen3000]

on the latest master branch, the problems present with multiparameter fluids seem to be gone:

[longemen3000]

On pharmaPCSAFT, there was a bug on the initial point for gases. You see, we use a virial approximation for initial guess of the gas phase (and its a pretty good guess), but it fails above the virial inflection point (B > 0).

The 0.6.13 release contains a lot of fixes related to this issue, feel free to test it.

[B-LIE]

Hm. I just upgraded to the latest release of Clapeyron.jl (v0.6.13). I still get problems with Gerg2008 and MultiFluid:

NOTES:

• I only provided an initial guess of p0, not the K-value (I don't know what the K value is, or how I should guess an initial value...).
• I used the same value of p0 in all cases. I realize that I could use the previous value of p0 as I gradually increase $\chi_\mathrm{H_2}$, but then I have to change my code a bit.

I tested the second bullet point, i.e., using the previously computed p0:

sol_mf = [vt_flash(mod_mf,V_c,T,[n_H2O_c,X*n_H2_c]; p0=3e6) for X in Xvec_H2]
#=
sol_mf = Vector(undef,nX)
_p0 = 3e6
for i in 1:nX
    sol_mf[i] = vt_flash(mod_mf,V_c,T,[n_H2O_c,Xvec_H2[i]*n_H2_c]; p0=_p0)
    _p0 = sol_mf[i][4].p
end
=#

The above code shows my original way to do it, where vt_flash computes virtually zero mole fraction H2 for low pressures/amounts of H2 (in the order of 1e-10; see plot below), and then suddenly jumps a little.

So I thought the second approach should work better, where I use the previous pressure as initial guess. Not good: this approach fails to find the solution at all when i>9.

So I am curious how you managed to get around the problem with Gerg2008 (and MultiFluid)?

---

A final question: which of the models do you think it is best? Consider:

NOTE that Gerg2008 and MultiFluid finds zero H2 in the liquid phase up to a little more than $\chi_\mathrm{H_2}=1$.

Compare this to the cubic EoS'es:

[golghanddashtih]

@longemen3000 Hi. I am trying to obtain fugacity coefficient for a mixture of water and hydrogen! The calculated fugacity coefficients are surprisingly big! Does that make sense? The code is given below with the given fugacity coefficients:

using Pkg
using Clapeyron
pr3_model_H2H2O = PR(["hydrogen", "water"]) 
P= 10e6
T=298
z=[0.5,0.5] 
fugacity_coefficient(pr3_model_H2H2O, P, T, z; phase=:v)

Output:
2-element Vector{Float64}:
 35.06811757640551
  0.0040527247770437294