Use one feed tangent for both phase stability trials

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

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

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()

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