Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
CZRMWMW opened this issue Jun 7, 2024 · 3 comments
Labels

Comments

@CZRMWMW
Copy link

CZRMWMW commented Jun 7, 2024

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

QQ图片20240607153529

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

@fwitte
Copy link
Member

fwitte commented Jun 8, 2024

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 (

def init_val0(self, c, key):
r"""
Set starting values for fluid properties.
The component classes provide generic starting values for their inlets
and outlets.
Parameters
----------
c : tespy.connections.connection.Connection
Connection to initialise.
"""
if np.isnan(c.get_attr(key).val0):
# starting value for mass flow is random between 1 and 2 kg/s
# (should be generated based on some hash maybe?)
if key == 'm':
c.get_attr(key).val0 = np.random.random() + 1
# generic starting values for pressure and enthalpy
else:
# retrieve starting values from component information
val_s = c.source.initialise_source(c, key)
val_t = c.target.initialise_target(c, key)
if val_s == 0 and val_t == 0:
if key == 'p':
c.get_attr(key).val0 = 1e5
elif key == 'h':
c.get_attr(key).val0 = 1e6
elif val_s == 0:
c.get_attr(key).val0 = val_t
elif val_t == 0:
c.get_attr(key).val0 = val_s
else:
c.get_attr(key).val0 = (val_s + val_t) / 2
# change value according to specified unit system
c.get_attr(key).val0 = hlp.convert_from_SI(
key, c.get_attr(key).val0, self.get_attr(key + '_unit')
)
). 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 #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.

@fwitte fwitte added the bug label Jun 8, 2024
@CZRMWMW
Copy link
Author

CZRMWMW commented Jun 8, 2024

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
Copy link
Member

fwitte commented Jun 8, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants