Pressure difference across CycleCloser for CO2 loop

[mvernacc]

First, thank you for the effort of creating TESpy -- it looks like a great tool.

I am trying model closed-cycle networks with gaseous CO2 as the working fluid. Per the advice in the documentation to "set up your system step by step", I've started with a simple model: A compressor circulates CO2 through a heat exchanger and the heat exchanger dissipates the compressor's pressure rise and rejects heat from the compressor work. The network also includes a CycleCloser.

When solving the model, I expect the pressures on either side of the CycleCloser to be equal. However, the solver fails to converge, and returns different pressures at the inlet and outlet of the CycleCloser.

Interestingly, this error only occurs with CO2, and goes away if I change the fluid to another gas.

I would like help figuring out what is wrong with this model. Have I not properly initialized the component and connection attributes? Am I using the CycleCloser improperly? Or is this a bug in TESpy?

Here is the script:

from tespy.components import CycleCloser, Compressor, HeatExchangerSimple
from tespy.connections import Connection
from tespy.networks import Network
from CoolProp.CoolProp import PropsSI

# Build the network
# -----------------
# Diagram of components and connections:
#
#             +------------+
#             ↑            |
#            ▷             |
#          / ↑             ↓
# compressor +-*--+--\/\/\-+
#           /   \  \   \    \
#          c1    \ c0   \    c2
#                 \     cooler
#              cycle closer
#

# The error only occurs if the fluid is CO2. With N2, NH3, He or Ar the solver converges, and
# pressure c0 == pressure c1  == 4.0 bar.
fluid = "CO2"
nw = Network(fluids=[fluid], p_unit="bar", T_unit="K", m_unit="g / s")

# Components
closer = CycleCloser("cycle closer")
compressor = Compressor("compressor")
cooler = HeatExchangerSimple("cooler")

# Connections
c0 = Connection(cooler, "out1", closer, "in1", label="0")
c1 = Connection(closer, "out1", compressor, "in1", label="1")
c2 = Connection(compressor, "out1", cooler, "in1", label="2")
connections = [c0, c1, c2]
nw.add_conns(*connections)

# Set attributes
# --------------
compressor.set_attr(eta_s=0.80)
c1.set_attr(fluid={fluid: 1.0})

# Cooler outlet conditions
p_cooler_bar = 4.0  # [bar] cooler outlet pressure
T_cooler = 273.15 + 35.0  # [K] cooler outlet temperature
h_cooler = PropsSI("H", "P", 1e5 * p_cooler_bar, "T", T_cooler, fluid)  # [J/kg] cooler outlet enthalpy
c1.set_attr(p=p_cooler_bar, h=h_cooler)

# Compressor outlet conditions
p_comp_bar = 10.0  # [bar] compressor outlet pressure
c2.set_attr(p=p_comp_bar)

cooler.set_attr(Q=-1e3)
# If I size the system by mass flow, instead of power, the same error still happens.
# c1.set_attr(m=1.0)

# Guess starting values to help the solver
# ----------------------------------------
# Setting initial values for the fluid around the loop does not help.
# for c in connections:
#     c.set_attr(fluid0={fluid: 1.0})

# Setting initial values for p,h on the other side of the CycleCloser does not help.
# c0.set_attr(p0=p_cooler_bar, h0=h_cooler)

# Solve
# -----
nw.solve("design")
nw.print_results()

And here is the output (using python 3.11.1, tespy 0.6.2):

iter    | residual | massflow | pressure | enthalpy | fluid
--------+----------+----------+----------+----------+---------
1       | 5.19e+05 | 3.91e-03 | 3.00e+05 | 4.29e+05 | 0.00e+00
2       | 1.18e+05 | 9.42e-03 | 1.18e+05 | 6.86e-11 | 0.00e+00
3       | 1.18e+05 | 3.73e-18 | 1.18e+05 | 3.64e-11 | 0.00e+00
4       | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
5       | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
6       | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
7       | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
8       | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
9       | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
10      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
11      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
12      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
13      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
14      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
15      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
16      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
17      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
18      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
19      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
20      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
21      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
22      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
23      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
24      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
25      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
26      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
27      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
28      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
29      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
30      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
31      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
32      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
33      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
34      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
35      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
36      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
37      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
38      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
39      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
40      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
41      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
42      | 1.18e+05 | 1.16e-19 | 1.18e+05 | 3.64e-11 | 0.00e+00
--------+----------+----------+----------+----------+---------
Total iterations: 42, Calculation time: 0.0 s, Iterations per second: 1429.71
WARNING:root:The solver does not seem to make any progress, aborting calculation. Residual value is 1.18e+05. This frequently happens, if the solver pushes the fluid properties out of their feasible range.
##### RESULTS (HeatExchangerSimple) #####
+--------+-----------+----------+----------+-----+-----+------+------+--------+
|        |         Q |       pr |     zeta |   D |   L |   ks |   kA |   Tamb |
|--------+-----------+----------+----------+-----+-----+------+------+--------|
| cooler | -1.00e+03 | 3.45e-01 | 1.76e+11 | nan | nan |  nan |  nan |    nan |
+--------+-----------+----------+----------+-----+-----+------+------+--------+
##### RESULTS (CycleCloser) #####
+--------------+------------------+-------------------+
|              |   mass_deviation |   fluid_deviation |
|--------------+------------------+-------------------|
| cycle closer |         0.00e+00 |          0.00e+00 |
+--------------+------------------+-------------------+
##### RESULTS (Compressor) #####
+------------+----------+----------+----------+--------+
|            |        P |    eta_s |       pr |   igva |
|------------+----------+----------+----------+--------|
| compressor | 1.00e+03 | 8.00e-01 | 3.75e+00 |    nan |
+------------+----------+----------+----------+--------+
##### RESULTS (Connection) #####
+----+-----------+-----------+-----------+-----------+
|    |         m |         p |         h |         T |
|----+-----------+-----------+-----------+-----------|
|  0 | 9.174e+00 | 5.180e+00 | 5.118e+05 | 3.093e+02 |
|  1 | 9.174e+00 | 4.000e+00 | 5.118e+05 | 3.081e+02 |
|  2 | 9.174e+00 | 1.500e+01 | 6.208e+05 | 4.310e+02 |
+----+-----------+-----------+-----------+-----------+

Note that the pressure at c0 is 5.180 bar, whereas the pressure at c1 is 4.0 bar. I would expect the pressure at both c0 and c1 to be 4.0 bar.

I noticed that the warning message says "This frequently happens, if the solver pushes the fluid properties out of their feasible range." However, the returned pressure and temperature range (4-15 bar, 308-431 K) are well within the "feasible range" for gaseous CO2.

[fwitte]

Dear @mvernacc,

thank you for reaching out and sharing your code snippet, I hope you can make use of the software!

To your issue: It looks like CoolProp only supports pressure higher than triple point pressure for CO2 (also see http://www.coolprop.org/fluid_properties/fluids/CarbonDioxide.html#fluid-carbondioxide and http://www.coolprop.org/_downloads/9a4fdd673eb02b49856d14485b27e6e9/CarbonDioxide.pdf):

>>> print(PropsSI("P_min", "CO2"))  # minimum possible pressure in CoolProp in Pa
517964.34344772575

If you specify a higher pressure than that value at the compressor's inlet, the problem can be solved. To also quickly answer some of your questions:

When solving the model, I expect the pressures on either side of the CycleCloser to be equal. However, the solver fails to converge, and returns different pressures at the inlet and outlet of the CycleCloser.

The outcome of the equations (not necessarily all of them) applied by the different components will be wrong, if the solver fails to converge. That is why you see the pressure differences. TESPy does not allow the pressure to be lower than the CoolProp provided minimum pressure. Since the pressure c0 will be a result of the equation of the CycleCloser and the specified pressure c1, the solver tries to find a valid value for c0, but is limited to the CoolProp value range.

I would like help figuring out what is wrong with this model. Have I not properly initialized the component and connection attributes? Am I using the CycleCloser improperly? Or is this a bug in TESpy?

Overall looks very good, one small hints for your models:

Starting values for fluid composition are rarely required. With single fluid networks you do not need any actually. In case some connections miss a starting value for the fluid composition, a warning should be printed in the preprocessing. Then you can make use of fluid0 to try and fix it.

Have a nice sunny weekend!

Best

Francesco