Question: substances available in various EoS:s?

[B-LIE]

I am playing around with Clapeyron with ModelingToolkit, and consider a chemical reaction: Glycol production:

Reaction:
C2H4O + H2O -> C2H6O2

("oxirane", "water", "ethylene glycol")

PR accepts these substances, but not PCSAFT or VTPR.

Since these substances constitute a non-ideal solution, it would have been nice to test PCSAFT. But that seems to not be possible?

[longemen3000]

Hello,

I checked on the parameters, and the PCPSAFT EoS seems to have the parameters for oxirane. Another option would be to use fit PCSAFT parameters using the i-PC-SAFT approach (there is a new paper that lines the steps to do this https://doi.org/10.1016/j.fluid.2025.114649).

Also, While VTPR does not have parameters for oxirane, UNIFAC and PSRK do ((oxirane, ["H2COCH2"=>1] in Clapeyron group notation), you could mix-and match a cubic that uses UNIFAC as an activity coefficient model to get something similar to VTPR. (UMRPR uses UNIFAC as their mixing rule, but their alpha correlation seems to be fitted on linear hydrocarbons, i would change that)

[B-LIE]

Thanks! I'll check out PCPSAFT first.

[B-LIE]

Follow-up question... I need to specify the reference in the form of enthalpy of formation + entropy + temperature + pressure + phase.

Google AI Studio suggests the follow statement:

Clapeyron.set_reference_state!(model; 
        T0=298.15, 
        P0=101325.0, 
        H0=Hf0_vec, 
        S0=S0_vec, 
        phase=:gas
    )

Here, Hf0_vec and S0_vec are vectors of real values for the formation enthalpy and the reference entropy.

The documentation says:

set_reference_state!(model, new_ref) where new_ref is a ReferenceState??? So does not list the form that Google AI Studio suggests... (AI Studio suggests keywords for T0, etc., which is very different from a ReferenceState argument.)

In the past, I have used (e.g.):

my_ref = ReferenceState(:volume; H0=Hf0_vec, S0=S0_vec, T0=298.15, P0=101325.0, phase=:gas)
model = PCPSAFT(["water", "ethanol"]; reference_state=my_ref)

Question: which of these methods, if any, are correct at the moment?

[longemen3000]

the last one is the correct method, but you need to provide an adequate ideal model for that, like ReidIdeal. the following code should work:

my_ref = ReferenceState(:volume; H0=Hf0_vec, S0=S0_vec, T0=298.15, P0=101325.0, phase=:gas)
model = PCPSAFT(["water", "ethanol"]; reference_state=my_ref,idealmodel = ReidIdeal)

if you don't have any Cp data, you could use MonomerIdeal, but keep in mind that MonomerIdeal assumes the same Cp (2.5R) for all substances, same as BasicIdeal (What MonomerIdeal does is defining their own reference state in an implicit manner)

I probably need a better error for the case.

Setting the reference state after the model construction is also valid with the correct ideal model:

model = PCPSAFT(["water", "ethanol"];idealmodel = ReidIdeal)
my_ref = ReferenceState(:volume; H0=Hf0_vec, S0=S0_vec, T0=298.15, P0=101325.0, phase=:gas)
Clapeyron.set_reference_state(model,my_ref)

[B-LIE]

Thanks. Doesn't the PCPSAFT model include Cp data? Alternatively, does the ReidIdeal contain adequate CP data? Sent from Outlook for Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Andrés Riedemann ***@***.***> Sent: Sunday, February 22, 2026 9:35:58 PM To: ClapeyronThermo/Clapeyron.jl ***@***.***> Cc: Bernt Lie ***@***.***>; Author ***@***.***> Subject: Re: [ClapeyronThermo/Clapeyron.jl] Question: substances available in various EoS:s? (Issue #548) [https://avatars.githubusercontent.com/u/38795484?s=20&v=4]longemen3000 left a comment (ClapeyronThermo/Clapeyron.jl#548)<#548 (comment)> the last one is the correct method, but you need to provide an adequate ideal model for that, like ReidIdeal. the following code should work: my_ref = ReferenceState(:volume; H0=Hf0_vec, S0=S0_vec, T0=298.15, P0=101325.0, phase=:gas) model = PCPSAFT(["water", "ethanol"]; reference_state=my_ref,idealmodel = ReidIdeal) if you don't have any Cp data, you could use MonomerIdeal, but keep in mind that MonomerIdeal assumes the same Cp (2.5R) for all substances, same as BasicIdeal (What MonomerIdeal does is defining their own reference state in an implicit manner) I probably need a better error for the case. — Reply to this email directly, view it on GitHub<#548 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AM5VM77FLLHJ3QSUPI5LKVD4NIHK5AVCNFSM6AAAAACV3JPTYOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTSNBRGY3DKMJZGA>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[longemen3000]

The PCPSAFT model does not include Cp data. For ReidIdeal we have data of 308 components that was taken mostly from Poling, 4th ed. Additional components could be fitted easily if needed (we did this for R-1234ZE(E), using the CoolProp ideal model as reference)

[B-LIE]

Ok. Hope the ReidIdeal contains CP for my substances... Water, oxirane, ethyl glycol, or that I can figure out how to add them Sent from Outlook for Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Andrés Riedemann ***@***.***> Sent: Sunday, February 22, 2026 9:46:49 PM To: ClapeyronThermo/Clapeyron.jl ***@***.***> Cc: Bernt Lie ***@***.***>; Author ***@***.***> Subject: Re: [ClapeyronThermo/Clapeyron.jl] Question: substances available in various EoS:s? (Issue #548) [https://avatars.githubusercontent.com/u/38795484?s=20&v=4]longemen3000 left a comment (ClapeyronThermo/Clapeyron.jl#548)<#548 (comment)> The PCPSAFT model does not include Cp data. For ReidIdeal we have data of 308 components that was taken mostly from Poling, 4th ed. Additional components could be fitted easily if needed (we did this for R-1234ZE(E), using the CoolProp ideal model as reference) — Reply to this email directly, view it on GitHub<#548 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AM5VM7ZM6PDX4APIZPKB2I34NIITTAVCNFSM6AAAAACV3JPTYOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTSNBRGY3TSMZWHA>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[B-LIE]

Seems like ReidIdeal does not contain data for the substances I planned to use...

julia> ReidIdeal(["water", "oxirane", "ethylene glycol"])
MissingException: Missing values exist ∈ single parameter c: ["oxirane", "ethylene glycol"].

Stacktrace:...

MonomerIdeal works, but as you say, that gives unrealistic values for Cp...

I do have the 6th edition of Properties of Gases and Liquids from 2023 (new authors). The sixth edition includes values for both "oxirane" and "ethylene glycol". The fifth edition does not.

Is there a way I can "specify" the values for Cp? I assume I should give a temperature function, right? I think the 6th edition provides (at least) two different correlations for Cp, but if these materials are given in a form that does not fit the standard in Clapeyron.jl, I can refit the expressions to the valid one.

[longemen3000]

i got access to the specific appendix. there are actually three types of Cp correlations:

• Polynomial
• Aly-Lee terms
• Plank-Einstein terms

For ethylene glycol and oxirane, you could fit 4th order polynomials and pass those to ReidIdeal in the following way:

components = ["component 1","component 2"]
a = [0.1,0.2]#a coefficients for comp1 and comp2
b = [0.2,0.2]#b coefficients for comp1 and comp2
c = [0.2,0.2]#c coefficients for comp1 and comp2
d = [0.2,0.2]#d coefficients for comp1 and comp2
e = [0.2,0.2]#e coefficients for comp1 and comp2

Mw = [16.0,32.0] 
model = ReidIdeal(components,userlocations = (a = a,b = b,c = c,d = d, e = e,,Mw = Mw))

For oxirane, i did a fit based on the coefficients in the NIST page, i got the coefficients (a,b,c,d,e) = (32.13903604109667, -0.030605030520446965, 0.0003832005115004867, -3.816262177445897e-7, 1.1345293825552631e-10):

The fit is not ideal, but it works up until around 1200 K:

[B-LIE]

Great! I tried something similar, with help or Google AI Studio (it is helpful, but extremely frustrating at times...). But your suggestion helps me understand things better! I will tinker a little on my code, and respond back. Sent from Outlook for Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Andrés Riedemann ***@***.***> Sent: Tuesday, February 24, 2026 5:14:06 AM To: ClapeyronThermo/Clapeyron.jl ***@***.***> Cc: Bernt Lie ***@***.***>; Author ***@***.***> Subject: Re: [ClapeyronThermo/Clapeyron.jl] Question: substances available in various EoS:s? (Issue #548) [https://avatars.githubusercontent.com/u/38795484?s=20&v=4]longemen3000 left a comment (ClapeyronThermo/Clapeyron.jl#548)<#548 (comment)> i got access to the specific appendix. there are actually three types of Cp correlations: * Polynomial * Aly-Lee terms * Plank-Einstein terms For ethylene glycol and oxirane, you could fit 4th order polynomials and pass those to ReidIdeal in the following way: components = ["component 1","component 2"] a = [0.1,0.2] #a coefficients for comp1 and comp2 b = [0.2,0.2]#b coefficients for comp1 and comp2 c = [0.2,0.2]#c coefficients for comp1 and comp2 e = [0.2,0.2]#d coefficients for comp1 and comp2 f = [0.2,0.2]#e coefficients for comp1 and comp2 Mw = [16.0,32.0] model = ReidIdeal(components,userlocations = (a = a,b = b,c = c,d = d, e = e,f = f,Mw = Mw)) For oxirane, i did a fit based on the coefficients in the NIST page, i got the coefficients (a,b,c,d,e,f) = (32.13903604109667, -0.030605030520446965, 0.0003832005115004867, -3.816262177445897e-7, 1.1345293825552631e-10): The fit is not ideal, but it works up until around 1200 K: image.png (view on web)<https://github.com/user-attachments/assets/881236cc-bf1d-4dbe-805f-3e2edd3af78b> — Reply to this email directly, view it on GitHub<#548 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AM5VM76NYXQFNWOIQL6ZV7T4NPFY5AVCNFSM6AAAAACV3JPTYOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTSNBYHA3DOOJTHA>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[B-LIE]

OK, I think my code works... [I have used Google AI Studio extensively. I could not have made it without that tool, but it is also at times frustrating to use...]

I have made two modules, one fairly general ThermoEoS module, and a specific Glycol_Chemistry module. The data are specified in Glycol_Chemistry:

module Glycol_Chemistry

using LinearAlgebra
using PhysicalConstants.CODATA2018: R
using Unitful: ustrip, @u_str

const substances = ["water", "oxirane", "ethylene glycol"]

# Absolute Thermo Data (Ideal Gas @ 298.15 K)
# [1: Water, 2: Oxirane, 3: Ethylene Glycol]
const Hf0 = [-241.826e3, -52.63e3, -394.4e3]
const S0  = [188.835, 242.4, 311.8]

# Heat Capacity Data: Properties of Gases and Liquids, 6th Ed. (2023)
const cp_raw = [
    # 1: Water (T range: 100 to 2273.1 K)
    (model=:AlyLee, Tmin=100.0, Tmax=2273.1, A=33.362, B=26.793, C=2610.5, D=8.896, E=1169.0),
    # 2: Oxirane (T range: 50 to 1500 K)
    (model=:AlyLee, Tmin=50.0,  Tmax=1500.0, A=6.388, B=125.41, C=911.2,  D=62.43, E=2201.0),
    # 3: Ethylene Glycol (T range: 20 to 1500 K)
    (model=:Einstein, Tmin=20.0, Tmax=1500.0, A=35.49, B=115.23, C=850.5, D=45.14, E=2105.0, F=0.0, G=0.0)
]

# Stoichiometry: 1 Water + 1 EO -> 1 EG
const nu = [-1.0, -1.0, 1.0]

# Kinetics (Liquid Phase)
const k0 = 4.5e8     # m3/(mol*s)
const Ea = 75000.0   # J/mol
const R_gas = ustrip(u"J/(mol*K)", R)

function production_rates(T, V, n)
    c = n ./ V 
    k = k0 * exp(-Ea / (R_gas * T))
    r = [k * c[1] * c[2]] 
    return (transpose(nu) * r)[:]
end

function get_thermo_package()
    return (names=substances, Hf0=Hf0, S0=S0, cp_raw=cp_raw, nu=nu)
end

end # module

OK -- here is your speciality: the ThermoEoS module -- the plan is to pass this on to the MTK components.

module ThermoEoS

using Clapeyron
using Polynomials
using ForwardDiff

const STATE = Ref{NamedTuple}()

# Helper: Evaluates the raw 6th Ed. equations for sampling
function evaluate_cp_ideal(data, T)
    if data.model == :Reid
        return data.A + data.B*T + data.C*T^2 + data.D*T^3
    elseif data.model == :AlyLee
        # Eq 2: A + B*(C/T / sinh(C/T))^2 + D*(E/T / cosh(E/T))^2
        u1 = data.C / T
        u2 = data.E / T
        term1 = iszero(u1) ? 1.0 : (u1 / sinh(u1))^2
        term2 = iszero(u2) ? 1.0 : (u2 / cosh(u2))^2
        return data.A + data.B * term1 + data.D * term2
    elseif data.model == :Einstein || data.model == :PlanckEinstein
        # Eq 3: A + Σ B_i * [ (C_i/T)^2 * exp(C_i/T) / (exp(C_i/T)-1)^2 ]
        cp = data.A
        u1 = data.C / T; e1 = exp(u1)
        cp += data.B * (u1^2 * e1) / (e1 - 1)^2
        u2 = data.E / T; e2 = exp(u2)
        cp += data.D * (u2^2 * e2) / (e2 - 1)^2
        if hasproperty(data, :F) && data.F != 0
            u3 = data.G / T; e3 = exp(u3)
            cp += data.F * (u3^2 * e3) / (e3 - 1)^2
        end
        return cp
    else
        error("Unknown Cp model type: $(data.model)")
    end
end

# Linear regression to convert complex forms to Reid polynomials
function regress_cp_to_reid(cp_raw_list)
    a_v, b_v, c_v, d_v = Float64[], Float64[], Float64[], Float64[]
    for data in cp_raw_list
        # Sample using the specific Tmin/Tmax from the handbook
        T_samples = collect(range(data.Tmin, data.Tmax, length=100))
        y = [evaluate_cp_ideal(data, T) for T in T_samples]
        p_fit = Polynomials.fit(T_samples, y, 3) 
        cf = Polynomials.coeffs(p_fit)
        # Ensure exactly 4 coefficients [a, b, c, d]
        push!(a_v, cf[1])
        push!(b_v, length(cf) >= 2 ? cf[2] : 0.0)
        push!(c_v, length(cf) >= 3 ? cf[3] : 0.0)
        push!(d_v, length(cf) >= 4 ? cf[4] : 0.0)
    end
    return Dict("a" => a_v, "b" => b_v, "c" => c_v, "d" => d_v)
end

function set_model(model_type, list::Vector{String}, cp_raw_list; phase=:liquid)
    # 1. Regress Cp data once
    reid_coeffs = regress_cp_to_reid(cp_raw_list)
    # 2. Instantiate base model with ReidIdeal
    model = model_type(list; idealmodel=ReidIdeal, ideal_userlocations=reid_coeffs)
    Mw = model.params.Mw.values
    STATE[] = (substances=list, model=model, Mw=Mw, phase=phase)
end

function set_absolute_reference(Hf0_vec, S0_vec)
    st = STATE[]
    my_ref = ReferenceState(:volume; 
        H0=Hf0_vec, S0=S0_vec, T0=298.15, P0=101325.0, phase=:gas)
    # Developer fix: returns a new wrapped model object
    new_model = Clapeyron.set_reference_state!(st.model, my_ref)
    STATE[] = (substances=st.substances, model=new_model, Mw=st.Mw, phase=st.phase)
end

function molar_mass()
    st = STATE[]
    return st.Mw ./ 1000.0 # kg/mol
end

function volume(p, T, m)
    st = STATE[]
    n = m ./ molar_mass()
    return Clapeyron.volume(st.model, p, T, n; phase=st.phase)
end

function specific_enthalpy_vt(V, T, m)
    st = STATE[]
    n = m ./ molar_mass()
    H_total = Clapeyron.VT_enthalpy(st.model, V, T, n)
    return H_total / sum(m)
end

function specific_enthalpy(p, T, m)
    v = volume(p, T, m)
    return specific_enthalpy_vt(v, T, m)
end

# Utility for checking the energy gap
function reaction_enthalpy(p,T, nu)
    st = STATE[]
    ns = length(st.substances)
    H_molar = [Clapeyron.enthalpy(st.model, p, T, [j==i ? 1.0 : 0.0 for j in 1:ns]) for i in 1:ns]
    return sum(nu[j] * H_molar[j] for j in 1:ns)
end

# Helpers for future @register_derivative use
function dV_dT(p, T, m)
    ForwardDiff.derivative(t -> volume(p, t, m), T)
end

function dV_dp(p, T, m)
    ForwardDiff.derivative(P -> volume(P, T, m), p)
end

function dV_dm(p, T, m)
    ForwardDiff.gradient(M -> volume(p, T, M), m)
end

end # module

HERE is how I create the model:

# 1. Fetch the process-specific data from the chemistry module
pkg = Glycol_Chemistry.get_thermo_package()

# 2. Configure the EoS with manual Cp regression (set to liquid phase)
ThermoEoS.set_model(ThermoEoS.Clapeyron.PCPSAFT, pkg.names, pkg.cp_raw, phase=:liquid)

# 3. Apply the ReferenceState (Enthalpy of Formation and Absolute Entropy)
ThermoEoS.set_absolute_reference(pkg.Hf0, pkg.S0)

The function for computing reaction_enthalpy is a utility function to check results; in the MTK code, I will only use functions volume and specific_enthalpy...

As an example, here is a plot of reaction enthalpy: