From 0989493b888e6aafa277521a9f84194272d63ce0 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Mon, 14 May 2018 06:45:34 +0200 Subject: [PATCH] added class for stoichometric combustion chamber, moved some lines in helpers --- tespy/components/components.py | 613 ++++++++++++++------------- tespy/helpers.py | 739 +++++++++++++++++---------------- 2 files changed, 722 insertions(+), 630 deletions(-) diff --git a/tespy/components/components.py b/tespy/components/components.py index 505d51d20..ee8d7f061 100644 --- a/tespy/components/components.py +++ b/tespy/components/components.py @@ -171,6 +171,10 @@ def set_attr(self, **kwargs): self.get_attr(key).set_attr(val=kwargs[key]) self.get_attr(key).set_attr(is_set=True) + elif isinstance(kwargs[key], dict): + self.get_attr(key).set_attr(val=kwargs[key]) + self.get_attr(key).set_attr(is_set=True) + elif kwargs[key] == 'var': self.get_attr(key).set_attr(val=1) self.get_attr(key).set_attr(is_set=True) @@ -2705,6 +2709,25 @@ class combustion_chamber(component): :align: center """ + def inlets(self): + return ['in1', 'in2'] + + def outlets(self): + return ['out1'] + + def attr(self): + return ['fuel', 'lamb', 'ti'] + + def attr_prop(self): + return {'fuel': dc_cp(), 'lamb': dc_cp(), 'ti': dc_cp()} + + def fuels(self): + return ['methane', 'ethane', 'propane', 'butane', + 'hydrogen'] + + def component(self): + return 'combustion chamber' + def comp_init(self, nw): if not self.fuel.is_set: @@ -2744,26 +2767,61 @@ def comp_init(self, nw): else: self.n[el] = 0 - self.hi = self.lhv() + self.lhv = self.calc_lhv() - def inlets(self): - return ['in1', 'in2'] + def calc_lhv(self): + r""" + calculates the lower heating value of the combustion chambers fuel - def outlets(self): - return ['out1'] + :returns: val (*float*) - lhv of the specified fuel - def attr(self): - return ['fuel', 'lamb', 'ti'] + **equation** - def attr_prop(self): - return {'fuel': dc_cp(), 'lamb': dc_cp(), 'ti': dc_cp()} + .. math:: + LHV = -\frac{\sum_i {\Delta H_f^0}_i - + \sum_j {\Delta H_f^0}_j } + {M_{fuel}}\\ + \forall i \in \text{reation products},\\ + \forall j \in \text{reation educts},\\ + \Delta H_f^0: \text{molar formation enthalpy} - def fuels(self): - return ['methane', 'ethane', 'propane', 'butane', - 'hydrogen'] + =============== ===================================== + substance :math:`\frac{\Delta H_f^0}{kJ/mol}` + =============== ===================================== + hydrogen 0 + methane -74.85 + ethane -84.68 + propane -103.8 + butane -124.51 + --------------- ------------------------------------- + oxygen 0 + carbondioxide -393.5 + water (g) -241.8 + =============== ===================================== - def component(self): - return 'combustion chamber' + """ + + hf = {} + hf['hydrogen'] = 0 + hf['methane'] = -74.85 + hf['ethane'] = -84.68 + hf['propane'] = -103.8 + hf['butane'] = -124.51 + hf[self.o2] = 0 + hf[self.co2] = -393.5 + # water (gaseous) + hf[self.h2o] = -241.8 + + key = set(list(hf.keys())).intersection( + set([a.replace(' ', '') + for a in CP.get_aliases(self.fuel.val)])) + + val = (-(self.n['H'] / 2 * hf[self.h2o] + self.n['C'] * hf[self.co2] - + ((self.n['C'] + self.n['H'] / 4) * hf[self.o2] + + hf[list(key)[0]])) / + molar_masses[self.fuel.val] * 1000) + + return val def equations(self, nw): r""" @@ -2895,60 +2953,6 @@ def derivatives(self, nw): return np.asarray(mat_deriv) - def lhv(self): - r""" - calculates the lower heating value of the combustion chambers fuel - - :returns: val (*float*) - lhv of the specified fuel - - **equation** - - .. math:: - LHV = -\frac{\sum_i {\Delta H_f^0}_i - - \sum_j {\Delta H_f^0}_j } - {M_{fuel}}\\ - \forall i \in \text{reation products},\\ - \forall j \in \text{reation educts},\\ - \Delta H_f^0: \text{molar formation enthalpy} - - =============== ===================================== - substance :math:`\frac{\Delta H_f^0}{kJ/mol}` - =============== ===================================== - hydrogen 0 - methane -74.85 - ethane -84.68 - propane -103.8 - butane -124.51 - --------------- ------------------------------------- - oxygen 0 - carbondioxide -393.5 - water (g) -241.8 - =============== ===================================== - - """ - - hf = {} - hf['hydrogen'] = 0 - hf['methane'] = -74.85 - hf['ethane'] = -84.68 - hf['propane'] = -103.8 - hf['butane'] = -124.51 - hf[self.o2] = 0 - hf[self.co2] = -393.5 - # water (gaseous) - hf[self.h2o] = -241.8 - - key = set(list(hf.keys())).intersection( - set([a.replace(' ', '') - for a in CP.get_aliases(self.fuel.val)])) - - val = (-(self.n['H'] / 2 * hf[self.h2o] + self.n['C'] * hf[self.co2] - - ((self.n['C'] + self.n['H'] / 4) * hf[self.o2] + - hf[list(key)[0]])) / - molar_masses[self.fuel.val] * 1000) - - return val - def reaction_balance(self, inl, outl, fluid): r""" calculates the reactions mass balance for one fluid @@ -3087,12 +3091,12 @@ def energy_balance(self, inl, outl): res += i.m.val_SI * (i.h.val_SI - h_mix_pT([i.m.val_SI, p_ref, i.h.val_SI, i.fluid.val], T_ref)) - res += i.m.val_SI * i.fluid.val[self.fuel.val] * self.hi + res += i.m.val_SI * i.fluid.val[self.fuel.val] * self.lhv for o in outl: res -= o.m.val_SI * (o.h.val_SI - h_mix_pT([o.m.val_SI, p_ref, o.h.val_SI, o.fluid.val], T_ref)) - res -= o.m.val_SI * o.fluid.val[self.fuel.val] * self.hi + res -= o.m.val_SI * o.fluid.val[self.fuel.val] * self.lhv return res @@ -3148,7 +3152,7 @@ def ti_func(self, inl, outl): for o in outl: m_fuel -= (o.m.val_SI * o.fluid.val[self.fuel.val]) - return (self.ti.val - m_fuel * self.hi) + return (self.ti.val - m_fuel * self.lhv) def drb_dx(self, inl, outl, dx, pos, fluid): r""" @@ -3348,7 +3352,7 @@ def calc_parameters(self, nw, mode): self.ti.val = 0 for i in inl: - self.ti.val += i.m.val_SI * i.fluid.val[self.fuel.val] * self.hi + self.ti.val += i.m.val_SI * i.fluid.val[self.fuel.val] * self.lhv n_fuel = 0 for i in inl: @@ -3362,13 +3366,13 @@ def calc_parameters(self, nw, mode): if mode == 'post': if not self.lamb.is_set: - self.lamb.val = n_oxygen / (n_fuel * - (self.n['C'] + self.n['H'] / 4)) + self.lamb.val = n_oxygen / ( + n_fuel * (self.n['C'] + self.n['H'] / 4)) if mode == 'pre': if 'lamb' in self.offdesign: - self.lamb.val = n_oxygen / (n_fuel * - (self.n['C'] + self.n['H'] / 4)) + self.lamb.val = n_oxygen / ( + n_fuel * (self.n['C'] + self.n['H'] / 4)) def print_parameters(self, nw): @@ -3386,212 +3390,257 @@ def print_parameters(self, nw): # %% -# maybe to come -#class combustion_chamber_stoich(combustion_chamber): -# """ -# .. note:: -# This component uses air as pseudo fluid and a mixture of air and -# stoichometric flue gas at the outlet. Therefore it implies some -# restrictions: -# -# - the fuel must always be connected to in1 -# - the fuel composition must be fixed (all fluid components have to be -# set) -# -# - the fluids of your network must then be: -# - 'air', -# - your fuel, e. g. 'CH4' and -# - your fuel_fg, e. g. 'CH4_fg'. -# -# In terms of this example, the fluid composition at the outlet of the -# combustion chamber is a mixture of 'air' and 'CH4_fg'. -# -# **available parameters** -# -# - fuel: fuel for combustion chamber -# - lamb: air to stoichiometric air ratio -# -# **available fuels** -# -# - methane -# - ethane -# - propane -# - butane -# - hydrogen -# -# **inlets and outlets** -# -# - in1, in2 -# - out1 -# -# .. image:: _images/combustion_chamber.svg -# :scale: 100 % -# :alt: alternative text -# :align: center -# """ -# -# def comp_init(self, nw): -# -# if not self.fuel_set: -# msg = 'Must specify fuel for combustion chamber.' -# raise MyComponentError(msg) -# -# if (len([x for x in nw.fluids if x in [a.replace(' ', '') for a in -# CP.get_aliases(self.fuel)]]) == 0): -# msg = ('The fuel you specified does not match the fuels available' -# ' within the network.') -# raise MyComponentError(msg) -# -# if (len([x for x in self.fuels() if x in [a.replace(' ', '') for a in -# CP.get_aliases(self.fuel)]])) == 0: -# msg = ('The fuel you specified is not available. Available fuels ' -# 'are: ' + str(self.fuels()) + '.') -# raise MyComponentError(msg) -# -# self.fuel = [x for x in nw.fluids if x in [a.replace(' ', '') for a in -# CP.get_aliases(self.fuel)]][0] -# -# self.o2 = [x for x in nw.fluids if x in -# [a.replace(' ', '') for a in CP.get_aliases('O2')]][0] -# self.co2 = [x for x in nw.fluids if x in -# [a.replace(' ', '') for a in CP.get_aliases('CO2')]][0] -# self.h2o = [x for x in nw.fluids if x in -# [a.replace(' ', '') for a in CP.get_aliases('H2O')]][0] -# self.n2 = [x for x in nw.fluids if x in -# [a.replace(' ', '') for a in CP.get_aliases('N2')]][0] -# -# structure = fluid_structure(self.fuel) -# -# self.n = {} -# for el in ['C', 'H', 'O']: -# if el in structure.keys(): -# self.n[el] = structure[el] -# else: -# self.n[el] = 0 -# -# self.hi = self.lhv() -# -# def inlets(self): -# return ['in1', 'in2'] -# -# def outlets(self): -# return ['out1'] -# -# def attr(self): -# return component.attr(self) + ['fuel', 'lamb'] -# -# def fuels(self): -# return ['methane', 'ethane', 'propane', 'butane', -# 'hydrogen'] -# -# def component(self): -# return 'combustion chamber stoichometric flue gas' -# -# def reaction_balance(self, inl, outl, fluid): -# r""" -# calculates the reactions mass balance for one fluid -# -# - determine molar mass flows of fuel and fresh air -# - calculate excess fuel -# - calculate residual value of the fluids balance -# -# :param inlets: the components connections at the inlets -# :type inlets: list -# :param outlets: the components connections at the outlets -# :type outlets: list -# :param fluid: fluid to calculate the reaction balance for -# :type fluid: str -# :returns: res (*float*) - residual value of mass balance -# -# **reaction balance equations** -# -# .. math:: -# res = \sum_i \left(x_{fluid,i} \cdot \dot{m}_{i}\right) - -# \sum_j \left(x_{fluid,j} \cdot \dot{m}_{j}\right) \; -# \forall i \in inlets, \; \forall j \in outlets -# -# \dot{m}_{fluid,m} = \sum_i \frac{x_{fluid,i} \cdot \dot{m}_{i}} -# {M_{fluid}} \; \forall i \in inlets\\ -# -# \lambda = \frac{\dot{m}_{f,m}}{\dot{m}_{O_2,m} \cdot -# \left(n_{C,fuel} + 0.25 \cdot n_{H,fuel}\right)} -# -# *fuel* -# -# .. math:: -# 0 = res - \left(\dot{m}_{f,m} - \dot{m}_{f,exc,m}\right) -# \cdot M_{fuel}\\ -# -# \dot{m}_{f,exc,m} = \begin{cases} -# 0 & \lambda \geq 1\\ -# \dot{m}_{f,m} - \frac{\dot{m}_{O_2,m}} -# {n_{C,fuel} + 0.25 \cdot n_{H,fuel}} & \lambda < 1 -# \end{cases} -# -# *oxygen* -# -# .. math:: -# 0 = res - \begin{cases} -# -\frac{\dot{m}_{O_2,m} \cdot M_{O_2}}{\lambda} & \lambda \geq 1\\ -# - \dot{m}_{O_2,m} \cdot M_{O_2} & \lambda < 1 -# \end{cases} -# -# *water* -# -# .. math:: -# 0 = res + \left( \dot{m}_{f,m} - \dot{m}_{f,exc,m} \right) -# \cdot 0.5 \cdot n_{H,fuel} \cdot M_{H_2O} -# -# *carbondioxide* -# -# .. math:: -# 0 = res + \left( \dot{m}_{f,m} - \dot{m}_{f,exc,m} \right) -# \cdot n_{C,fuel} \cdot M_{CO_2} -# -# *other* -# -# .. math:: -# 0 = res -# -# """ -# -# n_fuel = 0 -# for i in inl: -# n_fuel += i.m.val_SI * i.fluid.val[self.fuel] / molar_masses[self.fuel] -# -# n_oxygen = 0 -# for i in inl: -# n_oxygen += i.m.val_SI * i.fluid.val[self.o2] / molar_masses[self.o2] -# if not self.lamb_set: -# self.lamb = n_oxygen / (n_fuel * (self.n['C'] + self.n['H'] / 4)) -# -# n_fuel_exc = 0 -# if self.lamb < 1: -# n_fuel_exc = n_fuel - n_oxygen / (self.n['C'] + self.n['H'] / 4) -# -# if fluid == self.co2: -# dm = ((n_fuel - n_fuel_exc) * -# self.n['C'] * molar_masses[self.co2]) -# elif fluid == self.h2o: -# dm = ((n_fuel - n_fuel_exc) * -# self.n['H'] / 2 * molar_masses[self.h2o]) -# elif fluid == self.o2: -# if self.lamb < 1: -# dm = -n_oxygen * molar_masses[self.o2] -# else: -# dm = -n_oxygen / self.lamb * molar_masses[self.o2] -# elif fluid == self.fuel: -# dm = -(n_fuel - n_fuel_exc) * molar_masses[self.fuel] -# else: -# dm = 0 -# -# res = dm -# -# for i in inl: -# res += i.fluid.val[fluid] * i.m.val_SI -# for o in outl: -# res -= o.fluid.val[fluid] * o.m.val_SI -# return res + + +class combustion_chamber_stoich(combustion_chamber): + """r + + .. note:: + This combustion chamber uses fresh air and its fuel as the only + reactive gas components. Therefore note the following restrictions. You + are to + + - specify the oxigen mass fraction in the fresh air, + - fully define the fuel's fluid components, + - provide the aliases of the fresh air and the fuel and + - make sure, both of the aliases are part of the network fluid vector. + + The name of the flue gas will be: 'TESPy::yourfuelalias_fg'. It is also + possible to use fluid mixtures for the fuel, e. g. + :code:`fuel={CH4: 0.9, 'CO2': 0.1}`. If you specify a fluid mixture for + the fuel, TESPy will automatically create a custom fluid called: + 'TESPy::yourfuelalias'. For more information see the examples section. + + **available parameters** + + - fuel: fuel for combustion chamber + - lamb: air to stoichiometric air ratio + - ti: thermal input (:math:`{LHV \cdot \dot{m}_f}`) + + **equations** + + see :func:`tespy.components.components.combustion_chamber_stoich.equations` + + **available fuel gases** + + - methane + - ethane + - propane + - butane + - hydrogen + + **inlets and outlets** + + - in1, in2 + - out1 + + .. image:: _images/combustion_chamber.svg + :scale: 100 % + :alt: alternative text + :align: center + """ + + def inlets(self): + return ['in1', 'in2'] + + def outlets(self): + return ['out1'] + + def attr(self): + return ['fuel', 'fuel_alias', 'air_alias', 'oxy', 'lamb', 'ti'] + + def attr_prop(self): + return {'fuel': dc_cp(), 'fuel_alias': dc_cp(), + 'air_alias': dc_cp(), 'oxy': dc_cp(), + 'lamb': dc_cp(), 'ti': dc_cp()} + + def fuels(self): + return ['methane', 'ethane', 'propane', 'butane', + 'hydrogen'] + + def component(self): + return 'combustion chamber stoichometric flue gas' + + def comp_init(self, nw): + + if not self.fuel.is_set: + msg = 'Must specify fuel for combustion chamber.' + raise MyComponentError(msg) + + if not self.fuel_alias.is_set: + msg = 'Must specify fuel alias for combustion chamber.' + raise MyComponentError(msg) + + if not self.air_alias.is_set: + msg = 'Must specify air alias for combustion chamber.' + raise MyComponentError(msg) + + if not self.oxy.is_set: + msg = ('Must specify air oxygen mass fraction for' + 'combustion chamber.') + raise MyComponentError(msg) + + self.lhv = self.calc_lhv() + + def calc_lhv(self): + r""" + calculates the lower heating value of the combustion chambers fuel + + :returns: val (*float*) - lhv of the specified fuel + + **equation** + + .. math:: + LHV = -\frac{\sum_i {\Delta H_f^0}_i - + \sum_j {\Delta H_f^0}_j } + {M_{fuel}}\\ + \forall i \in \text{reation products},\\ + \forall j \in \text{reation educts},\\ + \Delta H_f^0: \text{molar formation enthalpy} + + =============== ===================================== + substance :math:`\frac{\Delta H_f^0}{kJ/mol}` + =============== ===================================== + hydrogen 0 + methane -74.85 + ethane -84.68 + propane -103.8 + butane -124.51 + --------------- ------------------------------------- + oxygen 0 + carbondioxide -393.5 + water (g) -241.8 + =============== ===================================== + + """ + + hf = {} + hf['hydrogen'] = 0 + hf['methane'] = -74.85 + hf['ethane'] = -84.68 + hf['propane'] = -103.8 + hf['butane'] = -124.51 + hf[self.o2] = 0 + hf[self.co2] = -393.5 + # water (gaseous) + hf[self.h2o] = -241.8 + + key = set(list(hf.keys())).intersection( + set([a.replace(' ', '') + for a in CP.get_aliases(self.fuel.val)])) + + val = (-(self.n['H'] / 2 * hf[self.h2o] + self.n['C'] * hf[self.co2] - + ((self.n['C'] + self.n['H'] / 4) * hf[self.o2] + + hf[list(key)[0]])) / + molar_masses[self.fuel.val] * 1000) + + return val + + def reaction_balance(self, inl, outl, fluid): + r""" + calculates the reactions mass balance for one fluid + + - determine molar mass flows of fuel and fresh air + - calculate excess fuel + - calculate residual value of the fluids balance + + :param inlets: the components connections at the inlets + :type inlets: list + :param outlets: the components connections at the outlets + :type outlets: list + :param fluid: fluid to calculate the reaction balance for + :type fluid: str + :returns: res (*float*) - residual value of mass balance + + **reaction balance equations** + + .. math:: + res = \sum_i \left(x_{fluid,i} \cdot \dot{m}_{i}\right) - + \sum_j \left(x_{fluid,j} \cdot \dot{m}_{j}\right) \; + \forall i \in inlets, \; \forall j \in outlets + + \dot{m}_{fluid,m} = \sum_i \frac{x_{fluid,i} \cdot \dot{m}_{i}} + {M_{fluid}} \; \forall i \in inlets\\ + + \lambda = \frac{\dot{m}_{f,m}}{\dot{m}_{O_2,m} \cdot + \left(n_{C,fuel} + 0.25 \cdot n_{H,fuel}\right)} + + *fuel* + + .. math:: + 0 = res - \left(\dot{m}_{f,m} - \dot{m}_{f,exc,m}\right) + \cdot M_{fuel}\\ + + \dot{m}_{f,exc,m} = \begin{cases} + 0 & \lambda \geq 1\\ + \dot{m}_{f,m} - \frac{\dot{m}_{O_2,m}} + {n_{C,fuel} + 0.25 \cdot n_{H,fuel}} & \lambda < 1 + \end{cases} + + *oxygen* + + .. math:: + 0 = res - \begin{cases} + -\frac{\dot{m}_{O_2,m} \cdot M_{O_2}}{\lambda} & \lambda \geq 1\\ + - \dot{m}_{O_2,m} \cdot M_{O_2} & \lambda < 1 + \end{cases} + + *water* + + .. math:: + 0 = res + \left( \dot{m}_{f,m} - \dot{m}_{f,exc,m} \right) + \cdot 0.5 \cdot n_{H,fuel} \cdot M_{H_2O} + + *carbondioxide* + + .. math:: + 0 = res + \left( \dot{m}_{f,m} - \dot{m}_{f,exc,m} \right) + \cdot n_{C,fuel} \cdot M_{CO_2} + + *other* + + .. math:: + 0 = res + + """ + + n_fuel = 0 + for i in inl: + n_fuel += i.m.val_SI * i.fluid.val[self.fuel] / molar_masses[self.fuel] + + n_oxygen = 0 + for i in inl: + n_oxygen += i.m.val_SI * i.fluid.val[self.o2] / molar_masses[self.o2] + if not self.lamb_set: + self.lamb = n_oxygen / (n_fuel * (self.n['C'] + self.n['H'] / 4)) + + n_fuel_exc = 0 + if self.lamb < 1: + n_fuel_exc = n_fuel - n_oxygen / (self.n['C'] + self.n['H'] / 4) + + if fluid == self.co2: + dm = ((n_fuel - n_fuel_exc) * + self.n['C'] * molar_masses[self.co2]) + elif fluid == self.h2o: + dm = ((n_fuel - n_fuel_exc) * + self.n['H'] / 2 * molar_masses[self.h2o]) + elif fluid == self.o2: + if self.lamb < 1: + dm = -n_oxygen * molar_masses[self.o2] + else: + dm = -n_oxygen / self.lamb * molar_masses[self.o2] + elif fluid == self.fuel: + dm = -(n_fuel - n_fuel_exc) * molar_masses[self.fuel] + else: + dm = 0 + + res = dm + + for i in inl: + res += i.fluid.val[fluid] * i.m.val_SI + for o in outl: + res -= o.fluid.val[fluid] * o.m.val_SI + return res # %% diff --git a/tespy/helpers.py b/tespy/helpers.py index f531bac4b..4c06b45ee 100644 --- a/tespy/helpers.py +++ b/tespy/helpers.py @@ -11,14 +11,10 @@ from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt -import collections - import math import numpy as np import sys -import time from scipy import interpolate -from scipy.optimize import fsolve import warnings warnings.simplefilter("ignore", RuntimeWarning) @@ -31,6 +27,215 @@ gas_constants = {} gas_constants['uni'] = 8.3144598 +# %% + + +class data_container: + """r + + The data container stores data on components and connections attributes. + There are subclasses for the following applications: + + - mass flow, pressure, enthalpy and temperature + - fluid + - component parameters + - component characteristics + + **allowed keywords** in kwargs: + + - see data_container.attr() + """ + + def __init__(self, **kwargs): + + invalid = [] + var = self.attr() + + # default values + for key in var.keys(): + self.__dict__.update({key: var[key]}) + + # specify values + for key in kwargs: + if key not in var.keys(): + invalid += [] + self.__dict__.update({key: kwargs[key]}) + + # print invalid keywords + if len(invalid) > 0: + print('The following keys are not available: ' + str(invalid)) + + def set_attr(self, **kwargs): + + invalid = [] + var = self.attr() + + # specify values + for key in kwargs: + if key not in var.keys(): + invalid += [] + self.__dict__.update({key: kwargs[key]}) + + # print invalid keywords + if len(invalid) > 0: + print('The following keys are not available: ' + str(invalid)) + + def get_attr(self, key): + if key in self.__dict__: + return self.__dict__[key] + else: + print('No attribute \"', key, '\" available!') + return None + + def attr(self): + return [] + + +class dc_prop(data_container): + """r + + data container for fluid properties + + **value specification** + + - val (*numeric*) - user specified value + - val0 (*numeric*) - user specified starting value + - val_SI (*numeric*) - value in SI unit + - val_set (*bool*) - is the specified value a parameter? + + **reference specification** + + - ref (*numeric*) - referenced connection + - ref_set (*bool*) - is the reference a parameter? + + **units** + + - unit (*str*) - unit + - unit_set (*bool*) - is the unit set for the corresponding value? if not, + network unit will be used in calculation (default) + """ + def attr(self): + return {'val': np.nan, 'val0': np.nan, 'val_SI': 0, 'val_set': False, + 'ref': None, 'ref_set': False, + 'unit': None, 'unit_set': False} + + +class dc_flu(data_container): + """r + + data container for fluid vector + + - val (*dict*) - user specified values + - val0 (*dict*) - user specified starting values + - val_set (*dict*) - which components of the fluid vector are set? + - balance (*bool*) - apply fluid balance equation? + """ + def attr(self): + return {'val': {}, 'val0': {}, 'val_set': {}, 'balance': False} + + +class dc_cp(data_container): + """r + + data container for component parameters + + - val (*numeric*) - user specified value + - val_set (*bool*) - is the specified value set? + - is_var (*bool*) - make this parameter a variable of the system? if so, + val will be used as starting value + """ + def attr(self): + return {'val': 0, 'val_SI': 0, 'is_set': False, 'is_var': False} + + +class dc_cc(data_container): + """r + + data container for component characteristics + + - func (*tespy.components.characteristics.characteristics object*) - + characteristic function to be applied + - func_set (*bool*) - is the characteristic function set? + + **using default characteristics** + + see tespy.components.characteristics module for default methods and + parameters, also see tespy.components.components module for available + parameters. + + - method (*str*) - which method of the characteristic function should be + applied? + - param (*str*) - to which parameter should the characteristic function be + applied? + + **using custom characteristics** + + linear interpolation will be applied, it is possible to use default + characteristics and overwrite x-values or y-values + + - x (*np.array*) - array for the x-values of the characteristic line + - y (*np.array*) - array for the y-values of the characteristic line + + """ + def attr(self): + return {'func': None, 'is_set': False, + 'method': 'default', 'param': None, + 'x': None, 'y': None} + +# %% + + +class MyNetworkError(Exception): + pass + + +class MyConnectionError(Exception): + pass + + +class MyComponentError(Exception): + pass + + +class MyConvergenceError(Exception): + pass + + +def query_yes_no(question, default='yes'): + """ + in prompt query + + :param question: question to ask in prompt + :type question: str + :param default: default answer + :type default: str + :returns: bool + """ + valid = {'yes': True, + 'y': True, + 'ye': True, + 'no': False, + 'n': False} + if default is None: + prompt = '[y / n]' + elif default == 'yes': + prompt = '[Y / n]' + elif default == 'no': + prompt = '[y / N]' + + while True: + sys.stdout.write(question + prompt) + choice = input().lower() + if default is not None and choice == '': + return valid[default] + elif choice in valid: + return valid[choice] + else: + sys.stdout.write('Please respond with \'yes\' or \'no\' ' + '(or \'y\' or \'n\').\n') + +# %% + class tespy_fluid: """r @@ -43,20 +248,22 @@ class tespy_fluid: - enthalpy, - entropy, - - density, - - viscoity and + - density and + - viscoity from pressure and temperature. Additionally molar mass and gas constant will be calculated. Inverse functions, e. g. entropy from pressure and enthalpy are calculated via newton algorithm from these tables. - :param alias: name of the fluid mixture will be "TESPY::alias" + :param alias: name of the fluid mixture will be "TESPy::alias" :type alias: str :param fluid: fluid vector for composition {fluid_i: mass fraction, ...} :type fluid: dict :param p_range: range of feasible pressures for newly created fluid + (provide in SI units) :type p_range: list :param T_range: range of feasible temperatures for newly created fluid + (provide in SI units) :type T_range: list :returns: no return value :raises: - :code:`TypeError`, if alias is not of type string @@ -67,7 +274,7 @@ class tespy_fluid: - plot (*bool*), plot the lookup table after creation """ - def __init__(self, alias, fluid, p_range, T_range, nw, **kwargs): + def __init__(self, alias, fluid, p_range, T_range, **kwargs): if not isinstance(alias, str): msg = 'Alias must be of type String.' @@ -78,13 +285,12 @@ def __init__(self, alias, fluid, p_range, T_range, nw, **kwargs): raise ValueError(msg) # process parameters - self.alias = 'TESPY::'+alias + self.alias = 'TESPy::'+alias self.fluid = fluid # adjust value ranges according to specified unit system - self.p_range = np.array(p_range) * nw.p[nw.p_unit] - self.T_range = ((np.array(T_range) + nw.T[nw.T_unit][0]) * - nw.T[nw.T_unit][1]) + self.p_range = np.array(p_range) + self.T_range = np.array(T_range) # set up grid self.p = np.linspace(self.p_range[0], self.p_range[1]) @@ -94,6 +300,10 @@ def __init__(self, alias, fluid, p_range, T_range, nw, **kwargs): self.plot = kwargs.get('plot', False) # calculate molar mass and gas constant + for f in self.fluid: + molar_masses[f] = CPPSI('M', f) + gas_constants[f] = CPPSI('GAS_CONSTANT', f) + molar_masses[self.alias] = 1 / molar_massflow(self.fluid) gas_constants[self.alias] = (gas_constants['uni'] / molar_masses[self.alias]) @@ -110,13 +320,11 @@ def __init__(self, alias, fluid, p_range, T_range, nw, **kwargs): for key in params.keys(): tespy_fluid.fluids[self.alias][key] = ( - self.create_lookup(params[key])) + self.create_lookup(key, params[key])) - p = 3e5 - s = s_mix_pT([0, p, 1e6, {'TESPY::test': 1, 'CH4': 0, 'CO2': 0}], 500) - T_mix_ps([0, p, 0, {'TESPY::test': 1, 'CH4': 0, 'CO2': 0}], s) + print('Successfully created LUTs for custom fluid ' + self.alias) - def create_lookup(self, func): + def create_lookup(self, name, func): """ create lookup table @@ -125,29 +333,30 @@ def create_lookup(self, func): :returns: y (scipy.interpolate.interp2d object) - lookup table """ + print('Creating LUT for ' + name) + x1 = self.p x2 = self.T y = np.empty((0, x1.shape[0]), float) # iterate - for i in self.p: + for p in x1: row = [] - for j in self.T: - row += [func([0, i, 0, self.fluid], j)] + for T in x2: + row += [func([0, p, 0, self.fluid], T)] y = np.append(y, [np.array(row)], axis=0) - self.T, self.p = np.meshgrid(self.T, self.p) - # plot table after creation? if self.plot: fig = plt.figure() ax = fig.add_subplot(111, projection='3d') - ax.plot_wireframe(self.p, self.T, y) - ax.set_xlabel('pressure') - ax.set_ylabel('temperature') - ax.set_zlabel('y') + ax.plot_wireframe(np.meshgrid(x2, x1)[0], + np.meshgrid(x2, x1)[1], y) + ax.set_xlabel('temperature') + ax.set_ylabel('pressure') + ax.set_zlabel(name) ax.view_init(10, 225) plt.show() @@ -184,9 +393,9 @@ def reverse_2d_deriv(params, y): """ d_u = 1 d_l = 1 - if k + d_u > params[0].y_max: + if y + d_u > params[0].y_max: d_u = 0 - if k - d_l < params[0].y_min: + if y - d_l < params[0].y_min: d_l = 0 return ((reverse_2d(params, y + d_u) - reverse_2d(params, y - d_l)) / (d_u + d_l)) @@ -195,6 +404,8 @@ def reverse_2d_deriv(params, y): # initialise the tespy_fluids.fluids container tespy_fluid.fluids = {} +# %% + class memorise: @@ -213,251 +424,50 @@ def __init__(self, fluids): memorise.visc_ph_f[fl] = [] memorise.s_ph[fl] = np.empty((0, num_fl + 3), float) memorise.s_ph_f[fl] = [] - memorise.count = 0 - - def del_memory(fluids): - - fl = tuple(fluids) - - mask = np.isin(memorise.T_ph[fl][:, -1], - memorise.T_ph_f[fl]) - memorise.T_ph[fl] = (memorise.T_ph[fl][mask]) - memorise.T_ph_f[fl] = [] - - mask = np.isin(memorise.T_ps[fl][:, -1], - memorise.T_ps_f[fl]) - memorise.T_ps[fl] = (memorise.T_ps[fl][mask]) - memorise.T_ps_f[fl] = [] - - mask = np.isin(memorise.v_ph[fl][:, -1], - memorise.v_ph_f[fl]) - memorise.v_ph[fl] = (memorise.v_ph[fl][mask]) - memorise.v_ph_f[fl] = [] - - mask = np.isin(memorise.visc_ph[fl][:, -1], - memorise.visc_ph_f[fl]) - memorise.visc_ph[fl] = (memorise.visc_ph[fl][mask]) - memorise.visc_ph_f[fl] = [] - - mask = np.isin(memorise.s_ph[fl][:, -1], - memorise.s_ph_f[fl]) - memorise.s_ph[fl] = (memorise.s_ph[fl][mask]) - memorise.s_ph_f[fl] = [] - - -memorise.T_ph = {} -memorise.T_ph_f = {} -memorise.T_ps = {} -memorise.T_ps_f = {} -memorise.v_ph = {} -memorise.v_ph_f = {} -memorise.visc_ph = {} -memorise.visc_ph_f = {} -memorise.s_ph = {} -memorise.s_ph_f = {} - - -class data_container: - """r - - The data container stores data on components and connections attributes. - There are subclasses for the following applications: - - - mass flow, pressure, enthalpy and temperature - - fluid - - component parameters - - component characteristics - - **allowed keywords** in kwargs: - - - see data_container.attr() - """ - - def __init__(self, **kwargs): - - invalid = [] - var = self.attr() - - # default values - for key in var.keys(): - self.__dict__.update({key: var[key]}) - - # specify values - for key in kwargs: - if key not in var.keys(): - invalid += [] - self.__dict__.update({key: kwargs[key]}) - - # print invalid keywords - if len(invalid) > 0: - print('The following keys are not available: ' + str(invalid)) - - def set_attr(self, **kwargs): - - invalid = [] - var = self.attr() - - # specify values - for key in kwargs: - if key not in var.keys(): - invalid += [] - self.__dict__.update({key: kwargs[key]}) - - # print invalid keywords - if len(invalid) > 0: - print('The following keys are not available: ' + str(invalid)) - - def get_attr(self, key): - if key in self.__dict__: - return self.__dict__[key] - else: - print('No attribute \"', key, '\" available!') - return None - - def attr(self): - return [] - - -class dc_prop(data_container): - """r - - data container for fluid properties - - **value specification** - - - val (*numeric*) - user specified value - - val0 (*numeric*) - user specified starting value - - val_SI (*numeric*) - value in SI unit - - val_set (*bool*) - is the specified value a parameter? - - **reference specification** - - - ref (*numeric*) - referenced connection - - ref_set (*bool*) - is the reference a parameter? - - **units** - - - unit (*str*) - unit - - unit_set (*bool*) - is the unit set for the corresponding value? if not, - network unit will be used in calculation (default) - """ - def attr(self): - return {'val': np.nan, 'val0': np.nan, 'val_SI': 0, 'val_set': False, - 'ref': None, 'ref_set': False, - 'unit': None, 'unit_set': False} - - -class dc_flu(data_container): - """r - - data container for fluid vector - - - val (*dict*) - user specified values - - val0 (*dict*) - user specified starting values - - val_set (*dict*) - which components of the fluid vector are set? - - balance (*bool*) - apply fluid balance equation? - """ - def attr(self): - return {'val': {}, 'val0': {}, 'val_set': {}, 'balance': False} - - -class dc_cp(data_container): - """r - - data container for component parameters - - - val (*numeric*) - user specified value - - val_set (*bool*) - is the specified value set? - - is_var (*bool*) - make this parameter a variable of the system? if so, - val will be used as starting value - """ - def attr(self): - return {'val': 0, 'val_SI': 0, 'is_set': False, 'is_var': False} - - -class dc_cc(data_container): - """r - - data container for component characteristics - - - func (*tespy.components.characteristics.characteristics object*) - - characteristic function to be applied - - func_set (*bool*) - is the characteristic function set? - - **using default characteristics** - - see tespy.components.characteristics module for default methods and - parameters, also see tespy.components.components module for available - parameters. - - - method (*str*) - which method of the characteristic function should be - applied? - - param (*str*) - to which parameter should the characteristic function be - applied? - - **using custom characteristics** - - linear interpolation will be applied, it is possible to use default - characteristics and overwrite x-values or y-values - - - x (*np.array*) - array for the x-values of the characteristic line - - y (*np.array*) - array for the y-values of the characteristic line - - """ - def attr(self): - return {'func': None, 'is_set': False, - 'method': 'default', 'param': None, - 'x': None, 'y': None} - - -class MyNetworkError(Exception): - pass + memorise.count = 0 + def del_memory(fluids): -class MyConnectionError(Exception): - pass + fl = tuple(fluids) + mask = np.isin(memorise.T_ph[fl][:, -1], + memorise.T_ph_f[fl]) + memorise.T_ph[fl] = (memorise.T_ph[fl][mask]) + memorise.T_ph_f[fl] = [] -class MyComponentError(Exception): - pass + mask = np.isin(memorise.T_ps[fl][:, -1], + memorise.T_ps_f[fl]) + memorise.T_ps[fl] = (memorise.T_ps[fl][mask]) + memorise.T_ps_f[fl] = [] + mask = np.isin(memorise.v_ph[fl][:, -1], + memorise.v_ph_f[fl]) + memorise.v_ph[fl] = (memorise.v_ph[fl][mask]) + memorise.v_ph_f[fl] = [] -class MyConvergenceError(Exception): - pass + mask = np.isin(memorise.visc_ph[fl][:, -1], + memorise.visc_ph_f[fl]) + memorise.visc_ph[fl] = (memorise.visc_ph[fl][mask]) + memorise.visc_ph_f[fl] = [] + mask = np.isin(memorise.s_ph[fl][:, -1], + memorise.s_ph_f[fl]) + memorise.s_ph[fl] = (memorise.s_ph[fl][mask]) + memorise.s_ph_f[fl] = [] -def query_yes_no(question, default='yes'): - """ - in prompt query - :param question: question to ask in prompt - :type question: str - :param default: default answer - :type default: str - :returns: bool - """ - valid = {'yes': True, - 'y': True, - 'ye': True, - 'no': False, - 'n': False} - if default is None: - prompt = '[y / n]' - elif default == 'yes': - prompt = '[Y / n]' - elif default == 'no': - prompt = '[y / N]' +memorise.T_ph = {} +memorise.T_ph_f = {} +memorise.T_ps = {} +memorise.T_ps_f = {} +memorise.v_ph = {} +memorise.v_ph_f = {} +memorise.visc_ph = {} +memorise.visc_ph_f = {} +memorise.s_ph = {} +memorise.s_ph_f = {} - while True: - sys.stdout.write(question + prompt) - choice = input().lower() - if default is not None and choice == '': - return valid[default] - elif choice in valid: - return valid[choice] - else: - sys.stdout.write('Please respond with \'yes\' or \'no\' ' - '(or \'y\' or \'n\').\n') +# %% def newton(func, deriv, params, k, **kwargs): @@ -513,6 +523,8 @@ def newton(func, deriv, params, k, **kwargs): return val +# %% + def T_mix_ph(flow): r""" @@ -575,65 +587,6 @@ def T_ph(p, h, fluid): return CPPSI('T', 'P', p, 'H', h, fluid) -def T_mix_ps(flow, s): - r""" - calculates the temperature from pressure and entropy, - uses CoolProp reverse functions for pure fluids, newton for mixtures - - :param flow: vector containing [mass flow, pressure, enthalpy, fluid] - :type flow: list - :param s: entropy in J / (kg * K) - :type s: numeric - :returns: T (float) - temperature in K - - **fluid mixtures** - - .. math:: - - T_{mix}\left(p,s\right) = T_{i}\left(pp_{i},s_{i}\right)\; - \forall i \in \text{fluid components}\\ - - s_{i} = s \left(pp_{i}, T_{mix} \right)\\ - pp: \text{partial pressure} - """ - fl = tuple(sorted(list(flow[3].keys()))) - a = memorise.T_ps[fl][:, 0:-1] - b = np.array([flow[1], flow[2]] + list(flow[3].values()) + [s]) - ix = np.where(np.all(abs(a - b) <= err**2, axis=1))[0] - if ix.size == 1: - T = memorise.T_ps[fl][ix, -1][0] - memorise.T_ps_f[fl] += [T] - return T - else: - if num_fluids(flow[3]) > 1: - val = newton(s_mix_pT, ds_mix_pdT, flow, s, - val0=300, valmin=70, valmax=3000, imax=10) - new = np.array([[flow[1], flow[2]] + list(flow[3].values()) + - [s, val]]) - memorise.T_ps[fl] = np.append(memorise.T_ps[fl], new, axis=0) - return val - else: - for fluid, x in flow[3].items(): - if x > err: - val = T_ps(flow[1], s, fluid) - new = np.array([[flow[1], flow[2]] + - list(flow[3].values()) + [s, val]]) - memorise.T_ps[fl] = np.append(memorise.T_ps[fl], - new, axis=0) - return val - - -def T_ps(p, s, fluid): - if 'IDGAS::' in fluid: - print('Ideal gas calculation not available by now.') - elif 'TESPY::' in fluid: - func = tespy_fluid.fluids[fluid]['s_pT'] - return newton(reverse_2d, reverse_2d_deriv, [func, p, s], 0, - valmin=func.y_min, valmax=func.y_max)[0] - else: - return CPPSI('T', 'P', p, 'S', s, fluid) - - def dT_mix_dph(flow): r""" calculates partial derivate of temperature to pressure at @@ -709,6 +662,69 @@ def dT_mix_ph_dfluid(flow): return np.asarray(vec_deriv) +# %% + + +def T_mix_ps(flow, s): + r""" + calculates the temperature from pressure and entropy, + uses CoolProp reverse functions for pure fluids, newton for mixtures + + :param flow: vector containing [mass flow, pressure, enthalpy, fluid] + :type flow: list + :param s: entropy in J / (kg * K) + :type s: numeric + :returns: T (float) - temperature in K + + **fluid mixtures** + + .. math:: + + T_{mix}\left(p,s\right) = T_{i}\left(pp_{i},s_{i}\right)\; + \forall i \in \text{fluid components}\\ + + s_{i} = s \left(pp_{i}, T_{mix} \right)\\ + pp: \text{partial pressure} + """ + fl = tuple(sorted(list(flow[3].keys()))) + a = memorise.T_ps[fl][:, 0:-1] + b = np.array([flow[1], flow[2]] + list(flow[3].values()) + [s]) + ix = np.where(np.all(abs(a - b) <= err**2, axis=1))[0] + if ix.size == 1: + T = memorise.T_ps[fl][ix, -1][0] + memorise.T_ps_f[fl] += [T] + return T + else: + if num_fluids(flow[3]) > 1: + val = newton(s_mix_pT, ds_mix_pdT, flow, s, + val0=300, valmin=70, valmax=3000, imax=10) + new = np.array([[flow[1], flow[2]] + list(flow[3].values()) + + [s, val]]) + memorise.T_ps[fl] = np.append(memorise.T_ps[fl], new, axis=0) + return val + else: + for fluid, x in flow[3].items(): + if x > err: + val = T_ps(flow[1], s, fluid) + new = np.array([[flow[1], flow[2]] + + list(flow[3].values()) + [s, val]]) + memorise.T_ps[fl] = np.append(memorise.T_ps[fl], + new, axis=0) + return val + + +def T_ps(p, s, fluid): + if 'IDGAS::' in fluid: + print('Ideal gas calculation not available by now.') + elif 'TESPY::' in fluid: + func = tespy_fluid.fluids[fluid]['s_pT'] + return newton(reverse_2d, reverse_2d_deriv, [func, p, s], 0, + valmin=func.y_min, valmax=func.y_max)[0] + else: + return CPPSI('T', 'P', p, 'S', s, fluid) + +# %% + def h_mix_pT(flow, T): r""" @@ -747,16 +763,25 @@ def h_pT(p, T, fluid): return CPPSI('H', 'P', p, 'T', T, fluid) -def h_ps(p, s, fluid): - if 'IDGAS::' in fluid: - print('Ideal gas calculation not available by now.') - elif 'TESPY::' in fluid: - func = tespy_fluid.fluids[fluid]['s_pT'] - T = newton(reverse_2d, reverse_2d_deriv, [func, p, s], 0, - valmin=func.y_min, valmax=func.y_max)[0] - return tespy_fluid.fluids[fluid]['h_pT'](p, T) - else: - return CPPSI('H', 'P', p, 'S', s, fluid) +def dh_mix_pdT(flow, T): + r""" + calculates partial derivate of enthalpy to temperature at constant pressure + + :param flow: vector containing [mass flow, pressure, enthalpy, fluid] + :type flow: list + :param T: temperature in K + :type T: numeric + :returns: dh / dT (float) - derivative in J / (kg * K) + + .. math:: + + \frac{\partial h_{mix}}{\partial T} = + \frac{h_{mix}(p,T+d)-h_{mix}(p,T-d)}{2 \cdot d} + """ + d = 2 + return (h_mix_pT(flow, T + d) - h_mix_pT(flow, T - d)) / (2 * d) + +# %% def h_mix_ps(flow, s): @@ -779,23 +804,18 @@ def h_mix_ps(flow, s): return h_mix_pT(flow, T_mix_ps(flow, s)) -def dh_mix_pdT(flow, T): - r""" - calculates partial derivate of enthalpy to temperature at constant pressure - - :param flow: vector containing [mass flow, pressure, enthalpy, fluid] - :type flow: list - :param T: temperature in K - :type T: numeric - :returns: dh / dT (float) - derivative in J / (kg * K) - - .. math:: +def h_ps(p, s, fluid): + if 'IDGAS::' in fluid: + print('Ideal gas calculation not available by now.') + elif 'TESPY::' in fluid: + func = tespy_fluid.fluids[fluid]['s_pT'] + T = newton(reverse_2d, reverse_2d_deriv, [func, p, s], 0, + valmin=func.y_min, valmax=func.y_max)[0] + return tespy_fluid.fluids[fluid]['h_pT'](p, T) + else: + return CPPSI('H', 'P', p, 'S', s, fluid) - \frac{\partial h_{mix}}{\partial T} = - \frac{h_{mix}(p,T+d)-h_{mix}(p,T-d)}{2 \cdot d} - """ - d = 2 - return (h_mix_pT(flow, T + d) - h_mix_pT(flow, T - d)) / (2 * d) +# %% def h_mix_pQ(flow, Q): @@ -854,6 +874,8 @@ def dh_mix_dpQ(flow, Q): l[1] -= d return (h_mix_pQ(u, Q) - h_mix_pQ(l, Q)) / (2 * d) +# %% + def v_mix_ph(flow): r""" @@ -907,6 +929,8 @@ def d_ph(p, h, fluid): else: return CPPSI('D', 'P', p, 'H', h, fluid) +# %% + def v_mix_pT(flow, T): r""" @@ -960,6 +984,8 @@ def d_pT(p, T, fluid): else: return CPPSI('D', 'P', p, 'T', T, fluid) +# %% + def visc_mix_ph(flow): r""" @@ -1013,6 +1039,8 @@ def visc_ph(p, h, fluid): else: return CPPSI('V', 'P', p, 'H', h, fluid) +# %% + def visc_mix_pT(flow, T): r""" @@ -1053,6 +1081,8 @@ def visc_pT(p, T, fluid): else: return CPPSI('V', 'P', p, 'T', T, fluid) +# %% + def s_mix_ph(flow): r""" @@ -1106,6 +1136,8 @@ def s_ph(p, h, fluid): else: return CPPSI('S', 'P', p, 'H', h, fluid) +# %% + def s_mix_pT(flow, T): r""" @@ -1166,6 +1198,8 @@ def ds_mix_pdT(flow, T): d = 2 return (s_mix_pT(flow, T + d) - s_mix_pT(flow, T - d)) / (2 * d) +# %% + def molar_massflow(flow): r""" @@ -1188,6 +1222,8 @@ def molar_massflow(flow): return mm +# %% + def num_fluids(fluids): r""" @@ -1212,6 +1248,8 @@ def num_fluids(fluids): return n +# %% + def single_fluid(fluids): r""" @@ -1229,6 +1267,8 @@ def single_fluid(fluids): else: return [] +# %% + def fluid_structure(fluid): """ @@ -1246,6 +1286,9 @@ def fluid_structure(fluid): return parts +# %% + + def lamb(re, ks, d): r""" calculates darcy friction factor