Merge pull request #36 from ypaul21/MixingRules

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,8 +40,8 @@ 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

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
@@ -37,12 +37,12 @@ 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

src/models/eos/cubic/alphas/BM.jl

@@ -0,0 +1,34 @@
+abstract type BMAlphaModel <: AlphaModel end
+
+struct BMAlphaParam <: EoSParam
+    acentricfactor::SingleParam{Float64}
+end
+
+@newmodelsimple BMAlpha BMAlphaModel BMAlphaParam
+
+export BMAlpha
+function BMAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
+    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
+    acentricfactor = SingleParam(params["w"],"acentric factor")
+    packagedparams = BMAlphaParam(acentricfactor)
+    model = BMAlpha(packagedparams, verbose=verbose)
+    return model
+end
+
+function α_function(model::RKModel,V,T,z,alpha_model::BMAlphaModel)
+    Tc = model.params.Tc.values
+    Tr = @. T/Tc
+    ω  = alpha_model.params.acentricfactor.values
+    m  = @. 0.480+1.547*ω-0.176*ω^2
+    α  = @. (Tr>1)*(exp((1-2/(2+m))*(1-Tr^(1+m/2))))^2+(Tr<=1)*(1+m*(1-√(Tr)))^2
+    return α
+end
+
+function α_function(model::PRModel,V,T,z,alpha_model::BMAlphaModel)
+    Tc = model.params.Tc.values
+    Tr = @. T/Tc
+    ω  = alpha_model.params.acentricfactor.values
+    m  = @. 0.37464+1.54226*ω-0.26992*ω^2
+    α  = @. (Tr>1)*(exp((1-2/(2+m))*(1-Tr^(1+m/2))))^2+(Tr<=1)*(1+m*(1-√(Tr)))^2
+    return α
+end

src/models/eos/cubic/alphas/PR78Alpha.jl

@@ -0,0 +1,24 @@
+abstract type PR78AlphaModel <: AlphaModel end
+
+struct PR78AlphaParam <: EoSParam
+    acentricfactor::SingleParam{Float64}
+end
+
+@newmodelsimple PR78Alpha PR78AlphaModel PR78AlphaParam
+
+export PR78Alpha
+function PR78Alpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
+    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
+    acentricfactor = SingleParam(params["w"],"acentric factor")
+    packagedparams = PR78AlphaParam(acentricfactor)
+    model = PR78Alpha(packagedparams, verbose=verbose)
+    return model
+end
+
+function α_function(model::CubicModel,V,T,z,alpha_model::PR78AlphaModel)
+    Tc = model.params.Tc.values
+    ω  = alpha_model.params.acentricfactor.values
+    m = @. (ω<=0.491)*(0.37464+1.54226*ω-0.26992*ω^2)+(ω>0.491)*(0.379642+1.487503*ω-0.164423*ω^2-0.016666*ω^3)
+    α = @. (1+m*(1-√(T/Tc)))^2
+    return α
+end

src/models/eos/cubic/alphas/PRAlpha.jl

@@ -0,0 +1,23 @@
+abstract type PRAlphaModel <: AlphaModel end
+
+struct PRAlphaParam <: EoSParam
+    acentricfactor::SingleParam{Float64}
+end
+
+@newmodelsimple PRAlpha PRAlphaModel PRAlphaParam
+
+export PRAlpha
+function PRAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
+    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
+    acentricfactor = SingleParam(params["w"],"acentric factor")
+    packagedparams = PRAlphaParam(acentricfactor)
+    model = PRAlpha(packagedparams, verbose=verbose)
+    return model
+end
+
+function α_function(model::CubicModel,V,T,z,alpha_model::PRAlphaModel)
+    Tc = model.params.Tc.values
+    ω  = alpha_model.params.acentricfactor.values
+    α = @. (1+(0.37464+1.54226*ω-0.26992*ω^2)*(1-√(T/Tc)))^2
+    return α
+end