The order of pipe input in the add_conns function leads to different calculation results.
CZRMWMW opened this issue · 3 comments
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
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 (
tespy/src/tespy/networks/network.py
Lines 1792 to 1832 in c90abe8
Best
Francesco
btw: Even if the simulation runs, the solver does not find a feasible solution. You might want to check your specifications.
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.
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