Ask for help with sCPA calculating the molar volume of water

[Benfzh]

I plan to realize the calculation of TP flash by sCPA state equation through matlab programming. However, I’ve been unable to obtain accurate results for the water component. Specifically, the molar volume obtained from Newton’s iteration for the state equation is inaccurate (no solution within the given interval). I’ve tried other compounds like triethylene glycol (associating compounds) and methane (non-associating compounds), both of which yield correct results. Therefore, I believe there are no issues with subsequent calculations involving fugacity coefficients and solving the nonlinear equation system.

Below is the part of my MATLAB program where I solve the water state equation (eq1 and eq2):

RDF = v / (v - 0.475 * bm_liq);
delta = RDF * (exp(epsilon_ij / R / T) - 1) * bij * beta_ij;
eq1 = -XA1 + v / (v + 2 * x(1) * XA1 * delta);
eq2 = -P + R * T / (v - bm_liq) - am_liq / (v * (v + bm_liq)) - R * T / v * (v / (v - 0.475 * bm_liq)) * (2 * x(1) * (1 - XA1));

The water properties in my code are as follows:

R = 8.314;
Tc = 647.29 ; %K
omega = 0.345; %ω
a0 = 0.12277; %Pa*m6/mol2
b = 1e-5 * 1.4515; %m3/mol
c1 = 0.67359;%
epsilon = 16655 ; %Pa*m3/mol
beta = 1e-3 * 69.2;

At 298.15 K and 1e5 Pa:

am_liq = 0.181664465033266;
bm_liq = 1.4515e-05;
epsilon_ij = 16655;
bij = 1.4515e-05;
beta_ij = 0.0692;

Using the Broyden method to solve eq1 and eq2, I got the following results:

XA1 = 0.551057164203367; % The correct value should be 0.1 to 0.2?
v = 1.130522838228e-3. % m^3/mol

Here are the correct results using Clapeyron's calculations:

using Clapeyron
model = sCPA(["water"])
V = volume(model, 1e5, 298.15, [1]; phase=:liquid)

v = 1.792601107455581e-5

This problem has bothered me for many days, forgive me for not reading the Clapeyron source code completely, my thermodynamics and programming foundation are not good, but I still hope that someone can help me solve this problem, if you are willing to take a short time for me, I would be very grateful.

[longemen3000]

hello,

at first glance, your delta is ok:

function delta_scpa(v,T)
    R = 8.314
    Tc = 647.29  #K
    omega = 0.345 #ω
    a0 = 0.12277 #Pa*m6/mol2
    b = 1e-5 * 1.4515 #m3/mol
    c1 = 0.67359
    bm_liq = b
    am_liq = a0*(1+c1*(1-sqrt(T/Tc)))^2
    c1 = 0.67359 # no units
    epsilon = 16655  #Pa*m3/mol
    beta = 1e-3 * 69.2;
    RDF = v / (v - 0.475 * bm_liq)
    delta = RDF * (exp(epsilon / R / T) - 1) * bm_liq * beta
    return delta
end

in the julia REPL:

julia> model = sCPA("water")
CPA{BasicIdeal, RK{BasicIdeal, Clapeyron.sCPAAlpha, NoTranslation, vdW1fRule}} with 1 component:
 "water"
RDF: Kontogeorgis (s-CPA)
Contains parameters: Mw, Tc, a, b, c1, epsilon_assoc, bondvol

julia> Clapeyron.delta_scpa(1.8e-5,373.15)
0.00034764314610167543

julia> Clapeyron.Δ(model,1.8e-5,373.15,[1.0]).values[1]*Clapeyron.N_A
0.00034764545404200635

So the problem is upstream. May i suggest using an exact solver for single components? for an associating model with just one association pair, the exact solution for the two sites that interact is the following:

function X_assoc_scpa(V,T)
    rho = 1/V
    #for water, the association scheme has two acceptor sites and two donor sites
    na = 2 
    nb = 2
    #for water, the  molar fraction of the donor site and the molar fraction of the acceptor site are the same
    zi = 1.0 
    zj = 1.0
    delta = delta_scpa(V,T)
    kia = na*zi*rho*delta
    kjb = nb*zj*rho*delta
    
    #solving a quadratic eq a*xx + b*x + c = 0
    a = kia
    b = 1 - kia + kjb
    c = -1
    denom = b + sqrt(b*b - 4*a*c)
    xia = -2*c/denom #equivalent to (-b + denom)/2a
    xk_ia = kia*xia
    xjb = (1- xk_ia)/(1 - xk_ia*xk_ia)
    #in this case, because na == nb, the nonbonding fractions have the same value. this is not always true.
    return xia,xjb 
end

what is x(1) in the code?

[Benfzh]

Thanks for the reply, for a single component, this x(1) I have always treated as 1, I checked this issue and found that when the correct v and xA1 values are entered, x(1) = 0.999935848398357 is not equal to 1.

Oddly enough, my program was able to get the correct results (I think the error is within the acceptable range) for Triethylene glycol (TEG), which is also a 4C-scheme.
In the calculation of one-component TEG (298.15 K, 1e5 Pa):

XA1 = 0.198240312753122  
v = 1.416272462825725e-04  # calc by my code
v = 1.416282e-04  # calc by aspen

In the TEG-CH4 system (298.15 K, 1e5 Pa, [0.5, 0.5]):

XA1 = 0.198257274068490  
v = 1.415728556685252e-04  # calc by my code
v = 1.415738e-04  # calc by aspen

and the calculation method you proposed for a single component looks cool (I don't fully understand it yet), but my main purpose is not to calculate a single component, to be precise, my code mainly wants to implement the flash calculation of the H2O-TEG-CH4 system, and the calculation implemented in the code is as follows:
Associating method:

The equation of state of the system is as follows:

site monomer fraction:

Associative Strength:

Enter my code to do the flash calculation:

z = [0.2,0.3,0.5] #molefrac of H2O-TEG-CH4
T = 298.15
P = 1e5

In my code I get the following result:

x = [0.396201838289065,0.603453370328999,0.000344791381935410] # liquid phase molefrac of H2O-TEG-CH4
y = [0.00603117682244540,2.97421988749969e-07,0.993968525755566] # vapor phase molefrac of H2O-TEG-CH4
X_assoc = [0.134668461505471,0.134668461505456,0.152715081690327,0.152715081690311]
v_liq = 9.36612614080473e-05 # liquid phase molar volume
v_gas = 0.024743465301345 # vapor phase molar volume

The v and x,y calculated in ASPEN are as follows:

x = [0.394397219917985,0.605254559321403,0.00034822076061161] # liquid phase molefrac of H2O-TEG-CH4
y = [0.00894929291534629,3.42475201093806E-07,0.991050364609453] # vapor phase molefrac of H2O-TEG-CH4
v_liq = 9.390943e-5 # liquid phase molar volume
v_gas = 0.02474414 # vapor phase molar volume

This seems to be only a partial error, but this error is much larger than when there is no water component, and it is fatal, and when I want to increase the molar fraction z of water, the program will error, indicating that the solved nonlinear equations (4eq of XAi, 1eq of CPA) have gone wrong, However, I could never find the cause of the error, leading me to think for a moment that there was an error in each of the water parameters (a, b, epsilon, beta), but there wasn't, and here is a look at all of the parameters in my H2O-TEG-CH4 calculations above:

 a = [0.181664465033266,8.129142482426758,0.182896028299481]
a0 = [0.12277,3.9126,0.23203]
aij = [0.181664465033266,1.572309120001429,0.181616621626639;
1.572309120001429,8.129142482426758,1.015429033482966;
0.181616621626639,1.015429033482966,0.182896028299481]
am_gas = 0.182881141152842
am_liq = 3.741109585759145
b = [1.4515e-05,1.321e-04,2.91e-05]
beta = [0.0692,0.0188,0]
beta_ij = [0.0692,0.036068823102508,0;
0.036068823102508,0.0188,0;
0,0,0]
bij = [1.4515e-05,7.33075e-05,2.18075e-05;
7.33075e-05,1.321e-04,8.061e-05;
2.18075e-05,8.061e-05,2.91e-05]
bm_gas = 2.901206592050947e-05
bm_liq = 8.547709333244092e-05
c1 = [0.67359,1.1692,0.4472]
epsilon = [16655,14337,0]
epsilon_ij = [16655,15496,8.3275e+03;
15496,14337,7.1685e+03;
8.3275e+03,7.1685e+03,0]

Thanks again for taking time out of your busy schedule to answer my questions.

[longemen3000]

hello,

at a quick glance and for earlier comments, you are solving association and volume simultaneously, right? as far as i know, association problems are really stiff, so if you don't have good initial points, the problem quickly converges to wrong values.

for CPA in particular, i recommend a nested iterative procedure for the solution of the volume at a set pressure:

1. solve the association fractions at the input volume
2. update volume

for the association solver, i reccomend a successive substitution approach, you start with non bonded fractions equal to 1 and solve iteratively using the relations for monomer fraction (with a damping factor of 0.5, that is XNi_new = 0.5*(XNi_old + v/(v + (...))) where N could be any of sites A,B,C,D. i saw in your flash code that you initialize X_assoc with zeros, maybe that is one of the causes for instability.

for the volume, the problem form of p-p0 is also not ideal, as it can venture into the metastable phase of the p-v diagram. in Clapeyron we use the log(v) = log(v_old) + (pset - p(v_old))/(v_old*dpdV(v_old), if you don't have an expression for dpdV, you can use a secant approximation instead with a damping to make it more stable (when using the exact derivative and starting from a point where p(v0) > p_set in liquid or p(v0) < p_set in gas, that iteration scheme always seems to converge)
I don't know about using the pure SRK values as initial points, while for the gas it seems fine, they can cause trouble for the liquid phase. this is because the volume calculated with association will be lower than the volume calculated without, so it can enter in the metastable phase. one example with sCPA and water:

julia> volume(model,1e5,250,phase = :l) #normal model
1.7360114307076252e-5

julia> volume(model,1e5,250,phase = :l) #here i manually deactivated the association part of CPA, same parameters otherwise
2.542812539096336e-5

On Clapeyron, due to how we solve volumes, we start (in CPA) with v0 = 1.25*b, on the gas phase, the SRK initial point seems fine. alternatively, if you dont want to solve the roots of an SRK model, you could use a virial approximation, and the second virial coefficient of a cubic equation of state is just B = b - a/RT.

[Benfzh]

The equation you mentioned in your last reply for volume calculationslog(v) = log(v_old) + (pset - p(v_old))/(v_old*dpdV(v_old)) (there's a bracket missing from the right side in your reply, right?) , I made an attempt but was unsuccessful. Because the calculated v is less than b, resulting in a negative value for RDF, a series of subsequent errors are reported.The attempted code is below:

v_liq = 1.25 * bm_liq;
        log_v = log(v_liq);
        while true
            X_assoc = 0.5*ones(1,4);
            v_liq_old = v_liq;
            log_v_old = log(v_liq);
            RDF = v_liq / (v_liq-0.475*bm_liq); %KG
            delta_A1B1 = RDF * (exp(epsilon_ij(1,1)/R/T)-1) * bij(1,1)*beta_ij(1,1);
            delta_A1B2 = RDF * (exp(epsilon_ij(1,2)/R/T)-1) * bij(1,2)*beta_ij(1,2);
            delta_A2B1 = RDF * (exp(epsilon_ij(2,1)/R/T)-1) * bij(2,1)*beta_ij(2,1);
            delta_A2B2 = RDF * (exp(epsilon_ij(2,2)/R/T)-1) * bij(2,2)*beta_ij(2,2);
            while true
                X_assoc_old=X_assoc;
                X_assoc(1) = 0.5*( X_assoc_old(1) + v_liq/( v_liq + (2*x(1)*X_assoc(2)*delta_A1B1 + 2*x(2)*X_assoc(4)*delta_A1B2) ) );
                X_assoc(2) = 0.5*( X_assoc_old(2) + v_liq/( v_liq + (2*x(1)*X_assoc(1)*delta_A1B1 + 2*x(2)*X_assoc(3)*delta_A2B1) ) );
                X_assoc(3) = 0.5*( X_assoc_old(3) + v_liq/( v_liq + (2*x(1)*X_assoc(2)*delta_A2B1 + 2*x(2)*X_assoc(4)*delta_A2B2) ) );
                X_assoc(4) = 0.5*( X_assoc_old(4) + v_liq/( v_liq + (2*x(1)*X_assoc(1)*delta_A1B2 + 2*x(2)*X_assoc(3)*delta_A2B2) ) );
                delta_X = sum(abs(X_assoc_old-X_assoc),'all');
                if delta_X < 1e-10
                    break
                end
            end
            assoc = ( 2*x(1)*(1-X_assoc(1)) + 2*x(1)*(1-X_assoc(2)) + 2*x(2)*(1-X_assoc(3)) + 2*x(2)*(1-X_assoc(4)) );
            P_calc = @(v_liq)R*T/(v_liq-bm_liq) - am_liq/(v_liq*(v_liq+bm_liq)) - 0.5*R*T/v_liq * (v_liq/(v_liq-0.475*bm_liq)) * assoc;
            dpdv_v_old = (P_calc(v_liq_old+1e-8)-P_calc(v_liq_old-1e-8))/2e-8;
            log_v = log_v_old + (P-P_calc(v_liq_old))/(v_liq_old*dpdv_v_old);
            v_liq = exp(log_v);
            delta_v = abs(v_liq_old-v_liq);
            if delta_v < 1e-10
                break
            end
        end

In my code, the gas phase is calculated entirely by the SRK equation, and since I don't think there is any associative in the gas phase, the CPA equation is simplified to the SRK equation.

Thanks again for providing guidance on the issues I raised!

[longemen3000]

Some comments:

• I recommend you to use inf-norm (max(X - X_old) < tol) for the termination criteria of the assoc loop
• normally pressure gradients are really steep, so finite differencing introduces a lot of error maybe using a complex step derivative could help here (imag(F(x+i*h))/h, where h = 1e-8*v)
• while association effects, are less pronounced in the vapour phase, they still exist, using SRK is a good aproximación though, just use a real SRK and not the sCPA SRK constants.

[Benfzh]

一些评论：

• 我建议您将 inf-norm （） 作为 assoc 循环的终止条件max(X - X_old) < tol
• 通常压力梯度非常陡峭，因此有限微分会引入很多误差，也许使用复阶阶梯导数可能会有所帮助（，其中imag(F(x+i*h))/h``h = 1e-8*v)
• 虽然在气相中关联效应不太明显，但它们仍然存在，但使用 SRK 是一个很好的近似方法，只需使用真正的 SRK 而不是 sCPA SRK 常数即可。

I have no more questions and request to close this issues. Thank you so much for your patience and help! Your guidance was invaluable in resolving my issue. I really appreciate the time and effort you took to answer my questions. Thanks again!