Basic questions on Clapeyron.jl

[B-LIE]

OK -- I'm not a thermodynamics expert, and perhaps my questions are "stupid". But I'm genuinely considering how I can use Clapeyron.jl in my work. So I hope the questions are not too out-of-place.

0. Preliminaries: two models

Water, Peng-Robinson, and a mixture of hydrogen, oxygen and water appearing in Fuel Cell models and Water electrolysis.

using Clapeyron
model_w_pr = PR(["water"])   # The familiar component
model_pem = PR(["hydrogen", "oxygen", "water"])    # mixture used in fuel cells, etc.

First tests:

julia> p_w_pr = pressure(model_w_pr,30e-3,373.15)
102502.26906925516
julia> volume(model_w_pr,p_w_pr,373.15)
2.2674765387029373e-5

Question: Why is the volume not found to be 30e-3???

OK -- what if I do:

julia> volume(model_w_pr,p_w_pr,373.15; phase=:liquid)
2.2674765387029373e-5
julia> volume(model_w_pr,p_w_pr,373.15; phase=:gas)
0.029999999999999995
julia> volume(model_w_pr,5*p_w_pr,373.15; phase=:liquid)
2.267239665085835e-5

Observations:

1. OK -- when I specify phase :liquid, the volume is the same as with the default phase.
2. When I specify phase :gas, the volume is the same (more or less) as the volume specified when computing the pressure
3. When I specify phase :liquid for a 5 times larger pressure, the (liquid) volume barely changes.

Questions: I'm confused. The "base case" is more or less water at 100C and 1 atm pressure.

• What does the volume even mean? Does this reflect that the water is in liquid form, and that the volume essentially is 2.27e-5 m^3 = 0.0227 L? Why don't I get back the volume 30e-3 = 30 L that was used to compute the pressure?
• What does it then mean when I specify phase :gas, where I do get back the volume that was used to compute the pressure?

1. Flash

In electrolysis of water to produce hydrogen, it is necessary to have two separation tanks: one for separating water and oxygen, and one for separating water and hydrogen. Both tanks will also be "cross-polluted", thus the tank with water and oxygen also contains traces of hydrogen, and the one with water and hydrogen contains traces of oxygen.

In a first attempt, I assume that the water is completely in water phase (temperature ca. 80C, pressures > 1 bar), and that there is no dissolution of the gases + they obey the ideal gas law. This allows me to compute the water volume, and then the gas volume in the tank. The partial pressures of the gases can be computed from the ideal gas law. I then "cheat" and add the saturation pressure of water ("cheat" in the sense that that probably means that there is water in the gas volume).

I need a more accurate model where I take into account dissolved gases in the liquid phase. At the outset, I know:

• separation tank volume (fixed)
• temperature (from energy balance, or similar)
• the number of moles of each component in the tank volume (from mole balances)

Question: How can I compute the flash using Clapeyron.jl so that I can find:

• tank pressure + partial pressures of each component in the gas phase
• liquid volume
• amount (mole) of each component in the liquid phase?

First ideas

• With given separator tank volume, temperature, and number of moles, compute the pressure(),
• Compute chemical potential in each phase for each component with the found pressure, given temperature, and amount (mole), and use some nonlinear solver to equate them?
• Or is there a suitable prebuilt Flash algorithm? If so, which is recommended?

2. Reaction equilibrium

I need to compute Gibb's (partial?/molar?) energy in order to find the so-called Nernst (equilibrium) voltage for a water Fuel Cell (or: water Electrolysis, which is the reverse).

Consider the reaction [preview doesn't work for LaTeX math...]:

H2 + 0.5 O2 -> H2O

Assuming water is the product, the stoichiometric vector is

nu = [-1, -1/2, +1]

where the components are ordered as ["hydrogen", "oxygen", "water"].

In Fuel Cell applications, hydrogen and oxygen are available at two different electrodes (anode and cathode, respectively) separated by a membrane. Water is produced at the oxygen side (cathode), and thus is present in a mixture of oxygen and water.

On the hydrogen side (anode), water vapor is also normally added to ensure that the membrane doesn't become dry. But water does not participate in the reaction on the anode side.

Gibb's (molar) energy of reaction is:

where in my notation, G with tilde on top and subscript j is partial (or: molar?) Gibbs' (free) energy, which equals chemical potential mu_j [when p and T are input arguments].

In electro-chemistry, it is common to write:

where circle superscript indicates standard state, and the left-hand side symbol is (partial/molar) Gibb's energy at standard state (i.e., pressure equals 1 atm or sometimes 1 bar), while a_j is activity.

For hydrogen and oxygen, it is common to set activity equal to partial pressure divided by standard pressure $p^\circ$, $a_j = p_j/p^\circ$ (assuming ideal gas). For water, it is common to set activity either equal to saturation pressure of water divided by standard pressure at the given temperature, $a_\mathrm{H_2O} = p_\mathrm{H_2O}^\mathrm{sat}/p^\circ$ [assuming water is produced as vapor], or equal to the model fraction of water $a_\mathrm{H_2O} = x_\mathrm{H_2O} \approx 1$ [assuming water is produced as liquid]. The choice of activity for water has ramification for how Gibb's energy for water is computed [choice of heat capacity...].

A simple expression for Gibb's energy at the standard state is:

where the two first term on the RHS are Gibbs energy of formation, and entropy of formation, respectively. These formation energies are critical in thermodynamics of reactions.

Questions:

1. Is the above more or less ok? I guess ideal gas has been assumed for the $\tilde{G}_j^\circ$?

• The above formalism is used extensively in publications for fuel cells and water electrolysis.
• Essentially, the above allows for computing Gibbs energy $\tilde{G}_j$ as a function of temperature, and partial pressure of hydrogen (on the anode side), partial pressure of oxygen (on the cathode side), and either saturation pressure of water (given by temperature) or mole fraction of water on the cathode side.

2. How can the above instead be done using Clapeyron.jl?

• Should independent models be used on the anode side and the cathode side? With, say, hydrogen + possible added water vapor to avoid drying up the membrane on the anode side, and oxygen + water + possibly other components of air on the cathode side?
• Then use the chemical_potential function to compute $\mu_j = \tilde{G}_j$ on each electrode, and then compute the weighted sum with stoichiometric coefficients??
• If as above, how can I specify partial pressures in the chemical_potential function? Sum partial pressures to find total pressure, and scale amount of the different components to unity?

Note: using chemical_potential requires that the Gibbs energy of formation and the entropy of formation are used in Clapeyron.jl.

3. Enthalpy flow rate?

In the energy balance used to compute the temperature, there will be enthalpy flows in and out of the system. Using decoration "dot" to indicate flow rate, enthalpy flow is $\dot{H}$.

In a simple description, assuming ideal solution, one can write $\dot{H} = \sum \dot{H}_j$ where $\dot{H}_j = \dot{n}_j \cdot \tilde{H}_j$ where $\tilde{H}_j$ is the molar enthalpy of pure component $j$.

Questions:

1. How can I compute $\dot{H}$ in a more formally correct way using Clapeyron.jl?

• Is the formally correct expression $\dot{H} = \sum \dot{n}_j \cdot \bar{H}_j$ where $\bar{H}_j$ is the partial molar enthalpy of component $j$, i.e., even for non-ideal solutions?
• If so, is there a Clapeyron.jl function for computing this partial molar enthalpy? (Similar to chemical_potential(), which I think is the partial molar Gibbs energy when $p$ and $T$ are specified...)

Note: just like for partial molar Gibbs energy for reactions, a partial molar enthalpy must include enthalpy of formation (either directly, or via Gibbs energy/entropy of formation) in hold necessary information about reaction enthalpy.

[pw0908]

Hi Bernt, happy to see you using Clapeyron! That's a lot of questions which I'll try to answer. One thing I'll mention is that we have more-detailed docs which can be found here: https://clapeyronthermo.github.io/Clapeyron.jl/dev/

0. As mentioned in the docs, volume(model,p,T,z) obtains the volume at a given pressure, temperature and composition of the most stable single phase. Internally, what Clapeyron does is look for the vapour and liquid roots, subsequently checking which has the lowest gibbs free energy, and then returns the corresponding volume. This is true for any bulk property written as prop(model,p,T,z). When you specify the phase, it ignores this check and just solves for the desired phase. In your case, at the conditions you specified, the liquid phase is more stable than the vapour phase. As such, when you computed p_w_pr = pressure(model_w_pr,30e-3,373.15), you were essentially solving for a metastable state.
1. At this current stage, we primarily support TP-Flash algorithms. However, @longemen3000 has been doing a lot of work implementing the other different flash algorithms and can probably give you more information on these. Our flash algorithms can all be found here: https://clapeyronthermo.github.io/Clapeyron.jl/dev/properties/multi/#TP-Flash.
2. In the main branch, Clapeyron does not support reactive systems yet. I've made some headway implementing a reactive solver which can be found here: Reactive Systems #284. This is still in early stages but I believe it could handle your specific system as all components are neutral. However, to do so with phase equilibrium will be very challenging. I'm not sure if the addition of electrolysis complicates the whole system (Nernst equations are basically just regular equilibrium constant equations but in different units). Unfortunately, I haven't had the time to test this yet.
3. Clapeyron handles equilibrium systems only. We do not consider fluctuations in time at all; this would require packages like ProcessSimulator.jl. Nevertheless, the equation you present for the enthalpy rate is correct even in the non-ideal case as you are using partial molar enthalpies. Unfortunately, we do not have methods built in which can compute partial molar enthalpies in Clapeyron (yet). It may be faster, if you had $\dot{n}_i$ to just do the following:

$$\dot{n}_T = \sum_i \dot{n}_i$$ $$z_i = \dot{n}_i/\dot{n}_T$$ $$\dot{H} = h(p,T,z)\times\dot{n}_T$$

where $h$ will just give you the total molar enthalpy.

Sorry these replies are brief. Let me know if there's anything you'd like to follow-up on.

[B-LIE]

It may be faster, if you had $\dot{n}_i$ to just do the following...

I'm not sure this will not work in the case of chemical reactions. OK -- it may work, but then I have to manually add a "heat of reaction".

I tried to look at the ProcessSimulator.jl web page, but it has not been published. I'll try to contact them directly.

[B-LIE]

@pw0908 and @longemen3000
OK -- I have tested the flash algorithm that @pw0908 suggested for a mixture of water and oxygen.

In practice, I know:

• Separator volume [physical equipment], say V_a = 1 (m3)
• Temperature [from energy balance], say T = 60 + 273.15 (K)
• Amounts in volume [from mole balance], say n_H2O_a = 18.5e3 and n_O2_a = 24.08 (mol)

First thing, I want to know the pressure p. Here is the system.

using Clapeyron
model_a_pr = PR(["water", "oxygen"])
V_a = 1 # m3
T = 273.15 + 60 # K
p_a = 1.2e5 # Pa, 1.2 bar
n_H2O_a = 1.85e4 # mol H2O
n_O2_a = 24.08 # mol O2

I have indicated an "intended pressure" of 1.2 bar.

Here is what I find using Clapeyron.jl:

julia> pressure(model_a_pr, V_a, T, [n_H2O_a, n_O2_a])
-1.237421762869517e8

This is, of course, a physically impossible pressure. My suspicion is that the data I gave are not realistic/compatible with the PR EoS. The amounts I have specified are based on steady state assumptions of a dynamic model where I have assumed that water is pure liquid and oxygen is an ideal gas.

If I reduce the amounts to 19.135% of what I specify above, I get more or less the assumed pressure...

OK -- if I don't care about computing the pressure (not doable in practice, but let me assume so). Then I get:

julia> sol_fl = tp_flash(model_a_pr, p_a, T, [n_H2O_a,n_O2_a]);
julia> sol_fl[1]
2×2 Matrix{Float64}:
 0.999994  5.93641e-6
 0.151313  0.848687

Here, I understand the first column is for the first component, i.e., water, and the second column is for oxygen.

Furthermore, the first row must be liquid phase, and the second must be vapor/gas phase?

Next,

julia> sol_fl[2]
2×2 Matrix{Float64}:
 18495.7       0.109799
     4.27365  23.9702
julia> sum(sol_fl[2], dims=1)
1×2 Matrix{Float64}:
 18500.0  24.08

OK -- summing over the rows gives back the amounts I provide, which is reassuring.

Question: Is there (or is it possible to create) a flash algorithm which uses:

• model + V + T + total amounts in the volume as input arguments, and
• returns pressure + mole fractions in various phases + mole numbers in various phases?

Would this be "VT flash"? See https://link.springer.com/chapter/10.1007/978-3-031-49432-1_5 for some work on this with Julia. See also https://www.sciencedirect.com/science/article/pii/S0021999122003370 .

[longemen3000]

Yes, this corresponds to a Volume-Temperature flash. I've been developing an algorithm for general flashes (any combination of P-V-T-H-S-U-gas fraction).

In particular, the flash in volume-temperature specification can be used via the vt_flash(model,V,T,z) function, but the current implementation in the latest release is not sufficiently stable to work with with example. In the master branch, however:

model_a_pr = PR(["water", "oxygen"])
V_a = 1 # m3
T = 273.15 + 60 # K
p_a = 1.2e5 # Pa, 1.2 bar
n_H2O_a = 1.85e4 # mol H2O
n_O2_a = 24.08 # mol O2
sol_fl = vt_flash(model_a_pr, V_a, T, [n_H2O_a, n_O2_a])
volume(sol_fl) - V_a #to test if the result is correct

that should result in the following:

julia> sol_fl = vt_flash(model_a_pr, V_a, T, [n_H2O_a, n_O2_a])
Flash result at T = 333.15, p = 1.29611e5 with 2 phases:
 (x = [0.999994, 6.4949e-6], β = 18496.2, v = 2.19268e-5)
 (x = [0.140186, 0.859814], β = 27.8663, v = 0.0213318)

julia> volume(sol_fl) - V_a
9.818812429784884e-13

julia> pressure(sol_fl)
129610.55274230446

[B-LIE]

Looks great! A question...

1. Any idea when this will be sufficiently tested to push out in the release version of Clapeyron?

Reason for my question: a PhD student will take a course in the spring of 2025 where one possibility is to use VT flash.

She will do it within a dynamic model, so she may have to develop a surrogate/proxy model to avoid costly implicit computations.

[B-LIE]

One more comment on partial molar properties... I've given it some thought. Essentially, you already compute the chemical potential, which is equal to the partial molar Gibb's energy.

What about the other ones? In wikipedia (https://en.wikipedia.org/wiki/Partial_molar_property) I found the following relations between partial molar properties, which are linear if $p$ and $T$ are known:

These three relations relate 6 partial molar properties: partial molar Gibbs energy $\bar{G}_j$ with internal energy $\bar{U}_j$, enthalpy $\bar{H}_j$, Helmholtz energy $\bar{A}_j$, entropy $\bar{S}_j$, and volume $\bar{V}_j$.

We already know partial molar Gibbs energy = chemical potential. So we need two more partial molar properties in order to compute all of them.

I found a book by Reinhard Hentschke where it is referred to a Gibbs-Helmholtz relation:

Similarly, the book suggests that:

Thus, to find partial molar enthalpy, if Clapeyron.jl can do automatic differentiation of the chemical potential (which is a function of p, T), then it should be straightforward to compute partial molar enthalpy??

Computing partial internal energy is perhaps more tricky, because that one assumes that the chemical potential is a function of V,T. But perhaps do-able?

If both partial molar enthalpy and partial molar internal energy can computed -- in addition to partial molar Gibbs energy = chemical potential, then all 6 partial molar properties should be available.

OK -- a big caveat: I am no thermo expert, so I may be wrong :-o .

[longemen3000]

Any idea when this will be sufficiently tested to push out in the release version of Clapeyron?

We can release a new version with partial properties and the fixed VT flash sometime next week. on testing, we always welcome tests submitted as issues with an MVE. for example, all the reported errors in #325 were added as tests.

On partial properties, in fact, we have the opposite problem, all our inner problems are in V-T space, so we have to transform all partial problems into a V-T formulation first. in particular, for a particular function F(V,T,n), one can obtain their partial property in the following way:

That relation is general (for example, one can calculate partial isobaric heat capacities), but takes more evaluations (two gradients + 1 dpdV + 1 dFdV), so alternative, faster forms are always good to add

[B-LIE]

This sounds great!!

Ah yes. You base the computations on the Helmholtz energy with variables V,T. [My understanding is that the reason for this is that the involved integrals are much simpler with V,T as variables -- e.g., with cubic EoS, the model is explicit in p with V,T as variables, while they are implicit cubic in V with p,T as variables.]

I look forward to the release!!! and will check it out [I have a research report due on January 31, 2025, so perhaps a little slow prior to that].

A few things I'm curious about:

• Is it realistic that the VT flash calculations can handle
  -> water + oxygen (as indicated) with small traces of hydrogen
  -> water + oxygen (as indicated) + nitrogen + argon (from air) with small traces of hydrogen
  -> water + hydrogen at a much higher pressure (30-50 bar) and small traces of oxygen
  -> hydrogen with traces of water at up to 750 bar
• Do the computation of enthalpy, Gibbs energy, etc. include energies of formation? [Necessary to get correct results with chemical reactions...]

OK -- I'll try to figure out how to use the functions with ModelingToolkit.jl, too. I'm trying to get in touch with the ProcessSimulator.jl people (whom I know a little bit), to learn how they are doing things.

Computation time is an issue, of course. Some 25 years ago, I did some work on water-ethanol distillation, with an implicit equilibrium model (azeotropic) which I solved for each time step using a Newton solver. I also tried to solve the equilibrium conditions off-line, and fitted an explicit polynomial function to the data.

Using the explicit polynomial to do the equilibrium computations speeded up the simulation time by a factor of 100. Conceptually, I could do something similar here, but it is much more demanding to fit an explicit function when multiple components (water, oxygen, etc.) are involved because the space is much larger.

[B-LIE]

@longemen3000
I made a very, very crude code for computing VT flash for a mixture of water and ideal gases; main case involves oxygen, hydrogen, nitrogen, argon, etc.

I'm curious as to how wrong results I get... NOTE: I assumed water is spread with substantial amounts in liquid and gas phases (Antoine's law), but that the ideal gases essentially are found in the gas phase. I would need to add some Henry's law constants to compute the (assumed) minor fraction of the ideal gases in the liquid phase...

Here is my code:

# Water saturation pressure
function p_H2O_sat(T)
    A_H2O = 5.11564
    B_H2O = 1687.537
    C_H2O = 230.17
    p_o = 1e5 # Pa
    p_sat = p_o*10^(A_H2O - B_H2O/(C_H2O + T - 273.15))
    return p_sat
end

#
function my_VT_flash(V,T,n)
    R = 8.31 # J/K/mol
    ρ_H2O_ℓ = 1e3 # kg/m3
    M_H2O = 18e-3 # kg/mol
    # 
    p_H2O = p_H2O_sat(T)
    n_H2O = n[1]
    n_ig = n[2:end]
    #
    V_ℓ = (n_H2O - p_H2O*V/(R*T))/(ρ_H2O_ℓ/M_H2O - p_H2O/(R*T))
    V_g = V - V_ℓ
    p_ig = n_ig*R*T/V_g
    p = p_H2O + sum(p_ig)
    n_H2O_ℓ = V_ℓ*ρ_H2O_ℓ/M_H2O
    n_H2O_g = n_H2O - n_H2O_ℓ
    return V_ℓ, V_g,  [p_H2O;p_ig], p, n_H2O_ℓ, n_H2O_g
end

My case is the following:

V = 1 # m3
T = 70 + 273.15 # K
n_H2O = 1.85e4
n_O2 = 24.08

When running the code, I get:

julia> V_L, V_G, p_vec, p, n_H2O_L, n_H2O_G = my_VT_flash(V, T, [n_H2O,n_O2]);
julia> V_L, V_G, V_L+V_G
(0.33286874943117706, 0.6671312505688229, 1.0)
julia> p_vec, sum(p_vec), p
([31167.53323759674, 102927.21568874586], 134094.74892634258, 134094.74892634258)
julia> n_H2O_L, n_H2O_G, n_H2O_L + n_H2O_G, n_H2O
(18492.70830173206, 7.291698267941683, 18500.0, 18500.0)
julia> p_vec/p
2-element Vector{Float64}:
 0.23242918523764772
 0.7675708147623524

Let me know when I can compare my results with a proper VT flash algorithm :-)

[longemen3000]

curiously, i got really close results if i use the GERG-2008 reference equation (for natural gas):

julia> model_b = GERG2008(["water","oxygen"])
MultiFluid{EmpiricAncillary, AsymmetricMixing, EmpiricDeparture} with 2 components:
 "water"
 "oxygen"
Contains parameters: Mw, Tc, Pc, Vc, Tr, Vr, acentricfactor, lb_volume, reference_state

julia> vt_flash(model_b,V,T,n)
Flash result at T = 343.15, p = 1.35194e5 with 2 phases:
 (x = [1.0, 3.68867e-7], β = 18492.8, v = 1.8424e-5)
 (x = [0.230743, 0.769257], β = 31.2941, v = 0.0210675)

[B-LIE]

Interesting. Looking forward to test your code more fully and comparing it to mine!! (I know mine is crude, but there is no iteration involved)

[longemen3000]

It does make sense, in those combinations of temperature-pressure, a really good approximation is supposing there is no dissolved gas.
But this approximation does not allow for calculating the fraction of dissolved gases in water, as you said before, you need additional information (that the EoS in this case provides).

[B-LIE]

In developing my code, I assumed partial pressure of water given by Antoine's law. All other gases were assumed to follow ideal gas law.

I then assumed that, by definition, partial pressures are given as mole fraction multiplied by total pressure (is this a definition, or does it assume ideal solution?).

To my surprise, it follows that water vapor also behaves like ideal gas... Is that because all other gases are ideal gases, or is it because my definition of partial pressure only is valid for ideal gas??

Just curious...

[B-LIE]

julia> vt_flash(model_b,V,T,n)
Flash result at T = 343.15, p = 1.35194e5 with 2 phases:
(x = [1.0, 3.68867e-7], β = 18492.8, v = 1.8424e-5)
(x = [0.230743, 0.769257], β = 31.2941, v = 0.0210675)

I'm trying to guess the interpretation of the result...

• x is the mole fractions in the order of components specified when defining the model, right?
• beta is the total amount in the given phase (mole), right?
• v is the molar volume of the phase, right? I.e., the volume is beta*v?

Apart from this, any particular recommendation about which EoS to use? Did you choose GERG2008 because that is the best?

[longemen3000]

Yes, you are correct. there are a lot of convenience functions to access the result (volume(result), pressure(result),Clapeyron.temperature(result), enthalpy(model,result)), and you can access the vectors of properties via result.volumes, result.compositions and result.fractions.

On the model selection, i selected GERG2008 just for their good reproduction of water properties (volume + saturation pressures) in conjunction with good properties for the gases + adequate extrapolation at moderate conditions (i could be wrong in this last one, but weak interactions between water + gases don't require additional parameters). the main disadvantage of the model is that it takes more time to solve than a cubic (around 10x). so there are tradeoff between model accuracy and computational time.

Of course, model selection is a whole area of research of its own. If i understand correctly, for this particular application, you require good reproductions of saturation pressures + volumes, so probably a SAFT could be a good option if you want predictive properties. In particular, pharmaPCSAFT (PCSAFT with water-specific temperature-dependent sigma), SAFTVRMie or PCPSAFT are good candidates for this particular combination of components. A translated cubic could also work, But in my experience, they overcorrect the water volume.

[B-LIE]

A translated cubic could also work,

I have tried model = VTPR(["water", "oxygen"]), but this doesn't work: I get an error message:

MissingException: Missing values exist ∈ single parameter groups: ["oxygen"].
Stacktrace:
  [1] SingleMissingError(param::Clapeyron.SingleParameter; all::Any)
    @ Clapeyron C:\Users\Bernt\.julia\packages\Clapeyron\2fEkA\src\database\database_rawparam.jl:360
  [2] SingleMissingError(param::Clapeyron.SingleParameter)
    @ Clapeyron C:\Users\Bernt\.julia\packages\Clapeyron\2fEkA\src\database\database_rawparam.jl:353
  [3] is_valid_param(param::Clapeyron.SingleParameter, options::Any)
....

It seems like VTPR doesn't support components such as oxygen, hydrogen, CO2, etc. But it works with water, ethanol, etc.

[longemen3000]

VTPR does have translation, but it also requires UNIFAC parameters. For this particular case, i was thinking along the lines of tcPR

[B-LIE]

SAFTVRMie

... gives error MissingException: Missing values exist ∈ single parameter sigma: ["oxygen"].

The two other SAFT models took much longer time to solve than the cubic equations.

I'll try to do the computations for various EoS over a range of temperatures tomorrow + compare them + compare to my crude algorithm tomorrow, and report back.

[B-LIE]

Some results...

Models I tested:

Here, mod_w is my coarse model (also referred to as simple).

It is interesting to see how long time it takes to solve flash the first time (including JIT compilation):

My simple function is by far the fastest initially (can be further improved). Cubics and GERG2008 are similar. The SAFT models are quite slow.

However, when I do flash for 72 temperatures between 15-95 C, things change and flash time is much more similar.

One observation: GERG2008 fails for most temperatures (either computing one phase, or getting NaN). Cubic EoS fails for some temperatures. Some results (I skipped GERG2008 because it failed a lot):

Note that PCSAFT and pharmaPCSAFT give identical (seemingly) results.

I'm curious as to what conclusions can be drawn wrt. which method is best... Not my simple algorithm, I think, although in many of the plots it gives a result quite similar to PCPSAFT. I can probably speed my simple algorithm up a little bit, but I'm not a super programmer...

[longemen3000]

Interesting observations indeed, some notes:

• SAFT models using water need to solve the association problem on each evaluation. We do have a sophisticated polyalgorithm to solve that problem, it takes a lot of compilation time (at least that time cost is only paid once).

• it is curious that GERG2008 fails this much, but i suppose that the EoS was not fitted with this particular application in mind, if we look at the range of compositions for the oxygen + water mixture that were used to fit the GERG2008 EoS, the conditions that are being tested are outside the composition range.

  

  Probably, the full multiparameter (IAPWS95 + Oxygen reference) could work in this case (using CoolProp; model = MultiFluid(["water","oxygen"]). if it doesn't then i will have to take another look at the solver, or implement a pressure-based algorithm (multiparameter EoS are really delicate around saturation lines and old ones (before 2010) normally show unphysical equilibria at low temperatures)

• pharmaPCSAFT and PCSAFT should be almost the same, but pharmaPCSAFT should give more accurate water volumes.

[B-LIE]

OK. I can test MultiFluid later today or tomorrow. Note: my computation times and the relative performances vary depending on how many temperatures I solve for. I simply used @timed to see how long it took to find the solutions for the 72 (or 9 or 18 or what not) temperature points. I should probably use package BenchmarkTools.jl instead and do a proper test. I was impressed at how fast the Clapeyron flash solver seemed to be.

For my own routine, I returned the results in an array -- I didn't spend time trying to figure out the data structure of Clapeyron -- essentially, I made the call as similar as possible to Clapeyron's vt_flash. In my case, the mod_w contains the function(s) for computing the Henry law volatilities.

[B-LIE]

OK... I added the MultiFluid. In addition, I changed the timing comparison from using @timed on building up the solution via [vt_flash(..., V, T,...) for T in Tvec], which includes timing for the for loop, into using ´@btimedonvt_flash(..., V, T, ...)` , which should measure the timing of the flash more directly (??).

Here is one result:

With this way of measuring computational time, my simple-minded algorithm really is much faster than iterative flash -- which is what I would have expected.

I don't know which of the methods is the faster one, but I could probably tune my model as follows:

• using pure water and flash for "the best method". I use Antoine's model for water saturation pressure (virtually no difference from Riedel's more accurate correlation in this temperature span), and Rackett's correlation for molar water volume as a function of temperature. I suspect Rackett's correlation can be re-tuned.
• I use the Wikipedia model for Henry's law parameters, which are taken from a recent paper by Sander. I suspect I can tune the Henry's law parameters and get quite good fit with, say, the MultiFluid model.
• Of course, in my simple model, I don't include combination effects if I add more components like Hydrogen, Nitrogen, etc.

So -- I'm curious of your thoughts. [Btw: I measure Gerg2008 to be 23_000 times slower than my simple-minded algorithm.]

[longemen3000]

In those conditions, i think the best model is a combination of:

• a good saturated liquid water correlation for volume (there is a chebyshev polynomial that reproduces the IAPWS95 results in the EoSSuperancillaries.jl repository. You can also use thr IF-97 standard, but any other accurate correlation for saturated water will suffice
• ideal or virial gas (you can extract virial coefficients from equations of state usingsecond_virial_coefficient(model, T, z), or you can use a correlation)
• an accurate correlation for vapour pressure (antoine or DIPPR101, there is also the IF-97 standard that provides a saturation pressure curve)
• henry coefficients for the fraction of gases dissolved on water (one claver way would be to extract the henry coefficients from an EoS model and use those. Another way would be to just suppose constant gas composition, T and P and directly find what is the liquid conposition, by solving the chemical potential equality system)
• 1 or 2 corrective steps accounting for the difference of pressure caused by the dissolved gases (if there is any difference)

[B-LIE]

Thanks a lot for quick responses and constructive inputs. I am super impressed by Clapeyron.jl and the work you guys are doing. I'm sure Clapeyron.jl could be a backbone for a Thermodynamics course, and make such courses much more concrete than they typically are. Also, for me -- even if I may end up using my super simple flash algorithm in my particular case (dynamic models solution via ModelingToolkit.jl), I could not have felt confident in it had it not been that Clapeyron.jl gives similar results for much more rigorous models. And I will, of course, also see if I can make the dynamic simulations work also with, say, one of the Clapeyron.jl routines.

One more question + a final comment for now.

1. Question: In my development, I assumed that partial pressure $p_j$ by definition is related to total pressure $p$ via the mole fraction $y_j$ as $p_j \triangleq y_j\cdot p$.

• Is it correct that this is a definition of partial pressure, or is there an underlying assumption of ideal gas (or ideal solution) here?

2. Comment: I will soon start to play around with Clapeyron.jl for computing with Gibbs (molar) energy and Enthalpy, incluiding with reactions.

[longemen3000]

Hello,

The sources i checked (Perry, Abbott & Van Ness) define partial pressure as yᵢ*p. This definition follows the criteria of sum(p_i) = p only for the case of ideal gases (dalton's law).

I couldn't find any working equation to define partial pressures on real gases using the definition of partial property (a partial property is defined at constant pressure, making dp/dn |p = p0 a constant with zero derivative). There are some alternative definitions, but none convince me for the specific context of partial properties.

Just to note, i'm looking at a partial pressure strictly as a definition of partial property, that is:

• sum(pᵢ) = p
• pᵢ is a function depending on the gradient of some scalar property.

For equilibria in general, the correct way is to work in terms of fugacities. partial pressures only appear when supposing that one of the phases is an ideal gas.