diff --git a/Project.toml b/Project.toml index 5a28e98..8d1e92c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MultiComponentFlash" uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff" authors = ["Olav Møyner "] -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"] diff --git a/docs/src/examples/basics.md b/docs/src/examples/basics.md index bf5bd43..3a8c4af 100644 --- a/docs/src/examples/basics.md +++ b/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,13 +85,13 @@ 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: @@ -98,8 +99,8 @@ From the chosen overall mole fractions `z`, and the flashed `K`-values together ```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, @@ -107,8 +108,8 @@ As expected, the liquid phase has more of the heavy component than in the overal ```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. diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl index 7e84721..6dc9182 100644 --- a/src/MultiComponentFlash.jl +++ b/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") diff --git a/src/eos_types.jl b/src/eos_types.jl index 2523e2f..00ad7ec 100644 --- a/src/eos_types.jl +++ b/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 \ No newline at end of file +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 diff --git a/src/kvalues_eos.jl b/src/kvalues_eos.jl new file mode 100644 index 0000000..f19d3ec --- /dev/null +++ b/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 diff --git a/src/mixture_types.jl b/src/mixture_types.jl index c81d92e..28e0ee1 100644 --- a/src/mixture_types.jl +++ b/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 diff --git a/src/rachford_rice.jl b/src/rachford_rice.jl index a55fd3a..9120e70 100644 --- a/src/rachford_rice.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index d14bf75..a6c3a9c 100644 --- a/test/runtests.jl +++ b/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