Spinodals

docs/src/properties/multi.md

@@ -63,6 +63,8 @@ Clapeyron.VT_mechanical_stability
 Clapeyron.VT_diffusive_stability
 Clapeyron.VT_chemical_stability
 Clapeyron.tpd
+Clapeyron.spinodal_pressure
+Clapeyron.spinodal_temperature
 ```
 
 ## TP Flash

ext/ClapeyronUnitfulExt.jl

@@ -207,6 +207,23 @@
     return uconvert(output, res)
 end
 
+function C.spinodal_pressure(model::EoSModel, T::Unitful.Temperature, x=SA[1.]; v0=nothing, phase=:unknown, output=(u"Pa",u"m^3"))
+    st = standardize(model,-1,T,x)
+    _,_T,_x = state_to_pt(model,st)
+    _v0 = v0===nothing ? nothing : total_volume(model, v0, x)
+    res = spinodal_pressure(model, _T, _x; v0=_v0, phase=phase).*(u"Pa",u"m^3")
+    return uconvert.(output, res)
+end
+
+function C.spinodal_temperature(model::EoSModel, p::Unitful.Pressure, x=SA[1.]; T0=nothing, v0=nothing, phase=:unknown, output=(u"K",u"m^3"))
+    st = standardize(model,p,-1,x)
+    _p,_,_x = state_to_pt(model,st)
+    _T0 = T0===nothing ? nothing : standardize(T0,1u"K")
+    _v0 = v0===nothing ? nothing : total_volume(model, v0, x)
+    res = spinodal_temperature(model, _p, _x; T0=_T0, v0=_v0, phase=phase).*(u"K",u"m^3")
+    return uconvert.(output, res)
+end
+
 # resolve ambiguity
 function C.chemical_potential(model::(C.SolidHfusModel), v::__VolumeKind, T::Unitful.Temperature, z=SA[1.]; output=u"J/mol")
     st = standardize(model,v,T,z)

src/methods/methods.jl

@@ -299,4 +299,4 @@ include("tpd.jl")
 include("stability.jl")
 include("pT.jl")
 include("property_solvers/Tproperty.jl")
-
+include("property_solvers/spinodal.jl")

src/methods/property_solvers/spinodal.jl

@@ -0,0 +1,86 @@
+
+"""
+    spinodal_pressure(model::EoSModel, T, x; v0, phase)
+
+Calculates the spinodal pressure and volume for a given temperature and composition. Returns a tuple, containing:
+- spinodal pressure [`Pa`]
+- spinodal volume [`m³`]    
+
+Calculates either the liquid or the vapor spinodal point depending on the given starting volume `v0` or the `phase`. The keyword `phase` is ignored if `v0` is given.
+"""
+function spinodal_pressure(model::EoSModel,T,x=SA[1.];v0=nothing,phase=:unknown)
+    x = x/sum(x)
+    model, idx_r = index_reduction(model,x)
+    x = x[idx_r]
+
+    # Determine initial guess (if not provided)
+    if isnothing(v0)
+        if is_liquid(phase)
+            v0 = bubble_pressure(model,T,x)[2]
+        elseif is_vapour(phase)
+            v0 = dew_pressure(model,T,x)[3]
+        else
+            error("Either `v0` or `phase` has to be specified!")
+        end
+    end
+
+    # Solve spinodal condition
+    f!(F,vz) = det_∂²A∂ϱᵢ²(model,F,exp10(vz[1]),T,x)
+    r = Solvers.nlsolve(f!,[log10(v0)],LineSearch(Newton()),NEqOptions(), ForwardDiff.Chunk{1}())
+    V_spin = exp10(Solvers.x_sol(r)[1])
+
+    if r.info.best_residual[1] < r.options.f_abstol # converged
+        return pressure(model,V_spin,T,x), V_spin
+    else                                            # not converged
+        return NaN, NaN
+    end
+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 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.
+"""
+function spinodal_temperature(model::EoSModel,p,x=SA[1.];T0=nothing,v0=nothing,phase=:unknown)
+    x = x/sum(x)
+    model, idx_r = index_reduction(model,x)
+    x = x[idx_r]
+
+    # Determine initial guess (if not provided)
+    if isnothing(T0) || isnothing(v0)
+        if is_liquid(phase)
+            Tv0 = bubble_temperature(model,p,x)[[1,2]]
+        elseif is_vapour(phase)
+            Tv0 = dew_temperature(model,p,x)[[1,3]]
+        else
+            error("Either `T0` and `v0` or `phase` have to be specified!")
+        end
+        T0 = isnothing(T0) ? Tv0[1] : T0
+        v0 = isnothing(v0) ? Tv0[2] : v0
+    end
+
+    # Solve spinodal condition
+    f!(F,Tz) = det_∂²A∂ϱᵢ²(model, F, volume(model, p, Tz[1], x; phase=phase, vol0=v0), Tz[1], x)
+    r = Solvers.nlsolve(f!,[T0],LineSearch(Newton()),NEqOptions(;f_abstol=1e-6), ForwardDiff.Chunk{1}())
+    T_spin = Solvers.x_sol(r)[1]
+
+    if all(r.info.best_residual .< r.options.f_abstol)  # converged
+        return T_spin, volume(model, p, T_spin, x; phase=phase, vol0=v0)
+    else                                                # not converged
+        return NaN, NaN
+    end
+end
+
+# Objective function for spinodal calculation -> det(∂²A/∂ϱᵢ) = 0
+function det_∂²A∂ϱᵢ²(model,F,v,T,x)
+    # calculates det(∂²A∂xᵢ² ⋅ ϱ) at V,T constant (see www.doi.org/10.1016/j.fluid.2017.04.009)
+    Av = ϱi -> eos(model, v, T, v.*ϱi)./v
+    F[1] = det(ForwardDiff.hessian(Av,x./v))
+    return F
+end
+
+export spinodal_pressure, spinodal_temperature

test/test_methods_api.jl

@@ -655,3 +655,18 @@ GC.gc()
     GC.gc()
 end
 
+@testset "spinodals" begin
+    # Example from Ref. https://doi.org/10.1016/j.fluid.2017.04.009
+    model = PCSAFT(["methane","ethane"])
+    T_spin = 223.
+    x_spin = [0.2,0.8]
+    (pl_spin, vl_spin) = spinodal_pressure(model,T_spin,x_spin;phase=:liquid)
+    (pv_spin, vv_spin) = spinodal_pressure(model,T_spin,x_spin;phase=:vapor)
+    @test vl_spin ≈ 7.218532167482202e-5 rtol = 1e-6
+    @test vv_spin ≈ 0.0004261109817247137 rtol = 1e-6
+
+    (Tl_spin_impl, xl_spin_impl) = spinodal_temperature(model,pl_spin,x_spin;T0=220.,v0=vl_spin)
+    (Tv_spin_impl, xv_spin_impl) = spinodal_temperature(model,pv_spin,x_spin;T0=225.,v0=vv_spin)
+    @test Tl_spin_impl ≈ T_spin rtol = 1e-6
+    @test Tv_spin_impl ≈ T_spin rtol = 1e-6
+end