Coverage for backend/idaes_service/solver/custom/thermal_utility_systems/header.py: 87%
231 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +0000
1# Pyomo core
2import pyomo.environ as pyo
3from pyomo.environ import (
4 Constraint,
5 Expression,
6 Param,
7 PositiveReals,
8 RangeSet,
9 Suffix,
10 Var,
11 value,
12 units as UNIT,
13)
14from pyomo.core.base.reference import Reference
15from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool
17# IDAES core
18from idaes.core import (
19 declare_process_block_class,
20 UnitModelBlockData,
21 useDefault,
22 StateBlock,
23)
24from idaes.core.util import scaling
25from idaes.core.util.config import is_physical_parameter_block
26from idaes.core.util.math import smooth_min, smooth_max
27from idaes.core.util.tables import create_stream_table_dataframe
28from idaes.core.solvers import get_solver
29from idaes.core.initialization import ModularInitializerBase
30from idaes.core.util.model_statistics import degrees_of_freedom
32# Logger
33import idaes.logger as idaeslog
35# Typing
36from typing import List
39__author__ = "Ahuora Centre for Smart Energy Systems, University of Waikato, New Zealand"
41# Set up logger
42_log = idaeslog.getLogger(__name__)
44class SimpleHeaderInitializer(ModularInitializerBase):
45 """Initialize a Header unit block with staged seeding and solves.
47 This routine performs a two-stage initialization:
48 1) Seed inlet and internal state variables, relax selected constraints, and
49 perform a first solve.
50 2) Reactivate/tighten constraints and perform a second solve.
52 Args:
53 blk: The Header unit model block to initialize.
54 **kwargs: Optional keyword arguments:
55 solver: A Pyomo/IDAES solver object. If not provided, uses ``get_solver()``.
56 solver_options (dict): Options to set on the solver, e.g. tolerances.
57 outlvl: IDAES log level (e.g., ``idaeslog.WARNING``).
59 Returns:
60 pyomo.opt.results.results_.SolverResults: The result object from the final solve.
62 Notes:
63 - Inlet state blocks are initialized via their own ``initialize`` if available.
64 - Mixed state is seeded from inlet totals/minimums (pressure) and average
65 enthalpy; works with temperature- or enthalpy-based property packages.
66 - Temporary seeds/relaxations are undone, leaving original DOF intact.
67 """
69 def initialize(self, blk, **kwargs):
70 # --- Solver setup
71 solver = kwargs.get("solver", None) or get_solver()
72 solver_options = kwargs.get("solver_options", {})
73 for k, v in solver_options.items(): 73 ↛ 74line 73 didn't jump to line 74 because the loop on line 73 never started
74 solver.options[k] = v
76 outlvl = kwargs.get("outlvl", idaeslog.WARNING)
77 log = idaeslog.getLogger(__name__)
79 # --- Time index
80 t0 = blk.flowsheet().time.first()
82 # --- 1) Initialize inlet state blocks
83 inlet_blocks = list(blk.inlet_blocks)
84 if len(inlet_blocks) < 1: 84 ↛ 85line 84 didn't jump to line 85 because the condition on line 84 was never true
85 raise ValueError("No inlet added to header.")
87 for sb in inlet_blocks:
88 if hasattr(sb, "initialize"): 88 ↛ 87line 88 didn't jump to line 87 because the condition on line 88 was always true
89 sb.initialize(outlvl=outlvl)
91 # --- 2) Aggregate inlet info for seeding mixed state block
92 F_mixed = sum(
93 value(sb[t0].flow_mol)
94 for sb in inlet_blocks
95 )
96 P_mixed = min(
97 value(sb[t0].pressure)
98 for sb in inlet_blocks
99 )
100 E_mixed = sum(
101 value(sb[t0].flow_mol * sb[t0].enth_mol, 0.0)
102 for sb in inlet_blocks
103 )
104 if F_mixed > 0:
105 h_mixed = E_mixed / F_mixed
106 else:
107 # Seed from the first inlet’s enthalpy (no double subscripting)
108 first_inlet = inlet_blocks[0]
109 h_mixed = value(first_inlet[t0].enth_mol)
111 # --- 3) Seed mixed_state: flow, pressure, enthalpy
112 ms = blk.mixed_state
113 ms[t0].flow_mol.set_value(
114 F_mixed
115 )
116 ms[t0].pressure.set_value(
117 P_mixed
118 )
119 ms[t0].enth_mol.set_value(
120 h_mixed
121 )
122 ms.initialize(outlvl=outlvl)
124 # --- 4) Seed outlet_states with pressure, enthalpy
125 flow_undefined = []
126 defined_flow = 0
127 for sb in blk.outlet_blocks:
128 sb[t0].pressure.set_value(
129 value(ms[t0].pressure)
130 )
131 sb[t0].enth_mol.set_value(
132 value(ms[t0].enth_mol)
133 )
134 if sb in [blk.outlet_condensate_state, blk.outlet_vent_state]:
135 sb[t0].flow_mol.set_value(
136 0.0
137 )
138 else:
139 if value(sb[t0].flow_mol, exception=False) is None: 139 ↛ 140line 139 didn't jump to line 140 because the condition on line 139 was never true
140 flow_undefined.append(sb)
141 else:
142 defined_flow += value(sb[t0].flow_mol)
144 tot_undefined_flow = max(sum(value(sb[t0].flow_mol) for sb in blk.inlet_blocks) - defined_flow, 0)
145 for sb in flow_undefined: 145 ↛ 146line 145 didn't jump to line 146 because the loop on line 145 never started
146 sb[t0].flow_mol.set_value(
147 tot_undefined_flow / len(flow_undefined)
148 )
150 for sb in blk.outlet_blocks:
151 sb.initialize(outlvl=outlvl)
153 res2 = solver.solve(blk, tee=False)
154 log.info(f"Header init status: {res2.solver.termination_condition}")
156 return res2
158def _make_config_block(config):
159 """Declare configuration options for the Header unit.
161 Declares property package references and integer counts for inlets and outlets.
163 Args:
164 config (ConfigBlock): The mutable configuration block to populate.
166 Raises:
167 ValueError: If invalid option values are provided by the caller (via IDAES).
168 """
170 config.declare(
171 "property_package",
172 ConfigValue(
173 default=useDefault,
174 domain=is_physical_parameter_block,
175 description="Property package to use for control volume",
176 ),
177 )
178 config.declare(
179 "property_package_args",
180 ConfigBlock(
181 implicit=True,
182 description="Arguments to use for constructing property packages",
183 ),
184 )
185 config.declare(
186 "num_inlets",
187 ConfigValue(
188 default=1,
189 domain=In(list(range(1, 100))),
190 description="Number of utility providers at inlets.",
191 ),
192 )
193 config.declare(
194 "num_outlets",
195 ConfigValue(
196 default=1,
197 domain=In(list(range(0, 100))),
198 description="Number of utility users at outlets." \
199 "Excludes outlets associated with condensate and vent flows.",
200 ),
201 )
202 config.declare(
203 "is_liquid_header",
204 ConfigValue(
205 default=False,
206 domain=Bool,
207 description="Flag for selecting liquid or vapour (including steam and other gases).",
208 ),
209 )
210@declare_process_block_class("simple_header")
211class SimpleHeaderData(UnitModelBlockData):
212 """Thermal utility header unit operation.
214 The Header aggregates multiple inlet providers and distributes utility to
215 multiple users, with optional venting, condensate removal (or liquid overflow), heat loss, and
216 pressure loss. A mixed (intermediate) state is used for balances and
217 pressure/enthalpy coupling across outlets.
219 Key features:
220 - Material, energy, and momentum balances with smooth min/max functions.
221 - Vapour/liquid equilibrium calculation for mixed state.
222 - Shared mixed enthalpy across outlets of the same phase.
223 - Computed excess flow from an overall flow balance.
224 - Optional heat and pressure losses.
226 Attributes:
227 inlet_list (list[str]): Names for inlet ports.
228 outlet_list (list[str]): Names for outlet ports (incl. condensate/ and vent).
229 inlet_blocks (list): StateBlocks for all inlets.
230 outlet_blocks (list): StateBlocks for all outlets.
231 mixed_state: Intermediate mixture StateBlock.
232 heat_loss (Var): Heat loss from the header (W).
233 pressure_loss (Var): Pressure drop from inlet minimum to mixed state (Pa).
234 makeup_flow_mol (Var): Required inlet makeup molar flow (mol/s).
235 """
237 default_initializer=SimpleHeaderInitializer
238 CONFIG = UnitModelBlockData.CONFIG()
239 _make_config_block(CONFIG)
241 def build(self) -> None:
242 # 1. Inherit standard UnitModelBlockData properties and functions
243 super().build()
245 # 2. Validate input parameters are valid
246 self._validate_model_config()
248 # 3. Create lists of ports with state blocks to add
249 self.inlet_list = self._create_inlet_port_name_list()
250 self.outlet_list = self._create_outlet_port_name_list()
252 # 4. Declare ports, state blocks and state property bounds
253 self.inlet_blocks = self._add_ports_with_state_blocks(
254 stream_list=self.inlet_list,
255 is_inlet=True,
256 has_phase_equilibrium=False,
257 is_defined_state=True,
258 )
259 self.outlet_blocks = self._add_ports_with_state_blocks(
260 stream_list=self.outlet_list,
261 is_inlet=False,
262 has_phase_equilibrium=False,
263 is_defined_state=False
264 )
265 self._internal_blocks = self._add_internal_state_blocks()
266 self._add_bounds_to_state_properties()
267 self._outlet_supply_blocks = self._create_custom_state_lists()
269 # 4. Declare references, variables and expressions for external and internal use
270 self._create_references()
271 self._create_variables()
272 self._create_expressions()
274 # 5. Set balance equations
275 self._add_material_balances()
276 self._add_energy_balances()
277 self._add_momentum_balances()
278 self._add_additional_constraints()
280 # 6. Other
281 self.scaling_factor = Suffix(direction=Suffix.EXPORT)
282 self.split_flow = self._create_flow_map_references()
284 def _validate_model_config(self) -> bool:
285 """Validate configuration for inlet and outlet counts.
287 Raises:
288 ValueError: If ``num_inlets < 1`` or ``num_outlets < 1``.
289 """
290 if self.config.num_inlets < 1: 290 ↛ 291line 290 didn't jump to line 291 because the condition on line 290 was never true
291 raise ValueError("Header requires at least one provider (num_inlets >= 1).")
292 if self.config.num_outlets < 1: 292 ↛ 293line 292 didn't jump to line 293 because the condition on line 292 was never true
293 raise ValueError("Header requires at least one user (num_outlets >= 1).")
294 return True
296 def _create_inlet_port_name_list(self) -> List[str]:
297 """Build ordered inlet port names.
299 Returns:
300 list[str]: Names ``["inlet_1", ..., "inlet_N"]`` based on ``num_inlets``.
301 """
302 return [
303 f"inlet_{i+1}" for i in range(self.config.num_inlets)
304 ]
306 def _create_outlet_port_name_list(self) -> List[str]:
307 """Build ordered outlet port names.
309 Returns:
310 list[str]: Names ``["outlet_1", ..., "outlet_n", "outlet_condensate", "outlet_vent"]``.
311 """
312 return [
313 f"outlet_{i+1}" for i in range(self.config.num_outlets)
314 ] + ["outlet_condensate"] + ["outlet_vent"]
316 def _add_ports_with_state_blocks(self,
317 stream_list: List[str],
318 is_inlet: List[str],
319 has_phase_equilibrium: bool=False,
320 is_defined_state: bool=None,
321 ) -> List[StateBlock]:
322 """Construct StateBlocks and expose them as ports.
324 Creates a StateBlock per named stream and attaches a corresponding inlet or
325 outlet Port. Inlet blocks are defined states; outlet blocks are calculated states.
327 Args:
328 stream_list (list[str]): Port/StateBlock base names to create.
329 is_inlet (bool): If True, create inlet ports with ``defined_state=True``;
330 otherwise create outlet ports with ``defined_state=False``.
331 has_phase_equilibrium (bool)
333 Returns:
334 list: The created StateBlocks, in the same order as ``stream_list``.
335 """
336 # Create empty list to hold StateBlocks for return
337 state_block_ls = []
339 # Setup StateBlock argument dict
340 tmp_dict = dict(**self.config.property_package_args)
341 tmp_dict["has_phase_equilibrium"] = has_phase_equilibrium
342 if is_defined_state == None: 342 ↛ 343line 342 didn't jump to line 343 because the condition on line 342 was never true
343 tmp_dict["defined_state"] = True if is_inlet else False
344 else:
345 tmp_dict["defined_state"] = is_defined_state
347 # Create an instance of StateBlock for all streams
348 for s in stream_list:
349 sb = self.config.property_package.build_state_block(
350 self.flowsheet().time, doc=f"Thermophysical properties at {s}", **tmp_dict
351 )
352 setattr(
353 self, s + "_state",
354 sb
355 )
356 state_block_ls.append(sb)
357 add_fn = self.add_inlet_port if is_inlet else self.add_outlet_port
358 add_fn(
359 name=s,
360 block=sb,
361 )
363 return state_block_ls
365 def _add_internal_state_blocks(self) -> List[StateBlock]:
366 """Create the intermediate (mixed) StateBlock.
368 The mixed state:
369 - Has phase equilibrium enabled.
370 - Is not a defined state (solved from balances).
371 """
372 tmp_dict = dict(**self.config.property_package_args)
373 tmp_dict["has_phase_equilibrium"] = True
374 tmp_dict["defined_state"] = False
376 self.mixed_state = self.config.property_package.build_state_block(
377 self.flowsheet().time,
378 doc=f"Thermophysical properties at intermediate mixed state.",
379 **tmp_dict
380 )
381 return [
382 self.mixed_state
383 ]
385 def _add_bounds_to_state_properties(self) -> None:
386 """Add lower and/or upper bounds to state properties.
388 - Set nonnegativity lower bounds on all inlet/outlet molar flows.
389 """
390 for sb in (self.inlet_blocks + self.outlet_blocks):
391 for t in sb:
392 sb[t].flow_mol.setlb(0.0)
394 def _create_custom_state_lists(self) -> List[StateBlock]:
395 """Partition outlet names into vapour outlets and capture their StateBlocks.
397 Populates:
398 - ``_outlet_supply_list``: Outlet names excluding condensate and vent.
399 - ``_outlet_supply_blocks``: Corresponding StateBlocks.
400 """
401 self._outlet_supply_list = [
402 v for v in self.outlet_list
403 if not v in ["outlet_condensate", "outlet_vent"]
404 ]
405 return [
406 getattr(self, n + "_state")
407 for n in self._outlet_supply_list
408 ]
410 def _create_references(self) -> None:
411 """Create convenient References.
413 Creates references to mixed_state properties:
414 - ``total_flow_mol``
415 - ``total_flow_mass``
416 - ``pressure``
417 - ``temperature``
418 - ``enth_mol``
419 - ``enth_mass``
420 - ``vapor_frac``
421 """
422 self.total_flow_mol = Reference(
423 self.mixed_state[:].flow_mol
424 )
425 self.total_flow_mass = Reference(
426 self.mixed_state[:].flow_mass
427 )
428 self.pressure = Reference(
429 self.mixed_state[:].pressure
430 )
431 self.temperature = Reference(
432 self.mixed_state[:].temperature
433 )
434 self.enth_mol = Reference(
435 self.mixed_state[:].enth_mol
436 )
437 self.enth_mass = Reference(
438 self.mixed_state[:].enth_mass
439 )
440 self.vapor_frac = Reference(
441 self.mixed_state[:].vapor_frac
442 )
444 def _create_variables(self) -> None:
445 """Create required variables.
447 Creates:
448 - ``heat_loss`` (W)
449 - ``pressure_loss`` (Pa)
450 """
451 self.heat_loss = Var(
452 self.flowsheet().time,
453 initialize=0.0,
454 doc="Heat loss",
455 units=UNIT.W
456 )
457 self.pressure_loss = Var(
458 self.flowsheet().time,
459 initialize=0.0,
460 doc="Pressure loss",
461 units=UNIT.Pa
462 )
464 def _create_expressions(self) -> None:
465 """Create convenient Expressions.
467 Creates:
468 - ``balance_flow_mol`` (mol/s)
469 - ``degree_of_superheat`` (K)
470 - ``makeup_flow_mol`` (mol/s)
471 - ``_partial_total_flow_mol`` (mol/s): used for scaling purposes in a material balance
472 """
473 self.degree_of_superheat = Expression(
474 self.flowsheet().time,
475 rule=lambda b, t: b.temperature[t] - b.outlet_condensate_state[t].temperature
476 )
477 self._partial_total_flow_mol = Expression(
478 self.flowsheet().time,
479 rule=lambda b, t: (
480 sum(
481 o[t].flow_mol
482 for o in (b.inlet_blocks + b._outlet_supply_blocks)
483 )
484 )
485 )
486 self.balance_flow_mol = Expression(
487 self.flowsheet().time,
488 rule=lambda b, t: (
489 sum(
490 i[t].flow_mol
491 for i in b.inlet_blocks
492 )
493 -
494 sum(
495 o[t].flow_mol
496 for o in (
497 b._outlet_supply_blocks +
498 [
499 b.outlet_vent_state
500 if self.config.is_liquid_header
501 else b.outlet_condensate_state
502 ]
503 )
504 )
505 )
506 )
507 self.makeup_flow_mol = Expression(
508 self.flowsheet().time,
509 rule=lambda b, t: (
510 (
511 b.outlet_condensate_state[t].flow_mol
512 if self.config.is_liquid_header
513 else b.outlet_vent_state[t].flow_mol
514 )
515 -
516 b.balance_flow_mol[t]
517 )
518 )
520 def _add_material_balances(self) -> None:
521 """Material balance equations summary.
523 Introduces:
524 - ``_partial_total_flow_mol``: Sum of known inlet and vapour outlet flows,
525 used for scaling a smooth vent calculation.
527 Constraints:
528 - ``mixed_state_material_balance``: Mixed flow equals total inlet flow.
529 - ``vent_flow_balance``: Depends on the header's primary phase: liquid vs gas
530 If gas header, smoothly enforces nonnegative vent flow.
531 If liquid header, determines flow from mixed-state vapour fraction
532 - ``condensate_flow_balance``: Depends on the header's primary phase: liquid vs gas
533 If gas header, determines flow from mixed-state vapour fraction
534 If liquid header, smoothly enforces nonnegative condensate flow.
535 """
537 @self.Constraint(
538 self.flowsheet().time,
539 doc="Mixed state material balance",
540 )
541 def mixed_state_material_balance(b, t):
542 return (
543 b.mixed_state[t].flow_mol
544 ==
545 sum(
546 i[t].flow_mol
547 for i in b.inlet_blocks
548 )
549 )
551 eps = 1e-5 # smoothing parameter; smaller = closer to exact max, larger = smoother
552 if self.config.is_liquid_header:
553 # Assigns excess liquid flow to outlet_condensate
554 @self.Constraint(
555 self.flowsheet().time,
556 doc="Condensate flow balance." \
557 "Determines the positive amount of excess flow that exits through outlet_condensate"
558 )
559 def condensate_flow_balance(b, t):
560 return (
561 b.outlet_condensate_state[t].flow_mol
562 ==
563 smooth_max(
564 b.balance_flow_mol[t] / (b._partial_total_flow_mol[t] + 1e-6),
565 0.0,
566 eps,
567 ) * (b._partial_total_flow_mol[t] + 1e-6)
568 )
570 # Removes any gas/vapour from a liquid header
571 @self.Constraint(
572 self.flowsheet().time,
573 doc="Vent balance."
574 )
575 def vent_flow_balance(b, t):
576 return b.outlet_vent_state[t].flow_mol == (
577 b.mixed_state[t].flow_mol * b.mixed_state[t].vapor_frac
578 )
579 else:
580 # Assigns excess steam/vapour flow to outlet_vent
581 @self.Constraint(
582 self.flowsheet().time,
583 doc="Vent flow balance." \
584 "Determines the positive amount of excess flow that exits through the vent"
585 )
586 def vent_flow_balance(b, t):
587 return (
588 b.outlet_vent_state[t].flow_mol
589 ==
590 smooth_max(
591 b.balance_flow_mol[t] / (b._partial_total_flow_mol[t] + 1e-6),
592 0.0,
593 eps,
594 ) * (b._partial_total_flow_mol[t] + 1e-6)
595 )
597 # Removes any condensate/liquid from a steam/gas header
598 @self.Constraint(
599 self.flowsheet().time,
600 doc="Condensate balance."
601 )
602 def condensate_flow_balance(b, t):
603 return (
604 b.outlet_condensate_state[t].flow_mol
605 ==
606 b.mixed_state[t].flow_mol * (1 - b.mixed_state[t].vapor_frac)
607 )
609 def _add_energy_balances(self) -> None:
610 """Energy balance equations summary.
612 Introduces:
613 - ``_liq_out_enth_mol``: Shared molar enthalpy for all liquid outlets,
614 including the condensate.
615 - ``_vap_out_enth_mol``: Shared molar enthalpy for all vapour outlets,
616 including the vent.
618 Constraints:
619 - ``inlets_to_mixed_state_energy_balance``: Inlet energy to mixed state (+ heat loss).
620 - ``mixed_state_to_outlets_energy_balance``: Mixed state to all outlets.
621 - ``molar_enthalpy_equality_eqn``: Common vapour enthalpy across vapour outlets and vent.
622 """
623 @self.Constraint(self.flowsheet().time, doc="Inlets to mixed state energy balance including heat loss")
624 def inlets_to_mixed_state_energy_balance(b, t):
625 return (
626 b.mixed_state[t].flow_mol * b.mixed_state[t].enth_mol
627 + b.heat_loss[t]
628 ==
629 sum(
630 i[t].flow_mol * i[t].enth_mol
631 for i in b.inlet_blocks
632 )
633 )
634 @self.Constraint(
635 self.flowsheet().time,
636 doc="Mixed state to outlets energy balance"
637 )
638 def mixed_state_to_outlets_energy_balance(b, t):
639 return (
640 b.mixed_state[t].enth_mol
641 *
642 sum(
643 o[t].flow_mol
644 for o in b.outlet_blocks
645 )
646 ==
647 sum(
648 o[t].flow_mol * o[t].enth_mol
649 for o in b.outlet_blocks
650 )
651 )
652 if self.config.is_liquid_header:
653 self._liq_out_enth_mol = Var(
654 self.flowsheet().time,
655 initialize=42.0 * 18,
656 doc="Molar enthalpy of the liquid outlets",
657 units=UNIT.J / UNIT.mol
658 )
659 @self.Constraint(
660 self.flowsheet().time,
661 self._outlet_supply_blocks + [self.outlet_condensate_state], # exclude vent outlet
662 doc="All liquid outlets (incl. condensate) share a common liquid enthalpy",
663 )
664 def molar_enthalpy_equality_eqn(b, t, o):
665 return (
666 o[t].enth_mol
667 ==
668 b._liq_out_enth_mol[t]
669 )
670 else:
671 self._vap_out_enth_mol = Var(
672 self.flowsheet().time,
673 initialize=2700.0 * 18,
674 doc="Molar enthalpy of the vapour outlets",
675 units=UNIT.J / UNIT.mol
676 )
677 @self.Constraint(
678 self.flowsheet().time,
679 self._outlet_supply_blocks + [self.outlet_vent_state], # exclude condensate outlet
680 doc="All vapour outlets (incl. vent) share a common vapour enthalpy",
681 )
682 def molar_enthalpy_equality_eqn(b, t, o):
683 return (
684 o[t].enth_mol
685 ==
686 b._vap_out_enth_mol[t]
687 )
689 def _add_momentum_balances(self) -> None:
690 """Momentum balance equations summary.
692 Computes the minimum inlet pressure via a sequential smooth minimum and
693 sets the mixed-state pressure to that minimum minus ``pressure_loss``,
694 then enforces equality to every outlet pressure.
696 Notes:
697 - Uses IDAES ``smooth_min`` for differentiable minimum pressure.
698 - ``_eps_pressure`` is a smoothing parameter (units of pressure).
699 """
700 inlet_idx = RangeSet(len(self.inlet_blocks))
701 # Get units metadata
702 units = self.mixed_state.params.get_metadata()
703 # Add variables
704 self._minimum_pressure = Var(
705 self.flowsheet().time,
706 inlet_idx,
707 doc="Variable for calculating minimum inlet pressure",
708 units=units.get_derived_units("pressure"),
709 )
710 self._eps_pressure = Param(
711 mutable=True,
712 initialize=1e-3,
713 domain=PositiveReals,
714 doc="Smoothing term for minimum inlet pressure",
715 units=units.get_derived_units("pressure"),
716 )
717 # Calculate minimum inlet pressure
718 @self.Constraint(
719 self.flowsheet().time,
720 inlet_idx,
721 doc="Calculation for minimum inlet pressure",
722 )
723 def minimum_pressure_constraint(b, t, i):
724 if i == inlet_idx.first():
725 return (
726 b._minimum_pressure[t, i]
727 ==
728 (b.inlet_blocks[i - 1][t].pressure)
729 )
730 else:
731 return (
732 b._minimum_pressure[t, i]
733 ==
734 smooth_min(
735 b._minimum_pressure[t, i - 1],
736 b.inlet_blocks[i - 1][t].pressure,
737 b._eps_pressure,
738 )
739 )
740 # Set mixed pressure to minimum inlet pressure minus any pressure loss
741 @self.Constraint(
742 self.flowsheet().time,
743 doc="Pressure equality constraint from minimum inlet to mixed state",
744 )
745 def mixture_pressure(b, t):
746 return (
747 b.mixed_state[t].pressure
748 ==
749 b._minimum_pressure[t, inlet_idx.last()] - b.pressure_loss[t]
750 )
751 # Set outlet pressures to mixed pressure
752 @self.Constraint(
753 self.flowsheet().time,
754 self.outlet_blocks,
755 doc="Pressure equality constraint from mixed state to outlets",
756 )
757 def pressure_equality_eqn(b, t, o):
758 return (
759 b.mixed_state[t].pressure
760 ==
761 o[t].pressure
762 )
764 def _add_additional_constraints(self) -> None:
765 """Add auxiliary constraints and bounds.
767 - Fix vent vapour fraction to near one (near 100% vapour).
768 OR
769 - Fix condensate vapour fraction to a small value (near 100% liquid).
770 """
771 if self.config.is_liquid_header:
772 @self.Constraint(self.flowsheet().time, doc="Vent vapour fraction.")
773 def vent_vapour_fraction(b, t):
774 return (
775 b.outlet_vent_state[t].vapor_frac
776 ==
777 1 #1 - 1e-6
778 )
779 else:
780 @self.Constraint(self.flowsheet().time, doc="Condensate vapour fraction.")
781 def condensate_vapour_fraction(b, t):
782 return (
783 b.outlet_condensate_state[t].vapor_frac
784 ==
785 0 #1e-6
786 )
788 def _create_flow_map_references(self):
789 """Create a two-key Reference for outlet flows over time and outlet name.
791 Builds a mapping ``(t, outlet_name) -> outlet_state[t].flow_mol`` and exposes it
792 as a Reference for compact access to outlet flow splits.
794 Returns:
795 pyomo.core.base.reference.Reference: A Reference indexed by ``(time, outlet)``.
796 """
797 self.outlet_idx = pyo.Set(initialize=self.outlet_list)
798 # Map each (t, o) to the outlet state's flow var
799 ref_map = {}
800 for o in self.outlet_list:
801 if o != "vent": 801 ↛ 800line 801 didn't jump to line 800 because the condition on line 801 was always true
802 outlet_state_block = getattr(self, f"{o}_state")
803 for t in self.flowsheet().time:
804 ref_map[(t, o)] = outlet_state_block[t].flow_mol
806 return Reference(ref_map)
808 def calculate_scaling_factors(self):
809 """Assign scaling factors to improve numerical conditioning.
811 Sets scaling factors for performance and auxiliary variables. If present,
812 also scales the shared vapour enthalpy variable ``_vap_out_enth_mol``.
813 """
814 super().calculate_scaling_factors()
815 scaling.set_scaling_factor(self.heat_loss, 1e-6)
816 scaling.set_scaling_factor(self.pressure_loss, 1e-6)
817 scaling.set_scaling_factor(self.balance_flow_mol, 1e-3)
818 scaling.set_scaling_factor(self._partial_total_flow_mol, 1e-3)
819 if hasattr(self, "_vap_out_enth_mol"):
820 scaling.set_scaling_factor(self._vap_out_enth_mol, 1e-6)
822 def _get_stream_table_contents(self, time_point=0):
823 """Create a stream table for all inlets and outlets.
825 Args:
826 time_point (int | float): Time index at which to extract stream data.
828 Returns:
829 pandas.DataFrame: A tabular view suitable for reporting via
830 ``create_stream_table_dataframe``.
831 """
832 io_dict = {}
834 for inlet_name in self.inlet_list:
835 io_dict[inlet_name] = getattr(self, inlet_name)
837 for outlet_name in self.outlet_list:
838 io_dict[outlet_name] = getattr(self, outlet_name)
840 return create_stream_table_dataframe(io_dict, time_point=time_point)
842 def _get_performance_contents(self, time_point=0, is_full_report=True):
843 """Collect performance results for reporting.
845 Args:
846 time_point (int | float): Time index at which to report values.
847 is_full_report (bool): Flag for full or partial performance report.
849 Returns:
850 dict: A report of internal unit model results.
851 """
852 return (
853 {
854 "vars": {
855 "Heat Loss": self.heat_loss[time_point],
856 "Pressure Drop": self.pressure_loss[time_point],
857 "Mass Flow": self.mixed_state[time_point].flow_mass,
858 "Molar Flow": self.mixed_state[time_point].flow_mol,
859 "Balance Flow": self.balance_flow_mol[time_point],
860 "Pressure": self.mixed_state[time_point].pressure,
861 "Temperature": self.mixed_state[time_point].temperature,
862 "Degree of Superheat": self.degree_of_superheat[time_point],
863 "Vapour Fraction": self.mixed_state[time_point].vapor_frac,
864 "Mass Specific Enthalpy": self.mixed_state[time_point].enth_mass,
865 "Molar Specific Enthalpy": self.mixed_state[time_point].enth_mol,
866 }
867 } if is_full_report else {
868 "vars": {
869 "Balance Flow": self.balance_flow_mol[time_point],
870 "Pressure": self.mixed_state[time_point].pressure,
871 "Temperature": self.mixed_state[time_point].temperature,
872 "Degree of Superheat": self.degree_of_superheat[time_point],
873 }
874 }
876 )
878 def initialize(self, *args, **kwargs):
879 """Initialize the Header unit using :class:`SimpleHeaderInitializer`.
881 Args:
882 *args: Forwarded to ``SimpleHeaderInitializer.initialize``.
883 **kwargs: Forwarded to ``SimpleHeaderInitializer.initialize`` (e.g., solver, options).
885 Returns:
886 pyomo.opt.results.results_.SolverResults: Results from the initializer's solve.
887 """
888 init = SimpleHeaderInitializer()
889 return init.initialize(self, *args, **kwargs)