tp_flash seems to work (with initial points)

Author: longemen3000

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

@@ -1,4 +1,4 @@
-function tp_flash_michelsen(model::ElectrolyteModel, p, T, z, method = MichelsenTPFlash(),reduced = false)
+function tp_flash_michelsen(model::ESElectrolyteModel, p, T, z, method = MichelsenTPFlash(),reduced = false)
 
     equilibrium = method.equilibrium
     K0 = method.K0
@@ -329,7 +329,7 @@ function rachfordrice(K, z, Z; β0=nothing, ψ0=nothing, non_inx=FillArrays.Fill
     end
 end
 
-function michelsen_optimization_of!(g,H,model::ElectrolyteModel,p,T,z,caches,ny_var_and_ψ,gz)
+function michelsen_optimization_of!(g,H,model::ESElectrolyteModel,p,T,z,caches,ny_var_and_ψ,gz)
     ny_var = @view ny_var_and_ψ[1:end-1]
     ψ = ny_var[end]
 

src/models/Electrolytes/ISElectrolyte.jl

@@ -53,18 +53,13 @@ Rgas(model::ISElectrolyteIdealWrapper) = Rgas(model.model)
 
 function eos_impl(model::ISElectrolyteWrapper,V,T,z)
     w = to_ion(model.salt,z)
-    ∑z = sum(z)
-    return ∑z*Rgas(model)*T*eos_impl(model,V,T,w)
+    return eos_impl(model.model,V,T,w)
 end
 
 function tp_flash_K0!(K,model::ISElectrolyteModel,p,T,z)
-    neutral = zeros(Bool,length(model))
-    r = eachrow(model.salt.mat)
-    for i in 1:length(model)
-        if count(!iszero,r[i]) == 1
-            neutral[i] = true
-        end
-    end
+    neutral = ones(Bool,length(model))
+    isalts = model.salt.isalts
+    neutral[isalts] .= false
     pures = split_model(model,neutral)
     psat = first.(extended_saturation_pressure.(pures,T))
     K .= 0
@@ -79,3 +74,27 @@ function each_split_model(model::ISElectrolyteWrapper,I_salt)
 end
 
 export ISElectrolyteWrapper
+
+function volume_impl(model::ISElectrolyteWrapper, p, T, z, phase, threaded, vol0)
+    w = to_ion(model.salt,z)
+    return volume(model.model, p, T, w, phase=phase, threaded=threaded, vol0=vol0)
+end
+
+function lb_volume(model::ISElectrolyteWrapper,T,z)
+    w = to_ion(model.salt,z)
+    return lb_volume(model.model,T,w)
+end
+
+function mw(model::ISElectrolyteWrapper)
+    return mw(model.neutralmodel)
+end
+
+function p_scale(model::ISElectrolyteWrapper,z)
+    w = to_ion(model.salt,z)
+    return p_scale(model.model,w)
+end
+
+function T_scale(model::ISElectrolyteWrapper,z)
+    w = to_ion(model.salt,z)
+    return T_scale(model.model,w)
+end

src/models/Electrolytes/SaltParam.jl

@@ -2,6 +2,7 @@ struct SaltParam <: ClapeyronParam
     explicit_solvent::Bool
     explicit_components::Vector{String} #neutrals + ions
     implicit_components::Vector{String} #neutrals + salts (ion)
+    isalts::Vector{Int} #Indices of the salts
     mat::Matrix{Float64} #used to calculate z(ion) -> z(salt)
     F::LU{Float64, Matrix{Float64}, Vector{Int64}} #used to calculate z(salt) -> z(ion)
 end
@@ -12,6 +13,7 @@ function explicit_salt_param(comps,salts,Z)
     nions = length(Z)
     implicit_components = Vector{String}(undef,nions - 1)
     mat = zeros(nions,nions)
+    isalts = Int[]
     nneutral = count(iszero,Z)
     #we suppose that first there are nneutral neutral components, followed by nions - nneutral ions
     for i in 1:nneutral
@@ -25,6 +27,7 @@ function explicit_salt_param(comps,salts,Z)
         salt = salts[k]
         salt_component = first(salt)
         implicit_components[i] = first(salt)
+        push!(isalts,i)
         pairings = last(salt)
         ri = rr[i]
         for ion_vals in pairings
@@ -40,8 +43,10 @@ function explicit_salt_param(comps,salts,Z)
     ∑ri = count(!iszero,ri)
     ri ./=  ∑ri
     end
-    rr[end] .= Z
-    return SaltParam(explicit_solvent,explicit_components,implicit_components,mat,lu(mat))
+    if nneutral < nions
+        rr[end] .= Z
+    end
+    return SaltParam(explicit_solvent,explicit_components,implicit_components,isalts,mat,lu(mat))
 end
 
 function SaltParam(model::ESElectrolyteModel)
@@ -64,7 +69,7 @@ function to_salt(m::SaltParam,z)
 end
 
 function to_ion(m::SaltParam,z)
-    if size(m.mat,1) == length(z)
+    if length(m.isalts) == 0
         zz = similar(z)
         zz .= z
     else
@@ -95,7 +100,16 @@ function IS_each_split_model(salt::SaltParam,I_salt)
         I_salt_plus_charge = vcat(I_salt,nions)
         mm = m[I_salt_plus_charge,I_ion_int]
     end
-    split_salt = SaltParam(false,salt.explicit_components[I_ion_int],salt.implicit_components[I_salt],mm,lu(mm))
+
+    isalts = Int[]
+    for i in 1:length(salt.isalts)
+        si = salt.isalts[i]
+        if salt.isalts[i] in I_salt
+            push!(isalts,si)
+        end
+    end
+
+    split_salt = SaltParam(false,salt.explicit_components[I_ion_int],salt.implicit_components[I_salt],isalts,mm,lu(mm))
     return split_salt,I_ion_int
 end
 

src/models/Electrolytes/electrolytes.jl

@@ -357,32 +357,30 @@ function a_res(ionmodel::IonModel, V, T, z, iondata, neutralmodel, neutral_data)
     return a_res(ionmodel, V, T, z, iondata)
 end
 
-function lb_volume(model::ElectrolyteModel,T,z)
+function lb_volume(model::ESElectrolyteModel,T,z)
     return lb_volume(model.neutralmodel,T,z)
 end
 
-function x0_volume_liquid(model::ElectrolyteModel,p,T,z)
+function x0_volume_liquid(model::ESElectrolyteModel,p,T,z)
     return x0_volume_liquid(model.neutralmodel,p,T,z)*1.15
 end
 
-function x0_volume_gas(model::ElectrolyteModel,p,T,z)
+function x0_volume_gas(model::ESElectrolyteModel,p,T,z)
     return x0_volume_gas(model.neutralmodel,p,T,z)
 end
 
-function mw(model::ElectrolyteModel)
+function mw(model::ESElectrolyteModel)
     return mw(model.neutralmodel)
 end
 
-function p_scale(model::ElectrolyteModel,z)
+function p_scale(model::ESElectrolyteModel,z)
     return p_scale(model.neutralmodel,z)
 end
 
-function T_scale(model::ElectrolyteModel,z)
+function T_scale(model::ESElectrolyteModel,z)
     return T_scale(model.neutralmodel,z)
 end
 
-
-
 function debye_length(model::ESElectrolyteModel,V,T,z,ϵ_r = @f(dielectric_constant),∑z = sum(z))
     Z = model.charge
     return debye_length(V,T,z,ϵ_r,Z)