Cubic Plus Association (CPA)

The combination of a cubic EOS with association with the association term. The sum of the terms goes like:

\[\alpha^{\rm r} = \alpha^{\rm r}_{\rm cub} + \alpha^{\rm r}_{\rm assoc}\]

Cubic part

The residual contribution to \(\alpha\) is expressed as the sum :

\[\alpha^{\rm r}_{\rm assoc} = \alpha^{\rm r}_{\rm cub,rep} +\alpha^{\rm r}_{\rm cub,att}\]

where the cubic parts come from two separate contributions. The repulsive part of the cubic EOS contribution is:

\[\alpha^{\rm r}_{\rm cub,rep} = -\ln(1 - b_{\rm mix}\rho)\]

The attractive part of the cubic EOS contribution:

\[\alpha^{\rm r}_{\rm cub,att} = -\frac{a_{\rm mix}}{RT}\dfrac{\ln\left(\frac{\Delta_1 b_{\rm mix}\rho + 1}{\Delta_2b_{\rm mix}\rho + 1}\right)}{b_{\rm mix}\cdot(\Delta_1 - \Delta_2)}\]

with the coefficients depending on the cubic type:

SRK: \(\Delta_1=1\), \(\Delta_2=0\)

PR: \(\Delta_1=1+\sqrt{2}\), \(\Delta_2=1-\sqrt{2}\)

The mixture models used for the \(a_{\rm mix}\) and \(b_{\rm mix}\) are the classical ones:

\[a_{\rm mix} = \sum_i\sum_jx_ix_j(1-k_{ij})a_{ij}(T)\]

with \(x\) the mole fraction, \(k_{ij}\) a weighting parameter

\[a_{ij}(T) = \sqrt{a_ia_j}\]

and

\[a_{i}(T) = a_{0i}\left[1+c_{1i}(1-\sqrt{T/T_{{\rm crit},i}})\right]^2\]

and for \(b\):

\[b_{\rm mix} = \sum_ix_ib_i\]

so there are three cubic parameters per fluid that need to be obtained though fitting: \(b_{i}\), \(a_{0i}\), \(c_{1i}\). The value of \(a_{ij}\) depends on temperature while \(b_{\rm mix}\) does not.

Association part

The residual contribution of the association to the Helmholtz energy \(\alpha^{\rm r}_{\rm assoc}\) is formulated as:

\[\alpha^{\rm r}_{\rm assoc} = \sum_{i} x_i \sum_{A_i} \left(\ln(Y_{A_i}) - \frac{Y_{A_i}}{2} \right) + \frac{M_i}{2},\]

where \(x_i\) is the mole fraction of molecule \(i\), \(A_i\) is a site of type \(A\) on molecule \(i\) that can interact with other site types (e.g., type \(B\)). \(Y_{A_i}\) then is the monomer mole fraction of site \(A_i\) and \(M_i\) is the count of association sites of molecule \(i\). The monomer mole fraction \(Y_{A_i}\) can be calculated with:

\[Y_{A_i} = \frac{1}{1+\rho N_{\rm A} \sum_j \sum_{B_j} x_j Y_{B_j} \Delta_{A_i B_j}}.\]

\(\rho\) is the molar density of the mixture, \(N_{\rm A}\) is the Avogadro constant, \(j\) is the index for a molecule that can be the same as molecule \(i\). \(B_j\) is a site of type \(B\) on molecule \(j\) and \(\Delta_{A_i B_j}\) is the interaction strength between the sites. The calculation of \(\Delta_{A_i B_j}\), \(Y_{A_i}\) and \(Y_{B_j}\) is explained in more detail in the documentation about association.

The important thing to understand is whether a site is an electron acceptor (positively charged) or electron donor (negatively charged) and how many sites are present on a molecule. For example, water is usually considered a molecule with two electron acceptor (labeled as “H” in teqp) and two electron donor (labeled as “e” in teqp) sites. The site types have to be specified.

[1]:

import teqp, numpy as np

[2]:

water = {
    "a0i / Pa m^6/mol^2": 0.12277,
    "bi / m^3/mol": 0.0000145,
    "c1": 0.6736,
    "Tc / K": 647.13,
    "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"]
}

jCPA = {
    "cubic": "SRK",
    "radial_dist": "CS",
#     "combining": "CR1", # No other option is implemented yet
    "pures": [water],
    "R_gas / J/mol/K": 8.31446261815324
}

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

T = 303.15 # K
rhomolar = 1000 # mol/m^3
molefracs = np.array([1.0])

# Note: passing data back and forth in JSON format is done for convenience and flexibility, not speed
res = model.get_assoc_calcs(T, rhomolar, molefracs)

print('D:', np.array(res['D']))
print('∆:', np.array(res['Delta']))
print('X_A:', np.array(res['X_A']))
print('siteid->(component, name):', res['to_CompSite'])
print('(component, name)->siteid:', res['to_siteid'])
print('multiplicities:', np.array(res['counts']))


print('pressure in Pa: ',T*rhomolar*model.get_R(molefracs) * ( 1+ model.get_Ar01(T,rhomolar,molefracs)))

D: [[0 2]
 [2 0]]
∆: [[1.24389806e-27 1.24389806e-27]
 [1.24389806e-27 1.24389806e-27]]
X_A: [0.54879025 0.54879025]
siteid->(component, name): [[0, [0, 'H']], [1, [0, 'e']]]
(component, name)->siteid: [[[0, 'H'], 0], [[0, 'e'], 1]]
multiplicities: [2 2]
pressure in Pa:  84414.03871039051