Tproperty bugs

[Sush1090]

Hello,
Thank you for adding Tproperty solver to Clapyeron.jl and making it robust.

There are a few bugs that are part of the release. I am mentioning them in the following with the possible fixes (only the ones I am aware of).

The following block of code will run as expected

using Clapeyron
model = cPR(["ethane"],idealmodel=ReidIdeal)
T0 = 300; p0 = 101325; z = [1] ;s0 = entropy(model,p0,T0,z);
p1 = 2*p0;
T1 = Tproperty(model,p1,s0,z,entropy) 
s1 = entropy(model,p1,T1,z)
@show s1-s0

The following will fail

using Clapeyron
model = cPR(["ethane"],idealmodel=ReidIdeal)
T0 = 300; p0 = 101325; z = [5] ;s0 = entropy(model,p0,T0,z);
p1 = 2*p0;
T1 = Tproperty(model,p1,s0,z,entropy) 
s1 = entropy(model,p1,T1,z)
@show s1-s0

This is because Tproperty_pure does not consider z to be anything apart form [1].

The other bug I am aware of is x0_Tproperty has an assert statement to keep the sum of moles to equal 1. I think this assertion can be removed.

The last one I know is that verbose = true does not show verbose in the REPL. This is because verbose is not passed to the underlying constructor of the function _Tproperty.

These are the ones that I know.
Thank you for implementing this functionality.

[Sush1090]

There is another small bug in FindEdge
I think it should be

    if (abs(∇f2) > abs(∇f1))
      FindEdge(f,c,b)
    else
      FindEdge(f,a,c)
    end

currently it is

    if (∇f2 > ∇f1)
      FindEdge(f,c,b)
    else
      FindEdge(f,a,c)
    end

The current version will only work when the value of property is increasing as temperature is increasing that is why it works for generic cases like entropy or enthalpy but for say mass_density it will not always work. Specifically when temperature lies between bubble and dew point

[Sush1090]

Hello,
As of now the implementation of Tproperty will return NaN when the pressure is super-critical.
Example:

using Clapeyron

# How to simulate super critical system?
fluids = ["R134a"]
model = cPR(fluids,idealmodel= ReidIdeal)

p_crit = crit_pure(model)[2]; T_crit = crit_pure(model)[1]

T1 = 300; p1 = saturation_pressure(model,T1)[1]+101325
s1 = entropy(model,p1,T1); h1 = enthalpy(model,p1,T1)

p2 = p_crit + 0.1; T2 = Tproperty(model,p2,s1,[1],entropy,verbose = true)
s2 = entropy(model,p2,T2);
@show s2 - s1 

This is because VT_identify_phase will return :unknown phase.

But a generic solution to this without using Tproperty exists:

using Roots
ff(T) = s1 - entropy(model,p2,T)
prob = ZeroProblem(ff,300.0)
sol = solve(prob)

@show entropy(model,p2,sol) - s1

So probably for this is it better to have something like a PT_identify_phase?

[Sush1090]

The fix I have for now is the following:

  if (p > Pc)
    verbose && @info "pressure is above critical pressure"
    return __Tproperty(model,p,prop,z,property,rootsolver,phase,abstol,reltol,threaded,Tc)
  end

Need to change the initialization of to critical temperature. This is only for Tproperty_pure. Need to add this after checking #case1.

Example:

fluids = ["R134A"]
model = cPR(fluids,idealmodel= ReidIdeal)

p_crit = crit_pure(model)[2]; T_crit = crit_pure(model)[1]

T1 = 300; p1 = saturation_pressure(model,T1)[1]+101325
s1 = entropy(model,p1,T1); h1 = enthalpy(model,p1,T1)

p2 = p_crit + 2*101325; T2 =  Tproperty(model,p2,s1,[1],entropy,verbose = true)
s2 = entropy(model,p2,T2);h2 = enthalpy(model,p2,T2)
@show s2 - s1 

[ Info: pressure is above critical pressure
s2 - s1 = 1.4210854715202004e-14
1.4210854715202004e-14

Probably will need to change it for mixtures as well.

[longemen3000]

hmm, maybe using the pressure extrapolation critical extrapolation makes sense in this particular case, but otherwise, the if case seems correct.

I also remembered that you wanted to calculate flashes with volume (or densities). I just recently added some special cases for those (using the volume obtained by the bubble/dew calculations as the property and using a linear interpolation of the temperatures when T (or p) is inside the saturation dome. With those changes, I did implement vt_flash

[Sush1090]

I think you are right. The "if" seems correct.
It is simply possible that for super critical pressure we need to initialize the root finder to critical temperature.
But one cannot be sure for all cases.
This might also need some thinking for mixtures.

For now this change is working fine for most cases with single components.

  if (p > Pc)
    verbose && @info "pressure is above critical pressure"
    return __Tproperty(model,p,prop,z,property,rootsolver,phase,abstol,reltol,threaded,Tc)
  end

Thank you for those implementations with vt_flash