From 67a8c512379986c9e9d38411d2f2ef1b25662927 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Thu, 19 Dec 2024 00:02:33 -0300 Subject: [PATCH] add a bunch of tests for flash (#320) --- src/methods/PH.jl | 7 +--- src/methods/PS.jl | 29 +++++++++------ .../property_solvers/multicomponent/flash.jl | 6 ++++ .../multicomponent/flash/QT.jl | 23 ++++++++---- .../multicomponent/flash/general_flash.jl | 18 +++++----- .../tp_flash/Michelsentp_flash.jl | 3 +- .../tp_flash/RachfordRicetp_flash.jl | 4 +++ .../multicomponent/tp_flash/multiphase.jl | 4 ++- test/test_methods_api.jl | 35 +++++++++++++++++++ 9 files changed, 95 insertions(+), 34 deletions(-) diff --git a/src/methods/PH.jl b/src/methods/PH.jl index 6d9f46086..9a4000159 100644 --- a/src/methods/PH.jl +++ b/src/methods/PH.jl @@ -9,7 +9,6 @@ function PH_property(model,p,h,z,f::F,phase,T0,threaded) where F if !is_unknown(phase) T,calc_phase = _Tproperty(model,p,h,z,T0 = T0,phase = phase,threaded = threaded) - #TODO: handle equilibria conditions if calc_phase != :eq && calc_phase != :failure return f(model,p,T,z;phase = calc_phase) elseif calc_phase == :eq && !supports_lever_rule(f) @@ -22,14 +21,10 @@ function PH_property(model,p,h,z,f::F,phase,T0,threaded) where F end else res = ph_flash(model,p,h,z,T0 = T0) - np = numphases(res) if f == temperature return temperature(res) - end - if np > 1 - return f(model,res) else - return f(model,p,T,z) + return f(model,res) end end end diff --git a/src/methods/PS.jl b/src/methods/PS.jl index 1b0e1a9f0..a2fdee369 100644 --- a/src/methods/PS.jl +++ b/src/methods/PS.jl @@ -6,17 +6,26 @@ function PS_property(model,p,s,z,f::F,phase,T0,threaded) where F if f == pressure return p end - T,calc_phase = _Tproperty(model,p,s,z,entropy,T0 = T0,phase = phase,threaded = threaded) - #TODO: handle equilibria conditions - if calc_phase != :eq || calc_phase != :failure - return f(model,p,T,z;phase = calc_phase) - elseif calc_phase == :eq && !supports_lever_rule(f) - thow(invalid_property_multiphase_error(f)) - elseif calc_phase == :eq && supports_lever_rule(f) && length(model) == 1 #TODO: lift single component condition - result = ps_flash(model,p,s,z,T0) - return f(model,result) + + if !is_unknown(phase) + T,calc_phase = _Tproperty(model,p,h,z,entropy,T0 = T0,phase = phase,threaded = threaded) + if calc_phase != :eq && calc_phase != :failure + return f(model,p,T,z;phase = calc_phase) + elseif calc_phase == :eq && !supports_lever_rule(f) + thow(invalid_property_multiphase_error(f)) + elseif calc_phase == :eq && supports_lever_rule(f) + result = ps_flash(model,p,s,z,T0) + return f(model,result) + else + return f(model,p,T,z;phase = phase) + end else - return f(model,p,T,z;phase = phase) + res = ps_flash(model,p,s,z,T0 = T0) + if f == temperature + return temperature(res) + else + return f(model,res) + end end end diff --git a/src/methods/property_solvers/multicomponent/flash.jl b/src/methods/property_solvers/multicomponent/flash.jl index e2a934c2b..e70d211f0 100644 --- a/src/methods/property_solvers/multicomponent/flash.jl +++ b/src/methods/property_solvers/multicomponent/flash.jl @@ -274,6 +274,12 @@ is_vle(method::FlashMethod) = is_vle(method.equilibrium) is_lle(method::FlashMethod) = is_lle(method.equilibrium) is_unknown(method::FlashMethod) = is_unknown(method.equilibrium) +@noinline function incorrect_np_flash_error(method,result) + np = numphases(result) + s = np == 1 ? "" : "s" + throw(ArgumentError("$method does not support an input with $np phase$s as an initial point. Got the following input: \n\n $result")) +end + include("flash/general_flash.jl") include("flash/PT.jl") include("flash/PH.jl") diff --git a/src/methods/property_solvers/multicomponent/flash/QT.jl b/src/methods/property_solvers/multicomponent/flash/QT.jl index 7777d6ea8..edde5c026 100644 --- a/src/methods/property_solvers/multicomponent/flash/QT.jl +++ b/src/methods/property_solvers/multicomponent/flash/QT.jl @@ -10,7 +10,7 @@ function qt_flash_x0(model,β,T,z,method::FlashMethod) x = z ./ sum(z) p,vl,vv,y = __x0_bubble_pressure(model,T,x) y ./= sum(y) - + return FlashResult(p,T,SA[x,y],SA[1.0-β,1.0*β],SA[vl,vv],sort = false) elseif 0.99 <= β <= 1.0 y = z ./ sum(z) @@ -20,12 +20,23 @@ function qt_flash_x0(model,β,T,z,method::FlashMethod) else pures = split_model(model) sat = extended_saturation_pressure.(pures,T) - ps = first.(sat) + ps = first.(sat) K = similar(ps) pmin,pmax = extrema(ps) - p0 = β*pmin + (1-β)*pmax fp(p) = qt_f0_p!(K,z,p,ps,β) - prob = Roots.ZeroProblem(fp,p0) + pm = β*pmin + (1-β)*pmax + pr1 = range(pmin,pm,5*length(model)) + pr2 = range(pm,pmax,5*length(model)) + δβ1 = abs.(fp.(pr1)) + δβ2 = abs.(fp.(pr2)) + δβ1_min,i1 = findmin(δβ1) + δβ2_min,i2 = findmin(δβ2) + if δβ1_min < δβ2_min + p00 = pr1[i1] + else + p00 = pr2[i2] + end + prob = Roots.ZeroProblem(fp,p00) p = Roots.solve(prob) end else @@ -40,7 +51,7 @@ function qt_flash(model::EoSModel,β,T,z;kwargs...) return qt_flash(model,β,T,z,method) end -function init_preferred_method(method::typeof(qt_flash),model::EoSModel,kwargs) +function init_preferred_method(method::typeof(qt_flash),model::EoSModel,kwargs) GeneralizedXYFlash(;kwargs...) end @@ -58,7 +69,7 @@ function qt_flash(model,β,T,z,method::FlashMethod) result1 = βflash_pure(model,temperature,T,βv,z) return index_expansion(result1,idx_r) end - + result = qt_flash_impl(model_r,β,T,z_r,method_r) if !issorted(result.volumes) #this is in case we catch a bad result. diff --git a/src/methods/property_solvers/multicomponent/flash/general_flash.jl b/src/methods/property_solvers/multicomponent/flash/general_flash.jl index 0b52d3512..b787cccc7 100644 --- a/src/methods/property_solvers/multicomponent/flash/general_flash.jl +++ b/src/methods/property_solvers/multicomponent/flash/general_flash.jl @@ -387,7 +387,7 @@ function xy_flash_neq(output,model,zbulk,np,input,state::F,μconfig) where F i_spec = spec1.k output[end-1] = β[i_spec] - val1 else #general caloric term, valid for enthalpy , entropy, internal energy, gibbs, helmholtz - output[end-1] = (val1 - ∑a - p_end*∑v - T*∑s)*RTinv + output[end-1] = (val1 - ∑a - p_end*∑v - T*∑s)/val1 end if spec2 == temperature @@ -400,7 +400,7 @@ function xy_flash_neq(output,model,zbulk,np,input,state::F,μconfig) where F i_spec = spec2.k output[end-1] = β[i_spec] - val2 else - output[end] = (val2 - ∑a - p_end*∑v - T*∑s)*RTinv + output[end] = (val2 - ∑a - p_end*∑v - T*∑s)/val1 end return output end @@ -431,7 +431,7 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,flash::FlashResult return xy_flash(model,spec,z,comps0,β0,volumes0,T0;rtol,atol,max_iters) end -function xy_flash(model::EoSModel,spec::FlashSpecifications,z,flash::FlashResult;rtol = 1e-12,atol = 1e-10,max_iters = 50) +function xy_flash(model::EoSModel,spec::FlashSpecifications,z,flash::FlashResult;rtol = 1e-14,atol = 1e-12,max_iters = 50) if numphases(flash) == 1 throw(ArgumentError("xy_flash cannot use single phase initial points as starting points.")) end @@ -542,7 +542,6 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,comps0,β0,volumes x_old .= x #bound positivity α0 = Solvers.positive_linesearch(x,s) - #backtrack linesearch, so the next result is strictly better than the last α,Θx = Solvers.backtracking_linesearch!(Θ,F,x_old,s,Θx,x,α0) Fnorm = sqrt(2*Θx) @@ -550,7 +549,6 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,comps0,β0,volumes snorm_old = snorm snorm = α*norm(s,2) δs = max(abs(snorm-snorm_old),norm((F[end-1],F[end]),Inf)) - Fnorm = norm(F,Inf) δs < rtol && (converged = true) xnorm = Solvers.dnorm(x,x_old,Inf) @@ -561,11 +559,10 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,comps0,β0,volumes isnan(δs) && (nan_converged = true) i == max_iters && (max_iters_reached = true) end - + max_iters_reached,converged,nan_converged if !converged || nan_converged || max_iters_reached x .= NaN end - if spec_norm > sqrt(rtol) x .= NaN end @@ -639,14 +636,15 @@ function GeneralizedXYFlash(;equilibrium = :unknown, x0 = nothing, y0 = nothing, v0 = nothing, - rtol = 1e-12, - atol = 1e-10, + rtol = 1e-14, + atol = 1e-12, max_iters = 100, flash_result = nothing) !(is_vle(equilibrium) | is_lle(equilibrium) | is_unknown(equilibrium)) && throw(error("invalid equilibrium specification for GeneralizedXYFlash")) if flash_result isa FlashResult comps,β,volumes = flash_result.compositions,flash_result.fractions,flash_result.volumes - @assert length(comps) == 2 + np = numphases(flash_result) + np != 2 && incorrect_np_flash_error(GeneralizedXYFlash,flash_result) w1,w2 = comps[1],comps[2] v = (volumes[1],volumes[2]) P00 = flash_result.data.p diff --git a/src/methods/property_solvers/multicomponent/tp_flash/Michelsentp_flash.jl b/src/methods/property_solvers/multicomponent/tp_flash/Michelsentp_flash.jl index fbe7c7a57..437ea9952 100755 --- a/src/methods/property_solvers/multicomponent/tp_flash/Michelsentp_flash.jl +++ b/src/methods/property_solvers/multicomponent/tp_flash/Michelsentp_flash.jl @@ -61,7 +61,8 @@ function MichelsenTPFlash(;equilibrium = :unknown, if flash_result isa FlashResult comps,β,volumes = flash_result.compositions,flash_result.fractions,flash_result.volumes - @assert length(comps) == 2 + np = numphases(flash_result) + np != 2 && incorrect_np_flash_error(MichelsenTPFlash,flash_result) w1,w2 = comps[1],comps[2] v = (volumes[1],volumes[2]) return Michelsentp_flash(;equilibrium,x0 = w1,y0 = w2,vol0 = v,K_tol,ss_iters,nacc,second_order,noncondensables,nonvolatiles) diff --git a/src/methods/property_solvers/multicomponent/tp_flash/RachfordRicetp_flash.jl b/src/methods/property_solvers/multicomponent/tp_flash/RachfordRicetp_flash.jl index 4658c354c..537ad929c 100755 --- a/src/methods/property_solvers/multicomponent/tp_flash/RachfordRicetp_flash.jl +++ b/src/methods/property_solvers/multicomponent/tp_flash/RachfordRicetp_flash.jl @@ -55,6 +55,10 @@ function RRTPFlash(;equilibrium = :unknown, nonvolatiles = nothing, flash_result = nothing) ss_iters = max_iters + if flash_result isa FlashResult + np = numphases(flash_result) + np != 2 && incorrect_np_flash_error(RRTPFlash,flash_result) + end #we call Michelsen to check if the arguments are correct. m = MichelsenTPFlash(;equilibrium,K0,x0,y0,v0,K_tol,ss_iters,nacc,noncondensables,nonvolatiles,flash_result) return RRTPFlash{eltype(m)}(m.equilibrium,m.K0,m.x0,m.y0,m.v0,m.K_tol,max_iters,m.nacc,m.noncondensables,m.nonvolatiles) diff --git a/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl b/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl index 1afa28e76..d4421a3e2 100644 --- a/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl +++ b/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl @@ -58,9 +58,11 @@ function MultiPhaseTPFlash(; if flash_result !== nothing comps = flash_result.compositions + np = numphases(flash_result) + np < 2 && incorrect_np_flash_error(MultiPhaseTPFlash,flash_result) ∑n = sum(flash_result.fractions) n00 = Vector{eltype(comps)}[] - for i in 1:numphases(flash_result) + for i in 1:np push!(n00,collect(comps[i])) end return MultiPhaseTPFlash(;K0 = nothing,x0 = nothing, y0 = nothing,n0 = n00,K_tol,ss_iters,nacc,second_order,max_phases,phase_iters) diff --git a/test/test_methods_api.jl b/test/test_methods_api.jl index 404536174..dd2d9a6d2 100644 --- a/test/test_methods_api.jl +++ b/test/test_methods_api.jl @@ -363,6 +363,41 @@ end end end +@testset "XY flash" begin + #1 phase (#320) + model = cPR(["ethane","methane"],idealmodel = ReidIdeal) + p = 101325.0 + z = [1.2,1.2] + T = 350.0 + h = enthalpy(model,p,T,z) + res0 = ph_flash(model,p,h,z) + @test Clapeyron.temperature(res0) ≈ T rtol = 1e-6 + @test enthalpy(model,res1) ≈ h rtol = 1e-6 + + #2 phases + h = -13831.0 + res1 = ph_flash(model,p,h,z) + @test enthalpy(model,res1) ≈ h rtol = 1e-6 + + + #examples for qt, qp flash (#314) + model = cPR(["ethane","propane"],idealmodel=ReidIdeal) + res2 = qt_flash(model,0.5,208.0,[0.5,0.5]) + @test Clapeyron.pressure(res2) ≈ 101634.82435966855 rtol = 1e-6 + res3 = qp_flash(model,0.5,120000.0,[0.5,0.5]) + @test Clapeyron.temperature(res3) ≈ 211.4972567716822 rtol = 1e-6 + + #1 phase input should error + model = PR(["IsoButane", "n-Butane", "n-Pentane", "n-Hexane"]) + z = [0.25, 0.25, 0.25, 0.25] + p = 1e5 + h = enthalpy(model, 1e5, 303.15, z) + r = Clapeyron.ph_flash(model, p, h, z) + @test_throws ArgumentError qt_flash(model,0.5,308,z,flash_result = r) + res4 = qp_flash(model,0.7,60000.0,z) + @test Clapeyron.temperature(res4) ≈ 289.6991395244328 rtol = 1e-6 +end + @testset "Saturation Methods" begin model = PR(["water"]) vdw = vdW(["water"])