Completing activity coefficient partial derivatives

[marcosfelt]

Hi! Thanks for creating this amazing package!

I want to get access to the partial derivatives of the activity coefficient with respect to the mole fraction. I see that this has been commented out as incorrect in the current implementation. Any pointers on how I might be able to contribute a fix?

https://github.com/CalebBell/thermo/blob/6df9cbea2c4d8ed63325041a733a87881f515102/thermo/activity.py#L938

[CalebBell]

Hi Kobi,

This is a tough problem. I thought I had the right formula but it turned out not to be correct, and after some more work I'm not less sure there is a proper general formula. My original hope and anticipation was that this could be calculated as a combination of the properties GE, dGE_dxs, and d2GE_dxixjs. I would be greatly appreciative if you could find the proper formula. I haven't found it published anywhere before.

You can check your work by taking the numerical derivative and comparing to your implementation. Right now, there is only one activity model with the correct analytical expression for dgammas_dxs - UNIFAC, which came from the shape of its equation. A sample check is shown below.

from thermo.unifac import UNIFAC, DOUFIP2006, DOUFSG
from fluids.numerics import *
import numpy as np

T = 373.15
xs = [0.2, 0.3, 0.1, 0.4]
chemgroups = [{9:6}, {78:6}, {1:1, 18:1}, {1:1, 2:1, 14:1}]
# m = Mixture(['benzene', 'cyclohexane', 'acetone', 'ethanol'], zs=xs, T=T, P=1e5)
GE = UNIFAC.from_subgroups(T=T, xs=xs, chemgroups=chemgroups, version=1,
                           interaction_data=DOUFIP2006, subgroups=DOUFSG)

def to_jac_gammas(xs):
    return GE.to_T_xs(T, xs).gammas()

dgammas_dxs_num = jacobian(to_jac_gammas, xs, scalar=False, perturbation=1e-6)
dgammas_dxs_num/np.array(GE.dgammas_dxs())

Sincerely, Caleb

[marcosfelt]

Hi Caleb,

I've been thinking about this for the past couple days. I stumbled across this resource from the people who created the CocoSimulator. There is also a more explicit version for difference activity coefficients from this older paper as well (can send individually if anyone doesn't have access).

I am still trying to wrap my head around their derivation, but I think the core idea is to use the relationship between constrained and unconstrained derivatives (equation 19 in the first link) to then find molar derivatives (equation 26).

[marcosfelt]

Spent some time reading in depth today and it seems like this is the correct equation:.

From what I understand you can arbitrarily choose m. This is the component which will be determined by all the others. In the paper, they choose the nth component. From the third page in this paper:

I am going to try putting this in code and cross-checking with UNIFAC as suggested.

[marcosfelt]

I gave this a try, and it didn't seem to work, so I'm going to close the issue for now. Code below for reference.

from thermo.unifac import UNIFAC, DOUFIP2006, DOUFSG
from thermo import GibbsExcess
from fluids.numerics import *
from fluids.constants import R_inv
import numpy as np

T = 373.15
xs = [0.2, 0.3, 0.1, 0.4]
chemgroups = [{9: 6}, {78: 6}, {1: 1, 18: 1}, {1: 1, 2: 1, 14: 1}]
# m = Mixture(['benzene', 'cyclohexane', 'acetone', 'ethanol'], zs=xs, T=T, P=1e5)
GE = UNIFAC.from_subgroups(
    T=T,
    xs=xs,
    chemgroups=chemgroups,
    version=1,
    interaction_data=DOUFIP2006,
    subgroups=DOUFSG,
)

def to_jac_gammas(xs):
    return GE.to_T_xs(T, xs).gammas()

def dgamma_dxs_constrained(ge: GibbsExcess, xs, N, T):
    q2 = ge.d2GE_dxixjs()
    gammas = ge.gammas()
    matrix = np.zeros((N, N))
    RT_inv = R_inv / T
    for i in range(N):
        for l in range(N):
            matrix[i, l] += q2[i][l]
            matrix[i, l] -= q2[i][N - 2]
            for k in range(N):
                matrix[i, l] += xs[k] * (q2[k][l] - q2[k][N - 2])
            matrix[i, l] *= gammas[i]

    return matrix * RT_inv

def dgamma_dxs_unconstrained(ge: GibbsExcess, xs, N, T):
    q2 = ge.d2GE_dxixjs()
    gammas = ge.gammas()
    matrix = np.zeros((N, N))
    RT_inv = R_inv / T
    for i in range(N):
        for l in range(N):
            matrix[i, l] += q2[i][l]
            for k in range(N):
                matrix[i, l] += xs[k] * q2[k][l]
            matrix[i, l] *= gammas[i]
    return matrix * RT_inv

# print(GE.d2GE_dxixjs())
# print(GE.gammas())
dgammas_dxs_num = jacobian(to_jac_gammas, xs, scalar=False, perturbation=1e-6)
unifac = np.array(GE.dgammas_dxs())
unconstrained = dgamma_dxs_unconstrained(GE, xs, 4, T)
constrained = dgamma_dxs_constrained(GE, xs, 4, T)
print("UNIFAC")
print(unifac)
print("My formula unconstrained")
print(unconstrained)
print("My formula constrained")
print(constrained)

print("Ratio unifac")
print(unifac / dgammas_dxs_num)
print("Ratio constrained")
print(constrained / dgammas_dxs_num)

print("Ratio unconstrained")
print(unconstrained / dgammas_dxs_num)

[CalebBell]

Hi Kobi, This is unfortunate. Not every formula published in the literature is accurate, for lots of reasons. Which activity model were you trying to obtain the derivatives with? I may still be able to offer help if you are more specific. With SymPy, I have always been able to obtain any derivatives for a specific model, although it can take some time.

Caleb

[marcosfelt]

Hi Caleb,

I was working with UNIFAC. I wrote out all the math here that I used to get the general formula I noted above: https://kobifelton.com/notes/composition-derivatives-of-activity-coefficients

Working on the specific activity model is actually a good idea. I'm using NRTL, so I'll try to work straight from that since there are some derivations in the literature. If I get it working, I'll submit a PR.