Add synthetic HXN regression tests; fix synthesis for non-equilibrium inlets, stale network ordering, and H_lim at the bubble point

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]

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]

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')
 

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 <yoelcortes@gmail.com>
+# Copyright (C) 2026-, Sarang Bhagwat <sarangbhagwat.developer@gmail.com>
+#
+# 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()

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.)