Skip to content

CPA KG giving same results of CPA CS for mixtures ? #185

Open
Open
CPA KG giving same results of CPA CS for mixtures ?#185
@Emanukka

Description

@Emanukka
Emanukka
opened on Feb 25, 2026

Hello!

I was validating my CPA implementation (https://github.com/Emanukka/reos) using teqp as a reference.
After validate cubics for pure and binary ( and everything was ok), I went to CPA
The CPA KG and CPA CS for pure water gave same results of teqp.

But when I try for a mixture of water and methanol, the KG model didn't match, only the CS.
To have other opnion, I simulated the same mixture in Clapeyron.jl, and the results matched my implemantion.

After some investigation, I've seen that teqp CPA KG and CPA CS gives the same results. This is strange, because I saw the issue ClapeyronThermo/Clapeyron.jl#239, and teqp results differ based on the radia dist. model (what is the expected):

reos CPA KG: {'A/R/T': -1.4579443886340535, 'Z': -0.9277464448546519, 'd(A/RT)_dT': np.float64(0.018614592198709968), 'd(A/RT)_dn': array([-2.60264098138663, -2.06026561164181]), 'X': array([0.59740598865315, 0.49566875655983, 0.59924543534437,
       0.50370256786259]), 'K': array([[0.              , 0.30206956317362, 0.              ,
        0.2082364510948 ],
       [0.30206956317362, 0.              , 0.2082364510948 ,
        0.              ],
       [0.              , 0.2082364510948 , 0.              ,
        0.12124917722626],
       [0.2082364510948 , 0.              , 0.12124917722626,
        0.              ]])}
teqp CPA KG : {'A/R/T': -1.4605958993455488, 'Z': -0.9313873373811081, 'd(A/RT)_dT': 0.018660061178390305, 'd(A/RT)_dn': array([-2.60853530634573, -2.06715513229804]), 'X': array([0.49505113890062, 0.59691154155244, 0.50308496586131,
       0.59875187893176]), 'K': array([[0.              , 0.30306714083355, 0.              ,
        0.20892414709904],
       [0.30306714083355, 0.              , 0.20892414709904,
        0.              ],
       [0.              , 0.20892414709904, 0.              ,
        0.1216495997952 ],
       [0.20892414709904, 0.              , 0.1216495997952 ,
        0.              ]])}
teqp CPA CS : {'A/R/T': -1.4605958993455488, 'Z': -0.9313873373811081, 'd(A/RT)_dT': 0.018660061178390305, 'd(A/RT)_dn': array([-2.60853530634573, -2.06715513229804]), 'X': array([0.49505113890062, 0.59691154155244, 0.50308496586131,
       0.59875187893176]), 'K': array([[0.              , 0.30306714083355, 0.              ,
        0.20892414709904],
       [0.30306714083355, 0.              , 0.20892414709904,
        0.              ],
       [0.              , 0.20892414709904, 0.              ,
        0.1216495997952 ],
       [0.20892414709904, 0.              , 0.1216495997952 ,
        0.              ]])}

@Comparing to Clapeyron.jl@

reos KG error : {'A/R/T': 6.834300182607953e-12, 'Z': 2.856637825256359e-11, 'd(A/RT)_dn': array([1.58373853653990e-11, 2.08394232495518e-01]), 'X': array([0.2052524611685 , 0.17029831322543, 0.18968112074682,
       0.15943862386366])}
teqp KG error : {'A/R/T': 0.0018186638270624462, 'Z': 0.003924447863158396, 'd(A/RT)_dn': array([0.00226474760059, 0.20574710570954]), 'X': array([0.00124602883642, 0.00082765692437, 0.00122612419599,
       0.00082362999508])}
teqp CS error : {'A/R/T': 3.2274546058031163e-13, 'Z': 4.5639679069040585e-12, 'd(A/RT)_dn': array([1.70176563486407e-12, 2.07541823462066e-01]), 'X': array([9.11645636512638e-12, 8.13540226922264e-12, 8.97275033479807e-12,
       8.10223850490108e-12])}

I don't kwon if I instantiating the teqp CS model in the wrong way :(

I will disponibilize the python and julia code used for validation!

Python

Imports

!pip install reos
!pip install teqp
!pip install si_units
from reos.eos import EquationOfState
from reos.state import State
from reos.consts import Consts
import teqp
import numpy as np
np.set_printoptions(precision=14)

from reos.cpa import CPAParameters, CPAPureRecord

R = Consts.ideal_gas_const()
from si_units import NAV, MOL

Functions used to calculate properties and compare

properties = ["A/R/T", "Z", "d(A/RT)_dn", "X","K"]

def get_properties(eos, t, rho, x = np.array([1.])):

    Z = eos.compressibility(t,rho,x)
    S = eos.entropy(t,rho,x) -  R * np.log(Z)
    A = eos.helmholtz(t, rho, x)

    dF_dT = (-S - A / t) / R / t

    X = eos.unbonded_sites_fraction(t, rho, x)
    K = eos.association_constants(t, rho, x)

    data = {
        "A/R/T" : A / R / t,
        "Z" : Z - 1 ,
        "d(A/RT)_dT" : dF_dT,
        "d(A/RT)_dn" : eos.chem_pot(t, rho, x) / R / t + np.log(eos.compressibility(t,rho,x)),
        "X" : X,
        "K" : K
    }

    return data

def teqp_get_properties(model, t, rho, x = np.array([1.])):

    alpha = model.get_Ar00(t, rho, x)

    rho_dalpha_drho = model.get_Ar01(t,rho,x) # =

    dalpha_dx = np.zeros(len(x))

    for j in range(len(x)):

        dalpha_dx[j] = (model.get_ATrhoXi(t, 0, rho, 0, x, j, 1))

    sum_j = dalpha_dx @ x


    dF_dn = np.zeros_like(dalpha_dx)

    for i in range(len(x)):

        dF_dn[i] = alpha + rho_dalpha_drho + dalpha_dx[i] - sum_j

    teqp_assoc_calcs = model.get_assoc_calcs(t, rho, x)
    X,K = get_X_and_K_from_teqp(teqp_assoc_calcs)
    teqpData = {
        "A/R/T": model.get_Ar00(t, rho, x),
        "Z": model.get_Ar01(t,rho,x),
        "d(A/RT)_dT": - model.get_Ar10(t, rho, x) / t,
        "d(A/RT)_dn": dF_dn,
        "X" : X,
        "K" : K

    }

    return teqpData

def comparison(eos, model_teqp, t, rho, x = np.array([1.])):

    err = {}

    reosData = get_properties(eos, t, rho ,x)
    teqpData = teqp_get_properties(model_teqp, t , rho, x)

    for prop in properties:

        err[prop] = abs((reosData[prop] - teqpData[prop]) / teqpData[prop])

    return err

def comparison_from_julia(reos_model, teqp_model, juliaData, t, rho, x = np.array([1.]), K_comparison = False):

    reos_err = {}
    teqp_err = {}

    reosData = get_properties(reos_model, t, rho ,x)
    teqpData = teqp_get_properties(teqp_model, t , rho, x)

    for prop in properties:

        if prop == "K" and K_comparison == False :
            
            continue
        else:
            reos_err[prop] = abs((reosData[prop] - juliaData[prop]) / juliaData[prop])
            teqp_err[prop] = abs((teqpData[prop] - juliaData[prop]) / juliaData[prop])

            
    return reos_err, teqp_err

def get_X_and_K_from_teqp(teqp_assoc_calcs):
    delta_non_phys = NAV * MOL  * np.array(teqp_assoc_calcs["Delta"])
    X_A = np.array(teqp_assoc_calcs["X_A"])

    map = teqp_assoc_calcs['to_CompSite']
    S = len(X_A)
    M = np.array(teqp_assoc_calcs["counts"])
    D = np.array(teqp_assoc_calcs["D"])

    m = np.zeros(S)
    mm = np.zeros((S,S))
    for j in range(S):
        i = map[j][1][0]
        m[j] = x[i]

        for l in range(S):

            k =map[l][1][0]
            mm[(j,l)] = x[i] * x[k]

    zerofy = D * 1
    for j in range(S):
        for l in range(S):

            if zerofy[(j,l)] != 1 and zerofy[(j,l)] != 0:
                zerofy[(j,l)] = 1

    delta_phys = delta_non_phys * zerofy

    K = rho * mm * delta_phys
    dQ_dX = m * (1 / X_A -1) - K@(X_A*M) # see https://doi.org/10.1021/ie060029x

    # print(dQ_dX,np.linalg.norm(dQ_dX))
    assert(np.linalg.norm(dQ_dX) < len(dQ_dX) * 1e-10 )

    return X_A, K

Models instantiating

_water = {
    "a0":0.12277,
    "b":0.000014515,
    "c1": 0.67359,
    "tc": 647.29,
    "epsilon": 16655.0,
    "kappa":  0.0692,
    "na": 2, # electron donor
    "nb": 2  # electron acceptor
}

_methanol = {
    "a0": 0.45897,
    "b":0.0000334,
    "c1": 1.0068,
    "tc": 512.64,
    "epsilon": 16070,
    "kappa": 0.0344,
    "na": 2, # electron donor
    "nb": 1  # electron acceptor
}

water = CPAPureRecord.new("water", molar_weight = 1.0, **_water)
methanol = CPAPureRecord.new("methanol", molar_weight = 2.0, **_methanol)

reos_parameters = CPAParameters.from_records([water,methanol])

reos_cpa_KG = EquationOfState.scpa(reos_parameters)
teqp_water = {
    "a0i / Pa m^6/mol^2": 0.12277,
    "bi / m^3/mol": 0.000014515,
    "c1": 0.67359,
    "Tc / K": 647.29,
    "epsABi / J/mol": 16655.0,
    "betaABi": 0.0692,
#   here the site types are declared. "e" for electron donor, "H" for electron acceptor
    "sites": ["e","e","H","H"]
}

teqp_methanol = {
    "a0i / Pa m^6/mol^2": 0.45897,
    "bi / m^3/mol": 0.0000334,
    "c1":1.0068,
    "Tc / K": 512.64,
    "epsABi / J/mol": 16070,
    "betaABi":0.0344,
#   here the site types are declared. "e" for electron donor, "H" for electron acceptor
    "sites": ["e","e","H"]
}

jCPA_KG = {
    "cubic": "SRK",
    "radial_dist": "KG",
#     "combining": "CR1", # No other option is implemented yet
    "pures": [teqp_water, teqp_methanol],
    "options": {"self_association_mask": [True, True]},
    "R_gas / J/mol/K": 8.31446261815324
}

jCPA_CS = {
    "cubic": "SRK",
    "radial_dist": "CS",
#     "combining": "CR1", # No other option is implemented yet
    "pures": [teqp_water, teqp_methanol],
    "options": {"self_association_mask": [True, True]},
    "R_gas / J/mol/K": 8.31446261815324
}

teqp_cpa_KG = teqp.make_model({"kind": "CPA", "model": jCPA_KG, "validate": False}, False)

teqp_cpa_CS = teqp.make_model({"kind": "CPA", "model": jCPA_CS, "validate": False}, False)

Comparison

t = 298.15 # K
rho = 1_000. # mol/m^3
x = np.array([0.6, 0.4])




reos_result_KG = get_properties(reos_cpa_KG, t, rho, x)
print(f"reos CPA KG: {reos_result_KG}")

teqp_result_KG = teqp_get_properties(teqp_cpa_KG, t, rho, x)
print(f"teqp CPA KG : {teqp_result_KG}")

teqp_result_CS = teqp_get_properties(teqp_cpa_CS, t, rho, x)
print(f"teqp CPA CS : {teqp_result_CS}\n")

juliaData_KG={
'Z':-0.9277464448281495,
'd(A/RT)_dT':0.018638018571771893,
'd(A/RT)_dn':-2.6026409814278524,
'X':[0.4956687564644897, 0.5974059887556956, 0.5037025677672323, 0.5992454354468074],
'A/R/T':-1.4579443886240895,
'K':[0.8390821199267238, 0.8676518795616637, 0.8676518795616637, 0.757807357664115],
}

juliaData_CS={
'Z':-0.9313873373768573,
'd(A/RT)_dT':0.018660061178296765,
'd(A/RT)_dn':-2.6085353063412935,
'X':[0.49505113889610836, 0.596911541557299, 0.5030849658567943, 0.5987518789366115],
'A/R/T':-1.4605958993450774,
'K':[0.8418531689820968, 0.8705172795793189, 0.8705172795793189, 0.7603099987200264],
}

reos_err, teqp_err_KG = comparison_from_julia(reos_cpa_KG,  teqp_cpa_KG, juliaData_KG,t, rho , x, K_comparison=False)
_, teqp_err_CS = comparison_from_julia(reos_cpa_KG,  teqp_cpa_CS, juliaData_CS, t, rho , x, K_comparison=False) # currently, CS in reos just in the Rust Package, but the CS match Clapeyron.jl calculations too

print("@Comparing to Clapeyron.jl@\n")

print(f"reos KG error : {reos_err}")
print(f"teqp KG error : {teqp_err_KG}")
print(f"teqp CS error : {teqp_err_CS}")

Julia

using Clapeyron
using StaticArrays
using Printf

R = Clapeyron.Rgas()
N_A = Clapeyron.N_A

function get_model(radial_dist)

    model = CPA(["Water","Methanol"], radial_dist = radial_dist; userlocations=(;
                a = [0.12277,0.45897],
                b = [0.000014515,0.0000334],
                c1 = [0.67359,1.0068],
                Mw = [1.,1.],
                Tc = [647.29, 512.64],
                Pc = [1e6,1e6],
                n_H = [2,1],
                n_e = [2,2],
                epsilon_assoc = Dict(
                    (("Water","H"),("Water","e"))=>16655.0 / R,
                    (("Methanol","H"),("Methanol","e"))=>16069.999999999998 / R,
                    (("Water","H"),("Methanol","e"))=> 0.5 * (16069.999999999998 + 16655.0 )/ R,
                    (("Water","e"),("Methanol","H"))=> 0.5 * (16069.999999999998 + 16655.0 )/ R,
                    
                    ),
                bondvol = Dict(
                    (("Water","H"),("Water","e"))=>0.0692, 
                    (("Methanol","H"),("Methanol","e"))=>0.0344,
                    (("Water","H"),("Methanol","e"))=>sqrt(0.0344 * 0.0692),
                    (("Water","e"),("Methanol","H"))=>sqrt(0.0344 * 0.0692),
                    
                    )
            ),
    )
    return model
end

function calc_properties(model, V, T, z = [1.]
    )
    A = Clapeyron.VT_helmholtz_energy_res(model, V, T, z)
    S = Clapeyron.VT_entropy_res(model, V, T, z)
    Z = Clapeyron.VT_compressibility_factor(model, V, T, z)
    dA_dni = Clapeyron.VT_chemical_potential_res(model, V, T, z)

    X = Clapeyron.assoc_fractions(model, V, T, z)
    Δ = Clapeyron.Δ(model, V, T, z)

    F = A / R / T
    dF_dT = (-S - A / T) / R / T

    d = Dict(
        "d(A/RT)_dn" => dA_dni[1] / R / T,
        "d(A/RT)_dT" => dF_dT,
        "Z" => -1 + Z,
        "A/R/T" => F, 
        "X" => X.v,
        "K" => N_A * Δ.values / V # this is wrong for mixtures, because I need the mjml matrix

    )

    return d

end

# println(model.sites)
# println(model.params.epsilon_assoc)
# println(model.params.bondvol)


V = 1/ 1000
T = 298.15
z = [0.6, 0.4]
# z = [1.]

# # println("R=",R)

# a_res =  Clapeyron.eos_res(mixture, V, T, z)  / R / T
# println("ARes/R/T = ",a_res) #ok

A = Clapeyron.VT_helmholtz_energy_res(model, V, T, z)
S = Clapeyron.VT_entropy_res(model, V, T, z)
Z = Clapeyron.VT_compressibility_factor(model, V, T, z)
dA_dni = Clapeyron.VT_chemical_potential_res(model, V, T, z)

F = A / R / T
dF_dT = (-S - A / T) / R / T

dF_dn = dA_dni / R / T


F
#%%

models = [get_model(:KG), get_model(:CS)]

model_name = ["KG", "CS"]
for i in 1:2
    
    d = calc_properties(models[i], V, T, z)
    
    @printf("juliaData_%s={\n", model_name[i] )
    
    for (k, v) in d
        @printf("'%s':%s,\n", k, v)
    end
    println("}")
end
#%%


# parse( dF_dn)
# partr
d = Dict(
    "d(A/RT)_dn" => dF_dn,
    "d(A/RT)_dT" => dF_dT,
    "Z" => -1 + Z,
    "A/R/T" => F, 

)

#%%
using Printf

println("{")
for (k, v) in d

    # if v == 
    @printf("'%s':%.17f,\n", k, v)

end
println("}")

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions