support TPD with activity coefficients

NEWS.md

@@ -17,6 +17,8 @@
 - calculation of volumes,saturation pressures and critical points of CPA models now defaults to the inner cubic model when there is no association present.
 - The default association implementation now uses a combination of accelerated successive substitution and newton optimization. While increasing allocations, the method is faster.
 - the default `volume` implementation now uses implicit AD to support derivatives. instead of propagating derivative information through the iterative procedure. This allows workloads of the type: `ForwardDiff.derivative(_p -> property(model,_p,T,z,phase = :l,vol0 = v0),p)` to be efficiently calculated.
+- `Clapeyron.tpd` code has been optimized. `tpd` has new keywords: `break_first`, that tries to return a negative tpd as early as possible, `lle` for only calculating TPD in liquid phases, `strategy`, that changes the search strategy between a K-value search (`:wilson`), a pure component search (`:pure`) or both strategies (`default`).
+- `Clapeyron.tpd` now supports activity models (if the keyword `lle` is set to `true`)
 - New EoS: modified Lee-Kesler-Plöcker with consistent parameters (`LKPmod`)
 - New EoS: Lee-Kesler-Plöker-equation of state, Sabozin-Jäger-Thol enhancement (`LKPSJT`, `enhancedLKP`)
 ## Bug fixes

src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl

@@ -205,7 +205,7 @@ function tp_flash_multi(model,p,T,z,options = MultiPhaseTPFlash())
         δn_add = true
         _result = (comps, βi, volumes, idx_vapour)
     elseif options.full_tpd #calculate full tpd.
-        comps,tpds,_,phase_w = tpd(model,p,T,z)
+        comps,tpds,_,phase_w = tpd(model,p,T,z,strategy = :pure)
         βi = initial_beta!(comps,z)
         volumes = similar(βi)
         n_phases = length(comps)

src/methods/tpd.jl

@@ -4,7 +4,6 @@ function mixture_fugacity(model,p,T,z;phase = :unknown, threaded = true, vol0 =
     return mixture_fugacity!(f,model,p,T,z;phase,threaded,vol0)
 end
 
-
 function mixture_fugacity!(f,model,p,T,z;phase = :unknown, threaded = true, vol0 = nothing)
     fugacity_coefficient!(f,model,p,T,z;phase,threaded,vol0)
     f .= f .* p .* z
@@ -133,6 +132,13 @@ end
 struct TPDKSolver end
 struct TPDPureSolver end
 
+function _tpd_and_v!(fxy,model,p,T,w,di,phase = :l)
+    lnϕw,v = lnϕ!(fxy,model,p,T,w,phase = phase)
+    tpd = @sum(w[i]*(lnϕw[i] + log(w[i]) - di[i])) - sum(w) + 1
+    return tpd,v
+end
+
+
 """
     wl,wv = tpd_solver(model,p,T,z,K0)
 
@@ -193,13 +199,12 @@ function tpd_solver(model,p,T,z,K0,
     fxy = ss_cache[3]
 
     if !stable_l
-        lnϕwl,vl = lnϕ!(fxy,model,p,T,wl,phase = :l)
-        tpd_l = @sum(wl[i]*(lnϕwl[i] + log(wl[i]) - di[i])) - sum(wl) + 1
+        tpd_l,vl = _tpd_and_v!(fxy,model,p,T,wl,di,:l)
         tpd_l < 0 && break_first && return wl,wv,tpd_l,tpd_v,vl,vv
     end
 
     if !stable_v
-        lnϕwv,vv = lnϕ!(fxy,model,p,T,wv,phase = :v)
+        tpd_v,vv = _tpd_and_v!(fxy,model,p,T,wv,di,:v)
         tpd_v = @sum(wv[i]*(lnϕwv[i] + log(wv[i]) - di[i])) - sum(wv) + 1
         tpd_v < 0 && break_first && return wl,wv,tpd_l,tpd_v,vl,vv
     end
@@ -269,18 +274,37 @@ function tpd_ss!(model,p,T,z,K0,is_liquid,solver = TPDKSolver(),cache = tpd_ss_c
     _tpd_ss!(model,p,T,z,K0,solver,is_liquid,cache,tol_equil,tol_trivial,maxiter)
 end
 
+function _tpd_fz_and_v!(solver::TPDKSolver,fxy,model,p,T,w,v0,liquid_overpressure = false,phase = :l)
+    v = volume(model,p,T,w,phase = phase,vol0 = v0)
+    if isnan(v) && liquid_overpressure && is_liquid(phase)
+        #michelsen recomendation: when doing tpd, sometimes, the liquid cannot be created at the 
+        #specified conditions. try elevating the pressure at the first iter.
+        v = volume(model,1.2p,T,w,phase = phase)
+    end
+    VT_mixture_fugacity!(fxy,model,v,T,w)
+    return fxy,v,true
+end
+
+function _tpd_fz_and_v!(solver::TPDPureSolver,fxy,model,p,T,w,v0,liquid_overpressure = false,phase = :l)
+    lnϕw,v = lnϕ!(fxy,model,p,T,w,phase = phase,vol0 = v0)
+    if isnan(v) && liquid_overpressure && is_liquid(phase)
+        #michelsen recomendation: when doing tpd, sometimes, the liquid cannot be created at the 
+        #specified conditions. try elevating the pressure at the first iter.
+        lnϕw,v = lnϕ!(fxy,model,1.2p,T,w,phase = phase)
+    end
+    return fxy,v,true
+end
 #tpd K-value solver
 #a port from MultiComponentFlash.jl
 function _tpd_ss!(model,p,T,z,K0,solver::TPDKSolver,is_liquid,cache,tol_equil, tol_trivial,maxiter)
-
     #phase of the trial composition
     phase = is_liquid ? :l : :v
-
     #is this a trivial solution?
     trivial = false
     S = 1.0
     iter = 0
     done = false
+    liquid_overpressure = false #if the liquid overpressure strategy was used
     K,fz,fw,wl,wv = cache
     w = is_liquid ? wl : wv
     K .= K0
@@ -294,8 +318,7 @@ function _tpd_ss!(model,p,T,z,K0,solver::TPDKSolver,is_liquid,cache,tol_equil, t
             S += w_i
         end
         @. w /= S
-        v = volume(model,p,T,w,phase = phase,vol0 = v)
-        VT_mixture_fugacity!(fw,model,v,T,w)
+        _,v,liquid_overpressure = _tpd_fz_and_v!(solver,fw,model,p,T,w,v,liquid_overpressure,phase)
         R_norm = zero(eltype(K))
         K_norm = zero(eltype(K))
         @inbounds for c in eachindex(K)
@@ -343,15 +366,7 @@ function _tpd_ss!(model,p,T,z,w0,solver::TPDPureSolver,_is_liquid,cache,tol_equi
     S = zero(eltype(w))
     while !done
         iter += 1
-        lnϕw, v = lnϕ!(lnϕw,model, p, T, w; phase=phase, vol0=v)
-        if isnan(v) && is_liquid(phase)
-            if !liquid_overpressure
-                #michelsen recomendation: when doing tpd, sometimes, the liquid cannot be created at the 
-                #specified conditions. try elevating the pressure at the first iter.
-                liquid_overpressure = true
-                lnϕw, v = lnϕ!(lnϕw,model, 1.2p, T, w; phase=phase, vol0=v)
-            end
-        end
+        _,v,liquid_overpressure = _tpd_fz_and_v!(solver,lnϕw,model,p,T,w,v,liquid_overpressure,phase)
         S = zero(eltype(w))
         for i in eachindex(w)
             wi = exp(di[i]-lnϕw[i])

src/models/Activity/equations.jl

@@ -223,4 +223,13 @@ function Obj_LLE(model::ActivityModel, F, T, x, xx)
     return F
 end
 
-export LLE
+export LLE
+
+function tpd(model::ActivityModel,p,T,z,cache = tpd_cache(model,p,T,z);reduced = false,break_first = false,lle = false,tol_trivial = 1e-5,strategy = :pure, di = nothing)
+    #TODO: support tpd with vle and activities?
+    if !lle
+        throw(ArgumentError("tpd only supports lle search with Activity Models. try using `tpd(model,p,T,z,lle = true)`"))
+    end
+    γϕmodel = __act_to_gammaphi(model,tpd,true)
+    return tpd(γϕmodel,p,T,z,cache;reduced,break_first,lle,tol_trivial,strategy,di)
+end

src/models/CompositeModel/GammaPhi.jl

@@ -210,4 +210,58 @@ function __eval_G_DETPFlash(wrapper::PTFlashWrapper{<:GammaPhi},p,T,x,equilibriu
     end
 end
 
+#=
+TPD support.
+
+TODO: support vle in TPD.
+=#
+
+function tpd_input_composition(model::GammaPhi,p,T,z,di,lle)
+    γ = activity_coefficient(model.activity,p,T,z)
+    #v = volume(model.fluid.model,p,T,z,phase = :l)
+    v = one(eltype(γ))
+    fz = γ .* p .* z
+    fz,:liquid,v
+end
+
+function _tpd_fz_and_v!(solver::TPDKSolver,fxy,model::GammaPhi,p,T,w,v0,liquid_overpressure = false,phase = :l)
+    #v = volume(model.fluid.model,p,T,w,phase = phase,vol0 = v0)
+    v = one(eltype(fxy))
+    fxy .= activity_coefficient(model.activity,p,T,w)
+    fxy .= fxy .* p .* w
+    return fxy,v,true
+end
+
+function _tpd_fz_and_v!(solver::TPDPureSolver,fxy,model::GammaPhi,p,T,w,v0,liquid_overpressure = false,phase = :l)
+    #v = volume(model.fluid.model,p,T,w,phase = phase,vol0 = v0)
+    fxy .= activity_coefficient(model.activity,p,T,w)
+    v = one(eltype(fxy))
+    fxy .= log.(fxy)
+    return fxy,v,true
+end
+
+function _tpd_and_v!(fxy,model::GammaPhi,p,T,w,di,phase = :l)
+    #v = volume(model.fluid.model,p,T,w,phase = phase)
+    v = one(eltype(fxy))
+    fxy .= activity_coefficient(model.activity,p,T,w)
+    fxy .= log.(fxy)
+    tpd = @sum(w[i]*(fxy[i] + log(w[i]) - di[i])) - sum(w) + 1
+    return tpd,v
+end
+
+function tpd_obj(model::GammaPhi, p, T, di, isliquid, cache = tpd_neq_cache(model,p,T,di,di), break_first = false)
+    vcache[] = one(eltype(di))
+    function f(α)
+        w = α .* α .* 0.25
+        w ./= sum(w)
+        γ = activity_coefficient(model.activity,p,T,w)
+        γ .= log(γ)
+        lnγw = γ
+        fx = @sum(w[i]*(lnγw[i] + log(w[i]) - di[i])) - sum(w) + 1
+    end
+
+    obj = Solvers.ADScalarObjective(f,di,ForwardDiff.Chunk{2}())
+    optprob = OptimizationProblem(obj = obj,inplace = true)
+end
+
 export GammaPhi

src/models/EmpiricHelmholtz/LKP/LKP.jl

@@ -117,7 +117,6 @@ function a_res(model::LKPModel,V,T,z = SA[1.0])
     δ = sum(z)*Vr/V
     τ = Tr/T
     δr = δ/Zr
-    @show δ,τ
     params_simple = lkp_params_simple(model)
     params_reference = lkp_params_reference(model)
     ω0,ωref = last(params_simple),last(params_reference)

src/modules/solvers/optimize.jl

@@ -1,5 +1,5 @@
 
-function ADScalarObjective(f,x0::AbstractArray,chunk = autochunk(x0),val::Val{N} = Val{2}()) where N
+function ADScalarObjective(f,x0::AbstractArray,chunk = autochunk(x0))
     Hres = DiffResults.HessianResult(x0)
     function _g(df,x,Gresult)
         ForwardDiff.gradient!(Gresult,f,x)

test/test_methods_api.jl

@@ -80,7 +80,7 @@ end
     fluid = PCSAFT(["water","methanol"]; assoc_options=AssocOptions(combining=:elliott))
     fluid.params.epsilon["water","methanol"] *= (1+0.18)
     v = volume(fluid, 1e5, 160.0, [0.5, 0.5],phase = :l)
-    @test Clapeyron.X(fluid,v,160.0[0.5,0.5]).v ≈ [0.0011693187791158642, 0.0011693187791158818, 0.0002916842981727242, 0.0002916842981727286] rtol = 1E-8
+    @test Clapeyron.X(fluid,v,160.0,[0.5,0.5]).v ≈ [0.0011693187791158642, 0.0011693187791158818, 0.0002916842981727242, 0.0002916842981727286] rtol = 1E-8
 end
 
 using EoSSuperancillaries
@@ -109,7 +109,11 @@ end
     T = 298.15
     p = 1e5
     phases,tpds,symz,symw = Clapeyron.tpd(system,p,T,[0.5,0.5])
-    @test tpds[1] ≈ -0.6081399681963373  rtol = 1e-6
+    @test tpds[1] ≈ -0.6081399681963373 rtol = 1e-6
+
+    act_system = UNIFAC(["water","cyclohexane"])
+    phases2,tpds2,symz2,symw2 = Clapeyron.tpd(act_system,p,T,[0.5,0.5],lle = true)
+    @test tpds2[1] ≈ -0.9412151812640561 rtol = 1e-6
     GC.gc()
 end