diff --git a/biosteam/facilities/hxn/_heat_exchanger_network.py b/biosteam/facilities/hxn/_heat_exchanger_network.py index 8e732a45..fbb3c447 100644 --- a/biosteam/facilities/hxn/_heat_exchanger_network.py +++ b/biosteam/facilities/hxn/_heat_exchanger_network.py @@ -10,9 +10,11 @@ @author: sarangbhagwat and yoelcp """ import biosteam as bst +import thermosteam as tmo import numpy as np from .hxn_synthesis import synthesize_network, StreamLifeCycle, plot_pinch_diagram from warnings import warn +import warnings __all__ = ('HeatExchangerNetwork',) @@ -225,7 +227,25 @@ def _cost(self): unit = i.unit if s_out: unit.ins[i.index] = s_out s_out = unit.outs[i.index] - self.HXN_sys = sys = bst.System(ID=None, path=all_units) + # Order the path by the rewired stream connections (and + # detect recycle loops) rather than using synthesis order: a + # hot-side exchanger is synthesized before the cold-side + # exchangers that feed it, and a single pass in synthesis + # order would leave it with stale inlets. HXprocess units are + # interaction units, which Network.from_units strips out (and + # disconnects) by default; keep them with interaction=False. + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter('always', RuntimeWarning) + network = tmo.Network.from_units(all_units, interaction=False) + for w in caught: + if 'network path could not be determined' in str(w.message): + warn('heat exchanger network path could not be fully ' + 'ordered from its stream connections; exchangers ' + 'fed by later ones in the path may be simulated ' + 'with stale inlets until convergence', RuntimeWarning) + else: + warn(w.message, w.category) + self.HXN_sys = sys = bst.System._from_network(None, network) sys.set_tolerance(method='fixedpoint', subsystems=True) original_purchase_costs = [hx.purchase_cost for hx in hxs] diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index 1041c205..c47dde15 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -368,16 +368,71 @@ def temperature_interval_pinch_analysis(hus, streams_quenched +def _end_state(stream_end, T_lo, T_hi): + """ + Return a copy of `stream_end` at equilibrium at its own enthalpy, or the + stream as given if that equilibrium state lies outside the stream's own + temperature range [T_lo, T_hi] (e.g. a non-condensable mislabelled as a + liquid, whose equilibrium state at the same enthalpy is a gas at an + absurd temperature). Either way the enthalpy is exactly `stream_end.H`. + """ + stream = stream_end.copy() + try: + stream.vle(H=stream_end.H, P=stream.P) + except Exception: + return stream_end.copy() + if T_lo <= stream.T <= T_hi: return stream + return stream_end.copy() + +def pinch_state(stream_in, stream_out, T_pinch): + """ + Return a copy of the stream in the state it has when it crosses the + pinch, with enthalpy guaranteed to lie within [min(H_in, H_out), + max(H_in, H_out)]. + + `stream_in` and `stream_out` are the stream's real end states (the + outlet quenched to equilibrium at its own enthalpy). For an interior + pinch the inlet copy is flashed at `T_pinch`; the result is used as is + when its enthalpy lies within the stream's own range. Otherwise the + stream never passes through that equilibrium state: a non-equilibrium + inlet (e.g. a superheated liquid from a non-rigorous HXutility) has + less enthalpy than the equilibrium fluid at the pinch, and the state + returned is instead the equilibrium state at the nearer end enthalpy. + The same end state is returned when `T_pinch` coincides with an end + temperature, because flashing a non-equilibrium inlet at its own + temperature does not reproduce `H_in` (and the result may even lie + inside the range). Using the *equilibrium* state at the end enthalpy, + rather than the stream as given, keeps the synthesizer consistent with + the problem table: the heat is offered at the temperature the + equilibrium model says it is available, not at a fictitious one; see + `_end_state` for the fallback when that state is unphysical. + + Either way the hot-side and cold-side loads split `|H_in - H_out|` + exactly and the transient stream used for matching never carries heat + the real stream does not have. This is the synthesizer's counterpart of + the clipping done by `_stream_H_at_boundaries` for the problem table. + """ + T_lo, T_hi = sorted((stream_in.T, stream_out.T)) + if T_pinch == stream_in.T: return _end_state(stream_in, T_lo, T_hi) + if T_pinch == stream_out.T: return _end_state(stream_out, T_lo, T_hi) + stream = stream_in.copy() + stream.vle(T=T_pinch, P=stream.P) + H_in = stream_in.H + H_out = stream_out.H + H_lo, H_hi = sorted((H_in, H_out)) + H = stream.H + if H_lo <= H <= H_hi: return stream + H_clipped = H_lo if H < H_lo else H_hi + return _end_state(stream_in if H_clipped == H_in else stream_out, T_lo, T_hi) + def load_duties(streams, streams_quenched, pinch_T_arr, T_out_arr, indices, is_cold, Q_hot_side, Q_cold_side): for index in indices: - stream = streams[index].copy() - H_in = stream.H - stream.vle(T = pinch_T_arr[index], P = stream.P) - H_pinch = stream.H + H_in = streams[index].H H_out = streams_quenched[index].H + H_pinch = pinch_state(streams[index], streams_quenched[index], pinch_T_arr[index]).H if not is_cold(index): - dH1 = abs(H_pinch - H_in) - dH2 = abs(H_out - H_pinch) + dH1 = H_in - H_pinch + dH2 = H_pinch - H_out if abs(dH1)<0.01: dH1 = 0 if abs(dH2)<0.01: dH2 = 0 Q_hot_side[index] = ['cool', dH1] @@ -418,14 +473,17 @@ def synthesize_network(hus, T_min_app=5., Qmin=1e-3, force_ideal_thermo=False, HXs_cold_side = [] streams_transient_cold_side = streams_inlet streams_transient_hot_side = [i.copy() for i in streams_inlet] + # Hot streams enter the cold-side design at their pinch state and cold + # streams enter the hot-side design at theirs; the enthalpy of that + # state is clipped to the stream's real range (see `pinch_state`). for i in hot_indices: s = streams_transient_cold_side[i] if s.T != pinch_T_arr[i]: - s.vle(T=pinch_T_arr[i], P=s.P) + streams_transient_cold_side[i] = pinch_state(s, streams_quenched[i], pinch_T_arr[i]) for i in cold_indices: s = streams_transient_hot_side[i] if s.T != pinch_T_arr[i]: - s.vle(T=pinch_T_arr[i], P=s.P) + streams_transient_hot_side[i] = pinch_state(s, streams_quenched[i], pinch_T_arr[i]) def get_stream_at_H_max(cold): s_cs = streams_transient_cold_side[cold] diff --git a/biosteam/units/design_tools/heat_transfer.py b/biosteam/units/design_tools/heat_transfer.py index 8ade3361..af53ebe9 100644 --- a/biosteam/units/design_tools/heat_transfer.py +++ b/biosteam/units/design_tools/heat_transfer.py @@ -77,7 +77,16 @@ def heat_exchange_to_condition(s_in, s_out, T=None, phase=None, else: if s_out.H < H_lim: s_out.H = H_lim else: + # At the bubble point: the most the stream can absorb (release) + # without leaving T is full vaporization (condensation). The + # enthalpy limit still applies; a limit short of the full phase + # change lands in the two-phase region, so solve it by VLE. s_out.phase = 'g' if heating else 'l' + if H_lim_given: + if heating: + if s_out.H > H_lim: s_out.vle(H=H_lim, P=s_out.P) + else: + if s_out.H < H_lim: s_out.vle(H=H_lim, P=s_out.P) else: s_out.vle(T=T, P=s_out.P) if H_lim_given: @@ -168,14 +177,14 @@ def counter_current_heat_exchange(s0_in, s1_in, s0_out, s1_out, if Q_hot_stream == Q_cold_stream == 0.: s0_out.copy_like(s0_in) - s1_in.copy_like(s1_out) + s1_out.copy_like(s1_in) return 0. if Q_hot_stream > 0 or Q_cold_stream < 0: # Sanity check if Q_hot_stream / s_hot_in.C < 0.1 or Q_cold_stream / s_cold_in.C > -0.1: s0_out.copy_like(s0_in) - s1_in.copy_like(s1_out) + s1_out.copy_like(s1_in) return 0. raise RuntimeError('inlet stream not in vapor-liquid equilibrium') diff --git a/tests/test_heat_exchange.py b/tests/test_heat_exchange.py new file mode 100644 index 00000000..5843498c --- /dev/null +++ b/tests/test_heat_exchange.py @@ -0,0 +1,82 @@ +# -*- coding: utf-8 -*- +# BioSTEAM: The Biorefinery Simulation and Techno-Economic Analysis Modules +# Copyright (C) 2020-, Yoel Cortes-Pena +# Copyright (C) 2026-, Sarang Bhagwat +# +# This module is under the UIUC open-source license. See +# github.com/BioSTEAMDevelopmentGroup/biosteam/blob/master/LICENSE.txt +# for license details. +""" +Tests for heat exchanger units and the counter-current heat exchange solver. +""" +import biosteam as bst +from numpy.testing import assert_allclose +from biosteam.units.design_tools.heat_transfer import heat_exchange_to_condition + +def test_heat_exchange_to_condition_respects_H_lim_at_bubble_point(): + # When the temperature limit coincides with the stream's bubble point + # (within the solver's 1e-3 K tolerance), the outlet is set to the + # saturated phase; the enthalpy limit must still be honored. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + s_in = bst.Stream('w_in', Water=2000., T=350., P=101325., phase='l', + units='kmol/hr') + T_bp = s_in.bubble_point_at_P().T + s_lim = s_in.copy(); s_lim.T = 360. + H_lim = s_lim.H + # heating: limit below full vaporization at the bubble point + s_out = s_in.copy() + Q = heat_exchange_to_condition(s_in, s_out, T=T_bp, H_lim=H_lim, heating=True) + assert_allclose(s_out.H, H_lim, rtol=1e-9) + assert_allclose(Q, H_lim - s_in.H, rtol=1e-9) + assert s_out.phase == 'l' and s_out.T < T_bp + # cooling: a saturated vapor whose limit is above full condensation + v_in = s_in.copy('w_vap'); v_in.phase = 'g'; v_in.T = T_bp + s_out = v_in.copy() + H_lim = v_in.H - 0.5 * (v_in.H - s_in.H) # half-way to liquid at 350 K + Q = heat_exchange_to_condition(v_in, s_out, T=T_bp, H_lim=H_lim, heating=False) + assert_allclose(s_out.H, H_lim, rtol=1e-9) + assert_allclose(Q, H_lim - v_in.H, rtol=1e-9) + +def test_HXprocess_H_lim_when_pinch_is_at_bubble_point(): + # Hot liquid water at exactly T_bp + dT against cold liquid water with an + # enthalpy limit below vaporization: the cold outlet may not exceed its + # enthalpy limit just because its temperature limit lands on the + # bubble point. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + cold = bst.Stream('cold', Water=2000., T=350., P=101325., phase='l', + units='kmol/hr') + T_bp = cold.bubble_point_at_P().T + hot = bst.Stream('hot', Water=1000., T=T_bp + 5., P=5e5, phase='l', + units='kmol/hr') + s_lim = cold.copy(); s_lim.T = 360. + H_lim = s_lim.H + hx = bst.HXprocess('hx', ins=(cold, hot), H_lim0=H_lim, T_lim1=355., dT=5.) + hx.simulate() + cold_out, hot_out = hx.outs + assert_allclose(cold_out.H, H_lim, rtol=1e-9) + assert_allclose(hx.Q, H_lim - cold.H, rtol=1e-9) + assert_allclose(hot_out.H, hot.H - hx.Q, rtol=1e-9) + assert hot_out.T > 355. # the hot stream was not the limiting side + +def test_HXprocess_never_modifies_inlets(): + # A superheated-liquid hot inlet (ethanol at 370 K, 1 atm; bp 351.4 K) + # against a two-phase cold inlet of the same fluid: the solver finds no + # feasible exchange. It must leave both inlet streams untouched rather + # than overwrite one with its (already modified) outlet. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + hot = bst.Stream('hot', Ethanol=700., T=370., P=101325., phase='l', + units='kmol/hr') + cold = bst.Stream('cold', Ethanol=400., T=348.6, P=101325., phase='l', + units='kmol/hr') + cold.vle(H=cold.H + 4e6, P=101325.) # two-phase at the bubble point + H_hot, H_cold = hot.H, cold.H + hx = bst.HXprocess('hx', ins=(cold, hot), dT=5.) + hx._run() + assert hx.Q == 0. + assert_allclose([hot.H, cold.H], [H_hot, H_cold], rtol=1e-12) + assert hot.phase == 'l' and hot.T == 370. + +if __name__ == '__main__': + test_heat_exchange_to_condition_respects_H_lim_at_bubble_point() + test_HXprocess_H_lim_when_pinch_is_at_bubble_point() + test_HXprocess_never_modifies_inlets() diff --git a/tests/test_hxn.py b/tests/test_hxn.py index 214372d5..4e7a147b 100644 --- a/tests/test_hxn.py +++ b/tests/test_hxn.py @@ -96,7 +96,7 @@ def test_energy_balance_error_contributions_ignored_none(): # --------------------------------------------------------------------------- from biosteam.facilities.hxn.hxn_synthesis import ( - temperature_interval_pinch_analysis, problem_table, + temperature_interval_pinch_analysis, problem_table, load_duties, pinch_state, ) def utility_hx(ID, T, P, phase, T_out, **flow): @@ -236,6 +236,87 @@ def test_problem_table_point_load_cannot_heat_above_itself(): assert table.pinch_T == 395. assert_allclose(table.cold_util_load, P_hot - (H_395 - H_392), rtol=1e-9) +def test_load_duties_non_equilibrium_inlet_conserves_energy(): + # A non-rigorous HXutility can carry a superheated liquid (ethanol at + # 370 K, 1 atm; bp 351.4 K). Flashing it at the 355 K pinch gives vapor + # with far more enthalpy than the inlet has; the pinch split must clip + # to the stream's real enthalpy range so that the hot-side and cold-side + # loads sum to the stream's duty (and are not phantom latent heat). + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + s_in = bst.Stream('ne_in', Ethanol=700., T=370., P=101325., phase='l', + units='kmol/hr') + s_out = bst.Stream('ne_out', Ethanol=700., T=310., P=101325., phase='l', + units='kmol/hr') + duty = s_in.H - s_out.H + flashed = s_in.copy(); flashed.vle(T=355., P=101325.) + assert flashed.H > s_in.H # the trap: equilibrium enthalpy above H_in + Q_hot_side, Q_cold_side = {}, {} + load_duties([s_in], [s_out], np.array([355.]), np.array([310.]), [0], + lambda i: False, Q_hot_side, Q_cold_side) + assert Q_hot_side[0][0] == Q_cold_side[0][0] == 'cool' + assert_allclose(Q_hot_side[0][1] + Q_cold_side[0][1], duty, rtol=1e-9) + # the stream cannot give any heat above the pinch: its whole enthalpy + # range lies below the equilibrium state at 355 K + assert_allclose(Q_hot_side[0][1], 0., atol=1e-6) + assert_allclose(Q_cold_side[0][1], duty, rtol=1e-9) + # and the transient state used for matching carries exactly H_in, at + # equilibrium (two-phase at the bubble point), not at the fictitious + # 370 K of the superheated liquid, so matching is consistent with the + # problem table + s = pinch_state(s_in, s_out, 355.) + assert_allclose(s.H, s_in.H, rtol=1e-9) + assert s.T < 355. and 'g' in s.phase and 'l' in s.phase + +def test_pinch_state_at_endpoints_uses_real_states(): + # pinch_T == T_in must put the whole duty on one side even when the + # inlet is not at equilibrium: flashing a superheated liquid at its own + # T does not reproduce H_in, and the resulting enthalpy can lie inside + # the stream's range so that clipping alone would not catch it. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + s_in = bst.Stream('sh_in', Water=100., T=380., P=101325., phase='l', + units='kmol/hr') + s_out = bst.Stream('sh_out', Water=100., T=400., P=101325., phase='g', + units='kmol/hr') + duty = s_out.H - s_in.H + flashed = s_in.copy(); flashed.vle(T=380., P=101325.) + assert s_in.H < flashed.H < s_out.H # the trap: in range, but not H_in + Q_hot_side, Q_cold_side = {}, {} + load_duties([s_in], [s_out], np.array([380.]), np.array([400.]), [0], + lambda i: True, Q_hot_side, Q_cold_side) + assert_allclose(Q_hot_side[0][1], duty, rtol=1e-9) + assert_allclose(Q_cold_side[0][1], 0., atol=1e-6) + assert_allclose(pinch_state(s_in, s_out, 380.).H, s_in.H, rtol=1e-9) + assert_allclose(pinch_state(s_in, s_out, 400.).H, s_out.H, rtol=1e-9) + # a phase-mislabelled non-condensable: clipping must return the real + # end state, never an equilibrium re-flash (N2 "liquid" at 400 K would + # otherwise come back as gas at thousands of K) + bst.settings.set_thermo(['Water', 'N2'], cache=True) + h_in = bst.Stream('n2_in', N2=100., T=400., P=101325., phase='l', + units='kmol/hr') + h_out = bst.Stream('n2_out', N2=100., T=320., P=101325., phase='l', + units='kmol/hr') + s = pinch_state(h_in, h_out, 355.) + assert h_out.H <= s.H <= h_in.H + assert 320. <= s.T <= 400. + +def test_unordered_network_path_warns_with_context(monkeypatch): + # thermosteam's Network.sort warns 'network path could not be determined' + # when its ordering heuristic does not settle; HXN re-raises that as its + # own warning so the user knows which facility it concerns. + import thermosteam as tmo + original = tmo.Network.from_units + def from_units_with_warning(*args, **kwargs): + network = original(*args, **kwargs) + warnings.warn('network path could not be determined', RuntimeWarning) + return network + monkeypatch.setattr(tmo.Network, 'from_units', from_units_with_warning) + units = synthetic_units() + HXN = bst.HeatExchangerNetwork('HXN', T_min_app=5.) + sys = bst.System.from_units('sys_unordered', units=[*units, HXN]) + with pytest.warns(RuntimeWarning, match='heat exchanger network path could not be fully ordered'): + sys.simulate() + assert abs(HXN.energy_balance_percent_error) < 1e-6 + def test_synthetic_network_reaches_MER(): units = synthetic_units() HXN = bst.HeatExchangerNetwork('HXN', T_min_app=5.) diff --git a/tests/test_hxn_regression.py b/tests/test_hxn_regression.py new file mode 100644 index 00000000..7c2d1864 --- /dev/null +++ b/tests/test_hxn_regression.py @@ -0,0 +1,228 @@ +# -*- coding: utf-8 -*- +# BioSTEAM: The Biorefinery Simulation and Techno-Economic Analysis Modules +# Copyright (C) 2020-, Yoel Cortes-Pena +# Copyright (C) 2026-, Sarang Bhagwat +# +# This module is under the UIUC open-source license. See +# github.com/BioSTEAMDevelopmentGroup/biosteam/blob/master/LICENSE.txt +# for license details. +""" +Regression tests for heat exchanger network synthesis on synthetic systems. + +Ten synthetic systems of increasing complexity (all with phase-changing +streams from case 3 on). For each, the synthesized network must + +(i) close its energy balance (|error| < 1e-6 %) without RuntimeWarnings, +(ii) never beat the minimum-energy-requirement (MER) targets of the problem + table computed on the same streams, and +(iii) recover at least as much heat as documented in ``CASES`` below, so that + no future change to ``hxn`` silently makes the synthesizer perform worse. + +The documented utility loads were recorded by running this file directly +(``python tests/test_hxn_regression.py`` prints them): cases 1-4 and 6-9 +at commit ``1ab689ff`` (branch ``hxn-pinch-diagram``); case 10 after the +synthesizer fixes on ``hxn-regression-tests`` (non-equilibrium inlets +clipped to the stream's enthalpy range; network path ordered by its +connections; H_lim honored at the bubble point); case 5 lowered after the +pinch state at an end temperature became the equilibrium state at that +end enthalpy (its 420 K, 5 bar vapor feed is below water's boiling point +there, a non-equilibrium inlet). Improvements leave slack; a +maintainer lowers the numbers deliberately when a better network is +intended. Never raise them to make a failing test pass. +""" +import warnings +import pytest +import biosteam as bst +from numpy.testing import assert_allclose +from biosteam.facilities.hxn.hxn_synthesis import problem_table + +EB_TOLERANCE = 1e-6 # percent; converged networks close to ~1e-10 % +MER_RTOL = 1e-3 # network may not beat the MER target by more than this +DOC_RTOL = 1e-3 # network may not be worse than documented by more than this + +def utility_hx(ID, T, P, phase, T_out, rigorous=None, **flow): + """A simulated HXutility acting as one process stream (kmol/hr flows).""" + s = bst.Stream(ID + '_in', T=T, P=P, phase=phase, units='kmol/hr', **flow) + if rigorous is None: rigorous = phase == 'g' + hx = bst.HXutility(ID, ins=s, T=T_out, rigorous=rigorous) + hx.simulate() + return hx + +def boiling_hx(ID, T, P, T_out, **flow): + """A cold liquid stream heated past its bubble point (rigorous VLE).""" + return utility_hx(ID, T, P, 'l', T_out, rigorous=True, **flow) + +def setup(name, chemicals=('Water', 'Ethanol')): + bst.settings.set_thermo(list(chemicals), cache=True) + bst.main_flowsheet.set_flowsheet('test_hxn_regression_' + name) + +# --------------------------------------------------------------------------- +# Cases +# --------------------------------------------------------------------------- + +def case_01_two_liquids(): + """Hot liquid, cold liquid; trivial counter-current match.""" + setup('01') + return [utility_hx('H1', 400., 5e5, 'l', 320., Water=1000.), + utility_hx('C1', 300., 101325., 'l', 360., Water=1000.)], 5. + +def case_02_pinch_limited(): + """Cold target above the hot inlet: part of the heating must be utility.""" + setup('02') + return [utility_hx('H1', 360., 5e5, 'l', 320., Water=1000.), + utility_hx('C1', 300., 5e5, 'l', 380., Water=800.)], 5. + +def case_03_condenser_two_colds(): + """Condensing ethanol vapor against two cold liquids.""" + setup('03') + return [utility_hx('H1', 355., 101325., 'g', 340., Ethanol=400.), + utility_hx('C1', 300., 101325., 'l', 345., Water=1500.), + utility_hx('C2', 310., 101325., 'l', 340., Water=300., Ethanol=300.)], 5. + +def case_04_report_case(): + """The 4-stream report case: 2 colds, condensing hot, hot liquid.""" + setup('04') + return [utility_hx('C1', 300., 101325., 'l', 390., Water=2000.), + utility_hx('C2', 310., 101325., 'l', 345., Water=500., Ethanol=500.), + utility_hx('H1', 352., 101325., 'g', 340., Ethanol=300.), + utility_hx('H2', 420., 5e5, 'l', 320., Water=800.)], 5. + +def case_05_boiling_cold(): + """A cold stream that boils (water -> steam) and a condensing hot stream.""" + setup('05') + return [boiling_hx('C1', 330., 101325., 380., Water=300.), + utility_hx('C2', 300., 101325., 'l', 350., Water=1000.), + utility_hx('H1', 420., 5e5, 'g', 330., Water=250.), + utility_hx('H2', 400., 5e5, 'l', 310., Water=1500.)], 5. + +def case_06_mixed_pressures(): + """5-bar condensing hot against 1-atm boiling cold; T_min_app = 10.""" + setup('06') + return [utility_hx('H1', 430., 5e5, 'g', 400., Water=200.), + utility_hx('H2', 380., 5e5, 'l', 320., Water=1200.), + boiling_hx('C1', 340., 101325., 375., Water=150.), + utility_hx('C2', 300., 101325., 'l', 360., Ethanol=800.), + utility_hx('C3', 320., 101325., 'l', 390., Water=600.)], 10. + +def case_07_threshold(): + """Threshold problem: heating dominates so the cold target is ~0; + one hot stream condenses and subcools.""" + setup('07') + return [utility_hx('H1', 375., 101325., 'g', 310., Water=80.), + utility_hx('H2', 360., 101325., 'l', 330., Water=200.), + utility_hx('C1', 300., 101325., 'l', 370., Water=1500.), + utility_hx('C2', 305., 101325., 'l', 350., Ethanol=800.), + boiling_hx('C3', 340., 101325., 355., Ethanol=200.), + utility_hx('C4', 320., 101325., 'l', 360., Water=700.)], 5. + +def case_08_two_condensers(): + """Two condensers at different temperatures (ethanol 1 atm, water 2 bar) + against three colds, one of which boils.""" + setup('08') + return [utility_hx('H1', 355., 101325., 'g', 335., Ethanol=300.), + utility_hx('H2', 400., 2e5, 'g', 360., Water=150.), + utility_hx('H3', 390., 5e5, 'l', 330., Water=900.), + utility_hx('C1', 300., 101325., 'l', 345., Water=1200.), + boiling_hx('C2', 330., 101325., 370., Ethanol=250.), + utility_hx('C3', 310., 101325., 'l', 380., Water=700.)], 5. + +def case_09_near_degenerate(): + """Eight streams with two near-equal pinch candidates, a partially + condensing hot stream (wet outlet) and a partially boiling cold stream.""" + setup('09') + return [utility_hx('H1', 380., 101325., 'g', 372., Water=120.), # partial condensation + utility_hx('H2', 365., 101325., 'g', 330., Ethanol=250.), + utility_hx('H3', 410., 5e5, 'l', 340., Water=700.), + utility_hx('H4', 345., 101325., 'l', 305., Ethanol=900.), + boiling_hx('C1', 350., 101325., 373.5, Water=200.), # partial boiling + utility_hx('C2', 300., 101325., 'l', 340., Water=1500.), + boiling_hx('C3', 320., 101325., 352., Ethanol=300.), + utility_hx('C4', 335., 101325., 'l', 395., Water=500.)], 5. + +def case_10_ten_streams(): + """Ten streams mixing liquids, condensers, boilers, and pressures.""" + setup('10') + return [utility_hx('H1', 355., 101325., 'g', 320., Ethanol=300.), + utility_hx('H2', 420., 5e5, 'g', 340., Water=150.), + utility_hx('H3', 395., 5e5, 'l', 330., Water=1000.), + utility_hx('H4', 370., 101325., 'l', 310., Ethanol=700.), + utility_hx('H5', 380., 101325., 'g', 372.5, Water=100.), # partial condensation + utility_hx('C1', 300., 101325., 'l', 360., Water=2000.), + boiling_hx('C2', 330., 101325., 380., Water=200.), + boiling_hx('C3', 320., 101325., 352., Ethanol=400.), + utility_hx('C4', 310., 101325., 'l', 345., Water=400., Ethanol=400.), + utility_hx('C5', 340., 101325., 'l', 390., Water=600.)], 5. + +# name -> (builder, documented hot utility load [kJ/hr], documented cold utility load [kJ/hr]) +# Documented values: see module docstring for provenance. +CASES = { + 'case_01_two_liquids': (case_01_two_liquids, 0, 1.53912e+06), + 'case_02_pinch_limited': (case_02_pinch_limited, 1.81522e+06, 0), + 'case_03_condenser_two_colds': (case_03_condenser_two_colds, 0, 9.49905e+06), + 'case_04_report_case': (case_04_report_case, 2.37319e+06, 3.56871e+06), + 'case_05_boiling_cold': (case_05_boiling_cold, 3.2224e+06, 7.81466e+06), + 'case_06_mixed_pressures': (case_06_mixed_pressures, 7.05541e+06, 4.93079e+06), + 'case_07_threshold': (case_07_threshold, 1.85977e+07, 0), + 'case_08_two_condensers': (case_08_two_condensers, 3.02237e+06, 7.36431e+06), + 'case_09_near_degenerate': (case_09_near_degenerate, 1.40965e+07, 9.66427e+06), + 'case_10_ten_streams': (case_10_ten_streams, 1.40742e+07, 8.06488e+06), +} + +# --------------------------------------------------------------------------- +# Harness +# --------------------------------------------------------------------------- + +def synthesize(builder): + units, T_min_app = builder() + HXN = bst.HeatExchangerNetwork('HXN', T_min_app=T_min_app) + sys = bst.System.from_units('sys', units=[*units, HXN]) + with warnings.catch_warnings(): + warnings.simplefilter('error', RuntimeWarning) + # thermosteam registry bookkeeping on temporary stream copies; not numerical + warnings.filterwarnings('ignore', message='.*has been replaced in registry', + category=RuntimeWarning) + sys.simulate() + return units, HXN, T_min_app + +def mer_targets(units, T_min_app): + hus = [hx.heat_utilities[0] for hx in units] + hus.sort(key=lambda hu: hu.duty) + streams_inlet = [hu.unit.ins[0].copy() for hu in hus] + streams_quenched = [hu.unit.outs[0].copy() for hu in hus] + for s in streams_quenched: s.vle(H=s.H, P=s.P) + is_hot = [hu.duty < 0 for hu in hus] + table = problem_table(streams_inlet, streams_quenched, is_hot, T_min_app) + return table.hot_util_load, table.cold_util_load + +def actual_loads(HXN): + hus = [hu for hx in HXN.new_HX_utils for hu in hx.heat_utilities] + heat = sum(hu.unit_duty for hu in hus if hu.unit_duty > 0) + cool = -sum(hu.unit_duty for hu in hus if hu.unit_duty < 0) + return heat, cool + +@pytest.mark.parametrize('name', list(CASES)) +def test_hxn_regression(name): + builder, doc_heat, doc_cool = CASES[name] + units, HXN, T_min_app = synthesize(builder) + # (i) energy balance + assert abs(HXN.energy_balance_percent_error) < EB_TOLERANCE, name + # (ii) MER targets are a lower bound; energy identity holds + heat, cool = actual_loads(HXN) + hot_target, cold_target = mer_targets(units, T_min_app) + net_duty = sum(hx.heat_utilities[0].unit_duty for hx in units) + assert heat >= hot_target * (1 - MER_RTOL), (name, heat, hot_target) + assert cool >= cold_target * (1 - MER_RTOL), (name, cool, cold_target) + assert_allclose(heat - cool, net_duty, rtol=1e-8, err_msg=name) + # (iii) never worse than documented + assert doc_heat is not None and doc_cool is not None, f'{name}: baseline not recorded' + assert heat <= doc_heat * (1 + DOC_RTOL) + 1e-9, (name, heat, doc_heat) + assert cool <= doc_cool * (1 + DOC_RTOL) + 1e-9, (name, cool, doc_cool) + +if __name__ == '__main__': + for name, (builder, *_) in CASES.items(): + units, HXN, T_min_app = synthesize(builder) + heat, cool = actual_loads(HXN) + hot_target, cold_target = mer_targets(units, T_min_app) + print(f"{name}: heat={heat:.6g} cool={cool:.6g} " + f"(MER hot={hot_target:.6g} cold={cold_target:.6g}; " + f"EB error={HXN.energy_balance_percent_error:.4f}%)")