fix #328

ClapeyronThermo · Jan 8, 2025 · 08539a5 · 08539a5
1 parent b8347e2
commit 08539a5
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 8 deletions.
diff --git a/src/methods/initial_guess.jl b/src/methods/initial_guess.jl
@@ -596,18 +596,18 @@ function x0_sat_pure_spinodal(model,T,v_lb,v_ub,B = second_virial_coefficient(mo
     pub = p(v_ub)
     psl = p(vsl)
     psv = p(vsv)
-    if plb < psv
+    pmid = 0.5*max(zero(psl),psl) + 0.5*psv
+
+    if plb <= pmid
         vsv_ub = volume(model,psv,T,phase = :l, vol0 = v_lb)
     else
         vsl_lb = one(psl)*v_lb
     end
-
-    if pub > psl && psl > 0
-        vsv_ub = volume_virial(B,psl,T)
+    if pub >= pmid
+        vsv_ub = Rgas(model)*T/pmid
     else
         vsv_ub = one(psv)*v_ub
     end
-
     return _x0_sat_pure_spinodal(model,T,vsl_lb,vsv_ub,vsl,vsv,B)
 end
 
@@ -618,8 +618,6 @@ function _x0_sat_pure_spinodal(model,T,vsl_lb,vsv_ub,vsl,vsv,B)
     psl_lb,dpsl_lb,d2psl_lb = Solvers.f∂f∂2f(p,vsl_lb)
     dpsl = zero(psl)
     poly_l = Solvers.hermite5_poly(vsl_lb,vsl,psl_lb,psl,dpsl_lb,dpsl,d2psl_lb,d2psl)
-
-
     ps_mid = 0.5*(psv + max(psl,zero(psl)))
     vl = volume_from_spinodal(ps_mid,poly_l,vsl_lb,0.5*(vsl_lb + vsl) - vsl_lb)
     if psl < 0
@@ -628,13 +626,34 @@ function _x0_sat_pure_spinodal(model,T,vsl_lb,vsv_ub,vsl,vsv,B)
     end
     psv_ub,dpsv_ub,d2psv_ub = Solvers.f∂f∂2f(p,vsv_ub)
     dpsv = zero(psl)
-    poly_v = Solvers.hermite5_poly(vsv,vsv_ub,psv,psv_ub,dpsv,dpsv_ub,d2psv,d2psv_ub)    
+    poly_v = Solvers.hermite5_poly(vsv,vsv_ub,psv,psv_ub,dpsv,dpsv_ub,d2psv,d2psv_ub)  
     vv = volume_from_spinodal(ps_mid,poly_v,vsv,(zero(vsv),vsv_ub - vsv))
     return ps_mid,vl,vv
 end
 
 function volume_from_spinodal(p,poly,vshift,v0)
     f(v) = p - evalpoly(v,poly)
+
+
+    if length(v0) == 2
+        v1,v2 = v0
+        f1,f2 = f(v1),f(v2)
+        if f1*f2 < 0
+            prob = Roots.ZeroProblem(f,v0)
+            return Roots.solve(prob) + vshift
+        else
+            #something really wrong happened with the bracketing,
+            #hopefully the hermite polynomial should reproduce the EoS in a vicinity of the
+            #interpolated region.
+            if abs(f1) < abs(f2)
+                prob = Roots.ZeroProblem(f,v1)
+                return Roots.solve(prob) + vshift
+            else
+                prob = Roots.ZeroProblem(f,v2)
+                return Roots.solve(prob) + vshift
+            end
+        end
+    end
     prob = Roots.ZeroProblem(f,v0)
     return Roots.solve(prob) + vshift
 end

diff --git a/test/test_methods_api_flash.jl b/test/test_methods_api_flash.jl
@@ -340,6 +340,9 @@ end
     @test Clapeyron.saturation_temperature(cPR("isobutane"),1.7855513185537157e6,crit_retry = false)[1] ≈ 366.52386488214876 rtol = 1e-6
     @test Clapeyron.saturation_temperature(cPR("propane"),2.1298218093361156e6,crit_retry = false)[1] ≈ 332.6046103831853 rtol = 1e-6
     @test Clapeyron.saturation_temperature(cPR("r134a"),2.201981727901889e6,crit_retry = false)[1] ≈ 344.6869001549851 rtol = 1e-6
+
+    #Issue 328
+    @test saturation_pressure(cPR("butane"),406.5487245045052)[1] ≈ 2.815259927796967e6 rtol = 1e-6
 end
 
 @testset "Tproperty" begin