From 6ae998f61bd2fe0a4507d8ef8c711f63ea94205f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 23 Oct 2026 20:52:50 +0200 Subject: [PATCH] Use a common feed tangent for both phase stability trials --- src/stability.jl | 33 +++++++++++++------ src/static.jl | 18 ++++++---- test/runtests.jl | 1 + test/stability_reference.jl | 65 +++++++++++++++++++++++++++++++++++++ 4 files changed, 100 insertions(+), 17 deletions(-) create mode 100644 test/stability_reference.jl diff --git a/src/stability.jl b/src/stability.jl index 3433a22..824dc9b 100644 --- a/src/stability.jl +++ b/src/stability.jl @@ -23,6 +23,20 @@ function stability_2ph!(storage, K, eos, c; kwarg...) return stability_2ph!(storage, K, eos, c, FlashConfig(); kwarg...) end +@inline function stability_reference_is_vapor(eos, cond, forces) + liquid = set_phase(cond, Val(:liquid)) + vapor = set_phase(cond, Val(:vapor)) + Z_l, scalars_l = prep(eos, liquid, forces) + Z_v, scalars_v = prep(eos, vapor, forces) + difference = zero(cond.p) + @inbounds for i in eachindex(cond.z) + lnphi_l = component_fugacity_coefficient(eos, liquid, i, Z_l, forces, scalars_l) + lnphi_v = component_fugacity_coefficient(eos, vapor, i, Z_v, forces, scalars_v) + difference += cond.z[i]*(lnphi_l - lnphi_v) + end + return difference > zero(difference) +end + function stability_2ph!(storage, K, eos, c, config::FlashConfig; verbose::Bool = false, extra_out::Bool = false, @@ -39,7 +53,12 @@ function stability_2ph!(storage, K, eos, c, config::FlashConfig; vapor = (p = p, T = T, z = y, phase = :vapor) current_as_liquid = (p = p, T = T, z = z, phase = :liquid) current_as_vapor = (p = p, T = T, z = z, phase = :vapor) - mixture_fugacities!(f_z, eos, current_as_vapor, forces) + # Both trial phases must be compared with ONE tangent plane. Choose the + # homogeneous feed with the lower Gibbs energy, including phase-dependent + # EOS coefficients, and retain its chemical potentials for both trials. + reference_is_vapor = stability_reference_is_vapor(eos, c, forces) + reference = reference_is_vapor ? current_as_vapor : current_as_liquid + mixture_fugacities!(f_z, eos, reference, forces) if check_vapor wilson_estimate!(K, eos, p, T) v = michelsen_test!(vapor, f_z, f_xy, vapor.z, z, K, eos, c, forces, Val(true), config; kwarg...) @@ -48,12 +67,6 @@ function stability_2ph!(storage, K, eos, c, config::FlashConfig; end stable_vapor, trivial_vapor, i_v = v if check_liquid - if forces_per_phase(eos) - # Need to recalculate fugacities for the liquid phase if the flash - # uses e.g. different bic coefficients for each phase. Otherwise, - # these are already ok. - mixture_fugacities!(f_z, eos, current_as_liquid, forces) - end wilson_estimate!(K, eos, p, T) l = michelsen_test!(liquid, f_z, f_xy, liquid.z, z, K, eos, c, forces, Val(false), config; kwarg...) else @@ -104,6 +117,7 @@ function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces, maxiter = 1000 ) trivial = false + converged = false S = 1.0 iter = 0 done = false @@ -138,12 +152,11 @@ function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces, ok = trivial || converged done = ok || iter == maxiter if done && !ok - trivial = true if print_output(config) - @warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K + @warn "Stability test failed to converge in $maxiter iterations." cond xy K_norm R_norm K end end end - stable = trivial || S <= 1.0 + tol_sat + stable = trivial || (converged && S <= 1.0 + tol_sat) return (stable, trivial, iter) end diff --git a/src/static.jl b/src/static.jl index dbc3181..777d348 100644 --- a/src/static.jl +++ b/src/static.jl @@ -526,21 +526,25 @@ end report = StabilityReport(true, true, true, true) return StaticStabilityResult(true, report, K, storage, true) end - vapor_phase = (p = cond.p, T = cond.T, z = cond.z, phase = Val(:vapor)) - f_z_vapor = static_fugacities(eos, vapor_phase, forces, F) + # A two-sided test uses two trial phases against the SAME feed tangent. + # Switching the reference with the trial can hide a cross-phase instability. + if stability_reference_is_vapor(eos, cond, forces) + reference = set_phase(cond, Val(:vapor)) + f_z = static_fugacities(eos, reference, forces, F) + else + reference = set_phase(cond, Val(:liquid)) + f_z = static_fugacities(eos, reference, forces, F) + end K_wilson = initial_guess_K(eos, cond, StaticConfig()) if check_vapor stable_vapor, trivial_vapor, i_v, K_vapor, y = static_michelsen_test( - f_z_vapor, cond.z, K_wilson, eos, cond, forces, Val(true); kwarg...) + f_z, cond.z, K_wilson, eos, cond, forces, Val(true); kwarg...) else stable_vapor, trivial_vapor, i_v, K_vapor, y = true, true, 0, K_wilson, cond.z end if check_liquid - liquid_phase = (p = cond.p, T = cond.T, z = cond.z, phase = Val(:liquid)) - f_z_liquid = forces_per_phase(eos) ? - static_fugacities(eos, liquid_phase, forces, F) : f_z_vapor stable_liquid, trivial_liquid, i_l, K_liquid, x = static_michelsen_test( - f_z_liquid, cond.z, K_wilson, eos, cond, forces, Val(false); kwarg...) + f_z, cond.z, K_wilson, eos, cond, forces, Val(false); kwarg...) else stable_liquid, trivial_liquid, i_l, K_liquid, x = true, true, 0, K_wilson, cond.z end diff --git a/test/runtests.jl b/test/runtests.jl index 7ab2dcd..dc256d0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ include("test_setup.jl") include("eos_equations.jl") include("implicit_derivatives.jl") +include("stability_reference.jl") @testset "Rachford-Rice" begin test_rachford_rice() diff --git a/test/stability_reference.jl b/test/stability_reference.jl new file mode 100644 index 0000000..9833325 --- /dev/null +++ b/test/stability_reference.jl @@ -0,0 +1,65 @@ +using Test, MultiComponentFlash, StaticArrays, LinearAlgebra + +@testset "One feed reference for both stability trials" begin + # A phase-dependent model isolates the reference error from any + # hydrogen correlation or reservoir input. Both self-trials used to report + # trivial stability, despite a negative cross-model tangent-plane distance. + mixture = MultiComponentMixture( + ["Water", "Methane", "Ethane"]; + A_ij = [0.0 0.8 -0.8; 0.8 0.0 0.0; -0.8 0.0 0.0] + ) + sw = SoreideWhitson( + mixture; A = (0.0, 0.0, 0.0), + A_mw = (0.0, 0.0, 0.0), alphas = (0.0, 0.0, 0.0) + ) + eos = GenericCubicEOS(mixture, sw) + immutable_eos = make_eos_immutable(eos) + cond = (p = 7.5e6, T = 333.15, z = @SVector [1 - 2.0e-6, 1.0e-6, 1.0e-6]) + static_result = stability_2ph_immutable(immutable_eos, cond; z_min = nothing) + mutable_cond = (p = cond.p, T = cond.T, z = collect(cond.z)) + mutable_stable, mutable_report = stability_2ph( + eos, mutable_cond; + extra_out = true + ) + @test !static_result.stable + @test !mutable_stable + @test mutable_report.liquid.stable == static_result.report.liquid.stable + @test mutable_report.vapor.stable == static_result.report.vapor.stable + @test isnan(static_result.storage.critical_distance) + + # Independently evaluate the resulting trial against one feed tangent. + MCF = MultiComponentFlash + forces = MCF.static_force_coefficients(immutable_eos, cond, Float64) + vapor_feed = (; cond..., phase = Val(:vapor)) + liquid_feed = (; cond..., phase = Val(:liquid)) + f_v = MCF.static_fugacities(immutable_eos, vapor_feed, forces, Float64) + f_l = MCF.static_fugacities(immutable_eos, liquid_feed, forces, Float64) + reference = dot(cond.z, log.(f_l ./ f_v)) > 0 ? f_v : f_l + K0 = initial_guess_K(immutable_eos, cond, MCF.StaticConfig()) + distances = map((Val(true), Val(false))) do inside_is_vapor + _, _, _, _, trial = MCF.static_michelsen_test( + reference, cond.z, K0, + immutable_eos, cond, forces, inside_is_vapor + ) + trial_cond = ( + p = cond.p, T = cond.T, z = trial, + phase = MCF.stability_phase(inside_is_vapor), + ) + trial_f = MCF.static_fugacities(immutable_eos, trial_cond, forces, Float64) + dot(trial, log.(trial_f ./ reference)) + end + @test minimum(distances) < -1.0e-10 +end + +@testset "An unfinished mutable stability trial is not stable" begin + eos = get_test_eos() + cond = test_conditions() + stable, report = stability_2ph( + eos, cond, initial_guess_K(eos, cond), + MultiComponentFlash.FlashConfig(print_output = false); + maxiter = 1, extra_out = true + ) + @test !stable + @test !report.liquid.trivial + @test !report.vapor.trivial +end