The order of pipe input in the add_conns function leads to different calculation results.

[CZRMWMW]

Dear Developer,I encountered a very peculiar issue. When I use nw.add_conns(L1, L2, L2_1, L2_2,L3_1,L3_2,c29,c30, c31, c32,L4,L5,L6,Z1,Z2) the returned result tells me that there are singular values in the calculation. 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. However, when I change the order to nw.add_conns(L1, L2, L2_1, L2_2,L3_1,L3_2,L4,L5,L6,Z1,Z2,c29,c30, c31, c32) it works just fine. I would like to know the reason for this. Thank you for your response. Here is the complete code, along with the corresponding data and process flow diagram.

import numpy as np
from tespy.tools.fluid_properties.wrappers import FluidPropertyWrapper
from tespy.tools.global_vars import gas_constants
from scipy.optimize import fsolve

COEF = {
    'N2': [6560.4875, -113.33376, 1.8053856, -0.0012660079, 1.3698985 * 10 ** -6, -7.1424530 * 10 ** -10,
           1.4957213 * 10 ** -13, -2.20796 * 10 ** 2],
    'O2': [-8900.9217, 125.9436, 0.29076177, 0.00055785819, -5.9211059 * 10 ** -8, -1.3143737 * 10 ** -10,
           5.3996454 * 10 ** -14, 8.72569 * 10 ** 2],
    'CO2': [9339.8089, -118.34503, 1.0016303, 0.00023651695, -1.3396751 * 10 ** -8, -3.6320845 * 10 ** -11,
            1.0767527 * 10 ** -14, -3.94483 * 10 ** 2],
    'H2O': [-18220.844, 265.64163, 0.4300414, 0.0016667329, -1.1295923 * 10 ** -6, 5.7171965 * 10 ** -10,
            -1.2340574 * 10 ** -13, 1.80771 * 10 ** 3],
    'Ar': [0, 0, 0.52033331, 0, 0, 0, 0, 150.225]
}

class KKHWrapper(FluidPropertyWrapper):
    print('KKHWrapper')
    def __init__(self, fluid, back_end=None,) -> None:
        super().__init__(fluid, back_end)
        if self.fluid not in COEF:
            msg = "Fluid not available in KKH database"
            raise KeyError(msg)
        self.coefficients = COEF[fluid]
        self._molar_mass = 1
        self._T_min = 100
        self._T_max = 2000
        self._p_min = 1000
        self._p_max = 10000000

    def h_pT(self, p, T):
        a = 6.744824499980115
        b = 6.779080644618959
        h = -(self.coefficients[0]) / T + (self.coefficients[1]) * a + self.coefficients[2] * T + (
        self.coefficients[3]) * T ** 2 + (self.coefficients[4]) * T ** 3 + (self.coefficients[5]) * T ** 4 + (
            self.coefficients[6]) * T ** 5 - (self.coefficients[7])
        return (h-0.09736829*(T- 273.15) + 56.52900916445292) *1000

    def h_from_T(self, T, target_h):
        a = 6.744824499980115
        T_k = T
        h = -(self.coefficients[0]) / T_k + (self.coefficients[1]) * a + self.coefficients[2] * T_k + \
            self.coefficients[3] * T_k ** 2 + self.coefficients[4] * T_k ** 3 + \
            self.coefficients[5] * T_k ** 4 + self.coefficients[6] * T_k ** 5 - \
            self.coefficients[7]
        return (h-0.09736829*(T- 273.15) + 56.52900916445292)*1000 - target_h
    def _h_pT(self, p, T):
        y = T * 1e-3
        return 1e6 * (
            self.coefficients[0]
            + self.coefficients[2] * y
            + self.coefficients[3] / 2 * y ** 2
            - self.coefficients[4] / y
            + self.coefficients[5] / 3 * y ** 3
        ) / self.coefficients[6]

    def T_ph(self, p, h):
        # Define a reasonable initial guess for temperature
        initial_guess = [400]  # Celsius
        # Solve for T that makes the enthalpy equation zero
        T_solution = fsolve(self.h_from_T,initial_guess, args=(h), xtol=1e-8, maxfev=1000)
        return T_solution[0]  # Return the first solution

    def _is_below_T_critical(self, T):
        pass

    def _make_p_subcritical(self, p):
        pass

from tespy.components import HeatExchanger, Sink, Source, Splitter, Drum,Pipe
from tespy.connections import Connection
from tespy.networks import Network
# from tespy.newfluid import KKHWrapper
import pandas as pd
df = pd.read_csv('E:/机理仿真/75工况热力性能实验表.csv',index_col='名称')
df = df[df.columns[1:]].loc[df.index.drop('烟气流向')].astype('float')

nw = Network(fluids=['water', 'O2', 'CO2', 'N2', 'Ar'])
nw.set_attr(p_unit='MPa', T_unit='C', h_unit='kJ / kg', m_unit="t / h", iterinfo=False)

w_source = Source('水源入口')  # 轴封蒸汽冷凝器
w_sink = Sink('水源出口')
w_sink2 = Sink('水源出口2')
w_sink3 = Sink('水源出口3')

g_source = Source('烟气入口')  # 轴封蒸汽冷凝器
g_sink = Sink('烟气出口')
dysmq1 = HeatExchanger('低压省煤器1')
dysmq2 = HeatExchanger('低压省煤器2')
pipe1 = Pipe('管道1')
dysmqflq = Splitter('低压省煤器分离器')
dyqb = Drum('低压汽包')
dyqz = HeatExchanger('低压汽蒸')
flq2 = Splitter('分离器2')

L1 = Connection(w_source, 'out1', dysmq1, 'in2')
L2 = Connection(dysmq1, 'out2', dysmqflq, 'in1')
L2_1 = Connection(dysmqflq, 'out1', dysmq2, 'in2')
L3_1 = Connection(dysmq2, 'out2', pipe1, 'in1')
L3_2 = Connection(pipe1, 'out1', flq2, 'in1')
L4 = Connection(flq2, 'out2', dyqb, 'in1')
L6 = Connection(dyqb, 'out2', w_sink, 'in1')

Z1 = Connection(dyqb, 'out1', dyqz, 'in2')
Z2 = Connection(dyqz, 'out2', dyqb, 'in2')
L5 = Connection(flq2, 'out1', w_sink3, 'in1')
L2_2 = Connection(dysmqflq, 'out2', w_sink2, 'in1')
c29 = Connection(g_source, 'out1', dyqz, 'in1')
c30 = Connection(dyqz, 'out1', dysmq2, 'in1')
c31 = Connection(dysmq2, 'out1', dysmq1, 'in1')
c32 = Connection(dysmq1, 'out1', g_sink, 'in1')

nw.add_conns(L1, L2, L2_1, L2_2,L3_1,L3_2,L4,L5,L6,Z1,Z2,c29,c30, c31, c32)
nw.check_network()
c29.set_attr(T=df['LP蒸发器']['烟气进口温度'])
c30.set_attr(T=df['LP省煤器2']['烟气进口温度'])
c31.set_attr(T=df['LP省煤器2']['烟气出口温度'])
c32.set_attr(m = df['LP省煤器2']['烟气流量'],p = 1,fluid={'water': 4.39 / 100, 'O2': 17.79 / 100, 'CO2': 3.34 / 100, 'N2': 73.17 / 100, 'Ar': 1.31 / 100},fluid_engines={"H2O": KKHWrapper,"O2": KKHWrapper,"CO2": KKHWrapper,"N2": KKHWrapper,"Ar": KKHWrapper})
L1.set_attr(fluid={'water': 1}, T=df['LP省煤器1']['工质进口温度'], p=0.101325 + df['LP省煤器1']['工质进口压力'], m=df['LP省煤器1']['工质流量'])
L2_1.set_attr(m=df['LP省煤器2']['工质流量'], T=df['LP省煤器2']['工质进口温度'], p=0.101325 + df['LP省煤器2']['工质进口压力'])
L3_1.set_attr(T = df['LP省煤器2']['工质出口温度'])
L3_2.set_attr(T = df['LP蒸发器']['工质进口温度'],p=0.101325 + df['LP蒸发器']['工质进口压力'])
L4.set_attr(m = df['LP蒸发器']['工质流量'])
dysmq1.set_attr(pr1=((df['LP省煤器1']['烟气设计压力']/1000000+ 0.101325) - df['LP省煤器1']['烟气压降']/1000000)/(df['LP省煤器1']['烟气设计压力']/1000000+ 0.101325))
dysmq2.set_attr(pr1=((df['LP省煤器2']['烟气设计压力']/1000000+ 0.101325) - df['LP省煤器2']['烟气压降']/1000000)/(df['LP省煤器2']['烟气设计压力']/1000000+ 0.101325))
dyqz.set_attr(pr1=1)
nw.solve("design")
nw.print_results()
# 60,706.21131

75工况热力性能实验表.csv

[fwitte]

Hi @ CZRMWMW,

thank you for reaching out and reporting the bug. I think, this is caused by the starting value generator. It will iterate through all connections and assign starting values according to the components attached to the inlet and to the outlet of the connection (https://github.com/oemof/tespy/blob/c90abe8468720f1ad846326dcaad45f73981c39a/src/tespy/networks/network.py#L1792-L1832). It sounds a little strange, that the order of connections has an effect here, but who knows. On top of that, there is also a little bit of randomness in the starting value assignment for mass flows, which could also be a reason (see https://github.com/oemof/tespy/issues/509). I'll have to investigate it further. I will update you, when I find the time to do so.

Best

Francesco

btw: Even if the simulation runs, the solver does not find a feasible solution. You might want to check your specifications.

[CZRMWMW]

Dear Francesco, thank you very much for your reply. I understand the issue now. If there is anything that I can assist with in the future, I'd be glad to help. I hope you have a pleasant day. Additionally, thank you very much for developing tespy. This tool has helped me solve many research problems.

[fwitte]

You can certainly do so, it is always appreciated. If you want, you can try to find, what causes the differences for this bug. For that, I'd suggest, you'd get a hand on the developer version of the software (see online documentation on how to install it) and then start checking, what the difference is for all SI_values (val_SI of the connection attributes) of the variables of the system variables_dict attribute of the Network class instance, after the initialization is complete and before the actual solution process starts. Have a nice weekend