From 3ab8ece54bd96cc9bae3f13b0ecc582c935e2466 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Tue, 30 Dec 2025 17:25:48 -0300 Subject: [PATCH 1/4] modify PatelTeja to use correlation first instead of experimental critical volume --- database/cubic/PatelTeja/PatelTeja_Vc_fit.csv | 3 ++ src/models/cubic/PatelTeja/PatelTeja.jl | 47 ++++++++++--------- 2 files changed, 28 insertions(+), 22 deletions(-) create mode 100644 database/cubic/PatelTeja/PatelTeja_Vc_fit.csv diff --git a/database/cubic/PatelTeja/PatelTeja_Vc_fit.csv b/database/cubic/PatelTeja/PatelTeja_Vc_fit.csv new file mode 100644 index 0000000000..f0a9dcac44 --- /dev/null +++ b/database/cubic/PatelTeja/PatelTeja_Vc_fit.csv @@ -0,0 +1,3 @@ +Clapeyron Database File +Patel-Teja Single fitted critical volumes +species,Vc_fit \ No newline at end of file diff --git a/src/models/cubic/PatelTeja/PatelTeja.jl b/src/models/cubic/PatelTeja/PatelTeja.jl index 0509cd036c..5ed7fbe93a 100755 --- a/src/models/cubic/PatelTeja/PatelTeja.jl +++ b/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,13 @@ 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") do - SingleParam("Vc",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n)) + Vc = get!(params,"Vc_fit") do + SingleParam("Vc_fit",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 + + acentricfactor = get!(params,"acentricfactor") do + SingleParam("acentricfactor",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n)) end return params @@ -68,15 +60,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 +97,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 +127,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 +158,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"]) model = CubicModel(PatelTeja,params,formatted_components; idealmodel,alpha,mixing,activity,translation, userlocations,ideal_userlocations,alpha_userlocations,activity_userlocations,mixing_userlocations,translation_userlocations, @@ -199,9 +193,18 @@ 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 + ω = model.params.acentricfactor + + for i in 1:length(model) + if _Vc.ismissingvalues[i] + w.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[i] / Pc[i] + end + end for i in 1:length(model) pci,Tci,Vci = _pc[i],_Tc[i],_Vc[i] @@ -220,7 +223,7 @@ end function c_premixing(model::PatelTejaModel) _Tc = model.params.Tc _pc = model.params.Pc - _Vc = model.params.Vc + _Vc = model.params.Vc_fit c = model.params.c _Zc = _pc .* _Vc ./ (R̄ .* _Tc) Ωc = @. 1 - 3*_Zc From 1f7c2fcc52a5995b76a22e3f80fb8bfa0a9e68ab Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Tue, 30 Dec 2025 18:02:58 -0300 Subject: [PATCH 2/4] move estimation of patel teja Zc into c_premixing, fix tests --- src/models/cubic/PatelTeja/PatelTeja.jl | 36 ++++++++++++++----------- src/models/cubic/equations.jl | 2 +- test/test_methods_eos.jl | 4 +-- test/test_models_cubic.jl | 2 +- 4 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/models/cubic/PatelTeja/PatelTeja.jl b/src/models/cubic/PatelTeja/PatelTeja.jl index 5ed7fbe93a..33ab4742b2 100755 --- a/src/models/cubic/PatelTeja/PatelTeja.jl +++ b/src/models/cubic/PatelTeja/PatelTeja.jl @@ -28,6 +28,10 @@ function transform_params(::Type{PatelTejaParam},params,components) SingleParam("acentricfactor",components,zeros(Base.promote_eltype(Pc,Tc),n),fill(true,n)) end + if all(Vc.ismissingvalues) && all(acentricfactor.ismissingvalues) + throw(MissingException("PatelTeja: cannot estimate Vc: missing acentricfactor parameter.")) + end + return params end @@ -158,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","cubic/PatelTeja/PatelTeja_Vc_fit.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, @@ -174,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) @@ -196,21 +200,11 @@ function ab_premixing(model::PatelTejaModel,mixing::MixingRule,k,l) _Vc = model.params.Vc_fit a = model.params.a b = model.params.b - ω = model.params.acentricfactor - - for i in 1:length(model) - if _Vc.ismissingvalues[i] - w.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[i] / Pc[i] - end - end 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 @@ -224,10 +218,20 @@ function c_premixing(model::PatelTejaModel) _Tc = model.params.Tc _pc = model.params.Pc _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 diff --git a/src/models/cubic/equations.jl b/src/models/cubic/equations.jl index e436d83de1..25bb5347c6 100644 --- a/src/models/cubic/equations.jl +++ b/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] diff --git a/test/test_methods_eos.jl b/test/test_methods_eos.jl index 7298eb8cb6..cc679383f6 100755 --- a/test/test_methods_eos.jl +++ b/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] diff --git a/test/test_models_cubic.jl b/test/test_models_cubic.jl index 4035770ed4..2a14bae53e 100644 --- a/test/test_models_cubic.jl +++ b/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 From ba8349a05fa05bcf9c3ccdab16d5dc305e97c71f Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Tue, 30 Dec 2025 18:09:23 -0300 Subject: [PATCH 3/4] fix test later --- test/test_solvers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_solvers.jl b/test/test_solvers.jl index 5b391ecaf2..9ee451ae2a 100755 --- a/test/test_solvers.jl +++ b/test/test_solvers.jl @@ -199,7 +199,7 @@ end @test x3 ≈ xs_1v x4 = SOL.solution(SOL.optimize(minlog,(1.5,2.5),SOL.BoundOptim1Var())) - @test x4 ≈ xs_1v + #@test x4 ≈ xs_1v end end @printline From e510050724300fb4d7d93b5c0bb82f46c3a90b35 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Tue, 30 Dec 2025 19:25:13 -0300 Subject: [PATCH 4/4] fix test --- test/test_solvers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_solvers.jl b/test/test_solvers.jl index 9ee451ae2a..5b391ecaf2 100755 --- a/test/test_solvers.jl +++ b/test/test_solvers.jl @@ -199,7 +199,7 @@ end @test x3 ≈ xs_1v x4 = SOL.solution(SOL.optimize(minlog,(1.5,2.5),SOL.BoundOptim1Var())) - #@test x4 ≈ xs_1v + @test x4 ≈ xs_1v end end @printline