Patel-Teja: use correlation first

database/cubic/PatelTeja/PatelTeja_Vc_fit.csv

@@ -0,0 +1,3 @@
+Clapeyron Database File
+Patel-Teja Single fitted critical volumes
+species,Vc_fit

src/models/cubic/PatelTeja/PatelTeja.jl

@@ -6,7 +6,8 @@ struct PatelTejaParam <: EoSParam
     c::SingleParam{Float64}
     Tc::SingleParam{Float64}
     Pc::SingleParam{Float64}
-    Vc::SingleParam{Float64}
+    Vc_fit::SingleParam{Float64}
+    acentricfactor::SingleParam{Float64}
     Mw::SingleParam{Float64}
 end
 
@@ -18,22 +19,17 @@ function transform_params(::Type{PatelTejaParam},params,components)
     c = get!(params,"c") do
         SingleParam("c",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n))
     end
-    Vc = params["Vc"]
 
+    Vc = get!(params,"Vc_fit") do
+        SingleParam("Vc_fit",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n))
+    end
 
-    Vc = get!(params,"Vc") do
-        SingleParam("Vc",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n))
+    acentricfactor = get!(params,"acentricfactor") do
+        SingleParam("acentricfactor",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n))
     end
-    ω = get(params,"acentricfactor",nothing)
-    if any(Vc.ismissingvalues)
-        isnothing(ω) && throw(MissingException("PatelTeja: cannot estimate Vc: missing acentricfactor parameter."))
-
-        for i in 1:n
-            if Vc.ismissingvalues[i]
-                ζc = evalpoly(ω[i],(0.329032,-0.076799,0.0211947))
-                Vc[i] = ζc * R̄ * Tc[i] / Pc[i]
-            end
-        end
+
+    if all(Vc.ismissingvalues) && all(acentricfactor.ismissingvalues)
+        throw(MissingException("PatelTeja: cannot estimate Vc: missing acentricfactor parameter."))
     end
 
     return params
@@ -68,15 +64,18 @@ end
 ## Input parameters
 - `Tc`: Single Parameter (`Float64`) - Critical Temperature `[K]`
 - `Pc`: Single Parameter (`Float64`) - Critical Pressure `[Pa]`
-- `Vc`: Single Parameter (`Float64`) - Critical Volume `[m³·mol⁻¹]`
+- `Vc_fit`: Single Parameter (`Float64`) - Fitted Critical Volume `[m³·mol⁻¹]`
 - `Mw`: Single Parameter (`Float64`) - Molecular Weight `[g·mol⁻¹]`
 - `k`: Pair Parameter (`Float64`) (optional)
 - `l`: Pair Parameter (`Float64`) (optional)
+- `acentricfactor`: Single Parameter (`Float64`) (optional) - acentric factor, used to fit the critical volume if no value is provided.
+
 
 ## Model Parameters
 - `Tc`: Single Parameter (`Float64`) - Critical Temperature `[K]`
 - `Pc`: Single Parameter (`Float64`) - Critical Pressure `[Pa]`
-- `Vc`: Single Parameter (`Float64`) - Critical Volume `[m³·mol⁻¹]`
+- `Vc_fit`: Single Parameter (`Float64`) - Fitted Critical Volume `[m³·mol⁻¹]`
+- `acentricfactor`: Single Parameter (`Float64`)
 - `Mw`: Single Parameter (`Float64`) - Molecular Weight `[g·mol⁻¹]`
 - `a`: Pair Parameter (`Float64`)
 - `b`: Pair Parameter (`Float64`)
@@ -102,7 +101,7 @@ Zcᵢ = Pcᵢ*Vcᵢ/(R*Tcᵢ)
 Δ₁ = -(ϵ + √δ)/2
 Δ₂ = -(ϵ - √δ)/2
 ```
-if `Vc` is not known, `Zc` can be estimated, using the acentric factor:
+if `Vc_fit` is not known, `Zc` can be estimated, using the acentric factor:
 
 ```
 Zc = 0.329032 - 0.076799ω + 0.0211947ω²
@@ -132,7 +131,6 @@ model = PatelTeja(["neon","hydrogen"]; userlocations = ["path/to/my/db","cubic/m
 model = PatelTeja(["neon","hydrogen"];
         userlocations = (;Tc = [44.492,33.19],
                         Pc = [2679000, 1296400],
-                        Vc = [4.25e-5, 6.43e-5],
                         Mw = [20.17, 2.],
                         acentricfactor = [-0.03,-0.21]
                         k = [0. 0.18; 0.18 0.], #k,l can be ommited in single-component models.
@@ -164,7 +162,7 @@ function PatelTeja(components;
     verbose = false)
 
     formatted_components = format_components(components)
-    params = getparams(formatted_components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations = userlocations, verbose = verbose,ignore_missing_singleparams = ["Vc"])
+    params = getparams(formatted_components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv","cubic/PatelTeja/PatelTeja_Vc_fit.csv"]; userlocations = userlocations, verbose = verbose,ignore_missing_singleparams = ["Vc_fit","acentricfactor"])
     model = CubicModel(PatelTeja,params,formatted_components;
                         idealmodel,alpha,mixing,activity,translation,
                         userlocations,ideal_userlocations,alpha_userlocations,activity_userlocations,mixing_userlocations,translation_userlocations,
@@ -180,7 +178,7 @@ end
 default_references(::Type{PatelTeja}) = ["10.1016/0009-2509(82)80099-7"]
 
 function PatelTeja_Ωb(ξc)
-    #poly = (-Zc^3,3Zc^2,2-3*Zc,1.0)
+    #poly = (-ξc^3,3ξc^2,2-3*ξc,1.0)
     #_,Ωb1,Ωb2,Ωb3 = Solvers.real_roots3(poly)
     #Ωb = max(Ωb1,Ωb2,Ωb3)
     Z̄c = complex(ξc)
@@ -199,15 +197,14 @@ end
 function ab_premixing(model::PatelTejaModel,mixing::MixingRule,k,l)
     _Tc = model.params.Tc
     _pc = model.params.Pc
-    _Vc = model.params.Vc
+    _Vc = model.params.Vc_fit
     a = model.params.a
     b = model.params.b
 
     for i in 1:length(model)
         pci,Tci,Vci = _pc[i],_Tc[i],_Vc[i]
         Zc = pci * Vci / (R̄ * Tci)
-        Z̄c = complex(Zc,zero(Zc))
-        d  = (36*Z̄c - 8 - 27*Z̄c^2 + 3*sqrt(3)*sqrt(27*Z̄c^4 -8*Z̄c^3))^(1/3)
+        #d  = (36*Z̄c - 8 - 27*Z̄c^2 + 3*sqrt(3)*sqrt(27*Z̄c^4 -8*Z̄c^3))^(1/3)
         Ωa,Ωb = PatelTeja_Ωab(Zc)
         a[i] = Ωa*R̄^2*Tci^2/pci
         b[i] = Ωb*R̄*Tci/pci
@@ -220,11 +217,21 @@ end
 function c_premixing(model::PatelTejaModel)
     _Tc = model.params.Tc
     _pc = model.params.Pc
-    _Vc = model.params.Vc
+    _Vc = model.params.Vc_fit
+    ω = model.params.acentricfactor
     c = model.params.c
-    _Zc = _pc .* _Vc ./ (R̄ .* _Tc)
-    Ωc = @. 1 - 3*_Zc
-    c .= Ωc .* R̄ .*_Tc ./ _pc
+
+    for i in 1:length(model)
+        Tc,Pc = _Tc[i],_pc[i]
+        if _Vc.ismissingvalues[i]
+            ω.ismissingvalues[i] && throw(MissingException("PatelTeja: cannot estimate Vc: missing acentricfactor parameter."))
+            ζc = evalpoly(ω[i],(0.329032,-0.076799,0.0211947))
+            _Vc[i] = ζc * R̄ * Tc / Pc
+        end
+        Zc = Pc*_Vc[i]/(R̄*Tc)
+        Ωc = 1 - 3*Zc
+        c[i] = Ωc*R̄*Tc/Pc
+    end
     return c
 end
 

src/models/cubic/equations.jl

@@ -248,7 +248,7 @@ function crit_pure(model::DeltaCubicModel)
     c = translation(model,Vc0,Tc,SA[1.0])
     Vc = Vc0 - c[1]
     #we know that in AB-cubics, the critical point is already determined.
-    model isa ABCubicModel && return (Tc,Pc,Vc)
+    #model isa ABCubicModel && return (Tc,Pc,Vc)
 
     #for a general cubic model, we check if the critical pressure corresponds to the calculated pressure
     a = model.params.a[1,1]

test/test_methods_eos.jl

@@ -197,10 +197,10 @@ GC.gc()
     p = 1e5
     T = 550.15
     @testset "Bulk properties" begin
-        @test Clapeyron.volume(system, p, T) ≈ 0.04557254632681239 rtol = 1e-6
+        @test Clapeyron.volume(system, p, T) ≈ 0.0456092301152745 rtol = 1e-6
     end
     @testset "VLE properties" begin
-        @test Clapeyron.saturation_pressure(system, T)[1] ≈ 4.51634223156497e6 rtol = 1E-6
+        @test Clapeyron.saturation_pressure(system, T)[1] ≈ 6.147730717654085e6 rtol = 1E-6
         Tc,Pc,Vc = Clapeyron.crit_pure(system)
         @test Tc == system.params.Tc.values[1]
         @test Pc == system.params.Pc.values[1]

test/test_models_cubic.jl

@@ -217,7 +217,7 @@
     @testset "Patel Teja Models" begin
         @testset "Patel Teja" begin
             system = PatelTeja(["ethane","undecane"])
-            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2284322450064429 rtol = 1e-6
+            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2326465478280517 rtol = 1e-6
             @test Clapeyron.cubic_p(system, V, T, z) ≈ Clapeyron.pressure(system, V, T, z) rtol = 1e-6
         end
         @testset "PTV" begin