Merge pull request #55 from ypaul21/cubic-optimization

Cubic (mixing) optimization

Project.toml

@@ -1,7 +1,7 @@
 name = "Clapeyron"
 uuid = "7c7805af-46cc-48c9-995b-ed0ed2dc909a"
 authors = ["Hon Wa Yew <yewhonwa@gmail.com>", "Pierre Walker <pwalker@mit.edu>", "Andrés Riedemann <andres.riedemann@gmail.com>"]
-version = "0.2.8"
+version = "0.2.9"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

src/Clapeyron.jl

@@ -56,6 +56,7 @@ include("models/eos/SAFT/BACKSAFT/BACKSAFT.jl")
 
 include("models/eos/SAFT/equations.jl")
 
+include("models/eos/cubic/equations.jl")
 include("models/eos/cubic/vdW.jl")
 include("models/eos/cubic/RK/RK.jl")
 include("models/eos/cubic/PR/PR.jl")
@@ -80,7 +81,6 @@ include("models/eos/cubic/PR/variants/PR78.jl")
 include("models/eos/cubic/PR/variants/VTPR.jl")
 include("models/eos/cubic/PR/variants/UMRPR.jl")
 
-include("models/eos/cubic/equations.jl")
 include("models/eos/Activity/equations.jl")
 
 include("models/eos/EmpiricHelmholtz/IAPWS95.jl")

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

@@ -38,29 +38,24 @@ function PR(components::Vector{String}; idealmodel=BasicIdeal,
      verbose=false)
     params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
     k  = params["k"]
-    _pc = params["pc"]
-    pc = _pc.values
+    pc = params["pc"]
     Mw = params["Mw"]
-    _Tc = params["Tc"]
-    Tc = _Tc.values
-    Ωa, Ωb = ab_consts(PR)
-    a = epsilon_LorentzBerthelot(SingleParam(params["pc"], @. Ωa*R̄^2*Tc^2/pc),k)
-    #check if this is correct in the general case.
-    b = sigma_LorentzBerthelot(SingleParam(params["pc"], @. Ωb*R̄*Tc/pc))
-
+    Tc = params["Tc"]
+    init_mixing = init_model(mixing,components,activity,mixing_userlocations,activity_userlocations,verbose)
+    a,b = ab_premixing(PR,init_mixing,Tc,pc,k)
     init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
     init_alpha = init_model(alpha,components,alpha_userlocations,verbose)
-    init_mixing = init_model(mixing,components,activity,mixing_userlocations,activity_userlocations,verbose)
     init_translation = init_model(translation,components,translation_userlocations,verbose)
-
     icomponents = 1:length(components)
-    packagedparams = PRParam(a, b, params["Tc"],_pc,Mw)
+    packagedparams = PRParam(a,b,Tc,pc,Mw)
     references = String[]
     model = PR(components,icomponents,init_alpha,init_mixing,init_translation,packagedparams,init_idealmodel,1e-12,references)
     return model
 end
 
 
+
+
 function ab_consts(::Type{<:PRModel})
     return 0.457235,0.077796
 end

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

@@ -37,23 +37,16 @@ function RK(components::Vector{String}; idealmodel=BasicIdeal,
      verbose=false)
     params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
     k  = params["k"]
-    _pc = params["pc"]
-    pc = _pc.values
+    pc = params["pc"]
     Mw = params["Mw"]
-    _Tc = params["Tc"]
-    Tc = _Tc.values
-    #T̄c = sum(sqrt.(Tc*Tc')) #is this term correctly calculated? sqrt(Tc*Tc') is a matrix sqrt
-    Ωa, Ωb = ab_consts(RK)
-    a = epsilon_LorentzBerthelot(SingleParam(params["pc"], @. Ωa*R̄^2*Tc^2/pc), k) 
-    b = sigma_LorentzBerthelot(SingleParam(params["pc"], @. Ωb*R̄*Tc/pc))
-
+    Tc = params["Tc"]
+    init_mixing = init_model(mixing,components,activity,mixing_userlocations,activity_userlocations,verbose)
+    a,b = ab_premixing(RK,init_mixing,Tc,pc,k)
     init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
     init_alpha = init_model(alpha,components,alpha_userlocations,verbose)
-    init_mixing = init_model(mixing,components,activity,mixing_userlocations,activity_userlocations,verbose)
     init_translation = init_model(translation,components,translation_userlocations,verbose)
-
     icomponents = 1:length(components)
-    packagedparams = RKParam(a, b, params["Tc"],_pc,Mw)
+    packagedparams = RKParam(a,b,Tc,pc,Mw)
     references = String[]
     model = RK(components,icomponents,init_alpha,init_mixing,init_translation,packagedparams,init_idealmodel,1e-12,references)
     return model

src/models/eos/cubic/equations.jl

@@ -1,3 +1,12 @@
+function ab_premixing(::Type{T},mixing,Tc,pc,kij) where T <: ABCubicModel
+    Ωa, Ωb = ab_consts(T)
+    _Tc = Tc.values
+    _pc = pc.values
+    a = epsilon_LorentzBerthelot(SingleParam(pc, @. Ωa*R̄^2*_Tc^2/_pc),kij)
+    b = sigma_LorentzBerthelot(SingleParam(pc, @. Ωb*R̄*_Tc/_pc))
+    return a,b
+end
+
 function cubic_ab(model::ABCubicModel,V,T,z=SA[1.0],n=sum(z))
     invn2 = (one(n)/n)^2
     a = model.params.a.values

src/models/eos/cubic/mixing/Kay.jl

@@ -17,10 +17,29 @@ end
 
 function mixing_rule(model::ABCubicModel,V,T,z,mixing_model::KayRuleModel,α,a,b,c)
     n = sum(z)
-    invn2 = (one(n)/n)^2
-    c̄ = dot(z,c)/n
-    b̄ = (dot(z,Symmetric(b.^(1/3)),z) * invn2).^3
-    ā = √(dot(z,Symmetric(a .* sqrt.(α*α') ./ b).^2,z)* invn2)  * b̄
+    invn = (one(n)/n)
+    invn2 = invn^2
+    #b̄ = (dot(z,Symmetric(b.^(1/3)),z) * invn2).^3
+    #ā = √(dot(z,Symmetric(a .* sqrt.(α*α') ./ b).^2,z)* invn2)  * b̄
+    b̄ = zero(eltype(z))
+    ab2 = zero(T+first(z))
+    for i in 1:length(z)
+        zi = z[i]
+        zi2 = zi*zi
+        αi = α[i]
+        b̄ += cbrt(b[i,i])*zi2
+        abij2 = (a[i,i]*αi/b[i,i])^2
+        ab2 += abij2*zi2
+        for j in 1:(i-1)
+            zij =zi*z[j]
+            abij2 = (a[i,j]*sqrt(αi*α[j])/b[i,j])^2
+            ab2 += 2*abij2*zij
+            b̄ += 2*cbrt(b[i,j])*zij
+        end
+    end    
+    b̄ = (b̄*invn2)^3
+    ā = sqrt(ab2*invn2)*b̄
+    c̄ = dot(z,c)*invn
     return ā,b̄,c̄
 end
 

src/models/eos/cubic/mixing/LCVM.jl

@@ -9,27 +9,26 @@ end
 export LCVMRule
 function LCVMRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
     init_activity = activity(components;userlocations = activity_userlocations,verbose)
-
     model = LCVMRule(components, init_activity)
     return model
 end
 
 
 function mixing_rule(model::PRModel,V,T,z,mixing_model::LCVMRuleModel,α,a,b,c)
     n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
+    #x = z./n
+    invn = (one(n)/n)
+    invn2 = invn^2
     g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
     b̄ = dot(z,Symmetric(b),z) * invn2
     c̄ = dot(z,c)/n
-
-    ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)
-
+    #ᾱ  = a.*sqrt.(α.*α')./(b*R̄*T)
+    Σxᾱ  = sum(α[i]*a[i,i]*z[i]/b[i,i] for i ∈ @comps)*invn/(R̄*T)
     λ  = 0.7
     AV = -0.52
     AM = -0.623
-    C1 = λ/AV+(1-λ)/AM
-
-    ā = b̄*R̄*T*(C1*(g_E/(R̄*T)-0.3*sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps))+sum(x[i]*ᾱ[i,i] for i ∈ @comps))
+    C1 = -1.8276947771329795 #λ/AV+(1-λ)/AM,is this ok?, that C1 is independent of the input conditions
+    Σlogb = sum(z[i]*log(b̄/b[i,i]) for i ∈ @comps)*invn #sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)
+    ā = b̄*R̄*T*(C1*(g_E/(R̄*T)-0.3*Σlogb)+Σxᾱ )
     return ā,b̄,c̄
 end

src/models/eos/cubic/mixing/MHV1.jl

@@ -14,24 +14,18 @@ function MHV1Rule(components::Vector{String}; activity = Wilson, userlocations::
     return model
 end
 
-function mixing_rule(model::RKModel,V,T,z,mixing_model::MHV1RuleModel,α,a,b,c)
+MHV1q(::PRModel) = 0.53
+MHV1q(::RKModel) = 0.593
+function mixing_rule(model::Union{RKModel,PRModel},V,T,z,mixing_model::MHV1RuleModel,α,a,b,c)
     n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
+    invn = (one(n)/n)
+    invn2 = invn^2
     g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
     b̄ = dot(z,Symmetric(b),z) * invn2
     c̄ = dot(z,c)/n
-    ā = b̄*R̄*T*(sum(x[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)-1/0.593*(g_E/(R̄*T)+sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)))
-    return ā,b̄,c̄
-end
-
-function mixing_rule(model::PRModel,V,T,z,mixing_model::MHV1RuleModel,α,a,b,c)
-    n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
-    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
-    b̄ = dot(z,Symmetric(b),z) * invn2
-    c̄ = dot(z,c)/n
-    ā = b̄*R̄*T*(sum(x[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)-1/0.53*(g_E/(R̄*T)+sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)))
+    q = MHV1q(model)
+    Σlogb = sum(z[i]*log(b̄/b[i,i]) for i ∈ @comps)*invn
+    Σab = invn*sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)
+    ā = b̄*R̄*T*(Σab-1/q*(g_E/(R̄*T)+Σlogb))
     return ā,b̄,c̄
 end

src/models/eos/cubic/mixing/MHV2.jl

@@ -14,24 +14,35 @@ function MHV2Rule(components::Vector{String}; activity = Wilson, userlocations::
     return model
 end
 
-function mixing_rule(model::RKModel,V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
+MHV2q(::PRModel) = (-0.4347,-0.003654)
+MHV2q(::RKModel) = (-0.4783,-0.0047)
+function mixing_rule(model::Union{PRModel,RKModel},V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
     n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
-    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
+    invn = (one(n)/n)
+    invn2 = invn^2
+    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)*invn
     b̄ = dot(z,Symmetric(b),z) * invn2
-    c̄ = dot(z,c)/n
-
-    ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)
-
-    q1 = -0.4783
-    q2 = -0.0047
-    c  = -q1*sum(x[i]*ᾱ[i,i] for i ∈ @comps)-q2*sum(x[i]*ᾱ[i,i]^2 for i ∈ @comps)-g_E/(R̄*T)-sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)
-
+    c̄ = dot(z,c)*invn
+    #ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)
+    q1,q2 = MHV2q(model)
+    Σlogb = zero(first(z))
+    Σab = zero(T+first(z))
+    Σab2 = Σab
+    for i ∈ @comps
+        ᾱ  = a[i,i]*α[i]/(b[i,i]*R̄*T)
+        zia = z[i]*ᾱ  
+        Σab += zia 
+        Σab2 += zia*ᾱ  
+        Σlogb += z[i]*log(b̄/b[i,i])
+    end
+    Σab *= invn
+    Σab2 *= invn
+    Σlogb *= invn
+    c  = -q1*Σab-q2*Σab2-g_E/(R̄*T)-Σlogb
     ā = b̄*R̄*T*(-q1-sqrt(q1^2-4*q2*c))/(2*q2)
     return ā,b̄,c̄
 end
-
+#=
 function mixing_rule(model::PRModel,V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
     n = sum(z)
     x = z./n
@@ -48,4 +59,4 @@ function mixing_rule(model::PRModel,V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
 
     ā = b̄*R̄*T*(-q1-sqrt(q1^2-4*q2*c))/(2*q2)
     return ā,b̄,c̄
-end
+end=#

src/models/eos/cubic/mixing/PSRK.jl

@@ -16,11 +16,13 @@ end
 
 function mixing_rule(model::RKModel,V,T,z,mixing_model::PSRKRuleModel,α,a,b,c)
     n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
+    invn = (one(n)/n)
+    invn2 = invn^2
     g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
     b̄ = dot(z,Symmetric(b),z) * invn2
     c̄ = dot(z,c)/n
-    ā = b̄*R̄*T*(sum(x[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)-1/0.64663*(g_E/(R̄*T)+sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)))
+    Σlogb = sum(z[i]*log(b̄/b[i,i]) for i ∈ @comps)*invn
+    Σab = sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)*invn
+    ā = b̄*R̄*T*(Σab-1/0.64663*(g_E/(R̄*T)+Σlogb))
     return ā,b̄,c̄
 end

src/models/eos/cubic/mixing/UMR.jl

@@ -8,23 +8,35 @@ end
 @registermodel UMRRule
 export UMRRule
 function UMRRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
-    init_activity = activity(components;userlocations = activity_userlocations,verbose)
-
+    init_activity = activity(components;userlocations = activity_userlocations,verbose)   
     model = UMRRule(components, init_activity)
     return model
 end
 
+function ab_premixing(::Type{PR},mixing::UMRRuleModel,Tc,pc,kij)
+    Ωa, Ωb = ab_consts(PR)
+    _Tc = Tc.values
+    _pc = pc.values
+    a = epsilon_LorentzBerthelot(SingleParam(pc, @. Ωa*R̄^2*_Tc^2/_pc),kij)
+    bi = @. Ωb*R̄*_Tc/_pc
+    bij = ((bi.^(1/2).+bi'.^(1/2))/2).^2
+    b = PairParam("b",Tc.components,bij)
+    return a,b
+end
+
 function mixing_rule(model::PRModel,V,T,z,mixing_model::UMRRuleModel,α,a,b,c)
     n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
+    #x = z./n
+    invn = (one(n)/n)
+    invn2 = invn^2
     lnγ_SG_  = lnγ_SG(mixing_model.activity,1e5,T,z)
     lnγ_res_ = lnγ_res(mixing_model.activity,1e5,T,z)
-    g_E = sum(x[i]*R̄*T*(lnγ_res_[i]+lnγ_SG_[i]) for i ∈ @comps)
-    b = Diagonal(b).diag
-    b = ((b.^(1/2).+b'.^(1/2))/2).^2
+    g_E = sum(z[i]*R̄*T*(lnγ_res_[i]+lnγ_SG_[i]) for i ∈ @comps)*invn
+    #b = Diagonal(b).diag
+    #b = ((b.^(1/2).+b'.^(1/2))/2).^2
     b̄ = dot(z,Symmetric(b),z) * invn2
-    c̄ = dot(z,c)/n
-    ā = b̄*R̄*T*(sum(x[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)-1/0.53*g_E/(R̄*T))
+    c̄ = dot(z,c)*invn
+    Σab = sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)*invn
+    ā = b̄*R̄*T*(Σab-1/0.53*g_E/(R̄*T))
     return ā,b̄,c̄
 end

src/models/eos/cubic/mixing/VTPR.jl

@@ -14,16 +14,26 @@ function VTPRRule(components::Vector{String}; activity = Wilson, userlocations::
     return model
 end
 
+function ab_premixing(::Type{PR},mixing::VTPRRule,Tc,pc,kij)
+    Ωa, Ωb = ab_consts(PR)
+    _Tc = Tc.values
+    _pc = pc.values
+    a = epsilon_LorentzBerthelot(SingleParam(pc, @. Ωa*R̄^2*_Tc^2/_pc),kij)
+    bi = @. Ωb*R̄*_Tc/_pc
+    bij = ((bi.^(3/4).+bi'.^(3/4))/2).^(4/3)
+    b = PairParam("b",Tc.components,bij)
+    return a,b
+end
+
 function mixing_rule(model::PRModel,V,T,z,mixing_model::VTPRRuleModel,α,a,b,c)
     n = sum(z)
-    x = z./n
-    invn2 = (one(n)/n)^2
+    invn = (one(n)/n)
+    invn2 = invn^2
     lnγ_res_ = lnγ_res(mixing_model.activity,1e5,T,z)
-    g_E_res = sum(x[i]*R̄*T*lnγ_res_[i] for i ∈ @comps)
-    b = Diagonal(b).diag
-    b = ((b.^(3/4).+b'.^(3/4))/2).^(4/3)
+    g_E_res = sum(z[i]*R̄*T*lnγ_res_[i] for i ∈ @comps)*invn
     b̄ = dot(z,Symmetric(b),z) * invn2
     c̄ = dot(z,c)/n
-    ā = b̄*R̄*T*(sum(x[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)-1/0.53087*(g_E_res/(R̄*T)))
+    Σab = invn*sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)
+    ā = b̄*R̄*T*(Σab-1/0.53087*(g_E_res/(R̄*T)))
     return ā,b̄,c̄
 end

src/models/eos/cubic/mixing/WS.jl

@@ -14,24 +14,18 @@ function WSRule(components::Vector{String}; activity = Wilson, userlocations::Ve
     return model
 end
 
-function mixing_rule(model::RKModel,V,T,z,mixing_model::WSRuleModel,α,a,b,c)
-    n = sum(z)
-    invn2 = (one(n)/n)^2
-    num = sum(z[i]*(b[i,i]-a[i,i]*α[i]/(R̄*T)) for i ∈ @comps)/n
-    den = 1 - (sum(z[i]*a[i,i]*α[i]/b[i,i]/(n*R̄*T) for i ∈ @comps)-excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)/(n*R̄*T)/log(2))
-    c̄ = dot(z,c)/n
-    b̄  = num/den
-    ā  = R̄*T*(b̄-num)
-    return ā,b̄,c̄
-end
+WS_λ(::PRModel) = 0.6232252401402305 #1/(2*√(2))*log((2+√(2))/(2-√(2)))
+WS_λ(::RKModel) = 0.6931471805599453#log(2)
 
-function mixing_rule(model::PRModel,V,T,z,mixing_model::WSRuleModel,α,a,b,c)
-    λ = 1/(2*√(2))*log((2+√(2))/(2-√(2)))
+function mixing_rule(model::Union{RKModel,PRModel},V,T,z,mixing_model::WSRuleModel,α,a,b,c)
+    λ = WS_λ(model)
     n = sum(z)
-    invn2 = (one(n)/n)^2
-    num = sum(z[i]*(b[i,i]-a[i,i]*α[i]/(R̄*T)) for i ∈ @comps)/n
-    den = 1 - (sum(z[i]*a[i,i]*α[i]/b[i,i]/(n*R̄*T) for i ∈ @comps)-excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)/(n*R̄*T)/λ)
-    c̄ = dot(z,c)/n
+    invn = (one(n)/n)
+    Σab = sum(z[i]*a[i,i]*α[i]/b[i,i]/(n*R̄*T) for i ∈ @comps)
+    num = sum(z[i]*(b[i,i]-a[i,i]*α[i]/(R̄*T)) for i ∈ @comps)*invn
+    gE = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)
+    den = 1 - (Σab-gE/(n*R̄*T)/λ)
+    c̄ = dot(z,c)*invn
     b̄  = num/den
     ā  = R̄*T*(b̄-num)
     return ā,b̄,c̄