Mixing rules

database/properties/critical.csv

@@ -12,7 +12,7 @@ n-Heptane,142825,540.2,2.72e6,0.346
 n-Octane,111659,568.7,2.47e6,0.396
 n-Nonane,111842,594.6,2.31e6,0.446
 n-Decane,124185,617.7,2.09e6,0.488
-n-Undecane,1120214,639,1.95e6,0.53
+undecane,1120214,639,1.95e6,0.53
 n-Dodecane,112403,658,1.82e6,0.577
 n-Tridecane,629505,675,1.68e6,0.617
 n-Tetradecane,629594,693,1.57e6,0.643

docs/src/api/parameters.md

@@ -27,4 +27,5 @@ Clapeyron.SiteParam
 ## Model Splitting
 ```@docs
 Clapeyron.split_model
+Clapeyron.is_splittable
 ```

src/Clapeyron.jl

@@ -48,11 +48,18 @@ include("models/eos/SAFT/CKSAFT/variants/sCKSAFT.jl")
 
 include("models/eos/SAFT/BACKSAFT/BACKSAFT.jl")
 
-include("models/eos/cubic/alphas/alphas.jl")
 
 include("models/eos/cubic/vdW.jl")
 include("models/eos/cubic/RK/RK.jl")
 include("models/eos/cubic/PR/PR.jl")
+
+include("models/eos/cubic/alphas/alphas.jl")
+
+include("models/eos/cubic/mixing/mixing.jl")
+
+include("models/eos/cubic/RK/variants/SRK.jl")
+include("models/eos/cubic/PR/variants/PR78.jl")
+
 include("models/eos/cubic/equations.jl")
 
 include("models/eos/EmpiricHelmholtz/IAPWS95.jl")

src/methods/differentials.jl

@@ -154,7 +154,7 @@ returns the second order volume (`V`) and temperature (`T`) derivatives of the t
 
 ```
 [ ∂²f/∂V²  ∂²f/∂V∂T
- ∂²f/∂V∂T  ∂²f/∂V²]
+ ∂²f/∂V∂T  ∂²f/∂T²]
  ```
 
 use this instead of the ∂2f if you only need second order information. ∂2f also gives zeroth and first order derivative information, but due to a bug in the used AD, it allocates more than necessary.

src/models/eos/SPUNG/SPUNG.jl

@@ -40,15 +40,16 @@ end
 function shape_factors(model::SPUNG{<:ABCubicModel},V,T,z=SA[1.0])
     n = sum(z)
     x = z * (1/n)
-    a,b = cubic_ab(model.shape_model,T,x)
-    a0,b0 = cubic_ab(model.shape_ref,T,SA[1.0])
+    a,b = cubic_ab(model.shape_model,V,T,x)
+    a0,b0 = cubic_ab(model.shape_ref,V,T,SA[1.0])
     h = b/b0
     fh = n*a/a0
     f = fh/h
     return f,h
 end
 
 mw(model::SPUNG) = mw(model.shape_model)
+molecular_weight(model::SPUNG,z=SA[1.0]) = molecular_weight(model.shape_model,z)
 
 function Base.show(io::IO,mime::MIME"text/plain",model::SPUNG)
     println(io,"Extended Corresponding States model")
@@ -74,14 +75,14 @@ function T_scale(model::SPUNG,z=SA[1.0])
     return T0*f
 end
 
-function p_scale(model::SPUNG,z=SA[1.0])
-    lb_v0 = lb_volume(model.model_ref)
-    T0 = T_scale(model.model_ref)
-    p0 = p_scale(model.model_ref)
-    f,h = shape_factors(model,lb_v0,T0,z) #h normaly should be independent of temperature
-    ps = p0*f/h
-    return ps
-end
+# function p_scale(model::SPUNG,z=SA[1.0])
+#     lb_v0 = lb_volume(model.model_ref)
+#     T0 = T_scale(model.model_ref)
+#     p0 = p_scale(model.model_ref)
+#     f,h = shape_factors(model,lb_v0,T0,z) #h normaly should be independent of temperature
+#     ps = p0*f/h
+#     return ps
+# end
 
 #=
 ideally we could perform SPUNG only providing x0, but i cant find the error here
@@ -115,14 +116,14 @@ function sat_pure_p(model::SPUNG,p::Real)
 end
 =#
 
-function split_model(model::SPUNG)
-    #only the shape model is splittable
-    shape_model = split_model(model.shape_model)
-    len = length(shape_model)
-    shape_ref = fill(model.shape_ref,len)
-    model_ref = fill(model.model_ref,len)
-    return SPUNG.(shape_model,shape_ref,model_ref)
-end
+# function split_model(model::SPUNG)
+#     #only the shape model is splittable
+#     shape_model = split_model(model.shape_model)
+#     len = length(shape_model)
+#     shape_ref = fill(model.shape_ref,len)
+#     model_ref = fill(model.model_ref,len)
+#     return SPUNG.(shape_model,shape_ref,model_ref)
+# end
 
 function shape_factors(model::SPUNG,V,T,z=SA[1.0])
     n = sum(z)

src/models/eos/cubic/PR/PR.jl

@@ -9,11 +9,11 @@ struct PRParam <: EoSParam
     Mw::SingleParam{Float64}
 end
 
-struct PR{T <: IdealModel,α,γ} <: RKModel
+struct PR{T <: IdealModel,α,γ} <:PRModel
     components::Array{String,1}
     icomponents::UnitRange{Int}
     alpha::α
-    activity::γ
+    mixing::γ
     params::PRParam
     idealmodel::T
     absolutetolerance::Float64
@@ -34,16 +34,16 @@ end
 
 Base.length(model::PR) = Base.length(model.icomponents)
 
-molecular_weight(model::PR,z=SA[1.0]) = group_molecular_weight(model.groups,mw(model),z)
+molecular_weight(model::PR,z=SA[1.0]) = comp_molecular_weight(mw(model),z)
 
 export PR
 function PR(components::Vector{String}; idealmodel=BasicIdeal,
-    alpha = SoaveAlpha,
-    activity = nothing,
+    alpha = PRAlpha,
+    mixing = vdW1fRule,
     userlocations=String[], 
     ideal_userlocations=String[],
     alpha_userlocations = String[],
-    activity_userlocations = String[],
+    mixing_userlocations = String[],
      verbose=false)
     params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
     k  = params["k"]
@@ -59,11 +59,11 @@ function PR(components::Vector{String}; idealmodel=BasicIdeal,
 
     init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
     init_alpha = init_model(alpha,components,alpha_userlocations,verbose)
-    init_activity = init_model(activity,components,activity_userlocations,verbose)
+    init_mixing = init_model(mixing,components,mixing_userlocations,verbose)
     icomponents = 1:length(components)
     packagedparams = PRParam(a, b, params["Tc"],_pc,Mw)
     references = String[]
-    model = PR(components,icomponents,init_alpha,init_activity,packagedparams,init_idealmodel,1e-12,references)
+    model = PR(components,icomponents,init_alpha,init_mixing,packagedparams,init_idealmodel,1e-12,references)
     return model
 end
 
@@ -72,39 +72,35 @@ function ab_consts(::Type{<:PRModel})
     return 0.457235,0.077796
 end
 
-function cubic_ab(model::PR{<:Any,SoaveAlpha},T,z=SA[1.0],n=sum(z))
+function cubic_ab(model::PR{<:Any,<:Any,<:Any},V,T,z=SA[1.0],n=sum(z))
+    invn2 = (one(n)/n)^2
     a = model.params.a.values
     b = model.params.b.values
-    α = model.alpha
-    ω = α.params.acentricfactor.values
-    Tc = model.params.Tc.values
-    invn = one(n)/n
-    αx = @. (1+(0.37464+1.54226*ω-0.26992*ω^2)*(1-√(T/Tc))) * z * invn 
-    āᾱ = dot(αx,Symmetric(a),αx)
-    b̄ = dot(z,Symmetric(b),z) * invn*invn
-    return āᾱ ,b̄
+    α = @f(α_function,model.alpha)
+    ā,b̄ = @f(mixing_rule,model.mixing,α,a,b)
+    return ā ,b̄
 end
 
-
 function cubic_abp(model::PRModel, V, T, z) 
     n = sum(z)
     v = V/n
-    āᾱ ,b̄ = cubic_ab(model,T,z,n)
+    āᾱ ,b̄ = cubic_ab(model,V,T,z,n)
     _1 = one(b̄)
     denom = evalpoly(v,(-b̄*b̄,2*b̄,_1))
     p = R̄*T/(v-b̄) - āᾱ /denom
     return āᾱ, b̄,p
 end
 
 function cubic_poly(model::PRModel,p,T,z)
-    a,b = cubic_ab(model,T,z)
+    a,b = cubic_ab(model,p,T,z)
     RT⁻¹ = 1/(R̄*T)
     A = a*p*RT⁻¹*RT⁻¹
     B = b*p*RT⁻¹
     k₀ = B*(B*(B+1.0)-A)
     k₁ = -B*(3*B+2.0) + A
     k₂ = B-1.0
     k₃ = 1.0
+    return [k₀,k₁,k₂,k₃]
 end
 #=
  (-B2-2(B2+B)+A)
@@ -113,7 +109,7 @@ end
 =#
 function a_res(model::PRModel, V, T, z)
     n = sum(z)
-    a,b = cubic_ab(model,T,z,n)
+    a,b = cubic_ab(model,V,T,z,n)
     Δ1 = 1+√2
     Δ2 = 1-√2
     ΔPRΔ = 2*√2

src/models/eos/cubic/PR/variants/PR78.jl

@@ -0,0 +1,19 @@
+#just a function, no struct
+function PR78(components::Vector{String}; idealmodel=BasicIdeal,
+    mixing = vdW1fRule,
+    userlocations=String[], 
+    ideal_userlocations=String[],
+    alpha_userlocations = String[],
+    mixing_userlocations = String[],
+     verbose=false)
+
+     return PR(components;
+     idealmodel = idealmodel,
+     alpha = PR78Alpha,
+     mixing=mixing,
+     ideal_userlocations = ideal_userlocations,
+     alpha_userlocations = alpha_userlocations,
+     mixing_userlocations = mixing_userlocations,
+     verbose = verbose)
+end
+export PR78

src/models/eos/cubic/RK/RK.jl

@@ -12,7 +12,7 @@ struct RK{T <: IdealModel,α,γ} <: RKModel
     components::Array{String,1}
     icomponents::UnitRange{Int}
     alpha::α
-    activity::γ
+    mixing::γ
     params::RKParam
     idealmodel::T
     absolutetolerance::Float64
@@ -33,16 +33,16 @@ end
 
 Base.length(model::RK) = Base.length(model.icomponents)
 
-molecular_weight(model::RK,z=SA[1.0]) = group_molecular_weight(model.groups,mw(model),z)
+molecular_weight(model::RK,z=SA[1.0]) = comp_molecular_weight(mw(model),z)
 
 export RK
 function RK(components::Vector{String}; idealmodel=BasicIdeal,
-    alpha = nothing,
-    activity = nothing,
+    alpha = RKAlpha,
+    mixing = vdW1fRule,
     userlocations=String[], 
     ideal_userlocations=String[],
     alpha_userlocations = String[],
-    activity_userlocations = String[],
+    mixing_userlocations = String[],
      verbose=false)
     params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
     k  = params["k"]
@@ -58,11 +58,11 @@ function RK(components::Vector{String}; idealmodel=BasicIdeal,
 
     init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
     init_alpha = init_model(alpha,components,alpha_userlocations,verbose)
-    init_activity = init_model(activity,components,activity_userlocations,verbose)
+    init_mixing = init_model(mixing,components,mixing_userlocations,verbose)
     icomponents = 1:length(components)
     packagedparams = RKParam(a, b, params["Tc"],_pc,Mw)
     references = String[]
-    model = RK(components,icomponents,init_alpha,init_activity,packagedparams,init_idealmodel,1e-12,references)
+    model = RK(components,icomponents,init_alpha,init_mixing,packagedparams,init_idealmodel,1e-12,references)
     return model
 end
 
@@ -71,28 +71,26 @@ function ab_consts(::Type{<:RKModel})
     Ωb = (2^(1/3)-1)/3
     return Ωa,Ωb
 end
-function cubic_ab(model::RK{<:Any,Nothing},T,z=SA[1.0],n=sum(z))
+function cubic_ab(model::RK{<:Any,<:Any,<:Any},V,T,z=SA[1.0],n=sum(z))
     invn2 = (one(n)/n)^2
     a = model.params.a.values
     b = model.params.b.values
-    Tc = model.params.Tc.values
-    T̄c = sum(sqrt.(Tc*Tc'))
-    āᾱ  = dot(z,Symmetric(a),z) * invn2*sqrt(T̄c/T)
-    b̄ = dot(z,Symmetric(b),z) * invn2
-    return āᾱ ,b̄
+    α = @f(α_function,model.alpha)
+    ā,b̄ = @f(mixing_rule,model.mixing,α,a,b)
+    return ā ,b̄
 end
 
 function cubic_abp(model::RKModel, V, T, z)
     n = sum(z)
-    a,b = cubic_ab(model,T,z,n)
+    a,b = cubic_ab(model,V,T,z,n)
     v = V/n
     p =  R̄*T/(v-b) - a/((v+b)*v)
     return a,b,p
 end
 
 function cubic_poly(model::RKModel,p,T,z)
     n = sum(z)
-    a,b = cubic_ab(model,T,z,n)
+    a,b = cubic_ab(model,p,T,z,n)
     RT⁻¹ = 1/(R̄*T)
     A = a*p* RT⁻¹* RT⁻¹
     B = b*p* RT⁻¹
@@ -103,7 +101,7 @@ end
 
 function a_res(model::RKModel, V, T, z)
     n = sum(z)
-    ā,b̄ = cubic_ab(model,T,z,n)
+    ā,b̄ = cubic_ab(model,V,T,z,n)
     ρ = n/V
     RT⁻¹ = 1/(R̄*T)
     return -log(1-b̄*ρ) - ā*RT⁻¹*log(b̄*ρ+1)/b̄
@@ -112,4 +110,4 @@ end
 
 cubic_zc(::RKModel) = 1/3
 
-include("variants/SRK.jl")
+# include("variants/SRK.jl")

src/models/eos/cubic/RK/variants/SRK.jl

@@ -1,32 +1,19 @@
-function cubic_ab(model::RK{<:Any,SoaveAlpha},T,z=SA[1.0],n=sum(z))
-    a = model.params.a.values
-    b = model.params.b.values
-    αmodel = model.alpha
-    ω = αmodel.params.acentricfactor.values
-    Tc = model.params.Tc.values
-    _1 = one(T+n)
-    invn = one(n)/n
-    αx = @. (_1+(0.480+1.547*ω-0.176*ω^2)*(1-√(T/Tc))) * z * invn 
-    āᾱ =dot(αx, Symmetric(a), αx)
-    b̄ = dot(z, Symmetric(b), z) * invn * invn
-    return āᾱ ,b̄
-end
 #just a function, no struct
 function SRK(components::Vector{String}; idealmodel=BasicIdeal,
-    activity = nothing,
+    mixing = vdW1fRule,
     userlocations=String[], 
     ideal_userlocations=String[],
     alpha_userlocations = String[],
-    activity_userlocations = String[],
+    mixing_userlocations = String[],
      verbose=false)
 
      return RK(components;
      idealmodel = idealmodel,
      alpha = SoaveAlpha,
-     activity=activity,
+     mixing=mixing,
      ideal_userlocations = ideal_userlocations,
      alpha_userlocations = alpha_userlocations,
-     activity_userlocations = activity_userlocations,
+     mixing_userlocations = mixing_userlocations,
      verbose = verbose)
 end
 export SRK