Søreide-Whitson EOS (#25)

Author: moyner

.gitignore

@@ -3,3 +3,4 @@
 *.jl.mem
 /docs/build/
 Manifest.toml
+*.json

docs/src/examples/advanced.md

@@ -31,7 +31,7 @@ If many flashes of the same mixture are to be performed at different conditions,
 m = SSIFlash()
 K = zeros(number_of_components(eos))
 S = flash_storage(eos, conditions, method = m)
-@allocated V = flash_2ph!(S, K, eos, conditions, method = m)
+@allocated flash_2ph!(S, K, eos, conditions, method = m)
 
 # output
 

src/MultiComponentFlash.jl

@@ -4,6 +4,8 @@ module MultiComponentFlash
     # Constants.
     const MINIMUM_COMPOSITION = 1e-10
     const IDEAL_GAS_CONSTANT = 8.3144598
+    @enum COMPONENT_ENUM COMPONENT_N2 COMPONENT_H2O COMPONENT_CO2 COMPONENT_H2S COMPONENT_OTHER_COMPONENT
+
     # Load types first
     include("mixture_types.jl")
     include("eos_types.jl")

src/eos.jl

@@ -8,6 +8,18 @@ Return number of components for the underlying mixture of the EOS.
 """
 number_of_components(e::AbstractEOS) = number_of_components(e.mixture)
 
+forces_per_phase(eos::GenericCubicEOS) = false
+
+function get_phase(cond)
+    return get(cond, :phase, :unknown)::Symbol
+end
+
+function set_phase(cond, phase::Symbol, throw::Bool = false)
+    if throw && haskey(cond, :phase) && cond.phase != :unknown
+        throw(ArgumentError("Phase state already set to $(cond.phase), cannot change to $phase."))
+    end
+    return (p = cond.p, T = cond.T, z = cond.z, phase = phase)
+end
 
 function static_coefficients(::AbstractGeneralizedCubic)
     return (0.4274802327, 0.08664035, 0.0, 1.0)
@@ -23,67 +35,32 @@ function weight_ai(eos::GenericCubicEOS{T, R}, cond, i) where {T<:AbstractGenera
     return eos.ω_a*T_r^(-1/2)
 end
 
-# PengRobinson specialization
-function static_coefficients(::AbstractPengRobinson)
-    return (0.457235529, 0.077796074, 1.0 + sqrt(2), 1.0 - sqrt(2))
-end
-
-function weight_ai(eos::GenericCubicEOS{T}, cond, i) where T<:AbstractPengRobinson
-    mix = eos.mixture
-    m = molecular_property(mix, i)
-    a = acentric_factor(m)
-    T_r = reduced_temperature(mix, cond, i)
-    tmp = (1.0 + (0.37464 + 1.54226*a - 0.26992*a*a)*(1-T_r^0.5))
-    return eos.ω_a*(tmp*tmp)
-end
+# Peng Robinson and variants
+include("eos/peng_robinson.jl")
+include("eos/peng_robinson_corrected.jl")
+include("eos/soreide_whitson.jl")
+# Soave Redlich-Kwong
+include("eos/soave_redlich_kwong.jl")
+# Zudkevitch-Joffe
+include("eos/zudkevitch_joffe.jl")
 
-function weight_ai(eos::GenericCubicEOS{PengRobinsonCorrected}, cond, i)
-    mix = eos.mixture
-    m = molecular_property(mix, i)
-    a = acentric_factor(m)
-    T_r = reduced_temperature(mix, cond, i)
-    aa = a*a
-    if a > 0.49
-        # Use alternate expression.
-        D = (0.379642 + 1.48503*a - 0.164423*aa + 0.016666*a*aa)
-    else
-        # Use standard expression.
-        D = (0.37464 + 1.54226*a - 0.26992*aa)
-    end
-    tmp = 1.0 + D*(1.0-T_r^0.5)
-    return eos.ω_a*(tmp*tmp);
+# Generic part
+Base.@propagate_inbounds function binary_interaction(eos::AbstractEOS, i::Int, j::Int, cond)
+    return binary_interaction(eos.mixture, i, j)
 end
 
-# ZudkevitchJoffe
-function weight_ai(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
-    zj = eos.type
-    mix = eos.mixture
-    T = cond.T
-    T_r = reduced_temperature(mix, cond, i)
-    return eos.ω_a*zj.F_a(T, i)*T_r^(-0.5)
+Base.@propagate_inbounds function binary_interaction(mixture::MultiComponentMixture{R}, i, j) where {R}
+    return binary_interaction(mixture.binary_interaction, i, j)::R
 end
 
-function weight_bi(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
-    zj = eos.type
-    T = cond.T
-    return eos.ω_b*zj.F_b(T, i)
+function binary_interaction(::Nothing, i, j)
+    return 0.0
 end
 
-# SoaveRedlichKwong
-function weight_ai(eos::GenericCubicEOS{SoaveRedlichKwong}, cond, i)
-    mix = eos.mixture
-    m = molecular_property(mix, i)
-    a = acentric_factor(m)
-    T_r = reduced_temperature(mix, cond, i)
-    return eos.ω_a*(1 + (0.48 + 1.574*a - 0.176*a^2)*(1-T_r^(1/2)))^2;
+Base.@propagate_inbounds function binary_interaction(B::AbstractMatrix, i, j)
+    return B[i, j]
 end
 
-# Generic part
-binary_interaction(eos::AbstractEOS, i, j) = binary_interaction(eos.mixture, i, j)
-binary_interaction(mixture::MultiComponentMixture{R}, i, j) where {R} = binary_interaction(mixture.binary_interaction, i, j)::R
-binary_interaction(::Nothing, i, j) = 0.0
-Base.@propagate_inbounds binary_interaction(B::AbstractMatrix, i, j) = B[i, j]
-
 """
     mixture_compressibility_factor(eos, cond, [forces, scalars, phase])
 
@@ -93,21 +70,31 @@ provided if they are already known.
 The compressibility factor adjusts the ideal gas law to account for non-linear behavior: ``pV = nRTZ``
 """
 
-function mixture_compressibility_factor(eos::AbstractCubicEOS, cond,
-                                            forces = force_coefficients(eos, cond),
-                                            scalars = force_scalars(eos, cond, forces),
-                                            phase = :unknown)
+function mixture_compressibility_factor(
+        eos::AbstractCubicEOS,
+        cond,
+        forces = force_coefficients(eos, cond),
+        scalars = force_scalars(eos, cond, forces),
+        phase = missing
+    )
+    if !ismissing(phase)
+        cond = set_phase(cond, phase)
+    end
     poly = eos_polynomial(eos, forces, scalars)
     roots = solve_roots(eos, poly)
-    r = pick_root(eos, roots, cond, forces, scalars, phase)
+    r = pick_root(eos, roots, cond, forces, scalars)
     return r
 end
 
 minimum_allowable_root(eos::AbstractCubicEOS, forces, scalars) = scalars.B
 minimum_allowable_root(eos, forces, scalars) = 1e-16
 
-@inline pick_root(eos, roots::Real, cond, forces, scalars, phase = :unknown) = roots
-function pick_root(eos, roots, cond, forces, scalars, phase = :unknown)
+@inline function pick_root(eos, roots::Real, cond, forces, scalars)
+    return roots
+end
+
+function pick_root(eos, roots, cond, forces, scalars)
+    phase = get_phase(cond)
     r_ϵ = minimum_allowable_root(eos, forces, scalars)
     max_r = maximum(roots)
     min_r = minimum((x) -> x > r_ϵ ? x : Inf, roots)
@@ -160,7 +147,23 @@ function force_coefficients(eos::AbstractCubicEOS, cond; static_size = false)
         B_i = zeros(eT, n)
     end
     coeff = (A_ij = A_ij, A_i = A_i, B_i = B_i)
-    force_coefficients!(coeff, eos, cond)
+    update_force_coefficients!(coeff, eos, cond)
+    return coeff
+end
+
+function get_force_coefficients(forces, eos, cond)
+    if forces_per_phase(eos)
+        phase = get_phase(cond)
+        if phase == :liquid
+            return forces.liquid
+        elseif phase == :vapor
+            return forces.vapor
+        else
+            error("Forces per phase are only supported for liquid and vapor phases, not $phase.")
+        end
+    else
+        return forces
+    end
 end
 
 """
@@ -170,10 +173,22 @@ In-place update of force coefficients.
 
 See also [`force_coefficients`](@ref)
 """
-function force_coefficients!(coeff, eos::AbstractCubicEOS, arg...)
-    update_attractive_linear!(coeff.A_i, eos, arg...)
-    update_attractive_quadratic!(coeff.A_ij, coeff.A_i, eos, arg...)
-    update_repulsive!(coeff.B_i, eos, arg...)
+function force_coefficients!(coeff, eos::AbstractCubicEOS, cond)
+    if forces_per_phase(eos)
+        cond = set_phase(cond, :liquid)
+        update_force_coefficients!(coeff.liquid, eos, cond)
+        cond = set_phase(cond, :vapor)
+        update_force_coefficients!(coeff.vapor, eos, cond)
+    else
+        update_force_coefficients!(coeff, eos, cond)
+    end
+    return coeff
+end
+
+function update_force_coefficients!(coeff, eos::AbstractCubicEOS, cond)
+    update_attractive_linear!(coeff.A_i, eos, cond)
+    update_attractive_quadratic!(coeff.A_ij, coeff.A_i, eos, cond)
+    update_repulsive!(coeff.B_i, eos, cond)
     return coeff
 end
 
@@ -201,7 +216,7 @@ function update_attractive_quadratic!(A_ij, A_i, eos::AbstractCubicEOS, cond)
     T = eltype(A_i)
     for i = 1:N
         @inbounds for j = i:N
-            a = sqrt(A_i[i]*A_i[j])*(one(T) - binary_interaction(eos, i, j))
+            a = sqrt(A_i[i]*A_i[j])*(one(T) - binary_interaction(eos, i, j, cond))
             a::T
             A_ij[i, j] = a
             A_ij[j, i] = a
@@ -227,7 +242,8 @@ end
 Compute EOS specific scalars for the current conditions based on the forces.
 """
 function force_scalars(eos::AbstractCubicEOS, cond, forces)
-    return cubic_scalars(forces.A_ij, forces.B_i, cond.z)
+    f = get_force_coefficients(forces, eos, cond)
+    return cubic_scalars(f.A_ij, f.B_i, cond.z)
 end
 
 function cubic_scalars(A_ij, Bv, z)
@@ -279,7 +295,8 @@ Base.@propagate_inbounds function component_fugacity_coefficient(eos::AbstractCu
     # NOTE: This returns ln(ψ), not ψ!
     m1 = eos.m_1
     m2 = eos.m_2
-    return component_fugacity_coefficient_cubic(m1, m2, cond.z, Z, scalars.A, scalars.B, forces.A_ij, forces.B_i, i)
+    f = get_force_coefficients(forces, eos, cond)
+    return component_fugacity_coefficient_cubic(m1, m2, cond.z, Z, scalars.A, scalars.B, f.A_ij, f.B_i, i)
 end
 
 Base.@propagate_inbounds function component_fugacity_coefficient_cubic(m1, m2, x, Z, A, B, A_mat, B_i, i)

src/eos/peng_robinson.jl

@@ -0,0 +1,14 @@
+
+# PengRobinson specialization
+function static_coefficients(::AbstractPengRobinson)
+    return (0.457235529, 0.077796074, 1.0 + sqrt(2), 1.0 - sqrt(2))
+end
+
+function weight_ai(eos::GenericCubicEOS{T}, cond, i) where T<:AbstractPengRobinson
+    mix = eos.mixture
+    m = molecular_property(mix, i)
+    a = acentric_factor(m)
+    T_r = reduced_temperature(mix, cond, i)
+    α_half = (1.0 + (0.37464 + 1.54226*a - 0.26992*a*a)*(1-T_r^0.5))
+    return eos.ω_a*(α_half*α_half)
+end

src/eos/peng_robinson_corrected.jl

@@ -0,0 +1,17 @@
+
+function weight_ai(eos::GenericCubicEOS{PengRobinsonCorrected}, cond, i)
+    mix = eos.mixture
+    m = molecular_property(mix, i)
+    a = acentric_factor(m)
+    T_r = reduced_temperature(mix, cond, i)
+    aa = a*a
+    if a > 0.49
+        # Use alternate expression.
+        D = (0.379642 + 1.48503*a - 0.164423*aa + 0.016666*a*aa)
+    else
+        # Use standard expression.
+        D = (0.37464 + 1.54226*a - 0.26992*aa)
+    end
+    tmp = 1.0 + D*(1.0-T_r^0.5)
+    return eos.ω_a*(tmp*tmp);
+end

src/eos/soave_redlich_kwong.jl

@@ -0,0 +1,8 @@
+# SoaveRedlichKwong
+function weight_ai(eos::GenericCubicEOS{SoaveRedlichKwong}, cond, i)
+    mix = eos.mixture
+    m = molecular_property(mix, i)
+    a = acentric_factor(m)
+    T_r = reduced_temperature(mix, cond, i)
+    return eos.ω_a*(1 + (0.48 + 1.574*a - 0.176*a^2)*(1-T_r^(1/2)))^2;
+end

src/eos/soreide_whitson.jl

@@ -0,0 +1,91 @@
+forces_per_phase(eos::GenericCubicEOS{T}) where T<:SoreideWhitson = true
+
+# Søreide-Whitson Peng-Robinson variant
+function weight_ai(eos::GenericCubicEOS{T}, cond, i) where T<:SoreideWhitson
+    sw = eos.type
+    mix = eos.mixture
+    m = molecular_property(mix, i)
+    a = acentric_factor(m)
+    T_r = reduced_temperature(mix, cond, i)
+    if sw.component_types[i] == COMPONENT_H2O
+        # Use the water-specific expression.
+        w1, w2, w3 = sw.water_coefficients
+        α_half = 1.0 + w1*(1.0 - w2*T_r) + w3*(T_r^(-3)-1.0)
+    else
+        α_half = 1.0 + (0.37464 + 1.54226*a - 0.26992*a*a)*(1-T_r^0.5)
+    end
+    return eos.ω_a*(α_half*α_half)
+end
+
+function binary_interaction(eos::GenericCubicEOS{T}, i::Int, j::Int, cond) where T<:SoreideWhitson
+    if i == j
+        bic = 0.0
+    else
+        phase = get_phase(cond)
+        sw = eos.type
+        if phase == :liquid
+            cat_i = eos.type.component_types[i]
+            cat_j = eos.type.component_types[j]
+            # Use eq 12 from paper
+            i_is_h2o = cat_i == COMPONENT_H2O
+            j_is_h2o = cat_j == COMPONENT_H2O
+            # in paper: i is hc, j is h2o
+            if i_is_h2o || j_is_h2o
+                # Special cases: H2O-H2S, H2O-N2 and H2O-CO2
+                if i_is_h2o
+                    cat_other = cat_j
+                    ix = j
+                else
+                    cat_other = cat_i
+                    ix = i
+                end
+                hc_property = molecular_property(eos.mixture, ix)
+                acf = acentric_factor(hc_property)
+                T_r = cond.T/critical_temperature(hc_property)
+                if cat_other == COMPONENT_CO2
+                    bic = soreide_whitson_co2_aqueous_bic(sw, acf, T_r)
+                elseif cat_other == COMPONENT_H2S
+                    bic = soreide_whitson_h2s_aqueous_bic(sw, acf, T_r)
+                elseif cat_other == COMPONENT_N2
+                    bic = soreide_whitson_n2_aqueous_bic(sw, acf, T_r)
+                else
+                    bic = soreide_whitson_hc_aqueous_bic(sw, acf, T_r)
+                end
+            else
+                # Use default.
+                bic = binary_interaction(eos.mixture, i, j)
+            end
+        elseif phase == :vapor
+            bic = binary_interaction(eos.mixture, i, j)
+        else
+            error("Søreide-Whitson requires the phase state to be specified for binary interactions (was: $phase).")
+        end
+    end
+    return bic
+end
+
+function soreide_whitson_n2_aqueous_bic(sw::SoreideWhitson, acf_i, T_ri)
+    c_sw = sw.molality
+    return -1.70235*(1 + 0.025587*c_sw^0.75) + 0.44338*(1 + 0.08126*c_sw^0.75)*T_ri
+end
+
+function soreide_whitson_h2s_aqueous_bic(sw::SoreideWhitson, acf_i, T_ri)
+    return -0.20441 + 0.234267*T_ri
+end
+
+function soreide_whitson_co2_aqueous_bic(sw::SoreideWhitson, acf_i, T_ri)
+    c_sw = sw.molality
+    return -0.31092*(1 + 0.15587*c_sw^0.7505) + 0.23580 * (1 + 0.17837*c_sw^0.979)*T_ri - 21.2566*exp(-6.7222*T_ri- c_sw)
+end
+
+function soreide_whitson_hc_aqueous_bic(sw::SoreideWhitson, acf_i, T_ri)
+    c_sw = sw.molality
+    a_0, a_1, a_2 = sw.A
+    a_mw_0, a_mw_1, a_mw_2 = sw.A_mw
+    α_0, α_1, α_2 = sw.alphas
+
+    A_0 = a_0 + a_mw_0*sign(acf_i)*abs(acf_i)^(-0.1)
+    A_1 = a_1 + a_mw_1*acf_i
+    A_2 = a_2 + a_mw_2*acf_i
+    return A_0*(1 + α_0*c_sw) + A_1*T_ri*(1 + α_1*c_sw) + A_2*T_ri*(1 + α_2*c_sw)
+end

src/eos/zudkevitch_joffe.jl

@@ -0,0 +1,14 @@
+# ZudkevitchJoffe
+function weight_ai(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
+    zj = eos.type
+    mix = eos.mixture
+    T = cond.T
+    T_r = reduced_temperature(mix, cond, i)
+    return eos.ω_a*zj.F_a(T, i)*T_r^(-0.5)
+end
+
+function weight_bi(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
+    zj = eos.type
+    T = cond.T
+    return eos.ω_b*zj.F_b(T, i)
+end