From 04573590d8cab6bd61559d5c57571efd7e3b8200 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 19 Sep 2026 15:01:07 +0200 Subject: [PATCH 1/5] Add negative flash option --- src/flash.jl | 95 ++++++++++++++----- src/rachford_rice.jl | 216 ++++++++++++++++++++++++++++--------------- src/static.jl | 107 ++++++++------------- test/runtests.jl | 113 ++++++++++++++++++++++ 4 files changed, 371 insertions(+), 160 deletions(-) diff --git a/src/flash.jl b/src/flash.jl index f9c98c9..5cc4be0 100644 --- a/src/flash.jl +++ b/src/flash.jl @@ -17,9 +17,14 @@ Two outcomes are possible: # Arguments - `eos`: the equation-of-state to be used for the flash -- `c`: conditions to flash the mixture at on the form `(p = 10e5, T = 303.15, z = [0.5, 0.3, 0.2])` -- `K`: optionally a buffer of length `number_of_components(eos)` used to hold K-values. Modified in-place. -- `V`: optionally the initial guess for V. If this value is not `NaN`, the stability check will be skipped. +- `c`: conditions to flash the mixture at on the form + `c = (p = 10e5, T = 303.15, z = [0.5, 0.3, 0.2])` +- `K`: optionally a buffer of length `number_of_components(eos)` used to hold + K-values. Modified in-place. +- `V=NaN`: optionally the initial guess for V. If this value is a finite value + other than `NaN`, the stability check will be skipped. Pass `Inf` to perform a + negative flash: initialize `V` from Rachford-Rice and allow solutions outside + `[0, 1]`. # Keyword arguments - `method = SSIFlash()`: Flash method to use. Can be `SSIFlash()`, `NewtonFlash()` or `SSINewtonFlash()`. @@ -103,6 +108,7 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig; if update_forces force_coefficients!(forces, eos, c) end + negative_flash = isinf(V) single_phase_init = isnan(V) || V == 1.0 || V == 0.0 if single_phase_init stable, stability_report = stability_2ph!(storage, K, eos, c, config; @@ -125,13 +131,19 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig; i = 0 else i = 1 - if isnan(V) - V = solve_rachford_rice(K, z, V) + if isnan(V) || negative_flash + V = solve_rachford_rice(K, z, NaN) end - while true - V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i) - converged = ϵ ≤ tolerance - if converged || i == maxiter + while isfinite(V) + V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i, + negative_flash) + # A negative flash has no admissible split if the updated K-values + # cease to straddle one. Do not evaluate fugacities at a RR pole. + isfinite(V) || break + residual_converged = ϵ ≤ tolerance + if residual_converged || i == maxiter + converged = residual_converged && + (!negative_flash || valid_negative_flash_solution(V, K, z)) if print_output(config) && verbose @info "Flash done in $i iterations." V K converged end @@ -146,6 +158,9 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig; end i += 1 end + if !isfinite(V) && print_output(config) && check && !negative_flash + error("No admissible Rachford-Rice root for flash") + end end return (V, K, (its = i, converged = converged, stability = stability_report)) end @@ -280,14 +295,30 @@ function update_value(v::T, newv::Real) where T<:ForwardDiff.Dual return T(newv, P) end -function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iteration) +function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, + iteration, negative_flash::Bool = false) z = cond.z x, y = storage.x, storage.y p, T = cond.p, cond.T - return ssi!(K, p, T, x, y, z, V, eos, forces) + return ssi!(K, p, T, x, y, z, V, eos, forces, negative_flash) +end + +@inline function valid_negative_flash_solution(V, K, z) + # The trivial K ≈ 1 fixed point has an indeterminate vapor fraction. + isfinite(V) || return false + V_lo, V_hi = positive_rachford_rice_bounds(K, z) + V_lo < V < V_hi || return false + nontrivial = false + @inbounds for K_i in K + isfinite(K_i) && K_i > zero(K_i) || return false + nontrivial |= abs(K_i - one(K_i)) > 1e-6 + end + return nontrivial && + rachford_rice_balance_error(V, objectiveRR(V, K, z)) <= 1e-8 end -function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real} +function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces, + negative_flash::Bool = false) where {F<:Real} # Initialize conditions for vapor and liquid phases based on K-values x = liquid_mole_fraction!(x, z, K, V) y = vapor_mole_fraction!(y, x, K) @@ -307,14 +338,17 @@ function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real} ϵ = max(ϵ, abs(1-r)) end V = solve_rachford_rice(K, z, V) - V = clamp(V, zero(V), one(V)) + if !negative_flash + V = clamp(V, zero(V), one(V)) + end return (V, ϵ)::Tuple{F, F} end cap_z(z) = min(max(z, MINIMUM_COMPOSITION), one(z)) cap_VL(v) = min(max(v, MINIMUM_COMPOSITION), 1 - MINIMUM_COMPOSITION) -function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration) +function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, + iteration, negative_flash::Bool = false) x, y = storage.x, storage.y z = cond.z x = liquid_mole_fraction!(x, z, K, V) @@ -322,11 +356,17 @@ function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iter # Newton part Δ = update_and_solve!(storage, eos, cond, forces, x, y, V) newton_dampen!(type.dMax, Δ) - V, ϵ = update_newton_from_increment!(K, x, y, V, Δ) + V, ϵ = update_newton_from_increment!(K, x, y, V, Δ, negative_flash) + if negative_flash + # The Newton composition update is clipped to positive values. Restore + # material balance with the RR root in the positive-composition window. + V = solve_rachford_rice(K, z, V) + end return (V, ϵ) end -function update_newton_from_increment!(K, x, y, V, Δ) +function update_newton_from_increment!(K, x, y, V, Δ, + negative_flash::Bool = false) ϵ = zero(eltype(K)) n = length(x) @inbounds for i in 1:n @@ -338,7 +378,10 @@ function update_newton_from_increment!(K, x, y, V, Δ) # Assign new value, overwriting old K[i] = K_next end - V = cap_VL(V - Δ[end]) + V -= Δ[end] + if !negative_flash + V = cap_VL(V) + end return (V, ϵ) end @@ -406,11 +449,14 @@ function newton_dampen!(dMax, Δ) @. Δ *= ω end -function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, iteration) +function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, + iteration, negative_flash::Bool = false) if iteration >= type.swap_iter - flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, V, iteration) + flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, + V, iteration, negative_flash) else - flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration) + flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration, + negative_flash) end end @@ -479,7 +525,9 @@ Base.@propagate_inbounds ∂(D, i) = D.partials[i] Compute liquid mole fractions from overall mole fraction `z`, vector with one K-value per component as`K` and vapor fraction `V`.""" -@inline liquid_mole_fraction(z, K, V) = z/(1 - V + V*K) +# This form avoids subtracting two large, nearly equal terms when |V| is +# large in a negative flash and K is close to one. +@inline liquid_mole_fraction(z, K, V) = z/muladd(V, K - 1, one(V)) """ y = vapor_mole_fraction(z, K, V) @@ -494,7 +542,10 @@ Compute vapor mole fractions from liquid mole fraction `x` and K-values `K`. """ @inline vapor_mole_fraction(x, K) = x*K -liquid_mole_fraction!(x, z, K, V) = begin x .= z ./ (1 .- V .+ V .* K);x end +liquid_mole_fraction!(x, z, K, V) = begin + @. x = z/muladd(V, K - 1, one(V)) + x +end vapor_mole_fraction!(y, x, K) = begin y .= x .* K;y end function vapor_mole_fraction!(y,z, K, V) x = y diff --git a/src/rachford_rice.jl b/src/rachford_rice.jl index 9120e70..750687f 100644 --- a/src/rachford_rice.jl +++ b/src/rachford_rice.jl @@ -2,10 +2,13 @@ solve_rachford_rice(K, z, [V]; ) Compute vapor mole fraction `V` for given equilibrium constants `K` and mole fractions `z`. +`V` may lie outside `[0, 1]` for a negative flash, but the corresponding phase +compositions must remain nonnegative. Return `NaN` when there is no such root. # Arguments `K` - Equal length to `z`, containing the equilibrium constants for each component. `z` - Mole fractions. Should sum up to unity. +`V` - Optional initial guess. `NaN` or `Inf` selects the midpoint of the admissible interval for iterative solves. # Keyword arguments @@ -20,95 +23,164 @@ 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 = true, verbose = false) - n = length(z) +function solve_rachford_rice(K, z, V = NaN; tol = 1e-12, maxiter = 1000, + ad = false, analytical = true, verbose = false) + V_lo, V_hi = positive_rachford_rice_bounds(K, z) + V_lo < V_hi || return oftype(K[1], NaN) if analytical - if n == 2 #exact solution, linear - z1,z2 = z - k1,k2 = K - b1,b2 = 1/(1-k1),1/(1-k2) - a1 = (z1*b2 + z2*b1) - a0 = z1+z2 - return a1/a0 - elseif n == 3 #exact solution, quadratic - z1,z2,z3 = z - k1,k2,k3 = K - b1,b2,b3 = 1/(1-k1),1/(1-k2),1/(1-k3) - a2 =(z1 + z2 + z3) - a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2) - a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1 - Δ = a1*a1 - 4*a0*a2 - inva2 = 1/(2*a2) - Δ2 = sqrt(Δ)*inva2 - x = -a1*inva2 - β1 = x-Δ2 - β2 = x+Δ2 - βmax = max(β1,β2) - βmin = min(β1,β2) - if 0 < β1 < 1 - return β1 - elseif 0 < β2 < 1 - return β2 - elseif isfinite(β1+β2) - kmin = min(k1,k2,k3) - kmax = max(k1,k2,k3) - kmin > 1 && return βmax - kmax < 1 && return βmin - else - return zero(Δ2)/zero(Δ2) + n = length(z) + if n == 2 + root = rachford_rice_analytic_2(K, z, V_lo, V_hi) + if isfinite(root) && + rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= tol + return root + end + elseif n == 3 + root = rachford_rice_analytic_3(K, z, V_lo, V_hi) + if isfinite(root) && + rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= tol + return root + end + end + end + return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi; + tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) +end + +@inline function rachford_rice_balance_error(V, residual) + # RR = sum(y) - sum(x). For normalized z, the two normalization errors + # are -V*RR and (1 - V)*RR. A small unscaled RR residual is insufficient + # when a negative flash has a very large |V|. + return abs(residual)*max(one(V), abs(V), abs(one(V) - V)) +end + +@inline function positive_rachford_rice_bounds(K, z) + # Whitson and Michelsen, Fluid Phase Equilibria 53 (1989), 51-71: + # the negative-flash window is the intersection of + # 1 + V*(K_i - 1) > 0 for every present component. Other intervals + # between poles can contain roots, but give negative phase compositions. + K_min = K_max = one(K[1]) + has_below = has_above = false + invalid = oftype(K[1], NaN) + @inbounds for i in eachindex(z) + z_i, K_i = z[i], K[i] + if !isfinite(z_i) || z_i < zero(z_i) || + (z_i > zero(z_i) && (!isfinite(K_i) || K_i <= zero(K_i))) + return invalid, invalid + end + if z_i > zero(z_i) + if K_i < one(K_i) + K_min = min(K_min, K_i) + has_below = true + elseif K_i > one(K_i) + K_max = max(K_max, K_i) + has_above = true end end end - # Lowest bound for solution - V_min = 1/(1 - maximum(K)) - # Largest bound for solution - V_max = 1/(1 - minimum(K)) - if isnan(V) - V = (V_min + V_max)/2 + has_below && has_above || return invalid, invalid + return inv(one(K_max) - K_max), inv(one(K_min) - K_min) +end + +@inline function rachford_rice_analytic_2(K, z, V_lo, V_hi) + k1, k2 = K + if k1 == one(k1) || k2 == one(k2) + return oftype(k1, NaN) + end + z1, z2 = z + b1, b2 = inv(one(k1) - k1), inv(one(k2) - k2) + root = (z1*b2 + z2*b1)/(z1 + z2) + return V_lo < root < V_hi ? root : oftype(root, NaN) +end + +@inline function rachford_rice_analytic_3(K, z, V_lo, V_hi) + k1, k2, k3 = K + if k1 == one(k1) || k2 == one(k2) || k3 == one(k3) + return oftype(k1, NaN) + end + z1, z2, z3 = z + b1, b2, b3 = inv(one(k1) - k1), inv(one(k2) - k2), inv(one(k3) - k3) + a2 = z1 + z2 + z3 + a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2) + a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1 + discriminant = a1*a1 - 4*a0*a2 + if discriminant >= zero(discriminant) + inv_2a2 = inv(2*a2) + offset = sqrt(discriminant)*inv_2a2 + center = -a1*inv_2a2 + root1, root2 = center - offset, center + offset + if V_lo < root1 < V_hi + return root1 + elseif V_lo < root2 < V_hi + return root2 + end end - if V_max < V_min - V_min, V_max = V_max, V_min + return oftype(k1, NaN) +end + +@inline function solve_rachford_rice_bounded(K, z, V, V_lo, V_hi; + tol = 1e-12, maxiter = 1000, ad = false, verbose = false) + # A supplied guess may be on a pole or outside the window after an SSI + # K-update. Never evaluate the Rachford-Rice function at that guess. + if !(V_lo < V < V_hi) || !isfinite(V) + V = V_lo/2 + V_hi/2 end - verbose && println("Solving Rachford-Rice for $n components.\nInitial guess V = $V\nBounds: [$V_min, $V_max]") - i = 1 - while i < maxiter + V += zero(V_lo) + best_V = V + best_error = oftype(V, Inf) + verbose && println("Solving Rachford-Rice in ($V_lo, $V_hi) from $V") + for iteration in 1:maxiter + residual = zero(V) + denominator = zero(V) + @inbounds for i in eachindex(z) + z_i = z[i] + iszero(z_i) && continue + delta_K = K[i] - one(K[i]) + term = muladd(V, delta_K, one(V)) + if !(term > zero(term)) || !isfinite(term) + return oftype(V, NaN) + end + residual += z_i*delta_K/term + denominator += z_i*delta_K^2/term^2 + end + if !isfinite(residual) || !isfinite(denominator) + return oftype(V, NaN) + end + balance_error = rachford_rice_balance_error(V, residual) + if balance_error < best_error + best_V, best_error = V, balance_error + end + balance_error <= tol && return V + denominator > zero(denominator) || return oftype(V, NaN) if ad - f(V) = objectiveRR(V, K, z) - r = f(V) - Jo = ForwardDiff.derivative(f, V) - V = V - Jo\r + denominator = -ForwardDiff.derivative(v -> objectiveRR(v, K, z), V) + end + if residual > zero(residual) + V_lo = V else - denum = 0.0 - r = 0.0 - @inbounds for i = 1:n - ΔK = K[i]-1.0 - a = z[i]*ΔK - b = (1.0 + V*ΔK) - r += a/b - denum += z[i]*ΔK^2/b^2 - end - V = V + r/denum + V_hi = V end - e = abs(r) - verbose && println("#$i V = $V, Δ=$e") - if e < tol - break + V_next = V + residual/denominator + if !(V_lo < V_next < V_hi) || !isfinite(V_next) || V_next == V + V_next = V_lo/2 + V_hi/2 end - if V < V_min || V > V_max - verbose && println("V = $V outside bounds [$V_min, $V_max], using bisection.") - obj = v -> objectiveRR(v, K, z) - V = find_zero(obj, (V_min, V_max), Bisection()) - break + if !(V_lo < V_next < V_hi) || V_next == V + # A root arbitrarily close to a pole can exhaust Float64 spacing + # before meeting the requested residual tolerance. Use the best + # representable point only if both phase sums remain accurate. + return best_error <= max(tol, 1e-8) ? best_V : oftype(V, NaN) end - i += 1 + verbose && println("#$iteration V = $V_next, residual = $residual") + V = V_next end - return V + return oftype(V, NaN) end function objectiveRR(V, K, z) eq = 0.0 for (i, k) in enumerate(K) - @inbounds eq = eq + ((k - 1.0)*z[i])/(1.0 + V*(k - 1.0)) + iszero(z[i]) && continue + @inbounds eq = eq + ((k - 1.0)*z[i])/muladd(V, k - 1.0, 1.0) end return eq end diff --git a/src/static.jl b/src/static.jl index 7bd8273..cf613ea 100644 --- a/src/static.jl +++ b/src/static.jl @@ -22,6 +22,9 @@ Run the immutable, accelerator-friendly two-phase flash implementation. scalar vapor fraction. When `storage` is omitted, a static storage marker is created automatically. +Pass `V0=Inf` to skip stability testing and perform a negative flash, allowing +the vapor fraction to converge outside `[0, 1]`. + Set `return_stability=true` to additionally return a [`StaticStabilityResult`](@ref). Its immutable `storage` can be supplied as `stability_storage` on the next call to enable the Michelsen bypass. The @@ -212,38 +215,25 @@ end end @inline function solve_rachford_rice(K::StaticVector{2}, z::StaticVector{2}, V = NaN) - z1, z2 = z - k1, k2 = K - b1, b2 = inv(1 - k1), inv(1 - k2) - return (z1*b2 + z2*b1)/(z1 + z2) + V_lo, V_hi = positive_rachford_rice_bounds(K, z) + V_lo < V_hi || return oftype(K[1], NaN) + root = rachford_rice_analytic_2(K, z, V_lo, V_hi) + if isfinite(root) && + rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= 1e-12 + return root + end + return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi) end @inline function solve_rachford_rice(K::StaticVector{3}, z::StaticVector{3}, V = NaN) - z1, z2, z3 = z - k1, k2, k3 = K - b1, b2, b3 = inv(1-k1), inv(1-k2), inv(1-k3) - a2 = z1 + z2 + z3 - a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2) - a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1 - discriminant = a1*a1 - 4*a0*a2 - if discriminant >= zero(discriminant) - inv_2a2 = inv(2*a2) - root_offset = sqrt(discriminant)*inv_2a2 - root_center = -a1*inv_2a2 - root1 = root_center - root_offset - root2 = root_center + root_offset - if zero(root1) < root1 < one(root1) - return root1 - elseif zero(root2) < root2 < one(root2) - return root2 - elseif isfinite(root1 + root2) - kmin = min(k1, k2, k3) - kmax = max(k1, k2, k3) - kmin > one(kmin) && return max(root1, root2) - kmax < one(kmax) && return min(root1, root2) - end + V_lo, V_hi = positive_rachford_rice_bounds(K, z) + V_lo < V_hi || return oftype(K[1], NaN) + root = rachford_rice_analytic_3(K, z, V_lo, V_hi) + if isfinite(root) && + rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= 1e-12 + return root end - return solve_rachford_rice_static_iterative(K, z, V) + return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi) end @inline solve_rachford_rice(K::StaticVector, z::StaticVector, V = NaN) = @@ -251,36 +241,10 @@ end @inline function solve_rachford_rice_static_iterative(K, z, V; tol = 1e-12, maxiter = 1000) - V_lo = inv(1 - maximum(K)) - V_hi = inv(1 - minimum(K)) - if V_hi < V_lo - V_lo, V_hi = V_hi, V_lo - end - if isnan(V) - V = (V_lo + V_hi)/2 - end - for _ in 1:maxiter - residual = zero(V) - denominator = zero(V) - @inbounds for i in eachindex(K) - delta_K = K[i] - one(K[i]) - term_denominator = one(V) + V*delta_K - residual += z[i]*delta_K/term_denominator - denominator += z[i]*delta_K^2/term_denominator^2 - end - abs(residual) < tol && break - if residual > zero(residual) - V_lo = V - else - V_hi = V - end - V_next = V + residual/denominator - if !(V_lo < V_next < V_hi) || !isfinite(V_next) - V_next = (V_lo + V_hi)/2 - end - V = V_next - end - return V + V_lo, V_hi = positive_rachford_rice_bounds(K, z) + V_lo < V_hi || return oftype(K[1], NaN) + return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi; + tol = tol, maxiter = maxiter) end @generated function static_fugacities(eos::GenericCubicEOS{E, R, N}, cond, forces, @@ -294,7 +258,7 @@ end end @inline function static_ssi(K::SVector{N, F}, p::F, T::F, z, V::F, - eos, forces) where {N, F<:Real} + eos, forces, negative_flash::Bool = false) where {N, F<:Real} x = SVector{N, F}(ntuple(i -> liquid_mole_fraction(z[i], K[i], V), Val(N))) y = SVector{N, F}(ntuple(i -> vapor_mole_fraction(x[i], K[i]), Val(N))) liquid = (p = p, T = T, z = x, phase = Val(:liquid)) @@ -308,7 +272,9 @@ end end K_next = SVector{N, F}(ntuple(i -> K[i]*ratios[i], Val(N))) V_next = solve_rachford_rice(K_next, z, V) - V_next = clamp(V_next, zero(V_next), one(V_next)) + if !negative_flash + V_next = clamp(V_next, zero(V_next), one(V_next)) + end return V_next, K_next, residual end @@ -346,6 +312,7 @@ end z = cond.z forces = static_force_coefficients(eos, cond, F) V = convert(F, V) + negative_flash = isinf(V) single_phase_init = isnan(V) || V == one(F) || V == zero(F) if single_phase_init stability_result = static_stability_2ph(K, eos, cond, forces; @@ -368,13 +335,21 @@ end iteration = 0 else iteration = 1 - if isnan(V) - V = solve_rachford_rice(K, z, V) + if isnan(V) || negative_flash + V = solve_rachford_rice(K, z, NaN) end - while true - V, K, residual = static_ssi(K, cond.p, cond.T, z, V, eos, forces) - converged = residual <= tolerance - (converged || iteration == maxiter) && break + while isfinite(V) + V, K, residual = static_ssi(K, cond.p, cond.T, z, V, eos, forces, negative_flash) + # A RR root outside the positive-composition window is not a + # negative flash. Stop before evaluating EOS fugacities at a pole. + isfinite(V) || break + residual_converged = residual <= tolerance + max_its_reached = iteration == maxiter + if residual_converged || max_its_reached + negative_converged = !negative_flash || valid_negative_flash_solution(V, K, z) + converged = residual_converged && negative_converged + break + end iteration += 1 end end diff --git a/test/runtests.jl b/test/runtests.jl index 0c5c80e..15631cc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,119 @@ end test_flash_inplace(SSIFlash()) end +@testset "Negative flash from V=Inf" begin + host_eos = get_test_eos() + static_eos = make_eos_immutable(host_eos) + z = @SVector [0.5, 0.3, 0.2] + static_config = MultiComponentFlash.StaticConfig() + + # These states are stable according to the ordinary flash, but their + # extrapolated equilibrium has V below zero or above one, respectively. + for (p, T, below_zero) in ((1e7, 250.0, true), (1e5, 400.0, false)) + cond = (p = p, T = T, z = collect(z)) + V_stable, _, stable_report = flash_2ph(host_eos, cond; + extra_out = true) + @test isnan(V_stable) + @test stable_report.stability.stable + + static_cond = (p = p, T = T, z = z) + K0 = initial_guess_K(static_eos, static_cond, static_config) + V_static, K_static, static_report = flash_2ph!(static_config, + K0, static_eos, static_cond, Inf; extra_out = true) + @test static_report.converged + @test !static_report.stability_result.stable + @test (below_zero ? V_static < 0 : V_static > 1) + V_immutable, K_immutable = flash_2ph_immutable( + static_eos, static_cond, Inf) + @test V_immutable ≈ V_static + @test K_immutable ≈ K_static + + for method in (SSIFlash(), NewtonFlash(), SSINewtonFlash()) + host_cond = (p = p, T = T, z = collect(z)) + V, K, report = flash_2ph(host_eos, host_cond, + initial_guess_K(host_eos, host_cond), Inf; + method = method, extra_out = true) + @test report.converged + @test !report.stability.stable + @test V ≈ V_static rtol = 1e-6 + @test K ≈ K_static rtol = 1e-6 + x = liquid_mole_fraction.(z, K, V) + y = vapor_mole_fraction.(x, K) + @test all(>(0), x) + @test all(>(0), y) + @test sum(x) ≈ 1 atol = 1e-7 + @test sum(y) ≈ 1 atol = 1e-7 + end + end + + # A trivial K ≈ 1 fixed point is not a converged negative flash. + trivial_cond = (p = 5e7, T = 250.0, z = collect(z)) + _, _, trivial_report = flash_2ph(host_eos, trivial_cond, + initial_guess_K(host_eos, trivial_cond), Inf; + extra_out = true, check = false, maxiter = 100) + @test !trivial_report.converged + + K4 = @SVector [0.2, 0.4, 2.0, 4.0] + z4 = @SVector [0.25, 0.25, 0.25, 0.25] + @test solve_rachford_rice(K4, z4, Inf) ≈ + solve_rachford_rice(K4, z4, NaN) + @test solve_rachford_rice(collect(K4), collect(z4), Inf) ≈ + solve_rachford_rice(collect(K4), collect(z4), NaN) + + # A guess at a pole must be reinitialized inside the positive-composition + # window, even when the correct negative-flash root lies outside [0, 1]. + for (z_negative, below_zero) in + ((@SVector([0.7, 0.2, 0.05, 0.05]), true), + (@SVector([0.05, 0.05, 0.2, 0.7]), false)) + for (K_test, z_test) in ((K4, z_negative), + (collect(K4), collect(z_negative))) + V_expected = solve_rachford_rice(K_test, z_test) + @test below_zero ? V_expected < 0 : V_expected > 1 + V_pole = 1/(1 - maximum(K_test)) + @test solve_rachford_rice(K_test, z_test, V_pole) ≈ V_expected + end + end + + # Regression for a phase-diagram failure: all K-values were above one, + # yet the old RR solve returned a root between two negative poles. That + # root gave negative/huge phase compositions and crashed the EOS. + K_runaway = @SVector [1.0496438050320456, 1.8228323947565352, + 3.6255787603685725, 8.673724738551629, + 23.930790867112876, 63.42734611565452] + z_runaway = @SVector [0.635, 0.115, 0.05, 0.1, 0.075, 0.025] + V_pole = -0.016018621040647378 + @test isnan(solve_rachford_rice(K_runaway, z_runaway, V_pole)) + @test isnan(solve_rachford_rice(collect(K_runaway), collect(z_runaway), V_pole)) + @test isnan(solve_rachford_rice(K_runaway, z_runaway)) + + # With K almost equal to one, a valid negative flash can have |V| >> 1. + # The RR stopping test and composition formula must both preserve the + # normalization of the extrapolated liquid and vapor phases. + K_near = @SVector [0.1, 0.5, 1.0 + 1e-14, 1.0 + 2e-14] + z_near = @SVector [0.01, 0.01, 0.49, 0.49] + for (K_test, z_test) in ((K_near, z_near), + (collect(K_near), collect(z_near))) + V_near = solve_rachford_rice(K_test, z_test) + x_near = liquid_mole_fraction.(z_test, K_test, V_near) + y_near = vapor_mole_fraction.(x_near, K_test) + @test V_near < 0 + @test all(>(0), x_near) + @test sum(x_near) ≈ 1 atol = 1e-10 + @test sum(y_near) ≈ 1 atol = 1e-10 + end + + K_no_split = @SVector [1.1, 2.0, 3.0] + V_bad, _, bad_report = flash_2ph!(static_config, K_no_split, + static_eos, (p = 1e7, T = 250.0, z = z), Inf; extra_out = true) + @test isnan(V_bad) + @test !bad_report.converged + V_bad_host, _, bad_host_report = flash_2ph(host_eos, + (p = 1e7, T = 250.0, z = collect(z)), collect(K_no_split), Inf; + extra_out = true) + @test isnan(V_bad_host) + @test !bad_host_report.converged +end + @testset "Static accelerator path" begin host_eos = get_test_eos() eos = make_eos_immutable(host_eos) From 98e9e10daa545c60db7d675c924cc53f7e98de5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 19 Sep 2026 15:40:44 +0200 Subject: [PATCH 2/5] Improve message --- src/flash.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/flash.jl b/src/flash.jl index 5cc4be0..cc901dd 100644 --- a/src/flash.jl +++ b/src/flash.jl @@ -159,7 +159,7 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig; i += 1 end if !isfinite(V) && print_output(config) && check && !negative_flash - error("No admissible Rachford-Rice root for flash") + error("No admissible Rachford-Rice root for flash for K = $K and cond = $c") end end return (V, K, (its = i, converged = converged, stability = stability_report)) From d5337be93118821ef37e8aa256671dc258d90f48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 19 Sep 2026 15:56:13 +0200 Subject: [PATCH 3/5] Handle split between negative flash and regular in terms of when to limit --- src/flash.jl | 24 +++++++++++++++++------- src/rachford_rice.jl | 24 ++++++++++++++++++++++++ src/static.jl | 11 +++++------ test/runtests.jl | 29 +++++++++++++++++++++++++++++ 4 files changed, 75 insertions(+), 13 deletions(-) diff --git a/src/flash.jl b/src/flash.jl index cc901dd..f840046 100644 --- a/src/flash.jl +++ b/src/flash.jl @@ -131,15 +131,27 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig; i = 0 else i = 1 + # if isnan(V) || negative_flash + # V = negative_flash ? solve_rachford_rice(K, z, NaN) : + # physical_vapor_fraction(K, z, NaN) + # end if isnan(V) || negative_flash - V = solve_rachford_rice(K, z, NaN) + if negative_flash + V = solve_rachford_rice(K, z, NaN) + else + V = physical_vapor_fraction(K, z, NaN) + end + V = negative_flash ? solve_rachford_rice(K, z, NaN) : + physical_vapor_fraction(K, z, NaN) end while isfinite(V) V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i, negative_flash) # A negative flash has no admissible split if the updated K-values - # cease to straddle one. Do not evaluate fugacities at a RR pole. - isfinite(V) || break + # cease to straddle one. An ordinary flash instead uses V = 0 or 1. + if !isfinite(V) + break + end residual_converged = ϵ ≤ tolerance if residual_converged || i == maxiter converged = residual_converged && @@ -337,10 +349,8 @@ function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces, K[c] *= r ϵ = max(ϵ, abs(1-r)) end - V = solve_rachford_rice(K, z, V) - if !negative_flash - V = clamp(V, zero(V), one(V)) - end + V = negative_flash ? solve_rachford_rice(K, z, V) : + physical_vapor_fraction(K, z, V) return (V, ϵ)::Tuple{F, F} end diff --git a/src/rachford_rice.jl b/src/rachford_rice.jl index 750687f..269b3c2 100644 --- a/src/rachford_rice.jl +++ b/src/rachford_rice.jl @@ -47,6 +47,30 @@ function solve_rachford_rice(K, z, V = NaN; tol = 1e-12, maxiter = 1000, tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) end +@inline function physical_vapor_fraction(K, z, V = NaN) + # In an ordinary flash, a single-phase RR iterate belongs on the nearest + # boundary of [0, 1]. Only negative flashes require an extrapolated root + # with two positive phase compositions. RR decreases on [0, 1] for K > 0. + r_liquid = r_vapor = zero(K[1]*z[1]) + @inbounds for i in eachindex(z) + K_i, z_i = K[i], z[i] + if !isfinite(K_i) || K_i <= zero(K_i) || + !isfinite(z_i) || z_i < zero(z_i) + return oftype(r_liquid, NaN) + end + delta_K = K_i - one(K_i) + r_liquid += z_i*delta_K + r_vapor += z_i*delta_K/K_i + end + if r_liquid <= zero(r_liquid) + return zero(r_liquid) + elseif r_vapor >= zero(r_vapor) + return one(r_vapor) + end + # Opposite endpoint signs imply a unique physical root in (0, 1). + return solve_rachford_rice(K, z, V) +end + @inline function rachford_rice_balance_error(V, residual) # RR = sum(y) - sum(x). For normalized z, the two normalization errors # are -V*RR and (1 - V)*RR. A small unscaled RR residual is insufficient diff --git a/src/static.jl b/src/static.jl index cf613ea..ef48c78 100644 --- a/src/static.jl +++ b/src/static.jl @@ -271,10 +271,8 @@ end residual = max(residual, abs(one(F) - ratios[i])) end K_next = SVector{N, F}(ntuple(i -> K[i]*ratios[i], Val(N))) - V_next = solve_rachford_rice(K_next, z, V) - if !negative_flash - V_next = clamp(V_next, zero(V_next), one(V_next)) - end + V_next = negative_flash ? solve_rachford_rice(K_next, z, V) : + physical_vapor_fraction(K_next, z, V) return V_next, K_next, residual end @@ -336,12 +334,13 @@ end else iteration = 1 if isnan(V) || negative_flash - V = solve_rachford_rice(K, z, NaN) + V = negative_flash ? solve_rachford_rice(K, z, NaN) : + physical_vapor_fraction(K, z, NaN) end while isfinite(V) V, K, residual = static_ssi(K, cond.p, cond.T, z, V, eos, forces, negative_flash) # A RR root outside the positive-composition window is not a - # negative flash. Stop before evaluating EOS fugacities at a pole. + # negative flash. Ordinary flashes use a single-phase boundary. isfinite(V) || break residual_converged = residual <= tolerance max_its_reached = iteration == maxiter diff --git a/test/runtests.jl b/test/runtests.jl index 15631cc..bd92fe1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -155,6 +155,35 @@ end @test !bad_host_report.converged end +@testset "Single-phase RR iterates in ordinary flashes" begin + # This K iterate was observed while sweeping the simple phase diagram. + # It has no negative-flash root, but an ordinary flash must remain on the + # vapor boundary and continue rather than turn V into NaN. + K_vapor = @SVector [1.2749359332142336, 1.0540125780535836, + 1.0710895839941093] + z = @SVector [0.3, 0.1, 0.6] + for (K_test, z_test) in ((K_vapor, z), + (collect(K_vapor), collect(z))) + @test isnan(solve_rachford_rice(K_test, z_test)) + @test MultiComponentFlash.physical_vapor_fraction(K_test, z_test) == 1 + @test MultiComponentFlash.physical_vapor_fraction(inv.(K_test), z_test) == 0 + end + + eos, _ = cubic_benchmark("simple") + cond = (p = 364529.05811623245, T = 274.15, z = collect(z)) + V, _, report = flash_2ph(eos, cond, collect(K_vapor), 0.5; + extra_out = true) + @test report.converged + @test 0 <= V <= 1 + + static_cond = (p = cond.p, T = cond.T, z = z) + V_static, _, static_report = flash_2ph!(MultiComponentFlash.StaticConfig(), + K_vapor, make_eos_immutable(eos), static_cond, 0.5; + extra_out = true) + @test static_report.converged + @test V_static ≈ V +end + @testset "Static accelerator path" begin host_eos = get_test_eos() eos = make_eos_immutable(host_eos) From 1c56dfe2b162961b520cc8e18a148cee5d8273b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 19 Sep 2026 17:46:35 +0200 Subject: [PATCH 4/5] Stability bypass fixes --- src/static.jl | 29 ++++++++++++++++++++--------- test/runtests.jl | 31 +++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 9 deletions(-) diff --git a/src/static.jl b/src/static.jl index ef48c78..4223471 100644 --- a/src/static.jl +++ b/src/static.jl @@ -30,6 +30,8 @@ Set `return_stability=true` to additionally return a `stability_storage` on the next call to enable the Michelsen bypass. The default `bypass_tolerance=10` uses the same conservative pressure, temperature and composition bounds as the mutable simulator integration. +Set `stability_bypass=false` to force a fresh stability test even when storage +is supplied; `return_stability=true` still records a new reference if eligible. The immutable path currently supports `SSIFlash` and generic cubic EOS values converted with [`make_eos_immutable`](@ref). @@ -62,7 +64,8 @@ end method = method, extra_out = true, stability_storage = stability_storage, - stability_bypass = stability_bypass || return_stability, + stability_bypass = stability_bypass, + update_bypass = stability_bypass || return_stability, kwarg...) return immutable_flash_output(V, K, report.stability_result, Val(return_stability)) @@ -171,9 +174,10 @@ end return false end reference = storage.reference - return maximum(abs, reference.z - cond.z) < b/tolerance && + stable = maximum(abs, reference.z - cond.z) < b/tolerance && abs(reference.p - cond.p) < b*abs(cond.p)/tolerance && abs(reference.T - cond.T) < b*tolerance + return stable end @inline static_minimum_eigenvalue(B::SMatrix) = @@ -301,6 +305,7 @@ end z_min = MINIMUM_COMPOSITION, stability_storage = nothing, stability_bypass::Bool = !isnothing(stability_storage), + update_bypass::Bool = stability_bypass, bypass_tolerance::Real = 10.0, kwarg... ) where {E, R, N} @@ -315,7 +320,8 @@ end if single_phase_init stability_result = static_stability_2ph(K, eos, cond, forces; storage = stability_storage_value(stability_storage), - update_bypass = stability_bypass, + use_bypass = stability_bypass, + update_bypass = update_bypass, bypass_tolerance = bypass_tolerance, maxiter = maxiter, kwarg...) @@ -394,6 +400,7 @@ end maxiter = 1000 ) where {N, F} trivial = false + converged = false S = one(F) iter = 0 xy = zero(SVector{N, F}) @@ -420,12 +427,13 @@ end converged = R_norm < tol_equil if trivial || converged break - elseif iter == maxiter - trivial = true + elseif iter >= maxiter + # An unfinished trial phase is not evidence of a trivial minimum. + # In particular, it must not create a single-phase bypass cache. break end end - stable = trivial || S <= one(F) + tol_sat + stable = trivial || (converged && S <= one(F) + tol_sat) return stable, trivial, iter, K, xy end @@ -434,10 +442,11 @@ end check_liquid::Bool = true, storage = nothing, update_bypass::Bool = false, + use_bypass::Bool = update_bypass, bypass_tolerance::Real = 10.0, kwarg... ) where {N, F} - if update_bypass && !isnothing(storage) && + if use_bypass && !isnothing(storage) && stability_bypass_available(storage, cond; tolerance = bypass_tolerance) report = StabilityReport(true, true, true, true) @@ -464,8 +473,10 @@ end report = StabilityReport(stable_liquid, trivial_liquid, stable_vapor, trivial_vapor) K_out = report.stable ? K_liquid : static_divide(y, x) - if update_bypass && report.stable && report.liquid.trivial && - report.vapor.trivial + # Only a complete, converged two-sided test can establish a new reference. + # Skipped trial phases are reported as stable above, but prove no such thing. + if update_bypass && check_liquid && check_vapor && report.stable && + report.liquid.trivial && report.vapor.trivial critical_distance = static_michelsen_critical_point_measure( eos, cond.p, cond.T, cond.z) next_storage = StaticStabilityStorage(cond, critical_distance) diff --git a/test/runtests.jl b/test/runtests.jl index bd92fe1..82c7338 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -271,6 +271,37 @@ end @test isnan(V_nearby) @test all(isfinite, K_nearby) @test nearby_stability.bypassed + _, _, forced_test = flash_2ph_immutable(eos, nearby; + stability_storage = flash_stability, + stability_bypass = false, + return_stability = true) + @test !forced_test.bypassed + @test forced_test.storage.reference == nearby + + # A trial that hits its iteration limit has not found a trivial + # stability minimum, even if the old implementation assumed stability. + unfinished = stability_2ph_immutable(eos, c; maxiter = 1) + @test !unfinished.stable + @test isnan(unfinished.storage.critical_distance) + near_unstable = (p = 1.01e6, T = c.T, z = c.z) + @test !stability_2ph_immutable(eos, near_unstable, unfinished).bypassed + @test !stability_2ph_immutable(eos, near_unstable).stable + + # A one-sided test may be useful, but cannot certify a reference for + # bypassing both trial phases at a later condition. + vapor_unstable = (p = 1e6, T = 200.0, z = c.z) + one_sided = stability_2ph_immutable(eos, vapor_unstable; + check_vapor = false) + @test one_sided.stable + @test isnan(one_sided.storage.critical_distance) + skipped_liquid = stability_2ph_immutable(eos, stable_cond; + check_liquid = false) + @test skipped_liquid.stable + @test isnan(skipped_liquid.storage.critical_distance) + near_vapor_unstable = (p = 1.01e6, T = 200.0, z = c.z) + @test !stability_2ph_immutable(eos, near_vapor_unstable, + one_sided).bypassed + @test !stability_2ph_immutable(eos, near_vapor_unstable).stable @test_throws ArgumentError stability_2ph_immutable(eos, nearby, stability.storage; bypass_tolerance = 0.0) end From c41f86e862fb3589b06c12657fda92da925694fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 19 Sep 2026 20:09:50 +0200 Subject: [PATCH 5/5] Clean up the API --- src/MultiComponentFlash.jl | 2 +- src/flash.jl | 19 +++++----- src/kvalues_eos.jl | 12 +++++-- src/rachford_rice.jl | 72 +++++++++++++++++++++++--------------- src/static.jl | 63 ++++++++++++++++++++++----------- test/runtests.jl | 55 +++++++++++++++++++++-------- 6 files changed, 147 insertions(+), 76 deletions(-) diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl index eb6a398..8794ec2 100644 --- a/src/MultiComponentFlash.jl +++ b/src/MultiComponentFlash.jl @@ -29,7 +29,7 @@ module MultiComponentFlash # K-values export wilson_estimate, wilson_estimate!, initial_guess_K, initial_guess_K! # Vapor-liquid equilibrium - export solve_rachford_rice + export solve_rachford_rice, solve_rachford_rice_unconstrained export liquid_mole_fraction, vapor_mole_fraction export component_fugacity, mixture_fugacities, mixture_compressibility_factor diff --git a/src/flash.jl b/src/flash.jl index f840046..0f85dc2 100644 --- a/src/flash.jl +++ b/src/flash.jl @@ -131,18 +131,12 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig; i = 0 else i = 1 - # if isnan(V) || negative_flash - # V = negative_flash ? solve_rachford_rice(K, z, NaN) : - # physical_vapor_fraction(K, z, NaN) - # end if isnan(V) || negative_flash if negative_flash - V = solve_rachford_rice(K, z, NaN) + V = solve_rachford_rice_unconstrained(K, z, NaN) else - V = physical_vapor_fraction(K, z, NaN) + V = solve_rachford_rice(K, z, NaN) end - V = negative_flash ? solve_rachford_rice(K, z, NaN) : - physical_vapor_fraction(K, z, NaN) end while isfinite(V) V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i, @@ -349,8 +343,11 @@ function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces, K[c] *= r ϵ = max(ϵ, abs(1-r)) end - V = negative_flash ? solve_rachford_rice(K, z, V) : - physical_vapor_fraction(K, z, V) + if negative_flash + V = solve_rachford_rice_unconstrained(K, z, V) + else + V = solve_rachford_rice(K, z, V) + end return (V, ϵ)::Tuple{F, F} end @@ -370,7 +367,7 @@ function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, if negative_flash # The Newton composition update is clipped to positive values. Restore # material balance with the RR root in the positive-composition window. - V = solve_rachford_rice(K, z, V) + V = solve_rachford_rice_unconstrained(K, z, V) end return (V, ϵ) end diff --git a/src/kvalues_eos.jl b/src/kvalues_eos.jl index 62307f3..d39ce7d 100644 --- a/src/kvalues_eos.jl +++ b/src/kvalues_eos.jl @@ -21,11 +21,19 @@ end flash_storage(eos::KValuesEOS, cond, method, config::FlashConfig) = nothing function flash_2ph!(storage, K, eos::KValuesEOS, cond, V = NaN; kwarg...) - return solve_rachford_rice(K, cond.z, V) + return kvalue_vapor_fraction(K, cond.z, V) end function flash_2ph!(storage, K, eos::KValuesEOS, cond, V, config::FlashConfig; kwarg...) - return solve_rachford_rice(K, cond.z, V) + return kvalue_vapor_fraction(K, cond.z, V) +end + +@inline function kvalue_vapor_fraction(K, z, V) + if isinf(V) + return solve_rachford_rice_unconstrained(K, z, V) + else + return solve_rachford_rice(K, z, V) + end end diff --git a/src/rachford_rice.jl b/src/rachford_rice.jl index 269b3c2..0a0cb6b 100644 --- a/src/rachford_rice.jl +++ b/src/rachford_rice.jl @@ -1,8 +1,48 @@ """ solve_rachford_rice(K, z, [V]; ) -Compute vapor mole fraction `V` for given equilibrium constants `K` and mole fractions `z`. -`V` may lie outside `[0, 1]` for a negative flash, but the corresponding phase +Compute the physical vapor mole fraction for equilibrium constants `K` and +overall mole fractions `z`. Return `0` for liquid-only conditions and `1` for +vapor-only conditions. Use [`solve_rachford_rice_unconstrained`](@ref) when a +negative flash requires a root outside `[0, 1]`. + +`V` is an optional initial guess for an interior two-phase root. Solver keyword +arguments are forwarded to `solve_rachford_rice_unconstrained` in that case. +""" +@inline function solve_rachford_rice(K, z, V = NaN; + tol = 1e-12, maxiter = 1000, ad = false, analytical = true, + verbose = false) + # RR decreases on [0, 1] for positive K. Its endpoint signs distinguish + # a physical split from a single-phase condition. + r_liquid = r_vapor = zero(K[1]*z[1]) + @inbounds for i in eachindex(z) + K_i, z_i = K[i], z[i] + if !isfinite(K_i) || K_i <= zero(K_i) || + !isfinite(z_i) || z_i < zero(z_i) + return oftype(r_liquid, NaN) + end + delta_K = K_i - one(K_i) + r_liquid += z_i*delta_K + r_vapor += z_i*delta_K/K_i + end + if r_liquid <= zero(r_liquid) + return zero(r_liquid) + elseif r_vapor >= zero(r_vapor) + return one(r_vapor) + end + return solve_rachford_rice_unconstrained(K, z, V; + tol = tol, maxiter = maxiter, ad = ad, analytical = analytical, + verbose = verbose) +end + +# Retain the previous internal spelling for callers that used it directly. +@inline physical_vapor_fraction(K, z, V = NaN) = solve_rachford_rice(K, z, V) + +""" + solve_rachford_rice_unconstrained(K, z, [V]; ) + +Compute a negative-flash vapor fraction for given equilibrium constants `K` +and mole fractions `z`. The root may lie outside `[0, 1]`, but both phase compositions must remain nonnegative. Return `NaN` when there is no such root. # Arguments @@ -19,11 +59,11 @@ compositions must remain nonnegative. Return `NaN` when there is no such root. # Examples ```julia-repl -julia> solve_rachford_rice([0.5, 1.5], [0.3, 0.7]) +julia> solve_rachford_rice_unconstrained([0.5, 1.5], [0.3, 0.7]) 0.8000000000000002 ``` """ -function solve_rachford_rice(K, z, V = NaN; tol = 1e-12, maxiter = 1000, +function solve_rachford_rice_unconstrained(K, z, V = NaN; tol = 1e-12, maxiter = 1000, ad = false, analytical = true, verbose = false) V_lo, V_hi = positive_rachford_rice_bounds(K, z) V_lo < V_hi || return oftype(K[1], NaN) @@ -47,30 +87,6 @@ function solve_rachford_rice(K, z, V = NaN; tol = 1e-12, maxiter = 1000, tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) end -@inline function physical_vapor_fraction(K, z, V = NaN) - # In an ordinary flash, a single-phase RR iterate belongs on the nearest - # boundary of [0, 1]. Only negative flashes require an extrapolated root - # with two positive phase compositions. RR decreases on [0, 1] for K > 0. - r_liquid = r_vapor = zero(K[1]*z[1]) - @inbounds for i in eachindex(z) - K_i, z_i = K[i], z[i] - if !isfinite(K_i) || K_i <= zero(K_i) || - !isfinite(z_i) || z_i < zero(z_i) - return oftype(r_liquid, NaN) - end - delta_K = K_i - one(K_i) - r_liquid += z_i*delta_K - r_vapor += z_i*delta_K/K_i - end - if r_liquid <= zero(r_liquid) - return zero(r_liquid) - elseif r_vapor >= zero(r_vapor) - return one(r_vapor) - end - # Opposite endpoint signs imply a unique physical root in (0, 1). - return solve_rachford_rice(K, z, V) -end - @inline function rachford_rice_balance_error(V, residual) # RR = sum(y) - sum(x). For normalized z, the two normalization errors # are -V*RR and (1 - V)*RR. A small unscaled RR residual is insufficient diff --git a/src/static.jl b/src/static.jl index 4223471..2811df1 100644 --- a/src/static.jl +++ b/src/static.jl @@ -218,37 +218,54 @@ end return static_minimum_eigenvalue(B) end -@inline function solve_rachford_rice(K::StaticVector{2}, z::StaticVector{2}, V = NaN) +@inline function solve_rachford_rice_unconstrained( + K::StaticVector{2}, z::StaticVector{2}, V = NaN; + tol = 1e-12, maxiter = 1000, ad = false, analytical = true, + verbose = false) V_lo, V_hi = positive_rachford_rice_bounds(K, z) V_lo < V_hi || return oftype(K[1], NaN) - root = rachford_rice_analytic_2(K, z, V_lo, V_hi) - if isfinite(root) && - rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= 1e-12 - return root + if analytical + root = rachford_rice_analytic_2(K, z, V_lo, V_hi) + if isfinite(root) && + rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= tol + return root + end end - return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi) + return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi; + tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) end -@inline function solve_rachford_rice(K::StaticVector{3}, z::StaticVector{3}, V = NaN) +@inline function solve_rachford_rice_unconstrained( + K::StaticVector{3}, z::StaticVector{3}, V = NaN; + tol = 1e-12, maxiter = 1000, ad = false, analytical = true, + verbose = false) V_lo, V_hi = positive_rachford_rice_bounds(K, z) V_lo < V_hi || return oftype(K[1], NaN) - root = rachford_rice_analytic_3(K, z, V_lo, V_hi) - if isfinite(root) && - rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= 1e-12 - return root + if analytical + root = rachford_rice_analytic_3(K, z, V_lo, V_hi) + if isfinite(root) && + rachford_rice_balance_error(root, objectiveRR(root, K, z)) <= tol + return root + end end - return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi) + return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi; + tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) end -@inline solve_rachford_rice(K::StaticVector, z::StaticVector, V = NaN) = - solve_rachford_rice_static_iterative(K, z, V) +@inline function solve_rachford_rice_unconstrained( + K::StaticVector, z::StaticVector, V = NaN; + tol = 1e-12, maxiter = 1000, ad = false, analytical = true, + verbose = false) + return solve_rachford_rice_static_iterative(K, z, V; + tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) +end @inline function solve_rachford_rice_static_iterative(K, z, V; - tol = 1e-12, maxiter = 1000) + tol = 1e-12, maxiter = 1000, ad = false, verbose = false) V_lo, V_hi = positive_rachford_rice_bounds(K, z) V_lo < V_hi || return oftype(K[1], NaN) return solve_rachford_rice_bounded(K, z, V, V_lo, V_hi; - tol = tol, maxiter = maxiter) + tol = tol, maxiter = maxiter, ad = ad, verbose = verbose) end @generated function static_fugacities(eos::GenericCubicEOS{E, R, N}, cond, forces, @@ -275,8 +292,11 @@ end residual = max(residual, abs(one(F) - ratios[i])) end K_next = SVector{N, F}(ntuple(i -> K[i]*ratios[i], Val(N))) - V_next = negative_flash ? solve_rachford_rice(K_next, z, V) : - physical_vapor_fraction(K_next, z, V) + if negative_flash + V_next = solve_rachford_rice_unconstrained(K_next, z, V) + else + V_next = solve_rachford_rice(K_next, z, V) + end return V_next, K_next, residual end @@ -340,8 +360,11 @@ end else iteration = 1 if isnan(V) || negative_flash - V = negative_flash ? solve_rachford_rice(K, z, NaN) : - physical_vapor_fraction(K, z, NaN) + if negative_flash + V = solve_rachford_rice_unconstrained(K, z, NaN) + else + V = solve_rachford_rice(K, z, NaN) + end end while isfinite(V) V, K, residual = static_ssi(K, cond.p, cond.T, z, V, eos, forces, negative_flash) diff --git a/test/runtests.jl b/test/runtests.jl index 82c7338..f280c99 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -96,10 +96,10 @@ end K4 = @SVector [0.2, 0.4, 2.0, 4.0] z4 = @SVector [0.25, 0.25, 0.25, 0.25] - @test solve_rachford_rice(K4, z4, Inf) ≈ - solve_rachford_rice(K4, z4, NaN) - @test solve_rachford_rice(collect(K4), collect(z4), Inf) ≈ - solve_rachford_rice(collect(K4), collect(z4), NaN) + @test solve_rachford_rice_unconstrained(K4, z4, Inf) ≈ + solve_rachford_rice_unconstrained(K4, z4, NaN) + @test solve_rachford_rice_unconstrained(collect(K4), collect(z4), Inf) ≈ + solve_rachford_rice_unconstrained(collect(K4), collect(z4), NaN) # A guess at a pole must be reinitialized inside the positive-composition # window, even when the correct negative-flash root lies outside [0, 1]. @@ -108,10 +108,17 @@ end (@SVector([0.05, 0.05, 0.2, 0.7]), false)) for (K_test, z_test) in ((K4, z_negative), (collect(K4), collect(z_negative))) - V_expected = solve_rachford_rice(K_test, z_test) - @test below_zero ? V_expected < 0 : V_expected > 1 + V_expected = solve_rachford_rice_unconstrained(K_test, z_test) + if below_zero + @test V_expected < 0 + @test solve_rachford_rice(K_test, z_test) == 0 + else + @test V_expected > 1 + @test solve_rachford_rice(K_test, z_test) == 1 + end V_pole = 1/(1 - maximum(K_test)) - @test solve_rachford_rice(K_test, z_test, V_pole) ≈ V_expected + @test solve_rachford_rice_unconstrained(K_test, z_test, + V_pole) ≈ V_expected end end @@ -123,9 +130,12 @@ end 23.930790867112876, 63.42734611565452] z_runaway = @SVector [0.635, 0.115, 0.05, 0.1, 0.075, 0.025] V_pole = -0.016018621040647378 - @test isnan(solve_rachford_rice(K_runaway, z_runaway, V_pole)) - @test isnan(solve_rachford_rice(collect(K_runaway), collect(z_runaway), V_pole)) - @test isnan(solve_rachford_rice(K_runaway, z_runaway)) + @test isnan(solve_rachford_rice_unconstrained(K_runaway, + z_runaway, V_pole)) + @test isnan(solve_rachford_rice_unconstrained( + collect(K_runaway), collect(z_runaway), V_pole)) + @test isnan(solve_rachford_rice_unconstrained(K_runaway, z_runaway)) + @test solve_rachford_rice(K_runaway, z_runaway) == 1 # With K almost equal to one, a valid negative flash can have |V| >> 1. # The RR stopping test and composition formula must both preserve the @@ -134,7 +144,7 @@ end z_near = @SVector [0.01, 0.01, 0.49, 0.49] for (K_test, z_test) in ((K_near, z_near), (collect(K_near), collect(z_near))) - V_near = solve_rachford_rice(K_test, z_test) + V_near = solve_rachford_rice_unconstrained(K_test, z_test) x_near = liquid_mole_fraction.(z_test, K_test, V_near) y_near = vapor_mole_fraction.(x_near, K_test) @test V_near < 0 @@ -164,9 +174,9 @@ end z = @SVector [0.3, 0.1, 0.6] for (K_test, z_test) in ((K_vapor, z), (collect(K_vapor), collect(z))) - @test isnan(solve_rachford_rice(K_test, z_test)) - @test MultiComponentFlash.physical_vapor_fraction(K_test, z_test) == 1 - @test MultiComponentFlash.physical_vapor_fraction(inv.(K_test), z_test) == 0 + @test isnan(solve_rachford_rice_unconstrained(K_test, z_test)) + @test solve_rachford_rice(K_test, z_test) == 1 + @test solve_rachford_rice(inv.(K_test), z_test) == 0 end eos, _ = cubic_benchmark("simple") @@ -379,6 +389,23 @@ end @test isbitstype(typeof(static_eos)) @test static_eos.K_values_evaluator isa SVector{2, Float64} @test round(flash_2ph(static_eos, static_cond), digits = 4) ≈ 0.8091 + + # This SPE11C state has no resolvable wide-bound root but is liquid-only. + K_near_pure = SVector(0.005779466034671553, 40.46848280887238) + near_pure_eos = KValuesEOS(K_near_pure, mixture) + near_pure_cond = (p = 2.20173869e7, T = 313.7125, + z = SVector(1.0, 1.0e-10)) + @test isnan(solve_rachford_rice_unconstrained( + K_near_pure, near_pure_cond.z)) + @test solve_rachford_rice(K_near_pure, near_pure_cond.z) == 0 + @test flash_2ph(near_pure_eos, near_pure_cond) == 0 + @test isnan(flash_2ph(near_pure_eos, near_pure_cond, + K_near_pure, Inf)) + + @test flash_2ph(KValuesEOS(SVector(1.5, 2.0), mixture), + static_cond) == 1 + @test flash_2ph(KValuesEOS(SVector(0.1, 0.5), mixture), + static_cond) == 0 end @testset "Static flashed mixture storage" begin