From 7c8e25577864bf07fa1f727542f7ac0cb80df74c Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Fri, 20 Dec 2024 15:46:58 -0300 Subject: [PATCH] fix error in #320 --- src/methods/property_solvers/Pproperty.jl | 9 +++++---- .../multicomponent/flash/general_flash.jl | 20 +++++++++---------- test/test_methods_api.jl | 7 +++++++ 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/methods/property_solvers/Pproperty.jl b/src/methods/property_solvers/Pproperty.jl index 5ef0f133b..d703f675f 100644 --- a/src/methods/property_solvers/Pproperty.jl +++ b/src/methods/property_solvers/Pproperty.jl @@ -196,6 +196,7 @@ function Pproperty_pure(model,T,prop,z,property::F,rootsolver,phase,abstol,relto Psat,vlsat,vvpat = saturation_pressure(model,T,crit = crit) if !is_unknown(phase) + return __Pproperty(model,T,prop,z,property,rootsolver,phase,abstol,reltol,threaded,Psat) end @@ -235,13 +236,13 @@ function __Pproperty(model,T,prop,z,property::F,rootsolver,phase,abstol,reltol,t end return __Pproperty(model,T,prop,z,property,rootsolver,new_phase,abstol,reltol,threaded,p0) end - f(p,prop) = property(model,p,T,z,phase = phase,threaded = threaded) - prop - prob = Roots.ZeroProblem(f,p0) + f(lnp,prop) = property(model,exp(lnp),T,z,phase = phase,threaded = threaded) - prop + prob = Roots.ZeroProblem(f,log(p0)) sol = Roots.solve(prob,rootsolver,p = prop,atol = abstol,rtol = reltol) - if !isfinite(sol) || sol < 0 + if isnan(sol) return sol,:failure end - return sol,phase + return exp(sol),phase end #= model = PCSAFT(["propane","dodecane"]) diff --git a/src/methods/property_solvers/multicomponent/flash/general_flash.jl b/src/methods/property_solvers/multicomponent/flash/general_flash.jl index d9e7cd0e9..f2f58586d 100644 --- a/src/methods/property_solvers/multicomponent/flash/general_flash.jl +++ b/src/methods/property_solvers/multicomponent/flash/general_flash.jl @@ -386,7 +386,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)/val1 + output[end-1] = (val1 - ∑a - p_end*∑v - T*∑s)/RTinv end if spec2 == temperature @@ -399,7 +399,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)/val1 + output[end] = (val2 - ∑a - p_end*∑v - T*∑s)/RTinv end return output end @@ -768,16 +768,16 @@ function px_flash_pure(model,p,x,z,spec::F,T0 = nothing) where F Ts,vl,vv = saturation_temperature(model,p) ∑z = sum(z) x1 = SA[1.0] - spec_to_vt(model,vl,T,x1,spec) - xl = ∑z*spec_to_vt(model,vl,T,x1,spec) - xv = ∑z*spec_to_vt(model,vv,T,x1,spec) + spec_to_vt(model,vl,Ts,x1,spec) + xl = ∑z*spec_to_vt(model,vl,Ts,x1,spec) + xv = ∑z*spec_to_vt(model,vv,Ts,x1,spec) βv = (x - xl)/(xv - xl) if !isfinite(βv) return FlashResultInvalid(1,βv) elseif βv < 0 || βv > 1 phase0 = βv < 0 ? :liquid : :vapour - T,_phase = _Tproperty(model,p,h/∑z,SA[1.0],spec,T0 = T0,phase = phase0) - return FlashResult(model,p,T,∑z,phase = _phase) + T,_phase = _Tproperty(model,p,x/∑z,SA[1.0],spec,T0 = T0,phase = phase0) + return FlashResult(model,p,T,SA[∑z*one(p)*one(T)],phase = _phase) else build_flash_result_pure(model,p,Ts,z,vl,vv,βv) end @@ -811,10 +811,10 @@ function tx_flash_pure(model,T,x,z,spec::F,P0 = nothing) where F return FlashResultInvalid(1,βv) elseif βv < 0 || βv > 1 phase0 = βv < 0 ? :liquid : :vapour - p,_phase = _Pproperty(model,T,x/∑z,SA[1.0],spec,p0 = P0,phase = phase0) - return FlashResult(model,p,T,∑z,phase = _phase) + p,_phase = _Pproperty(model,T,x/∑z,SA[1.0],spec,p0 = P0) + return FlashResult(model,p,T,SA[∑z*one(p)*one(T)],phase = _phase) else - build_flash_result_pure(model,p,Ts,z,vl,vv,βv) + build_flash_result_pure(model,p,T,z,vl,vv,βv) end end diff --git a/test/test_methods_api.jl b/test/test_methods_api.jl index ea5bfe7ae..8a2fee009 100644 --- a/test/test_methods_api.jl +++ b/test/test_methods_api.jl @@ -414,6 +414,13 @@ end result = xy_flash(model,spec,z,result0) #perform the flash @test Clapeyron.temperature(result) == 200.15 @test pressure(result) == 101325.0 + + #px_flash_pure/tx_flash_pure, 1 phase (#320) + model = cPR(["ethane"],idealmodel = ReidIdeal); + p = 101325;h = 100; z = Clapeyron.SA[1]; T = Clapeyron.PH.temperature(model,p,h,z) + @test enthalpy(model,p,T,z) ≈ h rtol = 1e-6 + res5 = Clapeyron.tx_flash_pure(model,T,h,z,enthalpy) + @test pressure(res5) ≈ p rtol = 1e-6 end @testset "Saturation Methods" begin