more work on MultiPhaseTPFlash

Author: longemen3000

src/methods/property_solvers/fugacity_coefficient.jl

@@ -333,18 +333,13 @@ function modified_lnϕ(model, p, T, z, cache; phase = :unknown, vol0 = nothing)
 end
 
 function modified_gibbs(model,p,T,w,phase = :unknown,vol = NaN)
-    TT = Base.promote_eltype(model,p,T,w)
-    RT = Rgas(model)*T
-    ∑w = sum(w)
-    iszero(∑w) && return zero(TT), zero(TT)
-    g_ideal = sum(xlogx,w) - xlogx(∑w)
     if isnan(vol)
         volw = volume(model,p,T,w,phase = phase)
     else
         volw = vol
     end
-    ∑zlogϕi,vv = ∑zlogϕ(model,p,T,w,phase = phase,vol = volw)
-    return ∑zlogϕi + g_ideal,vv
+    g =  VT_gibbs_energy(model,volw,T,w,p) #+ eos_g(BasicIdeal(),p,T,w)
+    return g,volw
 end
 
 function modified_∂lnϕ∂n(model, p, T, z, cache; phase = :unknown, vol0 = nothing)

src/methods/property_solvers/multicomponent/flash.jl

@@ -348,20 +348,12 @@ for prop in [:isochoric_heat_capacity, :isobaric_heat_capacity, :adiabatic_index
     end
 end
 
-function _multiphase_gibbs(model,p,T,result)
-    if result isa FlashResult
-        return _multiphase_gibbs(model,p,T,(result.compositions,result.fractions,result.volumes))
-    end
-    if model isa PTFlashWrapper && length(β) == 2 #TODO: in the future, PTFlashWrappers could be multiphase
-        if model.model isa RestrictedEquilibriaModel
-            return zero(eltype(β))
-        end
-    end
-    comps = result[1]
-    β = result[2]
-    volumes = result[3]
-    data = FlashResult(model,p,T,comps,β,volumes)
-    return Rgas(model)*T*data.data.g
+function _multiphase_gibbs(model,result,vapour_phase_index = 0)
+    gibbs_energy(model,result)/Rgas(model)/result.data.T
+end
+
+function _multiphase_gibbs(model::PTFlashWrapper,result,vapour_phase_index = 0)
+    modified_gibbs(model,result;vapour_phase_index)/Rgas(model)/result.data.T
 end
 
 function __mpflash_phase(vapour_phase_index,i)