diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index 5d5ff658..216485d8 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -10,11 +10,13 @@ @author: sarangbhagwat """ +from collections import namedtuple import numpy as np import biosteam as bst from warnings import warn -__all__ = ('StreamLifeCycle', 'synthesize_network') +__all__ = ('StreamLifeCycle', 'ProblemTable', 'problem_table', + 'synthesize_network') class LifeStage: @@ -124,8 +126,179 @@ def get_sorted_life_cycle(self): return self.life_cycle -def temperature_interval_pinch_analysis(hus, - T_min_app=10, +ProblemTable = namedtuple( + 'ProblemTable', + ['Ts', 'interval_H', 'point_H', 'residual', + 'hot_util_load', 'cold_util_load', 'pinch_T'] +) + +def _stream_H_at_boundaries(stream_in, H_in, H_out, T_lo, T_hi, Ts, shift, + stream_label): + """ + Enthalpies [kJ/hr] of one monotone stream at the grid boundaries + `Ts` (shifted scale, descending, all within [T_lo, T_hi]). + + Exact at the stream's own end points (H_in/H_out as given) by + *position*: `Ts[0]` and `Ts[-1]` are the stream's own T_hi/T_lo (every + monotone stream has `T_hi > T_lo` strictly, so `Ts` always has at least + these two entries) and are assigned H_in/H_out directly, without a float + comparison. In between, the inlet copy is flashed at the *real* + temperature `T + shift` and the result is clipped to + [min(H_in, H_out), max(H_in, H_out)] so that a non-equilibrium outlet + (e.g. a column reboiler/condenser product) can never inflate an + interval. A single copy is walked down the grid so each VLE is + warm-started from the previous boundary; `stream_label` (the inlet + stream's own ID) identifies the stream in the VLE-failure warning. + """ + assert Ts.size >= 2, ( + "boundary grid for a monotone stream must include both its own " + "end points" + ) + H_lo, H_hi = sorted((H_in, H_out)) + H_top, H_bottom = (H_in, H_out) if H_in > H_out else (H_out, H_in) + Hs = np.empty(Ts.size) + Hs[0] = H_top + Hs[-1] = H_bottom + stream = stream_in.copy() + for k in range(1, Ts.size - 1): + T = Ts[k] + T_real = T + shift + try: + stream.vle(T=T_real, P=stream.P) + H = stream.H + except Exception as error: + warn(f"could not solve VLE for stream {stream_label!r} at " + f"{T_real:.2f} K ({error!r}); interpolating enthalpy " + "linearly in temperature for the problem table", + RuntimeWarning) + # restart the warm start from a clean copy so the failed flash + # does not leave `stream` in a bad state for the next boundary + stream = stream_in.copy() + H = H_lo + (H_hi - H_lo) * (T - T_lo) / (T_hi - T_lo) + Hs[k] = min(max(H, H_lo), H_hi) + return Hs + +def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): + """ + Energy-consistent problem table (temperature-interval heat cascade). + + Parameters + ---------- + streams_inlet : list[Stream] + Inlet stream of each utility heat exchanger. + streams_quenched : list[Stream] + Corresponding outlet streams, re-flashed at their enthalpy. + is_hot : Sequence[bool] + True where the stream is cooled. + T_min_app : float + Minimum approach temperature [K]. + + Returns + ------- + ProblemTable + Grid temperatures `Ts` (shifted scale, descending), per-stream + `interval_H` (N x n-1) and `point_H` (N x n) contributions (+ for + hot, - for cold), the cascade `residual` (n) *leaving* each + boundary (i.e. after its point loads), `hot_util_load`, + `cold_util_load` and the shifted-scale `pinch_T`. + + Notes + ----- + Hot streams are shifted down by `T_min_app`; cold streams are not. For + monotone streams the contribution to interval (Ts[k], Ts[k+1]) is + sign * (H(Ts[k]) - H(Ts[k+1])) with H evaluated at the real temperature + and clipped to [H_in, H_out], so every stream's contributions telescope + exactly to sign * |H_out - H_in|. Isothermal streams, and streams whose + outlet temperature moves against their duty (a heated stream that exits + colder than it entered, e.g. a reboiler outlet at VLE), are point loads + at their outlet temperature. The cascade starting from zero hot utility + is residual[k] = sum(point_H[:, :k+1]) + sum(interval_H[:, :k]), the + heat *leaving* boundary Ts[k]. Feasibility must also hold for the heat + *arriving* at Ts[k] before its point loads are applied, + arriving[k] = residual[k] - sum(point_H[:, k]), because a source at + Ts[k] cannot serve a sink above Ts[k]. The minimum over both flows, + min(residual, arriving), fixes the hot utility target, + `residual[-1] + hot_util_load` the cold one, and its location the + pinch. With the per-stream identity above, hot_util_load - + cold_util_load equals the net heating demand. + + Examples + -------- + A threshold problem: 1000 kmol/hr of water cooled 400 -> 300 K supplies + every interval of 900 kmol/hr of water heated 300 -> 390 K, so no hot + utility is needed and the surplus leaves as cold utility. + + >>> import biosteam as bst + >>> from biosteam.facilities.hxn.hxn_synthesis import problem_table + >>> bst.settings.set_thermo(['Water']) + >>> hot_in = bst.Stream(Water=1000., T=400., P=5e5, phase='l', units='kmol/hr') + >>> hot_out = hot_in.copy(); hot_out.vle(T=300., P=5e5) + >>> cold_in = bst.Stream(Water=900., T=300., P=5e5, phase='l', units='kmol/hr') + >>> cold_out = cold_in.copy(); cold_out.vle(T=390., P=5e5) + >>> table = problem_table([hot_in, cold_in], [hot_out, cold_out], + ... [True, False], 5.) + >>> table.Ts + array([395., 390., 300., 295.]) + >>> round(table.hot_util_load, 3) + 0.0 + >>> round(table.cold_util_load, -1) + 1445550.0 + >>> table.pinch_T + 395.0 + """ + N = len(streams_inlet) + is_hot = np.asarray(is_hot, dtype=bool) + sign = np.where(is_hot, 1., -1.) + shift = np.where(is_hot, T_min_app, 0.) + T_in = np.array([s.T for s in streams_inlet]) + T_out = np.array([s.T for s in streams_quenched]) + H_in = np.array([s.H for s in streams_inlet]) + H_out = np.array([s.H for s in streams_quenched]) + monotone = (sign * (T_in - T_out)) > 0. + T_hi = np.where(monotone, np.maximum(T_in, T_out), T_out) - shift + T_lo = np.where(monotone, np.minimum(T_in, T_out), T_out) - shift + Ts = np.unique(np.concatenate([T_hi, T_lo]))[::-1] + n = Ts.size + interval_H = np.zeros((N, n - 1)) + point_H = np.zeros((N, n)) + for j in range(N): + if monotone[j]: + idx = np.flatnonzero((Ts <= T_hi[j]) & (Ts >= T_lo[j])) + Hs = _stream_H_at_boundaries(streams_inlet[j], H_in[j], H_out[j], + T_lo[j], T_hi[j], Ts[idx], shift[j], + streams_inlet[j].ID) + interval_H[j, idx[:-1]] = sign[j] * (Hs[:-1] - Hs[1:]) + else: + k = np.searchsorted(-Ts, -T_hi[j]) + point_H[j, k] = sign[j] * abs(H_out[j] - H_in[j]) + point_total = point_H.sum(axis=0) + residual = np.cumsum( + point_total + np.concatenate([[0.], interval_H.sum(axis=0)]) + ) + # heat arriving at each boundary, before that boundary's point loads: + # a point source at Ts[k] cannot serve sinks above Ts[k], so the cascade + # must be non-negative both before and after the point loads + arriving = residual - point_total + flow = np.minimum(residual, arriving) + k_pinch = int(np.argmin(flow)) + scale = np.abs(H_out - H_in).sum() + if -flow[k_pinch] <= 1e-9 * scale: # threshold problem: no hot utility + hot_util_load = 0. + k_pinch = 0 + else: + hot_util_load = -flow[k_pinch] + cold_util_load = residual[-1] + hot_util_load + if cold_util_load < 0.: + # only reachable in the threshold branch, by at most 1e-9 * scale: + # absorb the rounding into the hot utility so that + # hot_util_load - cold_util_load == sum(unit_duty) stays exact + hot_util_load -= cold_util_load + cold_util_load = 0. + return ProblemTable(Ts, interval_H, point_H, residual, + hot_util_load, cold_util_load, Ts[k_pinch]) + +def temperature_interval_pinch_analysis(hus, + T_min_app=10, force_ideal_thermo=False, sort_hus_by_T=False): hx_utils = hus @@ -150,77 +323,45 @@ def temperature_interval_pinch_analysis(hus, ID = 'Util_%s'%i stream.ID = 's_%s__%s'%(i,ID) N_heating = len(hus_heating) - is_cold_stream_index = lambda x: x < N_heating - T_in_arr = np.array([stream.T for stream in streams_inlet]) - T_out_arr = np.array([i.T for i in streams_quenched]) - adj_T_in_arr = T_in_arr.copy() - # adj_T_in_arr[:N_heating] -= T_min_app - adj_T_in_arr[N_heating:] -= T_min_app - adj_T_out_arr = T_out_arr.copy() - # adj_T_out_arr[:N_heating] -= T_min_app - adj_T_out_arr[N_heating:] -= T_min_app - T_changes_tuples = list(zip(adj_T_in_arr, adj_T_out_arr)) - all_Ts_descending = [*adj_T_in_arr, *adj_T_out_arr] - all_Ts_descending.sort(reverse=True) - stream_indices_for_T_intervals =\ - {(all_Ts_descending[i], all_Ts_descending[i+1]):[]\ - for i in range(len(all_Ts_descending)-1)} - H_for_T_intervals = dict.fromkeys(stream_indices_for_T_intervals, 0) cold_indices = list(range(N_heating)) hot_indices = list(range(N_heating, len(hxs))) indices = cold_indices + hot_indices - for i in range(len(all_Ts_descending)-1): - T_start = all_Ts_descending[i] - T_end = all_Ts_descending[i+1] - for stream_index in indices: - T1, T2 = T_changes_tuples[stream_index] - if (T1 >= T_start and T2 <= T_end) or (T2 >= T_start and T1 <= T_end): - multiplier = -1 if is_cold_stream_index(stream_index) else 1 - stream = streams_inlet[stream_index].copy() - if stream.T != T_start: stream.vle(T = T_start, P = stream.P) - H1 = stream.H - try: - stream.vle(T = T_end, P = stream.P) - except: - warn(f"could not solve VLE for {repr(stream)} at {repr(hxs[stream_index].owner)}", RuntimeWarning) - H2 = stream.H - H = multiplier*(H1 - H2) - H_for_T_intervals[(T_start, T_end)] += H - - res_H_vector = [] - prev_res_H = 0 - for interval, H in H_for_T_intervals.items(): - res_H_vector.append(prev_res_H + H) - prev_res_H = res_H_vector[len(res_H_vector)-1] - hot_util_load = - min(res_H_vector) - # assert hot_util_load>= 0, 'Hot utility load is negative' - if not hot_util_load>=0: - warn(f"Hot utility load is negative: {hot_util_load}", RuntimeWarning) - # print(hot_util_load) - # the lower temperature of the temperature interval for which the res_H is minimum - pinch_cold_stream_T = all_Ts_descending[res_H_vector.index(-hot_util_load)+1] + T_in_arr = np.array([stream.T for stream in streams_inlet]) + T_out_arr = np.array([i.T for i in streams_quenched]) + is_hot = np.zeros(len(hxs), dtype=bool) + is_hot[hot_indices] = True + table = problem_table(streams_inlet, streams_quenched, is_hot, T_min_app) + hot_util_load = table.hot_util_load + cold_util_load = table.cold_util_load + pinch_cold_stream_T = table.pinch_T pinch_hot_stream_T = pinch_cold_stream_T + T_min_app - cold_util_load = res_H_vector[len(res_H_vector)-1] + hot_util_load - # assert cold_util_load>=0, 'Cold utility load is negative' - if not cold_util_load>=0: - warn(f"Cold utility load is positive: {cold_util_load}", RuntimeWarning) + # Per-stream pinch temperature: where each stream is split between the + # hot-side and cold-side network designs. A stream already entirely on + # one side of the process pinch (T_in past pinch_cold_stream_T for a + # cold stream, or past pinch_hot_stream_T for a hot stream) is not + # split; its pinch_T is its own T_in. This clause also catches + # non-monotone streams (T_out on the wrong side of T_in for their duty, + # e.g. a cold stream whose VLE outlet ends up cooler than it entered): + # rather than split their problem_table point-load duty across the + # cascade, they get pinch_T = T_in too, so load_duties assigns their + # whole duty to a single side (Q_hot_side for a cold stream, + # Q_cold_side for a hot one). pinch_T_arr = [] for i in cold_indices: - if T_in_arr[i] > pinch_cold_stream_T: + if T_in_arr[i] > pinch_cold_stream_T or T_in_arr[i] > T_out_arr[i]: pinch_T_arr.append(T_in_arr[i]) elif T_out_arr[i] < pinch_cold_stream_T: pinch_T_arr.append(T_out_arr[i]) else: pinch_T_arr.append(pinch_cold_stream_T) for i in hot_indices: - if T_in_arr[i] < pinch_hot_stream_T: + if T_in_arr[i] < pinch_hot_stream_T or T_in_arr[i] < T_out_arr[i]: pinch_T_arr.append(T_in_arr[i]) elif T_out_arr[i] > pinch_hot_stream_T: pinch_T_arr.append(T_out_arr[i]) else: pinch_T_arr.append(pinch_hot_stream_T) pinch_T_arr = np.array(pinch_T_arr) - # print(pinch_T_arr, hot_util_load, cold_util_load,) return pinch_T_arr, hot_util_load, cold_util_load, T_in_arr, T_out_arr,\ hxs, hot_indices, cold_indices, indices, streams_inlet, hx_utils_rearranged, \ streams_quenched diff --git a/tests/test_hxn.py b/tests/test_hxn.py index 758a22e9..5520b6fe 100644 --- a/tests/test_hxn.py +++ b/tests/test_hxn.py @@ -10,6 +10,7 @@ Tests for the heat exchanger network facility. """ import warnings +import pytest import biosteam as bst import numpy as np from numpy.testing import assert_allclose @@ -90,8 +91,177 @@ def test_energy_balance_error_contributions_ignored_none(): assert len(errors) == N assert HXN.ignored is None +# --------------------------------------------------------------------------- +# Problem-table (pinch) analysis +# --------------------------------------------------------------------------- + +from biosteam.facilities.hxn.hxn_synthesis import ( + temperature_interval_pinch_analysis, problem_table, +) + +def utility_hx(ID, T, P, phase, T_out, **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) + hx = bst.HXutility(ID, ins=s, T=T_out, + rigorous=(phase == 'g')) # liquid streams stay liquid (report's case) + hx.simulate() + return hx + +def synthetic_units(): + """Report's 4-stream case: two cold, one condensing hot, one hot liquid.""" + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + bst.main_flowsheet.set_flowsheet('test_hxn_synthetic') + 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.)] + +def heat_utilities(units): + hus = [hx.heat_utilities[0] for hx in units] + hus.sort(key=lambda hu: hu.duty) + return hus + +def pinch_streams(hus): + """Inlet/quenched-outlet stream copies exactly as the synthesizer prepares them.""" + 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] + return streams_inlet, streams_quenched, is_hot + +def assert_energy_consistent(hus, T_min_app): + """Invariants of a correct problem table, independent of the synthesizer.""" + unit_duties = np.array([hu.unit_duty for hu in hus]) + table = problem_table(*pinch_streams(hus), T_min_app) + # (a) each stream's grid contributions telescope to its real duty + # (hot: +|dH|, cold: -dH) + per_stream = table.interval_H.sum(axis=1) + table.point_H.sum(axis=1) + assert_allclose(per_stream, -unit_duties, rtol=1e-9) + # (b) targets are non-negative and hot - cold is the net demand + assert table.hot_util_load >= 0. and table.cold_util_load >= 0. + assert_allclose(table.hot_util_load - table.cold_util_load, + unit_duties.sum(), rtol=1e-9) + # (c) a target can never exceed the un-integrated load + assert table.hot_util_load <= unit_duties[unit_duties > 0].sum() * (1 + 1e-9) + assert table.cold_util_load <= -unit_duties[unit_duties < 0].sum() * (1 + 1e-9) + # (d) the public wrapper reports the same targets + pinch_T_arr, hot, cold, *_ = temperature_interval_pinch_analysis(hus, T_min_app) + assert_allclose([hot, cold], [table.hot_util_load, table.cold_util_load], rtol=1e-12) + return table + +def test_problem_table_energy_consistency_doctest_system(): + sys, HXN, feed = build_system() + sys.simulate() + hus = HXN._get_original_heat_utilties() + hus.sort(key=lambda hu: hu.duty) + assert_energy_consistent(hus, 5.) + +@pytest.mark.parametrize('T_min_app', [5., 10., 20.]) +def test_problem_table_energy_consistency_synthetic(T_min_app): + hus = heat_utilities(synthetic_units()) + table = assert_energy_consistent(hus, T_min_app) + # the condensing ethanol stream (352 K in, ~351.4 K dew point) must keep + # its latent heat: hot target well below the un-integrated heating load + heating = sum(hu.unit_duty for hu in hus if hu.unit_duty > 0) + assert table.hot_util_load < 0.5 * heating + +def test_problem_table_two_streams_closed_form(): + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + bst.main_flowsheet.set_flowsheet('test_hxn_two_streams') + hot = utility_hx('Hw', 400., 5e5, 'l', 300., Water=1000.) + # Case 1 (threshold problem): the hot stream covers every interval of the + # smaller cold stream -> no hot utility, cold utility = net surplus. + cold = utility_hx('Cs', 300., 5e5, 'l', 390., Water=900.) + hus = heat_utilities([hot, cold]) + table = problem_table(*pinch_streams(hus), 5.) + Q_hot = -hot.heat_utilities[0].unit_duty + Q_cold = cold.heat_utilities[0].unit_duty + assert table.hot_util_load == 0. + assert_allclose(table.cold_util_load, Q_hot - Q_cold, rtol=1e-9) + assert table.pinch_T == table.Ts[0] # no pinch: everything sits below it + # Case 2: the larger cold stream is short of heat everywhere; the cascade + # minimum is at the cold inlet (300 K on the shifted scale) and the hot + # utility is the cold duty minus what the hot stream gives down to 305 K. + cold = utility_hx('Cb', 300., 5e5, 'l', 390., Water=1200.) + hus = heat_utilities([hot, cold]) + table = problem_table(*pinch_streams(hus), 5.) + s = hot.ins[0].copy(); s.vle(T=305., P=s.P) + Q_hot_above_pinch = hot.ins[0].H - s.H + Q_cold = cold.heat_utilities[0].unit_duty + assert_allclose(table.hot_util_load, Q_cold - Q_hot_above_pinch, rtol=1e-9) + assert table.pinch_T == 300. + # remaining hot-stream heat below the pinch leaves as cold utility + assert_allclose(table.cold_util_load, s.H - hot.outs[0].H, rtol=1e-9) + +def test_problem_table_non_monotone_stream_is_point_load(): + # A heated stream whose outlet is colder than its inlet (e.g. a column + # reboiler outlet at VLE): treated as an isothermal load at T_out. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + a = bst.Stream('a', Water=100., T=372., P=101325., phase='l', units='kmol/hr') + b = bst.Stream('b', Water=100., T=371., P=101325., phase='g', units='kmol/hr') + dH = b.H - a.H + assert dH > 0 + table = problem_table([a], [b], [False], 5.) + assert_allclose(table.point_H, [[-dH]]) + assert table.interval_H.shape == (1, 0) + assert_allclose([table.hot_util_load, table.cold_util_load, table.pinch_T], + [dH, 0., 371.]) + table = problem_table([b], [a], [True], 5.) + assert_allclose([table.hot_util_load, table.cold_util_load, table.pinch_T], + [0., dH, 372. - 5.]) # outlet T 372 K, shifted by T_min_app + +def test_problem_table_point_load_cannot_heat_above_itself(): + # A source at shifted temperature T cannot serve sinks above T: an + # isothermal condensing hot stream at 400 K (shifted to 395 K) must not + # cover the 395-398 K segment of a cold stream that runs 392 -> 398 K. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + hot_in = bst.Stream('hv', Water=100., T=400., P=101325., phase='g', + units='kmol/hr') + hot_out = bst.Stream('hl', Water=100., T=400., P=101325., phase='l', + units='kmol/hr') + P_hot = hot_in.H - hot_out.H + assert P_hot > 0 + cold_in = bst.Stream('cl', Water=1000., T=392., P=5e5, phase='l', + units='kmol/hr') + cold_out = cold_in.copy('cl_out'); cold_out.vle(T=398., P=5e5) + s = cold_in.copy(); s.vle(T=395., P=5e5) + H_392, H_395, H_398 = cold_in.H, s.H, cold_out.H + assert P_hot > H_398 - H_392 # more than enough heat overall + table = problem_table([hot_in, cold_in], [hot_out, cold_out], + [True, False], 5.) + assert_allclose(table.Ts, [398., 395., 392.]) + # heat arriving at 395 K, before the point load there, is short by the + # 395-398 K segment of the cold stream + assert_allclose(table.hot_util_load, H_398 - H_395, rtol=1e-9) + assert table.pinch_T == 395. + assert_allclose(table.cold_util_load, P_hot - (H_395 - H_392), rtol=1e-9) + +def test_synthetic_network_reaches_MER(): + units = synthetic_units() + HXN = bst.HeatExchangerNetwork('HXN', T_min_app=5.) + sys = bst.System.from_units('sys_synthetic', units=[*units, HXN]) + sys.simulate() + hus = heat_utilities(units) + table = problem_table(*pinch_streams(hus), 5.) + actual_heat = sum(hu.unit_duty for hx in HXN.new_HX_utils + for hu in hx.heat_utilities if hu.unit_duty > 0) + # the greedy heuristic reaches the (corrected) MER on this case; it can + # never legitimately beat it + assert_allclose(actual_heat, table.hot_util_load, rtol=1e-2) + # the lower bound is exact here only because every synthetic inlet is an + # equilibrium state (so the clipped table is exact) and the synthesizer + # respects T_min_app; with non-equilibrium inlets the table is conservative + assert actual_heat >= table.hot_util_load * (1 - 1e-3) + if __name__ == '__main__': test_cache_network_matches_fresh_synthesis() test_cache_network_perturbed_feed() test_cache_network_duplicate_IDs() test_energy_balance_error_contributions_ignored_none() + test_problem_table_energy_consistency_doctest_system() + for T_min_app in (5., 10., 20.): + test_problem_table_energy_consistency_synthetic(T_min_app) + test_problem_table_two_streams_closed_form() + test_problem_table_non_monotone_stream_is_point_load() + test_problem_table_point_load_cannot_heat_above_itself() + test_synthetic_network_reaches_MER()