fix saturation_temperature calculations (#290)

HISTORY.md

@@ -1,3 +1,21 @@
+# v0.6.4
+
+## New Features
+
+- New model: SAFT-VR-Mie with Gross-Vrabec quatrupolar contribution (`SAFTVRMieGV`)
+- New model: Co-Oriented Fluid Functional Equation for Electrostatic interactions (`COFFEE`)
+- Better support for evaluation of model properties at V == Inf (ideal gas limit)
+- New method: `adiabatic_index`, that calculates the ratio between the isobaric and isochoric heat capacities.
+- new API: `has_fast_crit_pure`, to indicate that models can calculate their pure critical point quickly. saturation initial guesses use the result of this function to decide if and when to call the `crit_pure` routine.
+- speed ups in some pressure routines
+-
+## Bug fixes
+
+- `MultiFluid` and `SingleFluid` models did not use the correct gas constant.
+- Fix mixing rule in `SAFTVRMie`.
+- `VT_identify_phase` now returns `:unknown` for an unstable state input.
+- Typos in `TProperty` for pure models.
+
 # v0.6.3
 
 ## New Features

NEWS.md

@@ -1,17 +1,16 @@
-# v0.6.4
+# v0.6.5
 
 ## New Features
 
-- New model: SAFT-VR-Mie with Gross-Vrabec quatrupolar contribution (`SAFTVRMieGV`)
-- New model: Co-Oriented Fluid Functional Equation for Electrostatic interactions (`COFFEE`)
-- Better support for evaluation of model properties at V == Inf (ideal gas limit)
-- New method: `adiabatic_index`, that calculates the ratio between the isobaric and isochoric heat capacities.
-- new API: `has_fast_crit_pure`, to indicate that models can calculate their pure critical point quickly. saturation initial guesses use the result of this function to decide if and when to call the `crit_pure` routine.
-- speed ups in some pressure routines
--
-## Bug fixes
+- Experimental: Bulk properties for Pressure-Enthalpy and Pressure-Entropy, the syntax is the following:
+  ```julia
+  using Clapeyron: PH
+  PH.entropy(model,p,h,z)
+  PH.adiabatic_index(model,p,h,z,T0 = T0) #suplying an initial point for the temperature
+  ```
+  The calculation is done via `Clapeyron.Tproperty`. there are also `PT` and `VT` functions for parity.
 
-- `MultiFluid` and `SingleFluid` models did not use the correct gas constant.
-- Fix mixing rule in `SAFTVRMie`.
-- `VT_identify_phase` now returns `:unknown` for an unstable state input.
-- Typos in `TProperty` for pure models.
+## Bug fixes
+- fixes in calculation of spinodal with cubics.
+- `MultiFluid` and `SingleFluid` errors when T_reducing != Tc.
+- fix `VT_identify_phase`.

src/methods/initial_guess.jl

@@ -578,7 +578,6 @@ function x0_sat_pure_spinodal(model,T,v_lb,v_ub,B = second_virial_coefficient(mo
     end
 
     psl = p(vsl)
-
     if isnan(vsl)
         return pressure(model,v_lb,T),v_lb,v_ub
     end
@@ -731,19 +730,29 @@ function x0_saturation_temperature end
 function x0_saturation_temperature(model,p)
     single_component_check(x0_saturation_temperature,model)
     coeffs = antoine_coef(model)
-    if coeffs !== nothing
-        return x0_saturation_temperature_antoine_coeff(model,p,coeffs)
-    end
+    #@show coeffs
+    #if coeffs !== nothing
+    #    return x0_saturation_temperature_antoine_coeff(model,p,coeffs)
+    #end
     #if obtaining the critical point is cheap, models can opt in by defining:
     #=
     x0_saturation_temperature(model::MyModel,p) = x0_saturation_temperature(model,p,crit_pure(model))
     =#
-    return x0_saturation_temperature(model,p,nothing)
+    if has_fast_crit_pure(model)
+
+        return x0_saturation_temperature_refine(model,p)
+    else
+        return x0_saturation_temperature_crit(model,p,crit_pure(model))
+    end
 end
 
 function x0_saturation_temperature(model::EoSModel,p,::Nothing)
     single_component_check(x0_saturation_temperature,model)
-    return x0_saturation_temperature_refine(model,p)
+    if has_fast_crit_pure(model)
+        return x0_saturation_temperature_crit(model,p,crit_pure(model))
+    else
+        return x0_saturation_temperature_refine(model,p)
+    end
 end
 
 function x0_saturation_temperature(model::EoSModel,p,crit::Tuple)
@@ -760,14 +769,7 @@ function x0_saturation_temperature_antoine_coeff(model,p,coeffs)
 end
 
 function x0_saturation_temperature_crit(model::EoSModel,p,crit)
-    Tc,Pc,Vc = crit
-    A,B,C = (6.668322465137264,6.098791871032391,-0.08318016317721941) #universal antoine constants (RK)
-    if Pc < p
-        nan = zero(p)/zero(p)
-        return (nan,nan,nan)
-    end
-    lnp̄ = log(p / Pc)
-    T0 = Tc*(B/(A-lnp̄)-C)
+    T0 = critical_tsat_extrapolation(model,p,crit)
     return x0_saturation_temperature_refine(model,p,T0)
 end
 
@@ -958,17 +960,17 @@ function critical_tsat_extrapolation(model,p,Tc,Pc,Vc)
         _0 = zero(Base.promote_eltype(model,p))
         return _0/_0
     end
-    _p(_T) = pressure(pure,Vc,_T)
+    _p(_T) = pressure(model,Vc,_T)
     dpdT = Solvers.derivative(_p,Tc)
     dTinvdlnp = -Pc/(dpdT*Tc*Tc)
     Δlnp = log(p/Pc)
     Tinv = 1/Tc + dTinvdlnp*Δlnp
     T = 1/Tinv
 end
 
-critical_tsat_extrapolation(model,T) = critical_tsat_extrapolation(model,p,crit_pure(model))
-critical_tsat_extrapolation(model,T,crit) = critical_tsat_extrapolation(model,p,crit[1],crit[2],crit[3])
-critical_tsat_extrapolation(model,T,Tc,Vc) = critical_tsat_extrapolation(model,p,Tc,pressure(model,Vc,Tc),Vc)
+critical_tsat_extrapolation(model,p) = critical_tsat_extrapolation(model,p,crit_pure(model))
+critical_tsat_extrapolation(model,p,crit) = critical_tsat_extrapolation(model,p,crit[1],crit[2],crit[3])
+critical_tsat_extrapolation(model,p,Tc,Vc) = critical_tsat_extrapolation(model,p,Tc,pressure(model,Vc,Tc),Vc)
 
 function dpdT_pure(model,v1,v2,T)
     #log(p/p0) = [-dpdT*T*T/p](p = p0,T = T0) * (1/T - 1/T0)

src/methods/property_solvers/multicomponent/flash.jl

@@ -0,0 +1,13 @@
+function flash(specifier,model,v1,v2,z,args...;kwargs...)
+
+end
+
+
+function ph_flash(model,p,h,z,T0 = Tproperty(model,p,h,z))
+
+end
+
+function ps_flash(model,p,s,z,T0 = Tproperty(model,p,s,z,entropy))
+
+end
+

src/methods/property_solvers/multicomponent/multicomponent.jl

@@ -262,7 +262,7 @@ include("krichevskii_parameter.jl")
 include("solids/sle_solubility.jl")
 include("solids/slle_solubility.jl")
 include("solids/eutectic_point.jl")
-
+include("flash.jl")
 export bubble_pressure_fug, bubble_temperature_fug, dew_temperature_fug, dew_pressure_fug
 export bubble_pressure,    dew_pressure,    LLE_pressure,    azeotrope_pressure, VLLE_pressure
 export bubble_temperature, dew_temperature, LLE_temperature, azeotrope_temperature, VLLE_temperature

src/methods/property_solvers/spinodal.jl

@@ -52,7 +52,7 @@ end
     spinodal_temperature(model::EoSModel, p, x; T0, v0, phase)
 
 Calculates the spinodal pressure and volume for a given pressure and composition. Returns a tuple, containing:
-- spinodal temperataure [`K`]
+- spinodal temperature [`K`]
 - spinodal volume [`m³`]    
 
 Calculates either the liquid or the vapor spinodal point depending on the given starting temperature `T0` and volume `v0` or the `phase`. The keyword `phase` is ignored if `T0` or `v0` is given.

src/models/cubic/equations.jl

@@ -312,17 +312,10 @@ function pure_spinodal(model::ABCubicModel,T::K,v_lb::K,v_ub::K,phase::Symbol,re
     Q1 = 2*b*b*(a*(1 - c1_c2) - bRT*c1c2*c1_c2)
     Q0 = b*b*b*(a*c1_c2 - bRT*c1c2*c1c2)
     dpoly = (Q0,Q1,Q2,Q3,Q4)
-    dfl = evalpoly(v_lb,dpoly)
-    dfv = evalpoly(v_ub,dpoly)
-    v_lbc = v_lb + c
-    v_ubc = v_ub + c
-    vx = ifelse(is_liquid(phase),v_lbc,v_ubc)
-    if dfl*dfv < 0
-        bracket = (v_lbc,v_ubc)
-        prob = Roots.ZeroProblem(Base.Fix2(evalpoly,dpoly),bracket)
-        vs = Roots.solve(prob)
-        return vs - c
-    end
+    v_lbc = v_lb + c #untranslated lb
+    v_ubc = v_ub + c #untranslated ub
+    dfl = evalpoly(v_lbc,dpoly)
+    dfv = evalpoly(v_ubc,dpoly)
     vx = ifelse(is_liquid(phase),v_lbc,v_ubc)
     vm = vcc #the critical volume is always between the two spinodals
     v_bracket = minmax(vx,vm)
@@ -358,7 +351,7 @@ function ab_consts(model::CubicModel)
     return ab_consts(typeof(model))
 end
 
-has_fast_crit_pure(model::ABCCubicModel) = true
+has_fast_crit_pure(model::ABCubicModel) = true
 
 function x0_saturation_temperature(model::ABCubicModel,p,::Nothing)
     crit = crit_pure(model)

test/test_methods_api.jl

@@ -410,6 +410,9 @@ end
 
     #Issue #290
     @test Clapeyron.saturation_temperature(cPR("R1233zde"),101325*20,crit_retry = false)[1] ≈ 405.98925205830335 rtol = 1e-6
+    @test Clapeyron.saturation_temperature(cPR("isobutane"),1.7855513185537157e6,crit_retry = false)[1] ≈ 366.52386488214876 rtol = 1e-6
+    @test Clapeyron.saturation_temperature(cPR("propane"),2.1298218093361156e6,crit_retry = false)[1] ≈ 332.6046103831853 rtol = 1e-6
+    @test Clapeyron.saturation_temperature(cPR("r134a"),2.201981727901889e6,crit_retry = false)[1] ≈ 344.6869001549851 rtol = 1e-6
 end
 
 @testset "Tproperty" begin