Heat exchanger in part load: Singularity in jacobian matrix

[YoannUniKS]

Dear all,

I use tespy to calculate the kA value of a gas/water heat exchanger in part load conditions. Here is a minimal code I'm using to define the heat exchanger at design conditions, following the tespy documentation:

from tespy.components import Sink, Source, HeatExchanger
from tespy.connections import Connection
from tespy.networks import Network

nw = Network(fluids=['water', 'methane'], T_unit='C', p_unit='bar', m_unit='t / h', v_unit='m3 / h', h_unit='kJ / kg', iterinfo=False)

# Defining components
pri_in = Source('Primary side inlet')
pri_out = Sink('Primary side outlet')
sec_in = Source('Secondary side inlet')
sec_out = Sink('Secondary side outlet')
he = HeatExchanger('heat exchanger')

# Defining connections
pri_he = Connection(pri_in, 'out1', he, 'in1')
he_pri = Connection(he, 'out1', pri_out, 'in1')
sec_he = Connection(sec_in, 'out1', he, 'in2')
he_sec = Connection(he, 'out2', sec_out, 'in1')

nw.add_conns(pri_he, he_pri, sec_he, he_sec)

#Gas side
pu = 80 #bar
Tu = 5 #°C
Td = 40 #°C
dpg = 0.3 # bar

#Water side
pin = 2.5 #bar
Tin = 55 #°C
Tout = 35 #°C
dpw = 0.4 #bar

Vdot = 1000

# set design attributes
he.set_attr(pr1=(pin-dpw)/pin, pr2=(pu-dpg)/pu, ttd_u=Tin-Td, design=['pr1', 'pr2', 'ttd_u'], offdesign=['zeta1', 'zeta2', 'kA_char'])

pri_he.set_attr(fluid={'water': 1, 'methane':0}, state='l', p=pin, T=Tin, design=['v'])
sec_he.set_attr(fluid={'water': 0, 'methane': 1}, state='g', T=Tu, v=Vdot, p=pu)
he_pri.set_attr(T=Tout, design=['T'])

nw.solve(mode='design')
nw.save('design_test')

print('Design: \nVdot (gas): {a:0.1f} m³/h'.format(a=sec_he.v.val), '\nT (gas), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=sec_he.T.val, b=he_sec.T.val), '\nVdot (water): {a:0.1f} m³/h'.format(a=pri_he.v.val), '\nT (water), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=pri_he.T.val, b=he_pri.T.val), '\nkA: {a:0.1f} W/K'.format(a=he.kA.val))

Design:
Vdot (gas): 1000.0 m³/h
T (gas), in: 5.0 °C, out: 40.0 °C
Vdot (water): 81.2 m³/h
T (water), in: 55.0 °C, out: 35.0 °C
kA: 85904.0 W/K

So far so good. The design point is defined. I now want to evaluate the kA value under part load (reduced flow rate on the gas side). Here the flow rate is reduced at 50% of its nominal value, other parameters are kept constant:

Vdot = 500

pri_he.set_attr(T=55)
he_sec.set_attr(T=40)
sec_he.set_attr(v=Vdot, T=5, p=80)

nw.solve(mode='offdesign', design_path='design_test')

print('\nOffdesign: \nVdot (gas): {a:0.1f} m³/h'.format(a=sec_he.v.val), '\nT (gas), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=sec_he.T.val, b=he_sec.T.val), '\nVdot (water): {a:0.1f} m³/h'.format(a=pri_he.v.val), '\nT (water), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=pri_he.T.val, b=he_pri.T.val), '\nkA: {a:0.1f} W/K'.format(a=he.kA.val))

Offdesign:
Vdot (gas): 500.0 m³/h
T (gas), in: 5.0 °C, out: 40.0 °C
Vdot (water): 26.9 m³/h
T (water), in: 55.0 °C, out: 24.9 °C
kA: 53547.7 W/K

Here again, no problem, the code works. However when the flow rate is too low, as here below, 16 % of the nominal flow rate, I get an error:

Vdot = 160

pri_he.set_attr(T=55)
he_sec.set_attr(T=40)
sec_he.set_attr(v=Vdot, T=5, p=80)

nw.solve(mode='offdesign', design_path='design_test')

print('\nOffdesign: \nVdot (gas): {a:0.1f} m³/h'.format(a=sec_he.v.val), '\nT (gas), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=sec_he.T.val, b=he_sec.T.val), '\nVdot (water): {a:0.1f} m³/h'.format(a=pri_he.v.val), '\nT (water), in: {a:0.1f} °C, out: {b:0.1f} °C'.format(a=pri_he.T.val, b=he_pri.T.val), '\nkA: {a:0.1f} W/K'.format(a=he.kA.val))

Singularity in jacobian matrix, calculation aborted! Make sure your network does not have any linear dependencies in the parametrisation. Other reasons might be
-> given temperature with given pressure in two phase region, try setting enthalpy instead or provide accurate starting value for pressure.
-> given logarithmic temperature differences or kA-values for heat exchangers,
-> support better starting values.
-> bad starting value for fuel mass flow of combustion chamber, provide small (near to zero, but not zero) starting value.

I have two questions regarding these results.
1/ It is unclear what is causing the error in the last case, is there a way to "help" the algorithm find a solution?
2/ I am actually calculating the kA value for several thousand operating points and I get the message above for maybe 10 % of them. From nw.solve() I'm getting the error message, but it is of no practical use when working with so many values. Is there an output variable from the function that I can use telling me that the calculation was not successful? I could not find that. I would use it in a first step to sort out the points where the calculation did not work.

Thanks in advance!

[fwitte]

Hi @YoannUniKS,

thank you for reaching out and reporting. I can have a closer look at this issue later today or tomorrow.

A general tip for improvement of numerical stability: change the parameters gradually from simulation to simulation. Have you tried that? Also, could you tell me more on the structure of your runs? How do you call the simulations, in some loops? Maybe I can give you more tips on that.

Finally, the attribute nw.lin_dep should be False, if a simulation did not run into any linear dependency and you can check if nw.res[-1] is smaller than 1e-3 after the simulation (residual value of the last iteration performed) to see, if the solver converged.

I think an attributed.converged would be a good addition to TESPy, I‘ll put it on the todos.

best

Francesco

[YoannUniKS]

Dear Francesco,

Thank you for your reply and willing to help!

I am using time series of the known varying parameters, i.e. on the gas side inlet flow rate, inlet pressure and inlet and outlet temperatures and on the water side inlet temperature. It is difficult to vary the parameters gradually as they all vary at the same time. Pressures and temperatures vary however in a relatively small range (+/-10 bar & +/- 15 K). Convergence problems occur when the gas flow rate is low compared to the design flow rate. I will try tomorrow to sort the rows by descending gas flow rate and see if there is an improvement.

The runs are not really optimised. All varying input parameters are in a dataframe (data), and I loop over the rows of the dataframe, as follow:

for i in range(x1,x2):
    he_sec.set_attr(T=data["Td"][i])
    sec_he.set_attr(v=data["Vdot"][i], T=data["Tu"][i], p=data["pu"][i])
    pri_he.set_attr(T=data["Tin"][i])
    nw.solve(mode='offdesign', design_path='design_test')

Then I pack the results from nw.solve() I need in another dataframe. It takes several minutes for a few thousands of runs. Any tip on how to improve that is of course welcome!

Thanks for mentioning both indicators nw.lin_dep and nw.res[-1]. From a quick check there seem to be a good correlation between nw.lin_dep=True and nw.res[-1] > 1e-3.

All the best.

[YoannUniKS]

Short update: I sorted the dataframe by descending gas flow rate and run the calculations again. It did not change the result. Convergence was not achieved for the exact same number of data points as without sorting the dataframe.

[fwitte]

Well, then there might be other issues. I cannot really test this without the data you have. You could join me tomorrow at 17:00 CET to discuss the issue, if you like: https://tespy.readthedocs.io/en/main/regular_meeting.html#online-stammtisch.

Other than that, only changing Vdot works for me:

import numpy as np

# [...]

for Vdot in np.arange(10, 500)[::-1]:
    sec_he.set_attr(v=Vdot)

    nw.solve(mode='offdesign', design_path='design_test')
    print(he.kA.val)

[YoannUniKS]

Ok, thanks, I will join the meeting this evening.

[YoannUniKS]

Thank you Francesco for your time.

So for anyone who has the same issue. The idea is to avoid "jumps" in the varied parameters (pressure, temperature, flow rate) and to start from a point where the solver converges and then reduce the varied parameters in small steps until the desired point. This helps the solver converge.