diff --git a/CHANGELOG.md b/CHANGELOG.md index cafb362..d587593 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,15 @@ - v0.1.0: - First implementation based on prior work - v0.1.1: - - Move plot_cycle function \ No newline at end of file + - Move plot_cycle function +- v0.2.0: + - Add LMTD MovingBoundary heat exchangers + - Add VDI atlas methods for plate heat exchanger + - Add scaling factor to TenCoefficientCompressor + - Add extrapolation option to TenCoefficientCompressor + - Improve warnings and parameters in TenCoefficientCompressor (eta_mech, T_sc) + - Add T_con_out as possible Input + - Automation is now with multiprocessing + - SDF saving is optional in automation + - Improve nominal design util + - Improve C10 regression generation, add .csv support for TenCoefficientCompressor diff --git a/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb b/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb index ea8dfc8..da02875 100644 --- a/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb +++ b/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb @@ -44,14 +44,14 @@ { "cell_type": "markdown", "metadata": {}, - "source": "As in the other example, we can specify save-paths,\nsolver settings and inputs to vary:\n" + "source": "As in the other example, we can specify save-paths,\nsolver settings and inputs to vary:\nNote that T_con can either be inlet or outlet, depending on the setting\nof `use_condenser_inlet`. Per default, we simulate the inlet, T_con_in\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": "save_path = r\"D:\\00_temp\\simple_heat_pump\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_in_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" + "source": "save_path = r\"D:\\00_temp\\simple_heat_pump\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" }, { "cell_type": "markdown", @@ -63,7 +63,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nsave_path_sdf, save_path_csv = utils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_in_ar=T_con_in_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" + "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nsave_path_sdf, save_path_csv = utils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_ar=T_con_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_condenser_inlet=use_condenser_inlet,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" }, { "cell_type": "markdown", @@ -99,19 +99,19 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "import matplotlib.pyplot as plt\nx_name = \"n in - (Relative compressor speed)\"\ny_name = \"COP in - (Coefficient of performance)\"\nplt.scatter(df[x_name], df[y_name])\nplt.ylabel(y_name)\nplt.xlabel(x_name)\nplt.show()\n" + "source": "import matplotlib.pyplot as plt\nx_name = \"n in - (Relative compressor speed)\"\ny_name = \"COP in - (Coefficient of performance)\"\nplt.scatter(df[x_name], df[y_name], s=20)\nplt.ylabel(y_name)\nplt.xlabel(x_name)\nplt.show()\n" }, { "cell_type": "markdown", "metadata": {}, - "source": "Looking at the results, we see that a higher frequency often leads to lower COP values.\nHowever, other inputs (temperatures) have a greater impact on the COP.\nWe can also use existing 3D-plotting scripts in vclibpy to analyze the\ndependencies. For this, we need the .sdf file. In the sdf, the field names are without\nthe unit and description, as those are accessible in the file-format in other columns.\n" + "source": "Looking at the results, we see that a higher frequency often leads to lower COP values.\nHowever, other inputs (temperatures) have a greater impact on the COP.\nWe can also use existing 3D-plotting scripts in vclibpy to analyze the\ndependencies. For this, we need the .sdf file. In the sdf, the field names are without\nthe unit and description, as those are accessible in the file-format in other columns.\nDepending on whether we varied the inlet or outlet, we have to specify the correct name.\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": "from vclibpy.utils.plotting import plot_sdf_map\nplot_sdf_map(\n filepath_sdf=save_path_sdf,\n nd_data=\"COP\",\n first_dimension=\"T_eva_in\",\n second_dimension=\"T_con_in\",\n fluids=[\"Propane\"],\n flowsheets=[\"Standard\"]\n)\n" + "source": "from vclibpy.utils.plotting import plot_sdf_map\nplot_sdf_map(\n filepath_sdf=save_path_sdf,\n nd_data=\"COP\",\n first_dimension=\"T_eva_in\",\n second_dimension=\"T_con_in\" if use_condenser_inlet else \"T_con_out\",\n fluids=[\"Propane\"],\n flowsheets=[\"Standard\"]\n)\n" } ], "metadata": { diff --git a/docs/jupyter_notebooks/e7_vapor_injection.ipynb b/docs/jupyter_notebooks/e7_vapor_injection.ipynb index f2599e1..2519f73 100644 --- a/docs/jupyter_notebooks/e7_vapor_injection.ipynb +++ b/docs/jupyter_notebooks/e7_vapor_injection.ipynb @@ -63,7 +63,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "save_path = r\"D:\\00_temp\\vapor_injection\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_in_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" + "source": "save_path = r\"D:\\00_temp\\vapor_injection\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" }, { "cell_type": "markdown", @@ -75,7 +75,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_in_ar=T_con_in_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" + "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_ar=T_con_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_condenser_inlet=True,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" }, { "cell_type": "markdown", @@ -123,7 +123,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "heat_pump = VaporInjectionEconomizer(\n evaporator=evaporator,\n condenser=condenser,\n fluid=\"Propane\",\n economizer=economizer,\n high_pressure_compressor=high_pressure_compressor,\n low_pressure_compressor=low_pressure_compressor,\n high_pressure_valve=high_pressure_valve,\n low_pressure_valve=low_pressure_valve\n)\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=r\"D:\\00_temp\\vapor_injection_economizer\",\n T_con_in_ar=T_con_in_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" + "source": "heat_pump = VaporInjectionEconomizer(\n evaporator=evaporator,\n condenser=condenser,\n fluid=\"Propane\",\n economizer=economizer,\n high_pressure_compressor=high_pressure_compressor,\n low_pressure_compressor=low_pressure_compressor,\n high_pressure_valve=high_pressure_valve,\n low_pressure_valve=low_pressure_valve\n)\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=r\"D:\\00_temp\\vapor_injection_economizer\",\n T_con_ar=T_con_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" } ], "metadata": { diff --git a/docs/source/examples/e6_simple_heat_pump.md b/docs/source/examples/e6_simple_heat_pump.md index 7b5c84e..988de47 100644 --- a/docs/source/examples/e6_simple_heat_pump.md +++ b/docs/source/examples/e6_simple_heat_pump.md @@ -71,11 +71,13 @@ heat_pump = StandardCycle( As in the other example, we can specify save-paths, solver settings and inputs to vary: +Note that T_con can either be inlet or outlet, depending on the setting +of `use_condenser_inlet`. Per default, we simulate the inlet, T_con_in ```python save_path = r"D:\00_temp\simple_heat_pump" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] -T_con_in_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] +T_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] n_ar = [0.3, 0.7, 1] ``` @@ -92,9 +94,10 @@ from vclibpy import utils save_path_sdf, save_path_csv = utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=use_condenser_inlet, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -127,7 +130,7 @@ You have to know the names, which are the column headers. import matplotlib.pyplot as plt x_name = "n in - (Relative compressor speed)" y_name = "COP in - (Coefficient of performance)" -plt.scatter(df[x_name], df[y_name]) +plt.scatter(df[x_name], df[y_name], s=20) plt.ylabel(y_name) plt.xlabel(x_name) plt.show() @@ -138,6 +141,7 @@ However, other inputs (temperatures) have a greater impact on the COP. We can also use existing 3D-plotting scripts in vclibpy to analyze the dependencies. For this, we need the .sdf file. In the sdf, the field names are without the unit and description, as those are accessible in the file-format in other columns. +Depending on whether we varied the inlet or outlet, we have to specify the correct name. ```python from vclibpy.utils.plotting import plot_sdf_map @@ -145,7 +149,7 @@ plot_sdf_map( filepath_sdf=save_path_sdf, nd_data="COP", first_dimension="T_eva_in", - second_dimension="T_con_in", + second_dimension="T_con_in" if use_condenser_inlet else "T_con_out", fluids=["Propane"], flowsheets=["Standard"] ) diff --git a/docs/source/examples/e7_vapor_injection.md b/docs/source/examples/e7_vapor_injection.md index 05e47cf..26d96c1 100644 --- a/docs/source/examples/e7_vapor_injection.md +++ b/docs/source/examples/e7_vapor_injection.md @@ -85,7 +85,7 @@ solver settings and inputs to vary: ```python save_path = r"D:\00_temp\vapor_injection" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] -T_con_in_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] +T_con_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] n_ar = [0.3, 0.7, 1] ``` @@ -102,9 +102,10 @@ from vclibpy import utils utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=True, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -161,7 +162,7 @@ heat_pump = VaporInjectionEconomizer( utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=r"D:\00_temp\vapor_injection_economizer", - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, use_multiprocessing=False, diff --git a/examples/e6_simple_heat_pump.py b/examples/e6_simple_heat_pump.py index e258898..48854a6 100644 --- a/examples/e6_simple_heat_pump.py +++ b/examples/e6_simple_heat_pump.py @@ -1,6 +1,6 @@ # # Example for a heat pump with a standard cycle -def main(): +def main(use_condenser_inlet: bool = True): # Let's start the complete cycle simulation with the # most basic flowsheet, the standard-cycle. As all flowsheets # contain a condenser and an evaporator, we defined a common BaseCycle @@ -61,9 +61,11 @@ def main(): ) # As in the other example, we can specify save-paths, # solver settings and inputs to vary: + # Note that T_con can either be inlet or outlet, depending on the setting + # of `use_condenser_inlet`. Per default, we simulate the inlet, T_con_in save_path = r"D:\00_temp\simple_heat_pump" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] - T_con_in_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] + T_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] n_ar = [0.3, 0.7, 1] # Now, we can generate the full-factorial performance map @@ -77,9 +79,10 @@ def main(): save_path_sdf, save_path_csv = utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=use_condenser_inlet, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -101,7 +104,7 @@ def main(): import matplotlib.pyplot as plt x_name = "n in - (Relative compressor speed)" y_name = "COP in - (Coefficient of performance)" - plt.scatter(df[x_name], df[y_name]) + plt.scatter(df[x_name], df[y_name], s=20) plt.ylabel(y_name) plt.xlabel(x_name) plt.show() @@ -110,16 +113,17 @@ def main(): # We can also use existing 3D-plotting scripts in vclibpy to analyze the # dependencies. For this, we need the .sdf file. In the sdf, the field names are without # the unit and description, as those are accessible in the file-format in other columns. + # Depending on whether we varied the inlet or outlet, we have to specify the correct name. from vclibpy.utils.plotting import plot_sdf_map plot_sdf_map( filepath_sdf=save_path_sdf, nd_data="COP", first_dimension="T_eva_in", - second_dimension="T_con_in", + second_dimension="T_con_in" if use_condenser_inlet else "T_con_out", fluids=["Propane"], flowsheets=["Standard"] ) if __name__ == "__main__": - main() + main(use_condenser_inlet=False) diff --git a/examples/e7_vapor_injection.py b/examples/e7_vapor_injection.py index 315b801..ca2a811 100644 --- a/examples/e7_vapor_injection.py +++ b/examples/e7_vapor_injection.py @@ -69,7 +69,7 @@ def main(): # solver settings and inputs to vary: save_path = r"D:\00_temp\vapor_injection" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] - T_con_in_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] + T_con_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] n_ar = [0.3, 0.7, 1] # Now, we can generate the full-factorial performance map @@ -83,9 +83,10 @@ def main(): utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=True, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -127,7 +128,7 @@ def main(): utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=r"D:\00_temp\vapor_injection_economizer", - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, use_multiprocessing=False, diff --git a/requirements.txt b/requirements.txt index 1692d98..5d3fc55 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy>=1.20.0 +numpy<=2.0 matplotlib>=3.3.4 coolprop>=6.4.1 sdf>=0.3.5 diff --git a/tests/test_examples.py b/tests/test_examples.py index cb8953d..1dde3ad 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -13,12 +13,12 @@ class TestExamples(unittest.TestCase): def setUp(self) -> None: self.timeout = 100 # Seconds which the script is allowed to run - def _run_example(self, example, timeout=None): + def _run_example(self, example, timeout=None, **kwargs): ex_py = pathlib.Path(__file__).absolute().parents[1].joinpath("examples", example) sys.path.insert(0, str(ex_py.parent)) module = importlib.import_module(ex_py.stem) example_main = getattr(module, "main") - example_main() + example_main(**kwargs) def test_e1_refrigerant_data(self): self._run_example(example="e1_refrigerant_data.py") @@ -38,6 +38,9 @@ def test_e5_expansion_valve(self): def test_e6_simple_heat_pump(self): self._run_example(example="e6_simple_heat_pump.py") + def test_e6_simple_heat_pump_outlet(self): + self._run_example(example="e6_simple_heat_pump.py", use_condenser_inlet=False) + def test_e7_vapor_injection(self): self._run_example(example="e7_vapor_injection.py") diff --git a/tests/test_regression.py b/tests/test_regression.py index 90196c5..9ab0c1f 100644 --- a/tests/test_regression.py +++ b/tests/test_regression.py @@ -131,7 +131,7 @@ def _regression_of_examples(self, flowsheet, fluid): # Just for quick study: Specify concrete points: T_eva_in_ar = [-10 + 273.15, 273.15] - T_con_in_ar = [30 + 273.15, 70 + 273.15] + T_con_ar = [30 + 273.15, 70 + 273.15] n_ar = [0.3, 1] os.makedirs(self.working_dir, exist_ok=True) @@ -146,7 +146,7 @@ def _regression_of_examples(self, flowsheet, fluid): _, path_csv = utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=self.working_dir, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, use_multiprocessing=False, @@ -155,6 +155,7 @@ def _regression_of_examples(self, flowsheet, fluid): m_flow_eva=0.9, dT_eva_superheating=5, dT_con_subcooling=0, + **kwargs ) path_csv_regression = pathlib.Path(__file__).parent.joinpath( "regression_data", "reference_results", f"{flowsheet}_{fluid}.csv" @@ -171,7 +172,7 @@ def test_vi_ps_propane(self): self._regression_of_examples("VaporInjectionPhaseSeparator", "Propane") def test_evi_propane(self): - self.skipTest("EVI works locally, only CI fails.") + #self.skipTest("EVI works locally, only CI fails.") self._regression_of_examples("VaporInjectionEconomizer", "Propane") @unittest.skip("not implemented") diff --git a/vclibpy/__init__.py b/vclibpy/__init__.py index ed9a16a..1be17e7 100644 --- a/vclibpy/__init__.py +++ b/vclibpy/__init__.py @@ -4,4 +4,4 @@ from vclibpy.datamodels import FlowsheetState, Inputs -__version__ = "0.1.2" +__version__ = "0.2.0" diff --git a/vclibpy/components/compressors/compressor.py b/vclibpy/components/compressors/compressor.py index 35d2a74..e4bb261 100644 --- a/vclibpy/components/compressors/compressor.py +++ b/vclibpy/components/compressors/compressor.py @@ -127,7 +127,7 @@ def calc_state_outlet(self, p_outlet: float, inputs: Inputs, fs_state: Flowsheet self.state_inlet.h + (state_outlet_isentropic.h - self.state_inlet.h) / eta_is ) - fs_state.set(name="eta_is", value=eta_is, unit="%", description="Isentropic efficiency") + fs_state.set(name="eta_is", value=eta_is, unit="-", description="Isentropic efficiency") self.state_outlet = self.med_prop.calc_state("PH", p_outlet, h_outlet) def calc_m_flow(self, inputs: Inputs, fs_state: FlowsheetState) -> float: @@ -148,7 +148,7 @@ def calc_m_flow(self, inputs: Inputs, fs_state: FlowsheetState) -> float: self.get_n_absolute(inputs.n) ) self.m_flow = self.state_inlet.d * V_flow_ref - fs_state.set(name="lambda_h", value=lambda_h, unit="%", description="Volumetric efficiency") + fs_state.set(name="lambda_h", value=lambda_h, unit="-", description="Volumetric efficiency") fs_state.set(name="V_flow_ref", value=V_flow_ref, unit="m3/s", description="Refrigerant volume flow rate") fs_state.set(name="m_flow_ref", value=self.m_flow, unit="kg/s", description="Refrigerant mass flow rate") return self.m_flow diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index cd09349..f6a87c5 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -1,11 +1,14 @@ -import warnings +import logging from abc import ABC +from typing import Union + import numpy as np import pandas as pd - from vclibpy.components.compressors.compressor import Compressor from vclibpy.datamodels import Inputs +logger = logging.getLogger(__name__) + def calc_ten_coefficients(T_eva, T_con, coef_list): """ @@ -55,40 +58,40 @@ class BaseTenCoefficientCompressor(Compressor, ABC): "m_flow": "Flow Rate(kg/h)", "capacity": "Capacity(W)", "input_power": "Input Power(W)", - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumentric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. + extrapolate (str, optional): + Method to handle extrapolation of data. + Default "hold" means no extrapolation """ def __init__(self, N_max, V_h, datasheet, **kwargs): """ Initialize the BaseTenCoefficientCompressor. - - Args: - N_max (float): Maximal rotations per second of the compressor. - V_h (float): Volume of the compressor in m^3. - datasheet (str): Path of the datasheet file. - parameter_names (dict, optional): Dictionary of parameter names. Defaults to None. - sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. """ super().__init__(N_max, V_h) sheet_name = kwargs.get('sheet_name', None) - self.md = pd.read_excel(datasheet, sheet_name=sheet_name) + if str(datasheet).endswith(".xlsx"): + self.md = pd.read_excel(datasheet, sheet_name=sheet_name) + else: + self.md = pd.read_csv(datasheet) parameter_names = kwargs.get('parameter_names', None) if parameter_names is None: self.parameter_names = { "m_flow": "Flow Rate(kg/h)", "capacity": "Capacity(W)", "input_power": "Input Power(W)", - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumentric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } else: self.parameter_names = parameter_names + self.extrapolate = kwargs.get("extrapolate", "hold") def get_parameter(self, T_eva, T_con, n, type_): """ @@ -117,7 +120,15 @@ def get_parameter(self, T_eva, T_con, n, type_): n_list.append(coefficients.pop(0)) param_list.append(calc_ten_coefficients(T_eva, T_con, coefficients)) - return np.interp(self.get_n_absolute(n), n_list, param_list) # linear interpolation + return self._interpolate(self.get_n_absolute(n), n_list, param_list) # linear interpolation + + def _interpolate(self, x_new, x, y): + if self.extrapolate == "hold": + # linear interpolation, no extrapolation + return np.interp(x_new, x, y) + if self.extrapolate == "linear": + return linear_interpolate_extrapolate(x_new, x, y) + raise KeyError(f"Given extrapolate option '{self.extrapolate}' is not supported!") class TenCoefficientCompressor(BaseTenCoefficientCompressor): @@ -151,7 +162,8 @@ class TenCoefficientCompressor(BaseTenCoefficientCompressor): T_sc (float): Subcooling according to datasheet in K. T_sh (float): Superheating according to datasheet in K. capacity_definition (str): Definition of "capacity" in the datasheet. "cooling" or "heating". - assumed_eta_mech (float): Assumed mechanical efficiency of the compressor (only needed if cooling). + assumed_eta_mech (float, callable): Assumed mechanical efficiency of the compressor (only needed if cooling). + If you pass a funtion, it must have this signature `eta_mech(self, p_outlet, inputs)` datasheet (str): Path of the modified datasheet. **kwargs: parameter_names (dict, optional): @@ -165,15 +177,32 @@ class TenCoefficientCompressor(BaseTenCoefficientCompressor): sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. """ - def __init__(self, N_max, V_h, T_sc, T_sh, capacity_definition, assumed_eta_mech, datasheet, **kwargs): + def __init__( + self, + N_max, + V_h, + T_sh, + capacity_definition, + datasheet, + T_sc: float = None, + assumed_eta_mech: Union[float, callable] = None, + scaling_factor: float = 1, + **kwargs + ): super().__init__(N_max=N_max, V_h=V_h, datasheet=datasheet, **kwargs) + if capacity_definition == "cooling" and assumed_eta_mech is None: + raise ValueError("capacity_definition cooling requires an assumption for eta_mech") + if capacity_definition == "heating" and T_sc is None: + raise ValueError("capacity_definition heating requires an assumption for T_sc") self.T_sc = T_sc self.T_sh = T_sh if capacity_definition not in ["cooling", "heating"]: raise ValueError("capacity_definition has to be either 'heating' or 'cooling'") self._capacity_definition = capacity_definition + # Don't use lambda function to cast float as a function, + # as local functions are not pickable for multiprocessing self.assumed_eta_mech = assumed_eta_mech - self.datasheet = datasheet + self.scaling_factor = scaling_factor def get_lambda_h(self, inputs: Inputs): """ @@ -192,18 +221,23 @@ def get_lambda_h(self, inputs: Inputs): T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - 273.15 # [°C] if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): - warnings.warn("The superheating of the given state is not " - "equal to the superheating of the datasheet. " - "State1.T_sh= " + str(round((self.state_inlet.T - T_eva - 273.15), 2)) + - ". Datasheet.T_sh = " + str(self.T_sh)) + logger.debug("The superheating of the given state is not " + "equal to the superheating of the datasheet. " + "State: T_sh=%s; Datasheet: T_sh=%s", + round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) + # The datasheet has a given superheating temperature which can # vary from the superheating of the real state 1 # which is given by the user. # Thus a new self.state_inlet_datasheet has to # be defined for all further calculations - state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 # [kg/s] + if self.T_sh != 0: + state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + else: + state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) + + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] lambda_h = m_flow / (n_abs * state_inlet_datasheet.d * self.V_h) return lambda_h @@ -223,19 +257,31 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): p_2=p_outlet, inputs=inputs ) - h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] h2s = self.med_prop.calc_state("PS", p_outlet, state_inlet_datasheet.s).h # [J/kg] + if callable(self.assumed_eta_mech): + eta_mech = self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) + else: + eta_mech = self.assumed_eta_mech + if self._capacity_definition == "heating": + if self.T_sc != 0: + h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] + else: + h3 = self.med_prop.calc_state("PQ", p_outlet, 0).h # [J/kg] h2 = h3 + capacity / m_flow # [J/kg] else: - h2 = h3 + (capacity + p_el * self.assumed_eta_mech) / m_flow # [J/kg] - - if h2s > h2: - raise ValueError("The calculated eta_s is above 1. You probably chose the wrong capacity_definition") - - eta_s = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) - return eta_s + h2 = state_inlet_datasheet.h + (p_el * eta_mech) / m_flow # [J/kg] + + eta_is = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) + if eta_is > 0.8: + logger.warning( + f"Calculated eta_is is {eta_is * 100} %, which is higher than " + f"typical maximal values of up to, e.g., 80 %. " + "You either chose the wrong capacity_definition, " + "or your assumed eta_mech is not realistic.", + ) + return eta_is def get_eta_mech(self, inputs: Inputs): """ @@ -250,13 +296,19 @@ def get_eta_mech(self, inputs: Inputs): p_outlet = self.get_p_outlet() if self._capacity_definition == "cooling": + if callable(self.assumed_eta_mech): + return self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) return self.assumed_eta_mech + # Else heating T_con, state_inlet_datasheet, m_flow, capacity, p_el = self._calculate_values( p_2=p_outlet, inputs=inputs ) + if self.T_sc != 0: + h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] + else: + h3 = self.med_prop.calc_state("PQ", p_outlet, 0).h # [J/kg] - h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] h2 = h3 + capacity / m_flow # [J/kg] eta_mech = (m_flow * (h2 - state_inlet_datasheet.h)) / p_el @@ -276,11 +328,14 @@ def _calculate_values(self, p_2: float, inputs: Inputs): T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] T_con = self.med_prop.calc_state("PQ", p_2, 0).T - 273.15 # [°C] - state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + if self.T_sh != 0: + state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + else: + state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 # [kg/s] - capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") # [W] - p_el = self.get_parameter(T_eva, T_con, inputs.n, "input_power") # [W] + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] + capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") * self.scaling_factor # [W] + p_el = self.get_parameter(T_eva, T_con, inputs.n, "input_power") * self.scaling_factor # [W] return T_con, state_inlet_datasheet, m_flow, capacity, p_el @@ -311,7 +366,7 @@ class DataSheetCompressor(BaseTenCoefficientCompressor): Dictionary to match internal parameter names (keys) to the names used in the table values. Default { - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumetric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } @@ -349,7 +404,7 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): """ T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 0).T T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - return self.get_parameter(T_eva, T_con, inputs.n, "eta_s") + return self.get_parameter(T_eva, T_con, inputs.n, "eta_is") def get_eta_mech(self, inputs: Inputs): """ @@ -365,3 +420,28 @@ def get_eta_mech(self, inputs: Inputs): T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 0).T T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T return self.get_parameter(T_eva, T_con, inputs.n, "eta_mech") + + +def linear_interpolate_extrapolate(x_new, x, y): + """ + Util function for linear 1D extrapolation or interpolation. + Used to avoid scipy requirement. + + x_new: points where to interpolate/extrapolate + x: known x values + y: known y values + """ + y_new = np.interp(x_new, x, y) + + # Handle left extrapolation + left_mask = x_new < x[0] + if x_new < x[0]: + slope = (y[1] - y[0]) / (x[1] - x[0]) + y_new = y[0] + slope * (x_new - x[0]) + + # Handle right extrapolation + if x_new > x[-1]: + slope = (y[-1] - y[-2]) / (x[-1] - x[-2]) + y_new = y[-1] + slope * (x_new - x[-1]) + + return y_new diff --git a/vclibpy/components/heat_exchangers/heat_exchanger.py b/vclibpy/components/heat_exchangers/heat_exchanger.py index 610ec8d..e56226a 100644 --- a/vclibpy/components/heat_exchangers/heat_exchanger.py +++ b/vclibpy/components/heat_exchangers/heat_exchanger.py @@ -37,11 +37,14 @@ def __init__( gas_heat_transfer: HeatTransfer, liquid_heat_transfer: HeatTransfer, two_phase_heat_transfer: TwoPhaseHeatTransfer, - secondary_medium: str + secondary_medium: str, + ratio_outer_to_inner_area: float = 1, + flow_type: str = "counter" ): super().__init__() self.A = A self.secondary_medium = secondary_medium.lower() + self.ratio_outer_to_inner_area = ratio_outer_to_inner_area self._wall_heat_transfer = wall_heat_transfer self._secondary_heat_transfer = secondary_heat_transfer @@ -61,7 +64,7 @@ def start_secondary_med_prop(self): # Set up the secondary_medium wrapper: med_prop_class, med_prop_kwargs = media.get_global_med_prop_and_kwargs() if self.secondary_medium == "air" and med_prop_class == media.RefProp: - fluid_name = "AIR.PPF" + fluid_name = "air.ppf" else: fluid_name = self.secondary_medium if self.med_prop_sec is not None: @@ -73,6 +76,7 @@ def start_secondary_med_prop(self): def terminate_secondary_med_prop(self): if self.med_prop_sec is not None: self.med_prop_sec.terminate() + self.med_prop_sec = None @abc.abstractmethod def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> Tuple[float, float]: @@ -164,6 +168,26 @@ def calc_alpha_secondary(self, transport_properties) -> float: m_flow=self.m_flow_secondary ) + def calc_k(self, alpha_pri: float, alpha_sec: float) -> float: + """ + Calculate the overall heat transfer coefficient (k) of the heat exchanger. + + Args: + alpha_pri (float): Heat transfer coefficient for the primary medium. + alpha_sec (float): Heat transfer coefficient for the secondary medium. + + Returns: + float: Overall heat transfer coefficient (k). + """ + k_wall = self.calc_wall_heat_transfer() + k = (1 / ( + (1 / alpha_pri) * self.ratio_outer_to_inner_area + + (1 / k_wall) * self.ratio_outer_to_inner_area + + (1 / alpha_sec) + ) + ) + return k + def calc_wall_heat_transfer(self) -> float: """ Calculate the heat transfer coefficient inside the wall. diff --git a/vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py b/vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py new file mode 100644 index 0000000..da28398 --- /dev/null +++ b/vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py @@ -0,0 +1,126 @@ +import logging +from dataclasses import dataclass + +import numpy as np + +from .heat_transfer import HeatTransfer +from vclibpy.media import TransportProperties + +logger = logging.getLogger(__name__) + + +@dataclass +class PlateHeatExchangerGeometry: + """ + Geometry of a fin and tube heat exchanger with two rows of pipes in a shifted arrangement. + + Source: VDI-Wärmeatlas, Druckverlust und Wärmeübergang in Plattenwärmeübertragern, 11. Auflage, S.1687 + + + """ + wall_thickness: float # thickness of a plate + lambda_w: float # thermal conductivity of plates + amplitude: float # amplitude of a wavy plate + wave_length: float # wavelength of a wavy plate + phi: float # embossing angle + width: float # width of a plate + height: float # height of a plate + n_plates: float # number of plates + + @property + def X(self) -> float: + """return the wave number""" + return 2 * np.pi * self.amplitude / self.wave_length + + @property + def enlargement_factor(self) -> float: + """return surface enlargement factor""" + return (1 + np.sqrt(1 + self.X ** 2) + 4 * np.sqrt(1 + self.X ** 2 / 2)) / 6 + + @property + def d_h(self) -> float: + """return hydraulic diameter""" + return 4 * self.amplitude / self.enlargement_factor + + @property + def A(self) -> float: + """return plate surface""" + return self.height * self.width * self.enlargement_factor * self.n_plates + + +class VDIAtlasPlateHeatTransfer(HeatTransfer): + """ + VDI-Atlas based heat transfer estimation. + Check `AirToWallTransfer` for further argument options. + + This class assumes two pipes with a shifted arrangement. + + Args: + A_cross (float): Free cross-sectional area. + characteristic_length (float): Characteristic length d_a. + geometry_parameters (AirSourceHeatExchangerGeometry): + Geometry parameters for heat exchanger according to VDI-Atlas + """ + + def __init__(self, + geometry_parameters: PlateHeatExchangerGeometry): + super().__init__() + self.geometry_parameters = geometry_parameters + + def calc(self, transport_properties: TransportProperties, m_flow: float) -> float: + """ + Calculate heat transfer coefficient from refrigerant to the wall of the heat exchanger. + + The flow is assumed to be always turbulent and is based on a calibrated + Nusselt correlation. + + Args: + transport_properties (TransportProperties): Transport properties of the fluid. + m_flow (float): Mass flow rate of the fluid. + + Returns: + float: Heat transfer coefficient from refrigerant to HE in W/(m^2*K). + """ + Re = self.calc_reynolds( + dynamic_viscosity=transport_properties.dyn_vis, + m_flow=m_flow + ) + Nu = self.calc_turbulent_plate_nusselt(Re, transport_properties.Pr) + return Nu * transport_properties.lam / self.geometry_parameters.d_h + + def calc_reynolds(self, dynamic_viscosity: float, m_flow: float) -> float: + """ + Calculate Reynolds number for flow between two plates. + + Args: + dynamic_viscosity (float): Dynamic viscosity of the fluid. + m_flow (float): Mass flow rate. + width (float): Characteristic length (e.g., width) of the plate. + enlargement_factor (float): surface enlargement factor. + + Returns: + float: Reynolds number. + + """ + + return 2 * m_flow / (self.geometry_parameters.enlargement_factor * self.geometry_parameters.width * dynamic_viscosity) + + def calc_turbulent_plate_nusselt(self, Re: float, Pr: float) -> float: + """ + Calculate Nusselt Number based on Nusselt correlation VDI Atlas. + + Args: + Re (float): Reynolds number. + Pr (float): Prandtl number. + Zeta (float): pressure drop coefficient. + eta_by_eta_w (float): dynamic viscosity temperature correction factor. + + Returns: + float: Apparent heat transfer coefficient. + """ + c_q = 0.122 + q = 0.374 + eta_by_eta_w = 1 + Zeta = 1 + + return c_q * Pr ** (1/3) * eta_by_eta_w ** (1/6) * (Zeta * Re ** 2 * np.sin(2 * self.geometry_parameters.phi * np.pi / 180)) ** q diff --git a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py new file mode 100644 index 0000000..d0db543 --- /dev/null +++ b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py @@ -0,0 +1,262 @@ +import logging + +import numpy as np + +from vclibpy.components.heat_exchangers import utils +from vclibpy.datamodels import FlowsheetState, Inputs +from vclibpy.components.heat_exchangers import HeatExchanger + +logger = logging.getLogger(__name__) + + +class MovingBoundaryTmCondenser(HeatExchanger): + """ + Condenser class which implements the actual `calc` method. + + Assumptions: + - No phase changes in secondary medium + - cp of secondary medium is constant over heat-exchanger + + See parent classes for arguments. + """ + + def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): + """ + Calculate the heat exchanger with the Tm-Method based on the given inputs. + + The flowsheet state can be used to save important variables + during calculation for later analysis. + + Both return values are used to check if the heat transfer is valid or not. + + Args: + inputs (Inputs): The inputs for the calculation. + fs_state (FlowsheetState): The flowsheet state to save important variables. + + Returns: + Tuple[float, float]: + error: Error in percentage between the required and calculated heat flow rates. + dT_min: Minimal temperature difference (can be negative). + """ + self.m_flow_secondary = inputs.m_flow_con # [kg/s] + self.calc_secondary_cp(T=inputs.T_con) + + # First we separate the flow: + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = utils.separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_inlet, + state_min=self.state_outlet, + p=self.state_inlet.p + ) + Q = Q_sc + Q_lat + Q_sh + + T_in, T_sc, T_sh, T_out = utils.get_condenser_phase_temperatures_and_alpha( + heat_exchanger=self, inputs=inputs, + Q_sc=Q_sc, Q_lat=Q_lat, Q_sh=Q_sh + ) + + tra_prop_med = self.calc_transport_properties_secondary_medium((T_in + T_out) / 2) + alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) + + # 1. Regime: Subcooling + Q_sc_Tm, A_sc, A_sc_available = 0, 0, 0 + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_outlet.T): + # Get transport properties: + tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(state_q0, self.state_outlet) + alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_con) + + # Only use still available area: + k_sc = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sc = utils.calc_mean_temperature( + T_hot_in=state_q0.T, T_hot_out=self.state_outlet.T, + T_cold_in=T_in, T_cold_out=T_sc + ) + A_sc = utils.calc_area(Q_sc, k_sc, T_m_sc) + A_sc_available = min(self.A, A_sc) + Q_sc_Tm = A_sc_available * k_sc * T_m_sc + + # 2. Regime: Latent heat exchange + Q_lat_Tm, A_lat, A_lat_available = 0, 0, 0 + if not np.isclose(Q_lat, 0): + # Get transport properties: + alpha_ref_wall = self.calc_alpha_two_phase( + state_q0=state_q0, + state_q1=state_q1, + fs_state=fs_state, + inputs=inputs + ) + k_lat = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_lat = utils.calc_mean_temperature( + T_hot_in=state_q1.T, T_hot_out=state_q0.T, + T_cold_in=T_sc, T_cold_out=T_sh + ) + A_lat = utils.calc_area(Q_lat, k_lat, T_m_lat) + # Only use still available area: + A_lat_available = min(max(self.A - A_sc_available, 0), A_lat) + Q_lat_Tm = A_lat_available * k_lat * T_m_lat + + logger.debug(f"con_lat: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + # 3. Regime: Superheat heat exchange + Q_sh_Tm, A_sh = 0, 0 + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_inlet.T, state_q1.T): + # Get transport properties: + tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(self.state_inlet, state_q1) + alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_con) + + k_sh = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sh = utils.calc_mean_temperature( + T_hot_in=self.state_inlet.T, T_hot_out=state_q1.T, + T_cold_in=T_sh, T_cold_out=T_out + ) + A_sh = utils.calc_area(Q_sh, k_sh, T_m_sh) + # Only use still available area: + A_sh_available = min(max(self.A - A_sc_available - A_lat_available, 0), A_sh) + Q_sh_Tm = A_sh_available * k_sh * T_m_sh + + logger.debug(f"con_sh: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + A_necessary = A_sh + A_lat + A_sc + Q_Tm = Q_sh_Tm + Q_sc_Tm + Q_lat_Tm + error_A = (1 - A_necessary / self.A) * 100 + error = (Q_Tm / Q - 1) * 100 + # Get possible dT_min: + dT_min_in = self.state_outlet.T - T_in + dT_min_out = self.state_inlet.T - T_out + dT_min_LatSH = state_q1.T - T_sh + + fs_state.set(name="A_con_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in condenser") + fs_state.set(name="A_con_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in condenser") + fs_state.set(name="A_con_sc", value=A_sc, unit="m2", description="Area for subcooling heat exchange in condenser") + fs_state.set(name="error_A", value=error_A, unit="%", description="Area-percentage error for heat exchange in condenser") + + return error, min(dT_min_in, + dT_min_LatSH, + dT_min_out) + + +class MovingBoundaryTmEvaporator(HeatExchanger): + """ + Evaporator class which implements the actual `calc` method. + + Assumptions: + - No phase changes in secondary medium + - cp of secondary medium is constant over heat-exchanger + + See parent classes for arguments. + """ + + def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): + """ + Calculate the heat exchanger with the Tm-Method based on the given inputs. + + The flowsheet state can be used to save important variables + during calculation for later analysis. + + Both return values are used to check if the heat transfer is valid or not. + + Args: + inputs (Inputs): The inputs for the calculation. + fs_state (FlowsheetState): The flowsheet state to save important variables. + + Returns: + Tuple[float, float]: + error: Error in percentage between the required and calculated heat flow rates. + dT_min: Minimal temperature difference (can be negative). + """ + self.m_flow_secondary = inputs.m_flow_eva # [kg/s] + self.calc_secondary_cp(T=inputs.T_eva_in) + + # First we separate the flow: + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = utils.separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_outlet, + state_min=self.state_inlet, + p=self.state_outlet.p + ) + + Q = Q_sc + Q_lat + Q_sh + + # Note: As Q_eva_Tm has to converge to Q_eva (m_ref*delta_h), we can safely + # calculate the output temperature. + T_mean = inputs.T_eva_in - Q / (self.m_flow_secondary_cp * 2) + tra_prop_med = self.calc_transport_properties_secondary_medium(T_mean) + alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) + # alpha_med_wall = 26 + + # Calculate secondary_medium side temperatures: + T_sh = inputs.T_eva_in - Q_sh / self.m_flow_secondary_cp + T_sc = T_sh - Q_lat / self.m_flow_secondary_cp + T_out = T_sc - Q_sc / self.m_flow_secondary_cp + + # 1. Regime: Superheating + Q_sh_Tm, A_sh, A_sh_available = 0, 0, 0 + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_outlet.T, state_q1.T): + # Get transport properties: + tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(self.state_outlet, state_q1) + alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_eva) + + k_sh = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sh = utils.calc_mean_temperature( + T_hot_in=inputs.T_eva_in, T_hot_out=T_sh, + T_cold_in=state_q1.T, T_cold_out=self.state_outlet.T + ) + # Only use still available area: + A_sh = utils.calc_area(Q_sh, k_sh, T_m_sh) + A_sh_available = min(self.A, A_sh) + Q_sh_Tm = A_sh_available * k_sh * T_m_sh + + logger.debug(f"eva_sh: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + # 2. Regime: Latent heat exchange + Q_lat_Tm, A_lat, A_lat_available = 0, 0, 0 + if not np.isclose(Q_lat, 0): + alpha_ref_wall = self.calc_alpha_two_phase( + state_q0=state_q0, + state_q1=state_q1, + fs_state=fs_state, + inputs=inputs + ) + k_lat = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_lat = utils.calc_mean_temperature( + T_hot_in=T_sh, T_hot_out=T_sc, + T_cold_in=state_q0.T, T_cold_out=state_q1.T + ) + A_lat = utils.calc_area(Q_lat, k_lat, T_m_lat) + # Only use still available area: + A_lat_available = min(max(self.A - A_sh_available, 0), A_lat) + Q_lat_Tm = A_lat_available * k_lat * T_m_lat + + logger.debug(f"eva_lat: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + # 3. Regime: Subcooling + Q_sc_Tm, A_sc = 0, 0 + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_inlet.T): + # Get transport properties: + tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(state_q0, self.state_inlet) + alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_eva) + + # Only use still available area: + k_sc = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sc = utils.calc_mean_temperature( + T_hot_in=T_sc, T_hot_out=T_out, + T_cold_in=self.state_inlet.T, T_cold_out=state_q0.T + ) + A_sc = utils.calc_area(Q_sc, k_sc, T_m_sc) + A_sc_available = min(max(self.A - A_sh_available - A_lat_available, 0), A_sc) + Q_sc_Tm = A_sc_available * k_sc * T_m_sc + + A_necessary = A_sc + A_lat + A_sh + Q_Tm = Q_sh_Tm + Q_sc_Tm + Q_lat_Tm + error_A = (A_necessary / self.A - 1) * 100 + error = (Q_Tm / Q - 1) * 100 + # Get dT_min + dT_min_in = inputs.T_eva_in - self.state_outlet.T + dT_min_out = T_out - self.state_inlet.T + + fs_state.set(name="A_eva_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in evaporator") + fs_state.set(name="A_eva_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in evaporator") + + return error, min(dT_min_out, dT_min_in) diff --git a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py index 19d29cd..e1cbe14 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py @@ -1,106 +1,69 @@ -import abc import logging import numpy as np from vclibpy.datamodels import FlowsheetState, Inputs from vclibpy.components.heat_exchangers.ntu import BasicNTU -from vclibpy.media import ThermodynamicState +from vclibpy.components.heat_exchangers.utils import separate_phases, get_condenser_phase_temperatures_and_alpha logger = logging.getLogger(__name__) -class MovingBoundaryNTU(BasicNTU, abc.ABC): +def iterate_area(heat_exchanger: BasicNTU, dT_max, alpha_pri, alpha_sec, Q) -> float: """ - Moving boundary NTU based heat exchanger. + Iteratively calculates the required area for the heat exchange. - See parent classe for arguments. - """ - - def separate_phases(self, state_max: ThermodynamicState, state_min: ThermodynamicState, p: float): - """ - Separates a flow with possible phase changes into three parts: - subcooling (sc), latent phase change (lat), and superheating (sh) - at the given pressure. - - Args: - state_max (ThermodynamicState): State with higher enthalpy. - state_min (ThermodynamicState): State with lower enthalpy. - p (float): Pressure of phase change. - - Returns: - Tuple[float, float, float, ThermodynamicState, ThermodynamicState]: - Q_sc: Heat for subcooling. - Q_lat: Heat for latent phase change. - Q_sh: Heat for superheating. - state_q0: State at vapor quality 0 and the given pressure. - state_q1: State at vapor quality 1 and the given pressure. - """ - # Get relevant states: - state_q0 = self.med_prop.calc_state("PQ", p, 0) - state_q1 = self.med_prop.calc_state("PQ", p, 1) - Q_sc = max(0.0, - min((state_q0.h - state_min.h), - (state_max.h - state_min.h))) * self.m_flow - Q_lat = max(0.0, - (min(state_max.h, state_q1.h) - - max(state_min.h, state_q0.h))) * self.m_flow - Q_sh = max(0.0, - min((state_max.h - state_q1.h), - (state_max.h - state_min.h))) * self.m_flow - return Q_sc, Q_lat, Q_sh, state_q0, state_q1 - - def iterate_area(self, dT_max, alpha_pri, alpha_sec, Q) -> float: - """ - Iteratively calculates the required area for the heat exchange. + Args: + heat_exchanger (BasicNTU): An instance of the BasicNTU or children classes + dT_max (float): Maximum temperature differential. + alpha_pri (float): Heat transfer coefficient for the primary medium. + alpha_sec (float): Heat transfer coefficient for the secondary medium. + Q (float): Heat flow rate. - Args: - dT_max (float): Maximum temperature differential. - alpha_pri (float): Heat transfer coefficient for the primary medium. - alpha_sec (float): Heat transfer coefficient for the secondary medium. - Q (float): Heat flow rate. - - Returns: - float: Required area for heat exchange. - """ - _accuracy = 1e-6 # square mm - _step = 1.0 - R = self.calc_R() - k = self.calc_k(alpha_pri, alpha_sec) - m_flow_cp_min = self.calc_m_flow_cp_min() - # First check if point is feasible at all - if dT_max <= 0: - return self.A - eps_necessary = Q / (m_flow_cp_min * dT_max) - - # Special cases: - # --------------- - # eps is equal or higher than 1, an infinite amount of area would be necessary. - if eps_necessary >= 1: - return self.A - # eps is lower or equal to zero: No Area required (Q<=0) - if eps_necessary <= 0: - return 0 - - area = 0.0 - while True: - NTU = self.calc_NTU(area, k, m_flow_cp_min) - eps = self.calc_eps(R, NTU) - if eps >= eps_necessary: - if _step <= _accuracy: - break - else: - # Go back - area -= _step - _step /= 10 - continue - if _step < _accuracy and area > self.A: + Returns: + float: Required area for heat exchange. + """ + _accuracy = 1e-6 # square mm + _step = 1.0 + R = heat_exchanger.calc_R() + k = heat_exchanger.calc_k(alpha_pri, alpha_sec) + m_flow_cp_min = heat_exchanger.calc_m_flow_cp_min() + # First check if point is feasible at all + if dT_max <= 0: + return heat_exchanger.A + eps_necessary = Q / (m_flow_cp_min * dT_max) + + # Special cases: + # --------------- + # eps is equal or higher than 1, an infinite amount of area would be necessary. + if eps_necessary >= 1: + return heat_exchanger.A + # eps is lower or equal to zero: No Area required (Q<=0) + if eps_necessary <= 0: + return 0 + + area = 0.0 + while True: + if heat_exchanger.flow_type == "cross" and area == 0.0: + eps = 0.0 + else: + NTU = heat_exchanger.calc_NTU(area, k, m_flow_cp_min) + eps = heat_exchanger.calc_eps(R, NTU) + if eps >= eps_necessary: + if _step <= _accuracy: break - area += _step + else: + # Go back + area -= _step + _step /= 10 + continue + if _step < _accuracy and area > heat_exchanger.A: + break + area += _step - return min(area, self.A) + return min(area, heat_exchanger.A) -class MovingBoundaryNTUCondenser(MovingBoundaryNTU): +class MovingBoundaryNTUCondenser(BasicNTU): """ Condenser class which implements the actual `calc` method. @@ -130,52 +93,50 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): dT_min: Minimal temperature difference (can be negative). """ self.m_flow_secondary = inputs.m_flow_con # [kg/s] - self.calc_secondary_cp(T=inputs.T_con_in) + self.calc_secondary_cp(T=inputs.T_con) # First we separate the flow: - Q_sc, Q_lat, Q_sh, state_q0, state_q1 = self.separate_phases( - self.state_inlet, - self.state_outlet, - self.state_inlet.p + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_inlet, + state_min=self.state_outlet, + p=self.state_inlet.p ) Q = Q_sc + Q_lat + Q_sh - # Note: As Q_con_ntu has to converge to Q_con (m_ref*delta_h), we can safely - # calculate the output temperature. + T_in, T_sc, T_sh, T_out = get_condenser_phase_temperatures_and_alpha( + heat_exchanger=self, inputs=inputs, + Q_sc=Q_sc, Q_lat=Q_lat, Q_sh=Q_sh + ) - T_mean = inputs.T_con_in + self.calc_secondary_Q_flow(Q) / (self.m_flow_secondary_cp * 2) - tra_prop_med = self.calc_transport_properties_secondary_medium(T_mean) + tra_prop_med = self.calc_transport_properties_secondary_medium((T_in + T_out) / 2) alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) - # Calculate secondary_medium side temperatures: - # Assumption loss is the same correlation for each regime - T_sc = inputs.T_con_in + self.calc_secondary_Q_flow(Q_sc) / self.m_flow_secondary_cp - T_sh = T_sc + self.calc_secondary_Q_flow(Q_lat) / self.m_flow_secondary_cp - T_out = T_sh + self.calc_secondary_Q_flow(Q_sh) / self.m_flow_secondary_cp - # 1. Regime: Subcooling Q_sc_ntu, A_sc = 0, 0 - if Q_sc > 0 and (state_q0.T != self.state_outlet.T): + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_outlet.T): self.set_primary_cp((state_q0.h - self.state_outlet.h) / (state_q0.T - self.state_outlet.T)) # Get transport properties: tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(state_q0, self.state_outlet) alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_con) # Only use still available area: - A_sc = self.iterate_area(dT_max=(state_q0.T - inputs.T_con_in), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_sc) + A_sc = iterate_area(heat_exchanger=self, + dT_max=(state_q0.T - T_in), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_sc) A_sc = min(self.A, A_sc) - Q_sc_ntu, k_sc = self.calc_Q_ntu(dT_max=(state_q0.T - inputs.T_con_in), + Q_sc_ntu, k_sc = self.calc_Q_ntu(dT_max=(state_q0.T - T_in), alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall, A=A_sc) # 2. Regime: Latent heat exchange Q_lat_ntu, A_lat = 0, 0 - if Q_lat > 0: + if not np.isclose(Q_lat, 0): self.set_primary_cp(np.inf) # Get transport properties: alpha_ref_wall = self.calc_alpha_two_phase( @@ -185,10 +146,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): inputs=inputs ) - A_lat = self.iterate_area(dT_max=(state_q1.T - T_sc), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_lat) + A_lat = iterate_area(heat_exchanger=self, + dT_max=(state_q1.T - T_sc), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_lat) # Only use still available area: A_lat = min(self.A - A_sc, A_lat) @@ -200,7 +162,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 3. Regime: Superheat heat exchange Q_sh_ntu, A_sh = 0, 0 - if Q_sh and (self.state_inlet.T != state_q1.T): + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_inlet.T, state_q1.T): self.set_primary_cp((self.state_inlet.h - state_q1.h) / (self.state_inlet.T - state_q1.T)) # Get transport properties: tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(self.state_inlet, state_q1) @@ -218,20 +180,22 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): Q_ntu = Q_sh_ntu + Q_sc_ntu + Q_lat_ntu error = (Q_ntu / Q - 1) * 100 # Get possible dT_min: - dT_min_in = self.state_outlet.T - inputs.T_con_in + dT_min_in = self.state_outlet.T - T_in dT_min_out = self.state_inlet.T - T_out dT_min_LatSH = state_q1.T - T_sh - fs_state.set(name="A_con_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in condenser") + fs_state.set(name="A_con_sh", value=A_sh, unit="m2", + description="Area for superheat heat exchange in condenser") fs_state.set(name="A_con_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in condenser") - fs_state.set(name="A_con_sc", value=A_sc, unit="m2", description="Area for subcooling heat exchange in condenser") + fs_state.set(name="A_con_sc", value=A_sc, unit="m2", + description="Area for subcooling heat exchange in condenser") return error, min(dT_min_in, dT_min_LatSH, dT_min_out) -class MovingBoundaryNTUEvaporator(MovingBoundaryNTU): +class MovingBoundaryNTUEvaporator(BasicNTU): """ Evaporator class which implements the actual `calc` method. @@ -264,10 +228,12 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): self.calc_secondary_cp(T=inputs.T_eva_in) # First we separate the flow: - Q_sc, Q_lat, Q_sh, state_q0, state_q1 = self.separate_phases( - self.state_outlet, - self.state_inlet, - self.state_inlet.p + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_outlet, + state_min=self.state_inlet, + p=self.state_outlet.p ) Q = Q_sc + Q_lat + Q_sh @@ -285,17 +251,18 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 1. Regime: Superheating Q_sh_ntu, A_sh = 0, 0 - if Q_sh and (self.state_outlet.T != state_q1.T): + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_outlet.T, state_q1.T): self.set_primary_cp((self.state_outlet.h - state_q1.h) / (self.state_outlet.T - state_q1.T)) # Get transport properties: tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(self.state_outlet, state_q1) alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_eva) if Q_lat > 0: - A_sh = self.iterate_area(dT_max=(inputs.T_eva_in - state_q1.T), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_sh) + A_sh = iterate_area(heat_exchanger=self, + dT_max=(inputs.T_eva_in - state_q1.T), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_sh) else: # if only sh is present --> full area: A_sh = self.A @@ -312,7 +279,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 2. Regime: Latent heat exchange Q_lat_ntu, A_lat = 0, 0 - if Q_lat > 0: + if not np.isclose(Q_lat, 0): self.set_primary_cp(np.inf) alpha_ref_wall = self.calc_alpha_two_phase( @@ -323,10 +290,10 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): ) if Q_sc > 0: - A_lat = self.iterate_area(dT_max=(T_sh - self.state_inlet.T), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_lat) + A_lat = iterate_area(heat_exchanger=self, dT_max=(T_sh - self.state_inlet.T), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_lat) else: A_lat = self.A - A_sh @@ -340,7 +307,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 3. Regime: Subcooling Q_sc_ntu, A_sc = 0, 0 - if Q_sc > 0 and (state_q0.T != self.state_inlet.T): + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_inlet.T): self.set_primary_cp((state_q0.h - self.state_inlet.h) / (state_q0.T - self.state_inlet.T)) # Get transport properties: tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(state_q0, self.state_inlet) @@ -360,7 +327,9 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): dT_min_in = inputs.T_eva_in - self.state_outlet.T dT_min_out = T_out - self.state_inlet.T - fs_state.set(name="A_eva_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in evaporator") - fs_state.set(name="A_eva_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in evaporator") + fs_state.set(name="A_eva_sh", value=A_sh, unit="m2", + description="Area for superheat heat exchange in evaporator") + fs_state.set(name="A_eva_lat", value=A_lat, unit="m2", + description="Area for latent heat exchange in evaporator") return error, min(dT_min_out, dT_min_in) diff --git a/vclibpy/components/heat_exchangers/ntu.py b/vclibpy/components/heat_exchangers/ntu.py index db7c92b..e1fe9de 100644 --- a/vclibpy/components/heat_exchangers/ntu.py +++ b/vclibpy/components/heat_exchangers/ntu.py @@ -28,7 +28,7 @@ class BasicNTU(HeatExchanger, abc.ABC): Typical values are around 20-30. """ - def __init__(self, flow_type: str, ratio_outer_to_inner_area: float, **kwargs): + def __init__(self, flow_type: str, **kwargs): """ Initializes BasicNTU. @@ -45,7 +45,6 @@ def __init__(self, flow_type: str, ratio_outer_to_inner_area: float, **kwargs): **kwargs: Additional keyword arguments passed to the parent class. """ super().__init__(**kwargs) - self.ratio_outer_to_inner_area = ratio_outer_to_inner_area # Set primary cp: self._primary_cp = None @@ -102,26 +101,6 @@ def calc_R(self) -> float: return self.m_flow_secondary_cp / m_flow_pri_cp return m_flow_pri_cp / self.m_flow_secondary_cp - def calc_k(self, alpha_pri: float, alpha_sec: float) -> float: - """ - Calculate the overall heat transfer coefficient (k) of the heat exchanger. - - Args: - alpha_pri (float): Heat transfer coefficient for the primary medium. - alpha_sec (float): Heat transfer coefficient for the secondary medium. - - Returns: - float: Overall heat transfer coefficient (k). - """ - k_wall = self.calc_wall_heat_transfer() - k = (1 / ( - (1 / alpha_pri) * self.ratio_outer_to_inner_area + - (1 / k_wall) * self.ratio_outer_to_inner_area + - (1 / alpha_sec) - ) - ) - return k - @staticmethod def calc_NTU(A: float, k: float, m_flow_cp: float) -> float: """ diff --git a/vclibpy/components/heat_exchangers/utils.py b/vclibpy/components/heat_exchangers/utils.py new file mode 100644 index 0000000..eb003f0 --- /dev/null +++ b/vclibpy/components/heat_exchangers/utils.py @@ -0,0 +1,93 @@ +import numpy as np + +from vclibpy import Inputs +from vclibpy.components.heat_exchangers import HeatExchanger + +from vclibpy.media import ThermodynamicState + + +def calc_area(Q, k, Tm): + if Tm * k == 0: + return 1e10 # A really large number, but not np.inf to still be able to calculate with it + else: + return Q / (k * Tm) + + +def calc_mean_temperature( + T_hot_in: float, T_hot_out: float, + T_cold_in: float, T_cold_out: float, + flow_type: str = "counter" +): + if flow_type == "counter": + dT_A = T_hot_in - T_cold_out + dT_B = T_hot_out - T_cold_in + elif flow_type == "parallel": + dT_A = T_hot_in - T_cold_in + dT_B = T_hot_out - T_cold_out + else: + raise TypeError("Given flow_type is not supported yet") + if dT_B < 0 or dT_A < 0: + return 0 # Heat can't be transferred + if np.isclose(dT_B, 0) or np.isclose(dT_A, 0) or np.isclose(dT_B, dT_A): + return abs((dT_A + dT_B) / 2) + return (dT_A - dT_B) / np.log(dT_A / dT_B) + + +def separate_phases(m_flow, med_prop, state_max: ThermodynamicState, state_min: ThermodynamicState, p: float): + """ + Separates a flow with possible phase changes into three parts: + subcooling (sc), latent phase change (lat), and superheating (sh) + at the given pressure. + + Args: + state_max (ThermodynamicState): State with higher enthalpy. + state_min (ThermodynamicState): State with lower enthalpy. + p (float): Pressure of phase change. + + Returns: + Tuple[float, float, float, ThermodynamicState, ThermodynamicState]: + Q_sc: Heat for subcooling. + Q_lat: Heat for latent phase change. + Q_sh: Heat for superheating. + state_q0: State at vapor quality 0 and the given pressure. + state_q1: State at vapor quality 1 and the given pressure. + """ + # Get relevant states: + state_q0 = med_prop.calc_state("PQ", p, 0) + state_q1 = med_prop.calc_state("PQ", p, 1) + Q_sc = max(0.0, + min((state_q0.h - state_min.h), + (state_max.h - state_min.h))) * m_flow + Q_lat = max(0.0, + (min(state_max.h, state_q1.h) - + max(state_min.h, state_q0.h))) * m_flow + Q_sh = max(0.0, + min((state_max.h - state_q1.h), + (state_max.h - state_min.h))) * m_flow + return Q_sc, Q_lat, Q_sh, state_q0, state_q1 + + +def get_condenser_phase_temperatures_and_alpha( + inputs: Inputs, + Q_sc: float, Q_lat: float, Q_sh: float, + heat_exchanger: HeatExchanger +): + cp = heat_exchanger.m_flow_secondary_cp + Q = Q_sc + Q_lat + Q_sh + + # Calculate secondary_medium side temperatures: + # Assumption loss is the same correlation for each regime + + # Calculate secondary_medium side temperatures: + # Assumption loss is the same correlation for each regime + if inputs.uses_condenser_inlet: + T_in = inputs.T_con_in + T_sc = inputs.T_con_in + heat_exchanger.calc_secondary_Q_flow(Q_sc) / cp + T_sh = T_sc + heat_exchanger.calc_secondary_Q_flow(Q_lat) / cp + T_out = T_sh + heat_exchanger.calc_secondary_Q_flow(Q_sh) / cp + else: + T_out = inputs.T_con_out + T_sh = T_out - heat_exchanger.calc_secondary_Q_flow(Q_sh) / cp + T_sc = T_sh - heat_exchanger.calc_secondary_Q_flow(Q_lat) / cp + T_in = T_sc - heat_exchanger.calc_secondary_Q_flow(Q_sc) / cp + return T_in, T_sc, T_sh, T_out diff --git a/vclibpy/datamodels.py b/vclibpy/datamodels.py index 969ac4c..a77ae7d 100644 --- a/vclibpy/datamodels.py +++ b/vclibpy/datamodels.py @@ -129,8 +129,17 @@ def convert_to_str_value_format(self, with_unit_and_description: bool) -> Dict[s dict: Containing all variables and values. """ if with_unit_and_description: - return {f"{k} in {v.unit} ({v.description})": v.value for k, v in self.items()} - return {k: v.value for k, v in self.items()} + return {f"{k} in {v.unit} ({v.description})": v.value for k, v in self.items_not_none()} + return {k: v.value for k, v in self.items_not_none()} + + def get_name(self): + """Get the name based on variable names and rounded values""" + return ";".join([k + "=" + str(round(v.value, 3)).replace(".", "_") + for k, v in self.items_not_none()]) + + def items_not_none(self): + """Get all items where the value is not None""" + return {k: v for k, v in self.items() if v.value is not None}.items() class FlowsheetState(VariableContainer): @@ -168,6 +177,7 @@ def __init__( n: float = None, T_eva_in: float = None, T_con_in: float = None, + T_con_out: float = None, m_flow_eva: float = None, m_flow_con: float = None, dT_eva_superheating: float = None, @@ -178,10 +188,17 @@ def __init__( Initializes an Inputs object with parameters representing external conditions for the vapor compression cycle. + Note that some inputs presuppose each other. For example, you can not + specify both `T_con_in`, `T_con_out`, and `m_flow_con`, + as one needs to be free for the simulation. + The current algorithms require `m_flow_con` + and either `T_con_in` or `T_con_out`. + Args: n (float): Relative compressor speed between 0 and 1 (unit: -). T_eva_in (float): Secondary side evaporator inlet temperature (unit: K). T_con_in (float): Secondary side condenser inlet temperature (unit: K). + T_con_out (float): Secondary side condenser inlet temperature (unit: K). m_flow_eva (float): Secondary side evaporator mass flow rate (unit: kg/s). m_flow_con (float): Secondary side condenser mass flow rate (unit: kg/s). dT_eva_superheating (float): Super-heating after evaporator (unit: K). @@ -189,6 +206,10 @@ def __init__( T_ambient (float): Ambient temperature of the machine (unit: K). """ super().__init__() + if T_con_out is not None and T_con_in is not None: + raise ValueError("You can either specify condenser inlet or " + "outlet temperature, not both.") + self.set( name="n", value=n, @@ -207,6 +228,12 @@ def __init__( unit="K", description="Secondary side condenser inlet temperature" ) + self.set( + name="T_con_out", + value=T_con_out, + unit="K", + description="Secondary side condenser outlet temperature" + ) self.set( name="m_flow_con", value=m_flow_con, @@ -239,3 +266,19 @@ def __init__( unit="K", description="Ambient temperature of machine" ) + + @property + def T_con(self): + """Returns either T_con_in or T_con_out, depending on what is specifiec.""" + if self.T_con_in is not None: + return self.T_con_in + return self.T_con_out + + @property + def uses_condenser_inlet(self): + """Indicates if condenser inlet or outlet is specified""" + return self.T_con_in is not None + + +if __name__ == '__main__': + print(Inputs().T_con_out) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index 9ab595b..20e6e6c 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -1,6 +1,7 @@ import logging from typing import List import numpy as np +from scipy.optimize import fsolve from abc import abstractmethod import matplotlib.pyplot as plt @@ -64,13 +65,52 @@ def setup_new_fluid(self, fluid): self.fluid = fluid def terminate(self): - self.med_prop.terminate() + if self.med_prop is not None: + self.med_prop.terminate() + self.med_prop = None for component in self.get_all_components(): component.terminate_secondary_med_prop() + component.med_prop = None def get_all_components(self) -> List[BaseComponent]: return [self.condenser, self.evaporator] + def get_start_condensing_pressure(self, inputs: Inputs): + if inputs.uses_condenser_inlet: + T_3_start = inputs.T_con_in + inputs.dT_con_subcooling + else: + try: + dT_con_start = inputs.dT_con_start + except AttributeError: + dT_con_start = 10 + T_3_start = inputs.T_con_out - dT_con_start + p_2_start = self.med_prop.calc_state("TQ", T_3_start, 0).p + return p_2_start + + def improve_first_condensing_guess(self, inputs: Inputs, m_flow_guess, p_2_guess, dT_pinch_assumption=0): + self.condenser.m_flow_secondary = inputs.m_flow_con # [kg/s] + self.condenser.calc_secondary_cp(T=inputs.T_con) + + def nonlinear_func(p_2, *args): + _flowsheet, _inputs, _m_flow_ref, _dT_pinch = args + state_q0 = _flowsheet.med_prop.calc_state("PQ", p_2, 0) + state_q1 = _flowsheet.med_prop.calc_state("PQ", p_2, 1) + state_3 = _flowsheet.med_prop.calc_state("PT", p_2, state_q0.T - _inputs.dT_con_subcooling) + Q_water_till_q1 = (state_q1.h - state_3.h) * _m_flow_ref + T_water_q1 = _inputs.T_con_in + Q_water_till_q1 / _flowsheet.condenser.m_flow_secondary_cp + return T_water_q1 + _dT_pinch - state_q1.T + + p_2_guess_optimized = fsolve( + func=nonlinear_func, + x0=p_2_guess, + args=(self, inputs, m_flow_guess, dT_pinch_assumption) + )[0] + return p_2_guess_optimized + + def get_start_evaporating_pressure(self, inputs: Inputs, dT_pinch: float = 0): + T_1_start = inputs.T_eva_in - inputs.dT_eva_superheating - dT_pinch + return self.med_prop.calc_state("TQ", T_1_start, 1).p + def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): """ Calculate the steady-state performance of a vapor compression cycle @@ -126,15 +166,22 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): # Settings min_iteration_step = kwargs.pop("min_iteration_step", 1) save_path_plots = kwargs.get("save_path_plots", None) - input_name = ";".join([k + "=" + str(np.round(v.value, 3)).replace(".", "_") - for k, v in inputs.get_variables().items()]) show_iteration = kwargs.get("show_iteration", False) use_quick_solver = kwargs.pop("use_quick_solver", True) err_ntu = kwargs.pop("max_err_ntu", 0.5) err_dT_min = kwargs.pop("max_err_dT_min", 0.1) max_num_iterations = kwargs.pop("max_num_iterations", 1e5) + dT_pinch_con_guess = kwargs.pop("dT_pinch_con_guess", 0) + dT_pinch_eva_guess = kwargs.pop("dT_pinch_eva_guess", 0) + improve_first_condensing_guess = kwargs.pop("improve_first_condensing_guess", False) p_1_history = [] p_2_history = [] + error_con_history = [] + error_eva_history = [] + dT_eva_history = [] + dT_con_history = [] + error_con, dT_min_eva, dT_min_con, error_eva = np.nan, np.nan, np.nan, np.nan + plot_last = -100 if use_quick_solver: step_p1 = kwargs.get("step_max", 10000) @@ -149,19 +196,15 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): self.setup_new_fluid(fluid) # First: Iterate with given conditions to get the 4 states and the mass flow rate: - T_1_start = inputs.T_eva_in - inputs.dT_eva_superheating - T_3_start = inputs.T_con_in + inputs.dT_con_subcooling - p_1_start = self.med_prop.calc_state("TQ", T_1_start, 1).p - p_2_start = self.med_prop.calc_state("TQ", T_3_start, 0).p - p_1_next = p_1_start - p_2_next = p_2_start + p_1_next = self.get_start_evaporating_pressure(inputs=inputs, dT_pinch=dT_pinch_eva_guess) + p_2_next = self.get_start_condensing_pressure(inputs=inputs) fs_state = FlowsheetState() # Always log what is happening in the whole flowsheet fs_state.set(name="Q_con", value=1, unit="W", description="Condenser heat flow rate") fs_state.set(name="COP", value=0, unit="-", description="Coefficient of performance") if show_iteration: - fig_iterations, ax_iterations = plt.subplots(2) + fig_iterations, ax_iterations = plt.subplots(3, 2, sharex=True) num_iterations = 0 @@ -178,13 +221,29 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): p_1 = p_1_next p_2 = p_2_next - p_1_history.append(p_1) - p_2_history.append(p_2) + p_1_history.append(p_1 / 1e5) + p_2_history.append(p_2 / 1e5) + error_con_history.append(error_con) + error_eva_history.append(error_eva) + dT_con_history.append(dT_min_con) + dT_eva_history.append(dT_min_eva) + if show_iteration: - ax_iterations[0].cla() - ax_iterations[1].cla() - ax_iterations[0].scatter(list(range(len(p_1_history))), p_1_history) - ax_iterations[1].scatter(list(range(len(p_2_history))), p_2_history) + for ax in ax_iterations.flatten(): + ax.clear() + iterations = list(range(len(p_1_history)))[plot_last:] + ax_iterations[0, 0].set_ylabel("error_eva in %") + ax_iterations[0, 1].set_ylabel("error_con in %") + ax_iterations[1, 0].set_ylabel("$\Delta T_\mathrm{Min}$ in K") + ax_iterations[1, 1].set_ylabel("$\Delta T_\mathrm{Min}$ in K") + ax_iterations[2, 0].set_ylabel("$p_1$ in bar") + ax_iterations[2, 1].set_ylabel("$p_2$ in bar") + ax_iterations[0, 0].scatter(iterations, error_eva_history[plot_last:]) + ax_iterations[0, 1].scatter(iterations, error_con_history[plot_last:]) + ax_iterations[1, 0].scatter(iterations, dT_eva_history[plot_last:]) + ax_iterations[1, 1].scatter(iterations, dT_con_history[plot_last:]) + ax_iterations[2, 0].scatter(iterations, p_1_history[plot_last:]) + ax_iterations[2, 1].scatter(iterations, p_2_history[plot_last:]) plt.draw() plt.pause(1e-5) @@ -206,18 +265,31 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): step_p1 /= 10 continue - # Calculate the states based on the given flowsheet try: - self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) + error_eva, dT_min_eva, error_con, dT_min_con = self.calculate_cycle_for_pressures( + p_1=p_1, p_2=p_2, inputs=inputs, fs_state=fs_state + ) except ValueError as err: logger.error("An error occurred while calculating states. " "Can't guess next pressures, thus, exiting: %s", err) return - if save_path_plots is not None and num_iterations == 1 and show_iteration: - self.plot_cycle(save_path=save_path_plots.joinpath(f"{input_name}_initialization.png"), inputs=inputs) - # Check heat exchangers: - error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) + if num_iterations == 1: + if improve_first_condensing_guess: + p_2_next = self.improve_first_condensing_guess( + inputs=inputs, + m_flow_guess=self.condenser.m_flow, + p_2_guess=p_2, + dT_pinch_assumption=dT_pinch_con_guess + ) + + if save_path_plots is not None and show_iteration: + input_name = inputs.get_name() + self.plot_cycle( + save_path=save_path_plots.joinpath(f"{input_name}_initialization.png"), + inputs=inputs + ) + if not isinstance(error_eva, float): print(error_eva) if error_eva < 0: @@ -233,7 +305,6 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): p_1_next = p_1 + step_p1 continue - error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) if error_con < 0: p_2_next = p_2 + step_p2 continue @@ -273,15 +344,99 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): if show_iteration: plt.close(fig_iterations) + return self.calculate_outputs_for_valid_pressures( + p_1=p_1, p_2=p_2, inputs=inputs, fs_state=fs_state, + save_path_plots=save_path_plots + ) + + def calc_steady_state_fsolve(self, inputs: Inputs, fluid: str = None, **kwargs): + # Settings + max_err = kwargs.pop("max_err_ntu", 0.5) + + save_path_plots = kwargs.get("save_path_plots", None) + dT_pinch_eva_guess = kwargs.pop("dT_pinch_eva_guess", 0) + + # Setup fluid: + if fluid is None: + fluid = self.fluid + self.setup_new_fluid(fluid) + + # First: Iterate with given conditions to get the 4 states and the mass flow rate: + p_1_next = self.get_start_evaporating_pressure(inputs=inputs, dT_pinch=dT_pinch_eva_guess) + p_2_next = self.get_start_condensing_pressure(inputs=inputs) + + def nonlinear_func(x, *args): + _flowsheet, _inputs, _fs_state, _max_err = args + _p_1, _p_2 = x + if not (_p_1 < _p_2 < self._p_max): + return 1000, 1000 + if not (self._p_min < _p_1 < _p_2): + return 1000, 1000 + _error_eva, dT_min_eva, _error_con, dT_min_con = _flowsheet.calculate_cycle_for_pressures( + p_1=x[0], p_2=x[1], inputs=_inputs, fs_state=_fs_state + ) + + if 0 <= _error_eva < _max_err: + _error_eva = 0 + if 0 <= _error_con < _max_err: + _error_con = 0 + if 0 > _error_eva: + _error_eva *= 5 + if 0 > _error_con: + _error_con *= 5 + return _error_eva, _error_con + + fs_state = FlowsheetState() # Always log what is happening in the whole flowsheet + try: + args = (self, inputs, fs_state, max_err) + p_optimized = fsolve( + func=nonlinear_func, + x0=np.array([p_1_next, p_2_next]), + args=args + ) + except Exception as err: + logger.error("An error occurred while calculating states using fsolve: %s", err) + return + error_con, error_eva = nonlinear_func(p_optimized, *args) + print(f"{error_con=}, {error_eva=}") + + p_1, p_2 = p_optimized + return self.calculate_outputs_for_valid_pressures( + p_1=p_1, p_2=p_2, inputs=inputs, fs_state=fs_state, + save_path_plots=save_path_plots + ) + + def calculate_cycle_for_pressures(self, p_1: float, p_2: float, inputs: Inputs, fs_state: FlowsheetState): + # Calculate the states based on the given flowsheet + self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) + # Check heat exchangers: + error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) + error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) + return error_eva, dT_min_eva, error_con, dT_min_con + + def calculate_outputs_for_valid_pressures( + self, + p_1, + p_2, + fs_state: FlowsheetState, + inputs: Inputs, + save_path_plots + ): + self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) # Calculate the heat flow rates for the selected states. Q_con = self.condenser.calc_Q_flow() Q_con_outer = self.condenser.calc_secondary_Q_flow(Q_con) Q_eva = self.evaporator.calc_Q_flow() Q_eva_outer = self.evaporator.calc_secondary_Q_flow(Q_eva) - self.evaporator.calc(inputs=inputs, fs_state=fs_state) - self.condenser.calc(inputs=inputs, fs_state=fs_state) + error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) + error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) P_el = self.calc_electrical_power(fs_state=fs_state, inputs=inputs) - T_con_out = inputs.T_con_in + Q_con_outer / self.condenser.m_flow_secondary_cp + if inputs.uses_condenser_inlet: + T_con_in = inputs.T_con_in + T_con_out = T_con_in + Q_con_outer / self.condenser.m_flow_secondary_cp + else: + T_con_out = inputs.T_con_out + T_con_in = T_con_out - Q_con_outer / self.condenser.m_flow_secondary_cp # COP based on P_el and Q_con: COP_inner = Q_con / P_el @@ -290,6 +445,17 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): COP_carnot = (T_con_out / (T_con_out - inputs.T_eva_in)) carnot_quality = COP_inner / COP_carnot # Calc return temperature: + if inputs.uses_condenser_inlet: + fs_state.set( + name="T_con_out", value=T_con_out, unit="K", + description="Condenser outlet temperature" + ) + else: + fs_state.set( + name="T_con_in", value=T_con_in, unit="K", + description="Condenser inlet temperature" + ) + fs_state.set( name="P_el", value=P_el, unit="W", description="Power consumption" @@ -319,8 +485,28 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): name="COP_outer", value=COP_outer, unit="-", description="Outer COP, including heat losses" ) - + fs_state.set( + name="error_con", value=error_con, + unit="%", description="Error in condenser heat exchanger model" + ) + fs_state.set( + name="error_eva", value=error_eva, + unit="%", description="Error in evaporator heat exchanger model" + ) + fs_state.set( + name="dT_min_eva", value=dT_min_eva, + unit="K", description="Evaporator pinch temperature" + ) + fs_state.set( + name="dT_min_con", value=dT_min_con, + unit="K", description="Condenser pinch temperature" + ) + fs_state.set( + name="eta_glob", value=fs_state.get("eta_is").value * fs_state.get("eta_mech").value, + unit="%", description="Global compressor efficiency" + ) if save_path_plots is not None: + input_name = inputs.get_name() self.plot_cycle(save_path=save_path_plots.joinpath(f"{input_name}_final_result.png"), inputs=inputs) return fs_state @@ -395,13 +581,19 @@ def _plot_secondary_heat_flow_rates(self, ax, inputs): self.evaporator.state_outlet.h * self.evaporator.m_flow - Q_eva ]) / self.evaporator.m_flow self.condenser.m_flow_secondary = inputs.m_flow_con - self.condenser.calc_secondary_cp(T=inputs.T_con_in) + self.condenser.calc_secondary_cp(T=inputs.T_con) self.evaporator.m_flow_secondary = inputs.m_flow_eva self.evaporator.calc_secondary_cp(T=inputs.T_eva_in) - ax.plot(delta_H_con / 1000, [ - inputs.T_con_in - 273.15, - inputs.T_con_in + Q_con / self.condenser.m_flow_secondary_cp - 273.15 - ], color="b") + if inputs.uses_condenser_inlet: + ax.plot(delta_H_con / 1000, [ + inputs.T_con_in - 273.15, + inputs.T_con_in + Q_con / self.condenser.m_flow_secondary_cp - 273.15 + ], color="b") + else: + ax.plot(delta_H_con / 1000, [ + inputs.T_con_out - Q_con / self.condenser.m_flow_secondary_cp - 273.15, + inputs.T_con_out - 273.15 + ], color="b") ax.plot(delta_H_eva / 1000, [ inputs.T_eva_in - 273.15, inputs.T_eva_in - Q_eva / self.evaporator.m_flow_secondary_cp - 273.15 diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index 2885b01..1c0b068 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -19,6 +19,9 @@ def calc_multiple_states( save_path: pathlib.Path, heat_pump: BaseCycle, inputs: List[Inputs], + use_multiprocessing: bool = False, + raise_errors: bool = False, + with_unit_and_description: bool = True, **kwargs): """ Function to calculate the flowsheet states for all given inputs. @@ -28,42 +31,55 @@ def calc_multiple_states( save_path (pathlib.Path): Location where to save the results as xlsx. heat_pump (BaseCycle): A valid flowsheet inputs (List[Inputs]): A list with all inputs to simulate + use_multiprocessing (bool): True to use all cores, default no multiprocessing **kwargs: Solver settings for the flowsheet """ rel_infos = [] - for i, single_inputs in enumerate(inputs): - fs_state = None - logger.info(f"Running combination {i+1}/{len(inputs)}.") - try: - fs_state = heat_pump.calc_steady_state(inputs=single_inputs, - **kwargs) - except Exception as e: - # Avoid loss of data if un-excepted errors occur. - logger.error(f"An error occurred: {e}") - if fs_state is None: - fs_state = FlowsheetState() + fs_states = [] + i = 0 + if use_multiprocessing: + mp_inputs = [[heat_pump, inputs_, kwargs, raise_errors] for inputs_ in inputs] + pool = multiprocessing.Pool(processes=min(multiprocessing.cpu_count(), len(inputs))) + for fs_state in pool.imap(_calc_single_hp_state, mp_inputs): + fs_states.append(fs_state) + i += 1 + logger.info(f"Ran {i} of {len(inputs)} points") + else: + for inputs_ in inputs: + fs_state = _calc_single_hp_state([heat_pump, inputs_, kwargs, raise_errors]) + fs_states.append(fs_state) + i += 1 + logger.info(f"Ran {i} of {len(inputs)} points") + + for fs_state, single_inputs in zip(fs_states, inputs): hp_state_dic = { - **single_inputs.convert_to_str_value_format(with_unit_and_description=True), - **fs_state.convert_to_str_value_format(with_unit_and_description=True) + **single_inputs.convert_to_str_value_format(with_unit_and_description), + **fs_state.convert_to_str_value_format(with_unit_and_description) } rel_infos.append(hp_state_dic) df = pd.DataFrame(rel_infos) df.index.name = "State Number" - df.to_excel(save_path.joinpath(f"{heat_pump}_{heat_pump.fluid}.xlsx"), sheet_name="HP_states", float_format="%.5f") + if os.path.isdir(save_path): + save_path = save_path.joinpath(f"{heat_pump}_{heat_pump.fluid}.xlsx") + df.to_excel(save_path, sheet_name="HP_states", float_format="%.5f") def full_factorial_map_generation( heat_pump: BaseCycle, T_eva_in_ar: Union[list, np.ndarray], - T_con_in_ar: Union[list, np.ndarray], + T_con_ar: Union[list, np.ndarray], n_ar: Union[list, np.ndarray], m_flow_con: float, m_flow_eva: float, save_path: Union[pathlib.Path, str], dT_eva_superheating=5, dT_con_subcooling=0, + use_condenser_inlet: bool = True, + dT_con_start: float = 10, use_multiprocessing: bool = False, save_plots: bool = False, + raise_errors: bool = False, + save_sdf: bool = True, **kwargs ) -> (pathlib.Path, pathlib.Path): """ @@ -71,24 +87,41 @@ def full_factorial_map_generation( used in other simulation tools like Modelica or to analyze the off-design of the flowsheet. The results are stored and returned as .sdf and .csv files. - Currently, only varying T_eva_in, T_con_in, and n is implemented. + Currently, only varying T_eva_in, T_con_in (or T_con_out), and n is implemented. However, changing this to more dimensions or other variables is not much work. In this case, please raise an issue. Args: heat_pump (BaseCycle): The flowsheet to use - T_eva_in_ar: Array with inputs for T_eva_in - T_con_in_ar: Array with inputs for T_con_in - n_ar: Array with inputs for n_ar - m_flow_con: Condenser mass flow rate - m_flow_eva: Evaporator mass flow rate - save_path: Where to save all results. - dT_eva_superheating: Evaporator superheating - dT_con_subcooling: Condenser subcooling + T_eva_in_ar (list): + Array with inputs for T_eva_in + T_con_ar (list): + Array with inputs for T_con_in or T_con_out, see `use_condenser_inlet` + n_ar (list): + Array with inputs for n_ar + m_flow_con (float): + Condenser mass flow rate + m_flow_eva (float): + Evaporator mass flow rate + save_path (Path): + Where to save all results. + dT_eva_superheating (float): + Evaporator superheating + dT_con_subcooling (float): + Condenser subcooling + use_condenser_inlet (bool): + True to consider T_con_ar as inlet, false for outlet. + dT_con_start (float): + Guess value for starting condenser pressure, + only relevant if `use_condenser_inlet=False` use_multiprocessing: True to use multiprocessing. May speed up the calculation. Default is False - save_plots: + save_plots (bool): True to save plots of each steady state point. Default is False + raise_errors (bool): + True to raise errors if they occur. + save_sdf (bool): + = False to not save sdf files. Default is True **kwargs: Solver settings for the flowsheet Returns: @@ -106,34 +139,47 @@ def full_factorial_map_generation( idx_for_access_later = [] for i_T_eva_in, T_eva_in in enumerate(T_eva_in_ar): for i_n, n in enumerate(n_ar): - for i_T_con_in, T_con_in in enumerate(T_con_in_ar): - idx_for_access_later.append([i_n, i_T_con_in, i_T_eva_in]) - inputs = Inputs(n=n, - T_eva_in=T_eva_in, - T_con_in=T_con_in, - m_flow_eva=m_flow_eva, - m_flow_con=m_flow_con, - dT_eva_superheating=dT_eva_superheating, - dT_con_subcooling=dT_con_subcooling) - list_mp_inputs.append([heat_pump, inputs, kwargs]) + for i_T_con, T_con in enumerate(T_con_ar): + idx_for_access_later.append([i_n, i_T_con, i_T_eva_in]) + base_inputs = dict( + n=n, + T_eva_in=T_eva_in, + m_flow_eva=m_flow_eva, + m_flow_con=m_flow_con, + dT_eva_superheating=dT_eva_superheating, + dT_con_subcooling=dT_con_subcooling + ) + if use_condenser_inlet: + inputs = Inputs( + T_con_in=T_con, + **base_inputs + ) + else: + inputs = Inputs( + T_con_out=T_con, + **base_inputs + ) + inputs.set(name="dT_con_start", value=dT_con_start, unit="K") + list_mp_inputs.append([heat_pump, inputs, kwargs, raise_errors]) list_inputs.append(inputs) fs_states = [] i = 0 if use_multiprocessing: - pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + # pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + pool = multiprocessing.Pool(processes=10) for fs_state in pool.imap(_calc_single_hp_state, list_mp_inputs): fs_states.append(fs_state) i += 1 logger.info(f"Ran {i} of {len(list_mp_inputs)} points") else: for inputs in list_inputs: - fs_state = _calc_single_hp_state([heat_pump, inputs, kwargs]) + fs_state = _calc_single_hp_state([heat_pump, inputs, kwargs, raise_errors]) fs_states.append(fs_state) i += 1 logger.info(f"Ran {i} of {len(list_mp_inputs)} points") # Save to sdf - result_shape = (len(n_ar), len(T_con_in_ar), len(T_eva_in_ar)) + result_shape = (len(n_ar), len(T_con_ar), len(T_eva_in_ar)) _dummy = np.zeros(result_shape) # Use a copy to avoid overwriting of values of any sort. _dummy[:] = np.nan # Get all possible values: @@ -144,21 +190,26 @@ def full_factorial_map_generation( all_variables.update({var: _dummy.copy() for var in fs_state.get_variable_names()}) all_variables_info.update({var: variable for var, variable in fs_state.get_variables().items()}) variables_to_excel.append({ - **inputs.convert_to_str_value_format(with_unit_and_description=True), - **fs_state.convert_to_str_value_format(with_unit_and_description=True), + **inputs.convert_to_str_value_format(with_unit_and_description=False), + **fs_state.convert_to_str_value_format(with_unit_and_description=False), }) # Save to excel save_path_sdf = save_path.joinpath(f"{heat_pump.flowsheet_name}_{heat_pump.fluid}.sdf") save_path_csv = save_path.joinpath(f"{heat_pump.flowsheet_name}_{heat_pump.fluid}.csv") pd.DataFrame(variables_to_excel).to_csv( - save_path_csv + save_path_csv, sep=";" ) + # Terminate heat pump med-props: + heat_pump.terminate() + if not save_sdf: + return save_path_csv + for fs_state, idx_triple in zip(fs_states, idx_for_access_later): - i_n, i_T_con_in, i_T_eva_in = idx_triple + i_n, i_T_con, i_T_eva_in = idx_triple for variable_name, variable in fs_state.get_variables().items(): - all_variables[variable_name][i_n][i_T_con_in][i_T_eva_in] = variable.value + all_variables[variable_name][i_n][i_T_con][i_T_eva_in] = variable.value _nd_data = {} for variable, nd_data in all_variables.items(): @@ -171,13 +222,13 @@ def full_factorial_map_generation( _scale_values = { "n": n_ar, - "T_con_in": T_con_in_ar, + "T_con_in" if use_condenser_inlet else "T_con_out": T_con_ar, "T_eva_in": T_eva_in_ar } inputs: Inputs = list_inputs[0] _parameters = {} for name, variable in inputs.items(): - if name not in _scale_values: + if name not in list(_scale_values.keys()) + ["T_con_out", "T_con_in"]: _parameters[name] = { "data": variable.value, "unit": variable.unit, @@ -199,20 +250,19 @@ def full_factorial_map_generation( } utils.save_to_sdf(data=sdf_data, save_path=save_path_sdf) - # Terminate heat pump med-props: - heat_pump.terminate() - return save_path_sdf, save_path_csv def _calc_single_hp_state(data): """Helper function for a single state to enable multiprocessing""" - heat_pump, inputs, kwargs = data + heat_pump, inputs, kwargs, raise_errors = data fs_state = None try: fs_state = heat_pump.calc_steady_state(inputs=inputs, **kwargs) except Exception as e: + if raise_errors: + raise e logger.error(f"An error occurred for input: {inputs.__dict__}: {e}") if fs_state is None: fs_state = FlowsheetState() diff --git a/vclibpy/utils/nominal_design.py b/vclibpy/utils/nominal_design.py index 370751b..6a373ab 100644 --- a/vclibpy/utils/nominal_design.py +++ b/vclibpy/utils/nominal_design.py @@ -45,34 +45,27 @@ def nominal_hp_design( m_flow_con_start = kwargs.get("m_flow_con_start", 0.2) m_flow_eva_start = kwargs.get("m_flow_eva_start", 1) accuracy = kwargs.get("accuracy", 0.001) - calculate_m_flow = dT_eva is not None and dT_con is not None - if calculate_m_flow: + + # We have to iterate to match the m_flows to the Q_cons: + m_flow_eva_next = m_flow_eva_start + m_flow_con_next = m_flow_con_start + while True: + # Set values + m_flow_eva = m_flow_eva_next + m_flow_con = m_flow_con_next + inputs.set("m_flow_con", m_flow_con) + inputs.set("m_flow_eva", m_flow_eva) # Get nominal value: fs_state = heat_pump.calc_steady_state(fluid=fluid, inputs=inputs) if fs_state is None: raise ValueError("Given configuration is infeasible at nominal point.") - - else: - # We have to iterate to match the m_flows to the Q_cons: - m_flow_eva_next = m_flow_eva_start - m_flow_con_next = m_flow_con_start - while True: - # Set values - m_flow_eva = m_flow_eva_next - m_flow_con = m_flow_con_next - inputs.set("m_flow_con", m_flow_con) - inputs.set("m_flow_eva", m_flow_eva) - # Get nominal value: - fs_state = heat_pump.calc_steady_state(fluid=fluid, inputs=inputs) - if fs_state is None: - raise ValueError("Given configuration is infeasible at nominal point.") - cp_eva = heat_pump.evaporator._secondary_cp - cp_con = heat_pump.condenser._secondary_cp - m_flow_con_next = fs_state.get("Q_con") / (dT_con * cp_con) - m_flow_eva_next = (fs_state.get("Q_con") * (1 - 1 / fs_state.get("COP"))) / (dT_eva * cp_eva) - # Check convergence: - if abs(m_flow_eva_next - m_flow_eva) < accuracy and abs(m_flow_con-m_flow_con_next) < accuracy: - break + cp_eva = heat_pump.evaporator._secondary_cp + cp_con = heat_pump.condenser._secondary_cp + m_flow_con_next = fs_state.get("Q_con").value / (dT_con * cp_con) + m_flow_eva_next = (fs_state.get("Q_con").value * (1 - 1 / fs_state.get("COP").value)) / (dT_eva * cp_eva) + # Check convergence: + if abs(m_flow_eva_next - m_flow_eva) < accuracy and abs(m_flow_con-m_flow_con_next) < accuracy: + break nominal_design_info = { **inputs.convert_to_str_value_format(with_unit_and_description=False), diff --git a/vclibpy/utils/sdf_.py b/vclibpy/utils/sdf_.py index 1b14d70..470dcfc 100644 --- a/vclibpy/utils/sdf_.py +++ b/vclibpy/utils/sdf_.py @@ -158,6 +158,6 @@ def _unpack_nd_data(data): if __name__ == '__main__': sdf_to_csv( - filepath=pathlib.Path(r"D:\00_temp\calibration_jmo\Optihorst_3D_vclibpy.sdf"), - save_path=pathlib.Path(r"D:\04_git\vclibpy\tests\regression_data") + filepath=pathlib.Path(r"D:\VCLibMap_old.sdf"), + save_path=pathlib.Path(r"D:\00_temp") ) diff --git a/vclibpy/utils/ten_coefficient_compressor_reqression.py b/vclibpy/utils/ten_coefficient_compressor_reqression.py index a3b484e..3722e68 100644 --- a/vclibpy/utils/ten_coefficient_compressor_reqression.py +++ b/vclibpy/utils/ten_coefficient_compressor_reqression.py @@ -1,6 +1,7 @@ -from typing import List +from typing import List, Union import csv import pandas as pd +from pathlib import Path from vclibpy.components.compressors import TenCoefficientCompressor from vclibpy import media, Inputs @@ -8,72 +9,58 @@ try: from sklearn.linear_model import LinearRegression from sklearn.preprocessing import PolynomialFeatures - from xlsxwriter.workbook import Workbook except ImportError as err: - raise ImportError("You have to install xlsxwriter and " - "sklearn to use this regression tool") + raise ImportError("You have to install sklearn to use this regression tool") def create_regression_data( - variables: List[str], - T_con: List[float], T_eva: List[float], n: List[float], - T_sh: float, T_sc: float, n_max: int, V_h: float, fluid: str, datasheet: str, - capacity_definition: str, assumed_eta_mech: int, - folder_path: str): + T_con: List[float], + T_eva: List[float], + n: List[float], + compressor: TenCoefficientCompressor, + save_path: Union[Path, str], + fluid: str = None, +): """ Performs multidimensional linear regression to create compressor learning data. Args: - variables: (List[str]): - Variable names to create regressions for. - Options are: eta_s, lambda_h, and eta_mech T_con (List[float]): Condensing temperatures in K T_eva (List[float]): Evaporation temperatures in K n (List[float]): Compressor speeds from 0 to 1 - T_sh (float): Superheat temperature. - T_sc (float): Subcooling temperature. - n_max (int): Maximum compressor speed. - V_h (float): Compressor volume. - fluid (str): Type of refrigerant. - datasheet (str): Path to the modified datasheet. - capacity_definition (str): Definition of compressor capacity (e.g., "cooling"). - assumed_eta_mech (int): Assumed mechanical efficiency. - folder_path (str): Path to the folder where the newly created table will be saved. + save_path (str): Path to the folder where the newly created table will be saved. + compressor (TenCoefficientCompressor): + Instance to create regression for. + If med_prop is not set, the given fluid will be used. + fluid (str): Type of refrigerant. Only used if not already specified in compressor. Returns: List[float]: A list containing the regression parameters [P0, P1, ..., P9]. Raises: ValueError: If any specified variable column is not present in the DataFrame. - - Example: - >>> create_regression_data(11.1, 8.3, 120, 20.9e-6, "PROPANE", "path/to/datasheet.xlsx", - ... "cooling", 1, 6, 5, 5, 5, 303.15, 243.15, "path/to/save", 3) - [intercept, P1, P2, P3, P4, P5, P6, P7, P8, P9] """ - # create RefProp, fluid & compressor instance - med_prop = media.CoolProp(fluid) - - compressor = TenCoefficientCompressor( - N_max=n_max, V_h=V_h, T_sc=T_sc, T_sh=T_sh, - capacity_definition=capacity_definition, - assumed_eta_mech=assumed_eta_mech, - datasheet=datasheet - ) + if compressor.med_prop is None: + med_prop_class, med_prop_kwargs = media.get_global_med_prop_and_kwargs() + med_prop = med_prop_class(fluid_name=fluid, **med_prop_kwargs) + compressor.med_prop = med_prop + else: + med_prop = compressor.med_prop - compressor.med_prop = med_prop keywords = { - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumentric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } + variables = ["eta_is", "eta_mech", "lambda_h"] + tuples_for_cols = [("", "n")] for _variable in variables: for _n in n: tuples_for_cols.append((keywords[_variable], compressor.get_n_absolute(_n))) # tuples_for_cols: - # eta_s eta_s eta_s lambda_h lambda_h lambda_h eta_mech ... + # eta_is eta_is eta_is lambda_h lambda_h lambda_h eta_mech ... # n 30 60 90 30 60 90 30 ... cols = pd.MultiIndex.from_tuples(tuples_for_cols) final_df = pd.DataFrame( @@ -83,7 +70,7 @@ def create_regression_data( # final_df: column names are tuples (tuples_for_cols). # First column is filled with P1, P2, ... - # for-loop for multiple types(eta_s, eta_mech, etc) + # for-loop for multiple types(eta_is, eta_mech, etc) for m, _variable in enumerate(variables): for k, _n in enumerate(n): # for-loop for multiple rotation speeds T_eva_list = [] @@ -95,7 +82,7 @@ def create_regression_data( for j in range(len(T_con)): if T_eva[i] < T_con[j]: p1 = med_prop.calc_state("TQ", T_eva[i], 1).p - state_1 = med_prop.calc_state("PT", p1, (T_eva[i] + T_sh)) + state_1 = med_prop.calc_state("PT", p1, (T_eva[i] + compressor.T_sh)) compressor.state_inlet = state_1 p2 = med_prop.calc_state("TQ", T_con[j], 1).p # The enthalpy and entropy of the outlet @@ -107,7 +94,7 @@ def create_regression_data( T_con_list.append(T_con[j]) inputs = Inputs(n=_n) - if _variable == "eta_s": + if _variable == "eta_is": result_list.append(compressor.get_eta_isentropic( p_outlet=p2, inputs=inputs) ) @@ -123,23 +110,11 @@ def create_regression_data( ) final_df[cols[m * len(n) + k + 1]] = create_regression_parameters(df, _variable) - - # dataframes with a double column header can't be saved as a - # .xlsx yet, if index = False (NotImplementedError). - # .xlsx format is necessary, because TenCoefficientCompressor.get_parameter() - # expects a .xlsx as an input - # --> workaround: save the dataframe as a .csv, read it again and save it as a .xlsx - # TODO: Revert once this feature is in pandas. - final_df.to_csv(folder_path + r"\regressions.csv", index=False) - - workbook = Workbook(folder_path + r"\regressions.xlsx") - worksheet = workbook.add_worksheet() - with open(folder_path + r"\regressions.csv", 'rt', encoding='utf8') as f: - reader = csv.reader(f) - for r, row in enumerate(reader): - for c, col in enumerate(row): - worksheet.write(r, c, col) - workbook.close() + df_new = final_df.copy() + df_new.columns = df_new.columns.get_level_values(0) # Use only first level + df_new.loc[-1] = final_df.columns.get_level_values(1) + final_df = df_new.sort_index().reset_index(drop=True) + final_df.to_csv(save_path, index=False) def create_regression_parameters(df: pd.DataFrame, variable: str):