critical Cross-over model

[longemen3000]

mockup for a general, critical cross-over model. what is missing is to actually implement the functions, But the recombine! structure is there.

[codecov]

Codecov Report

Attention: Patch coverage is 19.88304% with 137 lines in your changes missing coverage. Please review.

Project coverage is 72.30%. Comparing base (5a76873) to head (d2816b3).

Files with missing lines	Patch %	Lines
src/models/CrossOver/variants/Kiselev2000.jl	0.00%	69 Missing ⚠️
src/models/CrossOver/crossover.jl	0.00%	47 Missing ⚠️
src/utils/macros.jl	73.68%	10 Missing ⚠️
...models/SAFT/SAFTVRMie/variants/SAFTVRMieKiselev.jl	0.00%	8 Missing ⚠️
src/database/params/AssocOptions.jl	50.00%	2 Missing ⚠️
src/utils/split_model.jl	66.66%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #389      +/-   ##
==========================================
- Coverage   72.68%   72.30%   -0.38%     
==========================================
  Files         356      357       +1     
  Lines       29556    29623      +67     
==========================================
- Hits        21482    21419      -63     
- Misses       8074     8204     +130     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[branch171]

I can see there is a base code for crossover I will add and test, then if it looks the business I will do a pull request so you can get the implementation

[branch171]

Here is my code that I have developed so far seems to work as expected at least for pure CO2:

Note there are many Kiselev/Ely/McCabe papers all with slightly different implementations.
I have used what looks like the consenus. Also, for a mixture I think you would need to calculate the Tc0 and vc0 at the composition not for the pure components. But I had trouble with x0_crit_mix not converging for the example notebook.

For now I am using a very simple (likely wrong) linear mixing to get Tc0 and vc0 from the pure values.

parameters:

Clapeyron Database File,,,,,,,,,,
SAFTVRMie Like Parameters,,,,,,,,,,
species,DIPPR Number,Mw,Tc,vc,Tc0,vc0,d1,v1,Gi,segment,sigma,epsilon,lambda_r,lambda_a,n_H,n_e,source
carbon dioxide,,44.01,304.1282,9.41178357551188e-5,304.3330128243223,9.487219622412847e-5,1.277,0.0065,0.005,1.6936,3.085,222.89,26.408,5.255,1.0,0.0,10.1016/j.fluid.2015.12.044

struct SAFTVRMieKiselevParam{T} <: ParametricEoSParam{T}
    Mw::SingleParam{T}
    Tc::SingleParam{T}
    vc::SingleParam{T}
    Tc0::SingleParam{T} # this can be removed if calculated in the recombine/initiasation from crit_pure or better x0_crit_mix
    vc0::SingleParam{T} # this can be removed if calculated in the recombine/initiasation from crit_pure or better x0_crit_mix
    d1::SingleParam{T}
    v1::SingleParam{T}
    Gi::SingleParam{T}
    segment::SingleParam{T}
    sigma::PairParam{T}
    lambda_a::PairParam{T}
    lambda_r::PairParam{T}
    epsilon::PairParam{T}
    epsilon_assoc::AssocParam{T}
    bondvol::AssocParam{T} 
end

function SAFTVRMieKiselevParam(Mw,Tc,vc,Tc0,vc0,d1,v1,Gi,segment,sigma,lambda_a,lambda_r,epsilon,epsilon_assoc,bondvol)
    return build_parametric_param(SAFTVRMieKiselevParam,Mw,Tc,vc,Tc0,vc0,d1,v1,Gi,segment,sigma,lambda_a,lambda_r,epsilon,epsilon_assoc,bondvol) 
end

# base a_res without critical correction
function a_res_base(model ::SAFTVRMieKiselevModel, V, T, z, _data=@f(data))
    return @f(a_hs,_data) + @f(a_dispchain,_data) + @f(a_assoc,_data)
end

function Z_base(model ::SAFTVRMieKiselevModel, V, T, z)
    ares(x)  = a_res_base(model, x, T, z,)
    dares(x) = Solvers.derivative(ares,x) 
    return 1.0 - V*dares(V)
end

function a_res(model ::SAFTVRMieKiselevModel, V, T, z)
    ∑z = sum(z)
    Tc = model.params.Tc.values
    vc = model.params.vc.values
    Tc0 = model.params.Tc0.values
    vc0 = model.params.vc0.values
    d1 = model.params.d1.values
    v1 = model.params.v1.values
    Gi = model.params.Gi.values
    Δτ = zero(T+V+first(z))
    Δϕ = zero(T+V+first(z))
    _Tc0 = zero(T+V+first(z))
    _vc0 = zero(T+V+first(z))
    _d1 = zero(T+V+first(z))
    _v1 = zero(T+V+first(z))
    invGi = zero(T+V+first(z))
    for i ∈ @comps
        zᵢ,Tcᵢ,vcᵢ,Tc0ᵢ,vc0ᵢ,d1ᵢ,v1ᵢ,Giᵢ = z[i],Tc[i],vc[i],Tc0[i],vc0[i],d1[i],v1[i],Gi[i]
        Δτ += zᵢ*(Tcᵢ/Tc0ᵢ - 1.0)
        Δϕ += zᵢ*(vcᵢ/vc0ᵢ - 1.0)
        _Tc0 += zᵢ*Tc0ᵢ
        _vc0 += zᵢ*vc0ᵢ
        _d1 += zᵢ*d1ᵢ
        _v1 += zᵢ*v1ᵢ
        invGi += zᵢ/Giᵢ
    end
    Δτ /= ∑z
    Δϕ /= ∑z
    _Tc0 /= ∑z
    _vc0 /= ∑z
    _Tc = _Tc0 + Δτ*_Tc0
    _vc = _vc0 + Δϕ*_vc0
    _d1 /= ∑z
    _v1 /= ∑z
    _Gi = 1.0/invGi/∑z
    b = sqrt(1.359)
    m = 1.3
    beta = 0.325
    alpha = 0.11
    gamma = 2.0-alpha-2.0*beta
    D = 0.51  
    v = V/∑z
    τ = T/_Tc - 1.0
    ϕ = v/_vc - 1.0
    fnc = 4.0*(b/m*abs(ϕ*(1.0 + _v1*(ϕ^2)*exp(-8.5*ϕ)) + _d1*τ*(1.0 - 2.0*τ)))^(1.0/beta) + 2.0*τ
    q = sqrt((fnc + sqrt(fnc^2 + 12.0*τ^2))/6.0/_Gi)
    Y = (q/(1.0+q))^2
    τY = 0.0
    if Y > 0.0
        τY = τ*Y^(-alpha/2.0)
    end
    τs = τY + (1.0 + τ)*Δτ*Y^(2.0*(2.0 - alpha)/3.0)
    Ts = _Tc0*(τs + 1.0)
    ϕs = ϕ*Y^((gamma - 2.0*beta)/4.0) + (1.0 + ϕ)*Δϕ*Y^((2.0 - alpha)/2.0)
    vs = _vc0*(ϕs + 1.0)
    Δv = v/_vc0 - 1.0
    _a_res = a_res_base(model, vs, Ts, z) - a_res_base(model, _vc0, Ts, z) + a_res_base(model, _vc0, T, z) + ϕs*Z_base(model, _vc0, Ts, z) - Δv*Z_base(model, _vc0, T, z) + log((Δv + 1.0)/(ϕs + 1.0))
    return _a_res
end

[branch171]

Hi

my parameters wil be

struct Kiselev2000Param{T} <: ParametricEoSParam{T}
    Tc::SingleParam{T}
    vc::SingleParam{T}
    d1::SingleParam{T}
    v1::SingleParam{T}
    Gi::SingleParam{T}
end

How do I set up a crossover_userlocations or csv file for these?

,I made a new SAFTVRMieBase as the paramters when we use a crossover are note the same. The crossover impacts the liquid density. It also when fitted improves the speed of sound.

My call so far is:

modelbase = SAFTVRMieCrossOverBase(["carbon dioxide"];idealmodel=AlyLeeIdeal)
modelcrossover = Kiselev2000()
model = CrossOver(model;modelcrossover)

But Kiselev2000() needs parameters defined

[longemen3000]

I think the linear mixing rule is ok for the moment.

In this particular example, Kiselev and Mie parameters are bonded, so i would just make a SAFTVRMieKiselev model that returns a crossover model.

[longemen3000]

may i ask, are you using vanilla SAFTVRMie (LJ kernel) or the SAFTVRMie-15 version (with the mie kernel) ?

[longemen3000]

with the last PR, the model seems to work:

julia> model = SAFTVRMieKiselev("carbon dioxide")
Critical cross-over Model with 1 component:
 Base Model: SAFTVRMie15{BasicIdeal, Float64}("carbon dioxide")
 Liquid Model: Clapeyron.Kiselev2000("carbon dioxide")

julia> model.params.Tc0
SingleParam{Float64}("Tc0") with 1 component:
 "carbon dioxide" => 304.5806694166402 (missing)

julia> model.params.Vc0
SingleParam{Float64}("Vc0") with 1 component:
 "carbon dioxide" => 9.886971185386159e-5 (missing)

julia> model.params.Pc0
SingleParam{Float64}("Pc0") with 1 component:
 "carbon dioxide" => 7.402481658212735e6 (missing)
 
 julia> crit_pure(model)
(304.12821704404536, 7.322324594719085e6, 9.415265587182492e-5)

julia> crit_pure(model.basemodel)
(304.5806694166402, 7.402481658212735e6, 9.886971185386159e-5)

[branch171]

I used SAFTVRMie15 only for convenience. I don't know which is the better version. But I think SAFTVRMiwe might be a better starting choice.

I get the same results as you:
SAFTVRMieKiselev Tc, pc, vc = 304.12821704404536, 7.322324594719085e6, 9.415265587182492e-5

Note however, vc is not quite equal to vc_actual 9.41178357551188e-5. But maybe its just a numerical convergence thing

I think it should be like Tc is matching. The only difference is pc which is expected to be different as Kiselev does not match this unless the base model equals pc at Tc,vc. I will check if there is an issue.

I would be good to now see if we can get Metaheuristics via estimation to optimise the paramters. Then we can start a collect of components and test the mixture prediction. First one to test would be CO2, N2

[branch171]

I switched from SAFTVRMie15 to SAFTVRMie in any case as there is no CO2 in the database so its not good for comparisons.

I also added in whats called the Kernal term:

K = τ^2/2.0/(1.0 + τ^2)(_a20(Y^(-alpha/2.0/Δ) - 1.0) + _a21*(Y^(-(alpha - Δ)/Δ) - 1.0))
return _a_res - K

its only needed for CP and speed of sound which should go infinite and to zero respectively at the critical point.

I am happy with all the changes. I just need to try estimation now.

[branch171]

Here is the CO2 predictions and update files. I will try merge my new feature branch back into my master and then you should be able to see what I have as part of the HELD PR. you can then decide on what you want to take for HELD and CrossOver.

Adding the Kernal term seemed to fix the vc issue I now get:

GERG2008 Tc, pc, vc = 304.1282, 7.3773e6, 9.41178357551188e-5
SAFTVRMie Tc, pc, vc = 305.1913635309112, 7.495803823291793e6, 0.0001048302010667033
SAFTVRMieKiselev Tc, pc, vc = 304.1282035996962, 7.322322213172413e6, 9.411730538838142e-5

Which is close enough




SAFTVRMieKiselev_like.csv

[branch171]

function a_res_crossover(model::CrossOver,V,T,z,critmodel::Kiselev2000Model)
    ∑z = sum(z)
    Tc = critmodel.params.Tc.values
    vc = critmodel.params.Vc.values
    Tc0 = model.params.Tc0.values
    vc0 = model.params.Vc0.values
    d1 = critmodel.params.d1.values
    v1 = critmodel.params.v1.values
    Gi = critmodel.params.Gi.values
    a20 = critmodel.params.a20.values
    a21 = critmodel.params.a21.values
    basemodel = model.basemodel
    _0 = zero(Base.promote_eltype(critmodel,basemodel,V,T,z))
    Δτ = _0
    Δϕ = _0
    _Tc0 = _0
    _vc0 = _0
    _d1 = _0
    _v1 = _0
    invGi = _0
    _a20 = _0
    _a21 = _0
    for i ∈ @comps
        zᵢ,Tcᵢ,vcᵢ,Tc0ᵢ,vc0ᵢ,d1ᵢ,v1ᵢ,Giᵢ,a20ᵢ,a21ᵢ = z[i],Tc[i],vc[i],Tc0[i],vc0[i],d1[i],v1[i],Gi[i],a20[i],a21[i]
        Δτ += zᵢ*(Tcᵢ/Tc0ᵢ - 1.0)
        Δϕ += zᵢ*(vcᵢ/vc0ᵢ - 1.0)
        _Tc0 += zᵢ*Tc0ᵢ
        _vc0 += zᵢ*vc0ᵢ
        _d1 += zᵢ*d1ᵢ
        _v1 += zᵢ*v1ᵢ
        invGi += zᵢ/Giᵢ
        _a20 += zᵢ*a20ᵢ
        _a21 += zᵢ*a21ᵢ
    end
    Δτ /= ∑z
    Δϕ /= ∑z
    _Tc0 /= ∑z
    _vc0 /= ∑z
    _Tc = _Tc0 + Δτ*_Tc0
    _vc = _vc0 + Δϕ*_vc0
    _d1 /= ∑z
    _v1 /= ∑z
    _Gi = ∑z/invGi
    _a20 /= ∑z
    _a21 /= ∑z
    b = sqrt(1.359)
    m = 1.3
    beta = 0.325
    alpha = 0.11
    gamma = 2.0-alpha-2.0*beta
    Δ = 0.51
    v = V/∑z
    τ = T/_Tc - 1.0
    ϕ = v/_vc - 1.0
    fnc = 4.0*(b/m*abs(ϕ*(1.0 + _v1*(ϕ^2)*exp(-8.5*ϕ)) + _d1*τ*(1.0 - 2.0*τ)))^(1.0/beta) + 2.0*τ
    q = sqrt((fnc + sqrt(fnc^2 + 12.0*τ^2))/6.0/_Gi)
    Y = (q/(1.0+q))^2
    τs = τ*(Y+1e-6)^(-alpha/2.0) + (1.0 + τ)*Δτ*Y^(2.0*(2.0 - alpha)/3.0)
    Ts = _Tc0*(τs + 1.0)
    ϕs = ϕ*Y^((gamma - 2.0*beta)/4.0) + (1.0 + ϕ)*Δϕ*Y^((2.0 - alpha)/2.0)
    vs = _vc0*(ϕs + 1.0)
    Δv = v/_vc0 - 1.0
    _a_res = a_res(basemodel, vs, Ts, z) - a_res(basemodel, _vc0, Ts, z) + a_res(basemodel, _vc0, T, z) + ϕs*Z_base(basemodel, _vc0, Ts, z) - Δv*Z_base(basemodel, _vc0, T, z) + log((Δv + 1.0)/(ϕs + 1.0))
    K = τ^2/2.0/(1.0 + τ^2)*(_a20*((Y + 1e-6)^(-alpha/2.0/Δ) - 1.0) + _a21*(Y^(-(alpha - Δ)/Δ) - 1.0))
    return _a_res - K
end

[branch171]

Hi I have been working on this and tested using PatelTeja,see issues for an update on that eos. This eos was used by Kiselev and Ely in 2007 and seems to work well at least for pure components. I have tested against SARTVRMie and get a good match and better speed of sound predictions. I would like to apply to mixtures, particularly CO2/N2 as this is the most dominant mixture for CCSU. I need unlike pair parameters for the basemodel. In our frame work how do I specify SAFTVRMie unlike parameters. Do I simply add SAFTVRMie_unlike.cvs file with:
carbon dioxide,O=C(=O),nitrogen,,XXX,,

[pw0908]

Hi @branch171, in principle, yup, that's all you would need to do. There is also a manual way to build the model in REPL if you are interested.