Merge pull request #19 from moyner/kvalues

K value "eos"

Project.toml

@@ -1,7 +1,7 @@
 name = "MultiComponentFlash"
 uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "1.1.8"
+version = "1.1.9"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -18,6 +18,8 @@ julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [targets]
-test = ["Test"]
+test = ["Test", "ForwardDiff", "StaticArrays"]

docs/src/examples/basics.md

@@ -54,10 +54,11 @@ z = [0.4, 0.6] # 1 mole methane per 9 moles of decane
 conditions = (p = p, T = T, z = z)
 # Perform a flash to get the vapor fraction
 V = flash_2ph(eos, conditions)
+round(V, digits = 5)
 
 # output
 
-0.2660973082539688
+0.2661
 ```
 
 ## K-values and fractions
@@ -84,31 +85,31 @@ We can also get more output by turning on the `extra_out` flag. We can use this
 
 ```jldoctest
 V, K, report = flash_2ph(eos, conditions, extra_out = true);
-K
+round.(K, digits = 5)
 
 # output
 
 2-element Vector{Float64}:
- 4.1227246229737515
- 0.0004711582061535555
+ 4.12272
+ 0.00047
 ```
 
 From the chosen overall mole fractions `z`, and the flashed `K`-values together with the vapor fraction `V` we can get the phase mole fractions in the liquid phase:
 
 ```jldoctest
 julia> liquid_mole_fraction.(z, K, V)
 2-element Vector{Float64}:
- 0.24247146623483945
- 0.7575285337651605
+ 0.24247146623483776
+ 0.7575285337651623
 ```
 
 As expected, the liquid phase has more of the heavy component than in the overall mole fractions (0.75 relative to 0.6). If we compute the vapor fractions,
 
 ```jldoctest
 julia> vapor_mole_fraction.(z, K, V)
 2-element Vector{Float64}:
- 0.9996430842149212
- 0.0003569157850789261
+ 0.9996430842149211
+ 0.000356915785078925
 ```
 
 we see that the vapor phase is almost entirely made up of the lighter methane at the chosen conditions.

src/MultiComponentFlash.jl

@@ -13,6 +13,8 @@ module MultiComponentFlash
     export SoaveRedlichKwong, RedlichKwong, PengRobinson, PengRobinsonCorrected, ZudkevitchJoffe
     # The generic cubic form that supports the above
     export GenericCubicEOS
+    # KValue "fake" EOS
+    export KValuesEOS
     export number_of_components
     # Flash interfaces
     export flash_2ph, flash_2ph!, flash_storage
@@ -43,6 +45,7 @@ module MultiComponentFlash
 
     include("mixtures.jl")
     include("kvalues.jl")
+    include("kvalues_eos.jl")
     include("rachford_rice.jl")
     include("eos.jl")
     include("flash.jl")

src/eos_types.jl

@@ -94,7 +94,21 @@ end
 
 function GenericCubicEOS(setup::NamedTuple, mixture; volume_shift = nothing)
     if !isnothing(volume_shift)
-        @assert length(volume_shift) == number_of_components(mixture) "Volume shift must have one value per component."
+        length(volume_shift) == number_of_components(mixture) || throw(ArgumentError("Volume shift must have one value per component."))
     end
     return GenericCubicEOS(setup.type, mixture, setup.m_1, setup.m_2, setup.ω_a, setup.ω_b, volume_shift)
-end
+end
+
+struct KValuesEOS{T, R, N, V} <: AbstractEOS
+    "Callable on the form `cond -> V` or a set of constants (Tuple/AbstractVector)"
+    K_values_evaluator::T
+    mixture::MultiComponentMixture{R, N}
+    volume_shift::V
+end
+
+function KValuesEOS(K, mixture; volume_shift = nothing)
+    if !isnothing(volume_shift)
+        length(volume_shift) == number_of_components(mixture) || throw(ArgumentError("Volume shift must have one value per component."))
+    end
+    return KValuesEOS(K, mixture, volume_shift)
+end

src/kvalues_eos.jl

@@ -0,0 +1,23 @@
+function evaluate_K(eos::KValuesEOS, cond)
+    return eos.K_values_evaluator(cond)
+end
+
+function evaluate_K(eos::KValuesEOS{T, <:Any, <:Any}, cond) where T<:Tuple
+    return eos.K_values_evaluator
+end
+
+function evaluate_K(eos::KValuesEOS{T, <:Any, <:Any}, cond) where T<:AbstractVector
+    return eos.K_values_evaluator
+end
+
+function initial_guess_K(eos::KValuesEOS, cond)
+    return evaluate_K(eos, cond)
+end
+
+function flash_storage(eos::KValuesEOS, cond = missing; kwarg...)
+    return nothing
+end
+
+function flash_2ph!(storage, K, eos::KValuesEOS, cond, V = NaN; kwarg...)
+    return solve_rachford_rice(K, cond.z, V)
+end

src/mixture_types.jl

@@ -38,17 +38,21 @@ end
 """
     MolecularProperty("Name")
 
-Convenience constructor that looks up molecular properties from a table in `tabulated_properties`.
+Convenience constructor that looks up molecular properties from a table in
+`tabulated_properties`.
 
-The properties are taken from the wonderful MIT-licensed [CoolProp](http://coolprop.org). Please note that
-the equations of state included in this module may not be approprioate for all the available fluids, especially for mixtures!
+The properties are taken from the wonderful MIT-licensed
+[CoolProp](http://coolprop.org). Please note that the equations of state
+included in this module may not be appropriate for all the available species
+included in CoolProp, especially for mixtures!
 
-See list of species [at CoolProp website](http://coolprop.org/fluid_properties/PurePseudoPure.html#list-of-fluids).
+See list of species [at CoolProp
+website](http://coolprop.org/fluid_properties/PurePseudoPure.html#list-of-fluids).
 
 """
 function MolecularProperty(s::String)
     t = tabulated_properties();
-    @assert haskey(t, s)
+    haskey(t, s) || throw(ArgumentError("$s not found in `tabulated_properties()`."))
     return t[s]
 end
 
@@ -74,3 +78,13 @@ struct MultiComponentMixture{R, N}
         new{realtype, n}(name, names, properties, A_ij)
     end
 end
+
+"""
+    MultiComponentMixture(names::Vector{String}; A_ij = nothing, name = "UnnamedMixture")
+
+Create a multicomponent mixture using name lookup for species.
+"""
+function MultiComponentMixture(names::Vector{String}; kwarg...)
+    props = map(MolecularProperty, names)
+    return MultiComponentMixture(props; names = names, kwarg...)
+end

src/rachford_rice.jl

@@ -20,7 +20,7 @@ julia> solve_rachford_rice([0.5, 1.5], [0.3, 0.7])
 0.8000000000000002
 ```
 """
-function solve_rachford_rice(K, z, V = NaN; tol=1e-12, maxiter=1000, ad = false, analytical = false, verbose = false)
+function solve_rachford_rice(K, z, V = NaN; tol=1e-12, maxiter=1000, ad = false, analytical = true, verbose = false)
     n = length(z)
     if analytical
         if n == 2 #exact solution, linear
@@ -66,6 +66,9 @@ function solve_rachford_rice(K, z, V = NaN; tol=1e-12, maxiter=1000, ad = false,
     if isnan(V)
         V = (V_min + V_max)/2
     end
+    if V_max < V_min
+        V_min, V_max = V_max, V_min
+    end
     verbose && println("Solving Rachford-Rice for $n components.\nInitial guess V = $V\nBounds: [$V_min, $V_max]")
     i = 1
     while i < maxiter
@@ -109,3 +112,21 @@ function objectiveRR(V, K, z)
     end
     return eq
 end
+
+function objectiveRR_dV(V, K, z)
+    RR_dv = 0.0
+    for i in eachindex(K)
+        z_i = z[i]
+        K_i = K[i]
+        RR_dv -= (z_i*(K_i - 1)^2)/(1+V*(K_i-1))^2
+    end
+    return RR_dv
+end
+
+function objectiveRR_dK(V, K, z, i)
+    return z[i]/(1.0 + V*(K[i]-1))^2
+end
+
+function objectiveRR_dz(V, K, z, i)
+    return (K[i] - 1.0)/(1.0 + V*(K[i]-1))
+end

test/runtests.jl

@@ -58,3 +58,44 @@ end
     m2 = MolecularProperty(mw = 0.0440, p_c = 7.38e6, T_c = 304.1, V_c = 9.412e-5, acentric_factor = 0.224)
     @test m1 == m2
 end
+
+@testset "K-value EOS" begin
+    mixture = MultiComponentMixture(["CarbonDioxide", "Water"])
+    eos = KValuesEOS([0.01, 2.0], mixture)
+    cond = (p = 1e5, T = 273.15, z = (0.1, 0.9))
+    @test round(flash_2ph(eos, cond), digits = 4) ≈ 0.8091
+    @test number_of_components(eos) == 2
+    eos2 = KValuesEOS(cond -> [0.01, 2.0], mixture)
+    @test round(flash_2ph(eos2, cond), digits = 4) ≈ 0.8091
+end
+
+using ForwardDiff
+@testset "Rachford-Rice derivatives" begin
+    N = 25
+    for z_light in range(0.0, 1.0, length = N)
+        z = [z_light, 1.0 - z_light]
+        for K1 in range(0.001, 1000.0, length = N)
+            for K2 in range(0.001, 1000.0, length = N)
+                K = [K1, K2]
+                for V in range(0, 1, length = N)
+                    f_v(V) = MultiComponentFlash.objectiveRR(V, K, z)
+                    dv_ad = ForwardDiff.derivative(f_v, V)
+                    dv_a = MultiComponentFlash.objectiveRR_dV(V, K, z)
+                    @test dv_a ≈ dv_ad
+                    f_K(K) = MultiComponentFlash.objectiveRR(V, K, z)
+                    dK_ad = ForwardDiff.gradient(f_K, K)
+                    for i in eachindex(K)
+                        dK_a_i = MultiComponentFlash.objectiveRR_dK(V, K, z, i)
+                        @test dK_a_i ≈ dK_ad[i]
+                    end
+                    f_z(z) = MultiComponentFlash.objectiveRR(V, K, z)
+                    dz_ad = ForwardDiff.gradient(f_z, z)
+                    for i in eachindex(z)
+                        dz_a_i = MultiComponentFlash.objectiveRR_dz(V, K, z, i)
+                        @test dz_a_i ≈ dz_ad[i]
+                    end
+                end
+            end
+        end
+    end
+end