diff --git a/HISTORY.md b/HISTORY.md index d81e99266..cc99a4ec9 100644 --- a/HISTORY.md +++ b/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 diff --git a/NEWS.md b/NEWS.md index 81ed237ce..44f2b3a72 100644 --- a/NEWS.md +++ b/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. \ No newline at end of file +## Bug fixes +- fixes in calculation of spinodal with cubics. +- `MultiFluid` and `SingleFluid` errors when T_reducing != Tc. +- fix `VT_identify_phase`. diff --git a/src/methods/initial_guess.jl b/src/methods/initial_guess.jl index 5fab90812..ba640aba3 100755 --- a/src/methods/initial_guess.jl +++ b/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,7 +960,7 @@ 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) @@ -966,9 +968,9 @@ function critical_tsat_extrapolation(model,p,Tc,Pc,Vc) 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) diff --git a/src/methods/property_solvers/multicomponent/flash.jl b/src/methods/property_solvers/multicomponent/flash.jl new file mode 100644 index 000000000..8b85a5a7e --- /dev/null +++ b/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 + diff --git a/src/methods/property_solvers/multicomponent/multicomponent.jl b/src/methods/property_solvers/multicomponent/multicomponent.jl index 04acf2768..d0eb0bcf4 100755 --- a/src/methods/property_solvers/multicomponent/multicomponent.jl +++ b/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 diff --git a/src/methods/property_solvers/spinodal.jl b/src/methods/property_solvers/spinodal.jl index 10544334b..3878c4adb 100644 --- a/src/methods/property_solvers/spinodal.jl +++ b/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. diff --git a/src/models/cubic/equations.jl b/src/models/cubic/equations.jl index 92a4d5b0f..363b52088 100644 --- a/src/models/cubic/equations.jl +++ b/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) diff --git a/test/test_methods_api.jl b/test/test_methods_api.jl index 337803034..e181dcad8 100644 --- a/test/test_methods_api.jl +++ b/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