ClapeyronThermo / Clapeyron.jl

Clapeyron provides a framework for the development and use of fluid-thermodynamic models, including SAFT, cubic, activity, multi-parameter, and COSMO-SAC.
MIT License
194 stars 49 forks source link

dew temperature abnormality with water? #290

Open Sush1090 opened 2 weeks ago

Sush1090 commented 2 weeks ago

Hello,

I am interested in computing the bubble and dew temperatures for a given mixture with one component being water. I have the following code block

using Clapeyron, Plots
model = cPR(["water","ethanol"],idealmodel=ReidIdeal)
x = range(0,0.99,length=100)
X = Clapeyron.FractionVector.(x)
P = 101325

T_bubble = zeros(size(X,1))
T_dew = zeros(size(X,1))

for i in eachindex(X)
    T_bubble[i] = bubble_temperature(model,P,X[i],FugBubbleTemperature())[1]
    T_dew[i] = dew_temperature(model,P,X[i],FugDewTemperature())[1]
end 

fig1 = plot(x,T_dew,label = "Dew temperature")
plot!(x,T_bubble,label = "Bubble temperature")
xlabel!("Fraction 1st fluid")
ylabel!("Temperature (K)")

There seems like an abnormality in the output image

For most cases without water I have decent answer. I know that in some cases bubble_temperature and dew_temperature return NaN values but that is due to initialization of the non-linear solver.

But for such abnormality I am not sure what might be causing it.

Thank you

pw0908 commented 2 weeks ago

Hi @Sush1090,

This is indeed very odd. I quickly checked with an auto-plotter and I think the issue lies in FugBubbleTemperature / FugDewTemperature. For some reason, they are returning the initial guess for this mixture which leads to this discontinuity (@longemen3000 Im not sure why it's doing that). When I use our default method, I get a plot that looks like this: newplot (16)

While this plot looks worse, those 'loops' at lower temperatures are actually indicative of the presence of VLLE. This is very typical in a water+ethanol system when you don't include a fitted binary interaction parameter. You'll need to be careful when trying to make a 'clean' version of this plot.

Btw, you don't need to include ReidIdeal when doing phase equilibrium calculations.

Hope this helps!

Best regards,

Pierre

Sush1090 commented 1 week ago

Hello, Thank you for the insight. This helps me.

I have found another possible error with dew and bubble point calculations and maybe it is related to the above mentioned abnormality (Its just a guess -- I am not sure).

Here is my code block (same as before but different fluids)

using Clapeyron, Plots
model = cPR(["r134a","R1233zde"])
x = range(0,0.99,length=100)
X = Clapeyron.FractionVector.(x)
P = 101325*20

T_bubble = zeros(size(X,1))
T_dew = zeros(size(X,1))

for i in eachindex(X)
    T_bubble[i] = bubble_temperature(model,P,X[i],FugBubbleTemperature())[1]
    T_dew[i] = dew_temperature(model,P,X[i],FugDewTemperature())[1]
end 

fig1 = plot(x,T_dew,label = "Dew temperature")
plot!(x,T_bubble,label = "Bubble temperature")
xlabel!("Fraction 1st fluid")
ylabel!("Temperature (K)")

I get the following error:

ERROR: LoadError: saturation temperataure for PR{BasicIdeal, TwuAlpha, NoTranslation, vdW1fRule}("R1233zde") not found at p = 2.0265e6
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] fix_sat_ti!(sat::Vector{Tuple{…}}, pure::Vector{PR{…}}, crit::Vector{Tuple{…}}, p::Float64, in_media::Vector{Bool})
    @ Clapeyron C:\Users\sush9\.julia\packages\Clapeyron\vnjPG\src\methods\property_solvers\multicomponent\bubble_point.jl:108
  [3] __x0_bubble_temperature(model::PR{…}, p::Float64, x::Vector{…}, Tx0::Nothing, volatiles::Vector{…}, pure::Vector{…}, crit::Nothing)
    @ Clapeyron C:\Users\sush9\.julia\packages\Clapeyron\vnjPG\src\methods\property_solvers\multicomponent\bubble_point.jl:259
  [4] __x0_bubble_temperature
    @ C:\Users\sush9\.julia\packages\Clapeyron\vnjPG\src\methods\property_solvers\multicomponent\bubble_point.jl:253 [inlined]
  [5] bubble_temperature_init
    @ C:\Users\sush9\.julia\packages\Clapeyron\vnjPG\src\methods\property_solvers\multicomponent\bubble_point.jl:315 [inlined]
  [6] bubble_temperature_impl(model::PR{BasicIdeal, TwuAlpha, NoTranslation, vdW1fRule}, p::Float64, x::Vector{Float64}, method::FugBubbleTemperature{Nothing})
    @ Clapeyron C:\Users\sush9\.julia\packages\Clapeyron\vnjPG\src\methods\property_solvers\multicomponent\bubble_point\bubble_fugacity.jl:397
  [7] bubble_temperature(model::PR{…}, p::Int64, x::Clapeyron.Fractions.FractionVector{…}, method::FugBubbleTemperature{…})
    @ Clapeyron C:\Users\sush9\.julia\packages\Clapeyron\vnjPG\src\methods\property_solvers\multicomponent\bubble_point.jl:372
  [8] top-level scope
    @ C:\Users\sush9\.julia\dev\NonLinearSolverClapeyron\test.jl:11
  [9] include(fname::String)
    @ Base.MainInclude .\client.jl:489
 [10] top-level scope
    @ REPL[3]:1
in expression starting at C:\Users\sush9\.julia\dev\NonLinearSolverClapeyron\test.jl:10
Some type information was truncated. Use `show(err)` to see complete types.

If I change the pressure P = 101325*30 the solution seems accurate.

There is something different happening in fix_sat_ti! in the file bubble_point.jl. Computation of saturation_temperature(pure_i,p,crit = crit_i)[1] outputs NaN but if computed separately in the terminal I get an accurate value.

It is possible these these two are not related, in that case I'll make another issue.

Thanking you, Sushrut

pw0908 commented 1 week ago

Hi Sushrut,

This bug is unrelated and is honestly weird / frustrating. As I'm sure you've figured out, our saturation_temperature solver is failing for R1233zde at this specific pressure, but doesn't fail at a higher one. As such, it is an unrelated bug. If you want to close this issue and re-open a new one you can, but it seems @longemen3000 has already started to fix the bug in reference to this issue.

@longemen3000 Is it already fixed with these commits? The tests seem to be failing before reaching the new bug.

Best,

Pierre

longemen3000 commented 1 week ago

The saturation temperature bug is fixed, but there are other things failing right now, I plan to merge that branch as soon as possible