ValueError: Unable to find gaseous density

[Luis Jimenez]

Hello Everyone,

I'm modelling components for Kalina cycle system (a power system that uses the mixture ammonia+water as work fluid) as part of my thesis, with Peng-Robinson backend, using the python wrapper. Also I need to do some plots.

I have CoolProp version 6.4.1.

I've attached the code here:

import CoolProp
backend = 'PR'
fluid = 'Ammonia&Water'
AS = CoolProp.AbstractState(backend, fluid) 

# set mole fraction of NH3
x = 0.5
AS.set_mole_fractions([x, 1 - x])

# set binary interaction parameter
AS.set_binary_interaction_double(0, 1, "kij", -0.2401)

# pressure in Pa and temperature in K
P = 3e6
T = 390

# solve for P,T inputs
AS.update(CoolProp.PT_INPUTS, P, T)

# print results
print("T = %.2f K" % AS.T())
print("P = %.1f kPa" % (AS.p()/1000))
print("rho = %.5f kg/m^3" % (AS.rhomass()))
print("H = %.2f kJ/kg" % (AS.hmass()/1000))
print("S = %.5f kJ/kg*K" % (AS.smass()/1000))
print("phase = %s" % (AS.phase()))

However this return the following error

ValueError                                Traceback (most recent call last)

<ipython-input-18-5a46a4779b09> in <module>()
     17 
     18 # solve for P,T inputs
---> 19 AS.update(CoolProp.PT_INPUTS, P, T)
     20 
     21 # print results

CoolProp/AbstractState.pyx in CoolProp.CoolProp.AbstractState.update()

CoolProp/AbstractState.pyx in CoolProp.CoolProp.AbstractState.update()

ValueError: Unable to find gaseous density for T: 390 K, p: 3e+06 Pa

In fact, return the same ValueError for any mole fraction in the interval 0.47 <= x <= 0.62.

But if we modify the composition, for example to x = 0.4, this return the following results

T = 390.00 K
P = 3000.0 kPa
rho = 667.73139 kg/m^3
H = 474.96 kJ/kg
S = 2.07423 kJ/kg*K
phase = 5

will there be any solution for this problem?

 

Thanks

Regards,

Luis Jiménez

[Friederike Boehm]

Hi Luis,

I am pretty sure you are having the same problem I had a couple months ago: https://groups.google.com/g/coolprop-users/c/dE5aK0GWHhA/m/mcs6efxKAwAJ

The inputs and concentrations land in the two phase region (see attached picture for x=0.5). For x=0.46 you land in the liquid phase, for x=0.63 you are still in the two-phase region. PT inputs do not properly work for mixtures when you are in the two phase region. For some reason, CoolProp does calculate something for higher qualities (>0.3) but as far as I can tell, this will not be correct. What you need to do is figure out, whether your PT inputs land in the two-phase region and if they do, use PQ inputs and rootfinding to figure out, which quality correlates with your temperature.

For the rootfinding part you can have a look at the Jupyter Notebook in the thread linked above (look for the text saying "PQ inputs and root finding T works") or at this issue on GitHub.

I only implemented the "phase check" for a given PT input in Mathematica so you'll need to write this yourself. I have listed the steps below. For me one problem was, that once I built the phase envelope, some update calls (PQ and QT) would lead to Assertion Failure. I get around that by having two Abstract States, one for building and reading the phase envelope, one for the actual calculating part.

Create functions for boiling and dew temperatures:

• create two AbstractState instances with the same components, composition, interaction parameters etc.
• build and read the phase envelope for one of them. Use setting "veryfine"
  
• split the pressure and temperature data at the maximum pressure (should be the critical point)
• this will give you a list each for the dew and the boiling points
• interpolate both lists separately to have functions TBoil(p) and TDew(p)

Use these functions to determine the phase:

• determine whether the pressure input is below the critical pressure (if not you can use PT inputs)
  
• calculate the dew temperature and the boiling temperatures
• check if temperature input is between these temperatures
  • if it is, specify the phase to be "twophase" and use rootfinding and PQ inputs so the calculated T will match the input temperature
    
  • if not, specify the phase to "liquid" or "gas" and use your PT inputs

One last thing about using Peng-Robinson for your calculations: there was a coding error in the ideal gas contribution (see issue on GitHub) which I fixed a couple months ago. This fix is not in the official release yet (I think) so you might need to compile the python wrapper yourself. The flash calculation should work without the fix, but all properties derived from the ideal gas contribution will be wrong.

Hope, this helps :)

Best regards,

Friederike