Coverage for backend/idaes_service/solver/custom/thermal_utility_systems/steam_header.py: 90%
129 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +0000
1from pyomo.environ import Suffix, Var, Expression, Constraint
2from pyomo.common.config import ConfigBlock, ConfigValue, In
3from pyomo.network import Arc, SequentialDecomposition
4from pyomo.core.base.reference import Reference
5import pyomo.environ as pyo
7# Import IDAES cores
8from idaes.core import (
9 declare_process_block_class,
10 UnitModelBlockData,
11 useDefault,
12)
13from idaes.core.util.config import is_physical_parameter_block
14import idaes.logger as idaeslog
15from idaes.core.util.tables import create_stream_table_dataframe
16from idaes.core.util.math import smooth_max
17from idaes.models.unit_models.separator import SplittingType, EnergySplittingType
18from idaes.models.unit_models.mixer import Mixer
19from idaes.models.unit_models.heater import Heater
20from idaes_service.solver.custom.custom_separator import CustomSeparator
21from idaes_service.solver.custom.simple_separator import SimpleSeparator
22from ..inverted import add_inverted, initialise_inverted
23__author__ = "Team Ahuora"
25# Set up logger
26_log = idaeslog.getLogger(__name__)
29@declare_process_block_class("SteamHeader")
30class SteamHeaderData(UnitModelBlockData):
31 """
32 Steam Header unit operation:
33 Mixer -> Cooler -> Phase Separator -> Simple Separator
34 Separates 100% liquid to condensate_outlet and 100% vapor to splitter outlets.
35 Uses Sequential Decomposition for optimal initialization order.
36 """
38 CONFIG = UnitModelBlockData.CONFIG()
39 CONFIG.declare(
40 "property_package",
41 ConfigValue(
42 default=useDefault,
43 domain=is_physical_parameter_block,
44 description="Property package to use for control volume",
45 ),
46 )
47 CONFIG.declare(
48 "property_package_args",
49 ConfigBlock(
50 implicit=True,
51 description="Arguments to use for constructing property packages",
52 ),
53 )
54 CONFIG.declare(
55 "num_inlets",
56 ConfigValue(
57 default=2,
58 domain=int,
59 description="Number of inlets to add" "Index [-1]: Steam makeup",
60 ),
61 )
62 CONFIG.declare(
63 "num_outlets",
64 ConfigValue(
65 default=2,
66 domain=int,
67 description="Number of outlets to add"
68 "Index [-1]: Vent"
69 "Index [-2]: Condensate",
70 ),
71 )
73 def build(self):
74 super().build()
75 self.scaling_factor = Suffix(direction=Suffix.EXPORT)
77 self.inlet_list = [f"inlet_{i+1}" for i in range(self.config.num_inlets)]
78 self.outlet_list = [
79 f"outlet_{i+1}" for i in range(self.config.num_outlets)
80 ] + ["vent"]
82 # Create internal units
83 self.mixer = Mixer(
84 property_package=self.config.property_package,
85 property_package_args=self.config.property_package_args,
86 num_inlets=self.config.num_inlets,
87 inlet_list=self.inlet_list,
88 )
89 self.cooler = Heater(
90 property_package=self.config.property_package,
91 property_package_args=self.config.property_package_args,
92 has_pressure_change=True,
93 dynamic=self.config.dynamic,
94 has_holdup=self.config.has_holdup,
95 )
96 self.phase_separator = CustomSeparator(
97 property_package=self.config.property_package,
98 property_package_args=self.config.property_package_args,
99 outlet_list=["vapor_outlet", "condensate_outlet"],
100 split_basis=SplittingType.phaseFlow,
101 energy_split_basis=EnergySplittingType.enthalpy_split,
102 )
103 self.splitter = SimpleSeparator(
104 property_package=self.config.property_package,
105 property_package_args=self.config.property_package_args,
106 outlet_list=self.outlet_list,
107 )
108 self.unit_ops = [self.mixer, self.cooler, self.phase_separator, self.splitter]
110 # Updated internal arcs
111 self.mixer_to_cooler_arc = Arc(
112 source=self.mixer.outlet, destination=self.cooler.inlet
113 )
114 self.cooler_to_separator_arc = Arc(
115 source=self.cooler.outlet, destination=self.phase_separator.inlet
116 )
117 self.separator_to_splitter_arc = Arc(
118 source=self.phase_separator.vapor_outlet, destination=self.splitter.inlet
119 )
121 self.inlet_blocks = [getattr(self.mixer, f"{i}_state") for i in self.inlet_list]
122 self.outlet_blocks = [getattr(self.splitter, f"{o}_state") for o in self.outlet_list]
124 # Declare slack variables for internal use
125 self.balance_flow_mol = Var(
126 self.flowsheet().time, initialize=0.0, doc="Balance molar flow (negative means makeup is required, positive means venting)",
127 units=pyo.units.mol / pyo.units.s
128 )
130 # This is only used for scaling the smooth_min/smooth_max value.
131 self._partial_total_flow_mol = Var(
132 self.flowsheet().time, initialize=0.0, doc="Partial total fixed molar flow",
133 units=pyo.units.mol / pyo.units.s
134 )
136 self.makeup_flow_mol = Var(
137 self.flowsheet().time, initialize=0.0, doc="Makeup molar flow",
138 units=pyo.units.mol / pyo.units.s
139 )
142 # Add inverted transformers to heat_duty and deltaP
143 # (so that positive values correspond to heat loss and pressure drop)
145 add_inverted(self.cooler, "heat_duty")
146 add_inverted(self.cooler, "deltaP")
147 # Declare additional variables and aliases to expose to the user
148 self.heat_duty_inverted = Reference(self.cooler.heat_duty_inverted)
149 self.deltaP_inverted = Reference(self.cooler.deltaP_inverted)
150 self.heat_duty = Reference(self.cooler.heat_duty)
151 self.deltaP = Reference(self.cooler.deltaP)
153 self.total_flow_mol = Reference(
154 self.cooler.control_volume.properties_out[:].flow_mol
155 )
156 self.total_flow_mass = Reference(
157 self.cooler.control_volume.properties_out[:].flow_mass
158 )
159 self.pressure = Reference(
160 self.cooler.control_volume.properties_out[:].pressure
161 )
162 self.temperature = Reference(
163 self.cooler.control_volume.properties_out[:].temperature
164 )
165 self.enth_mol = Reference(
166 self.cooler.control_volume.properties_out[:].enth_mol
167 )
168 self.enth_mass = Reference(
169 self.cooler.control_volume.properties_out[:].enth_mass
170 )
171 self.vapor_frac = Reference(
172 self.cooler.control_volume.properties_out[:].vapor_frac
173 )
174 self.split_flow = self._create_split_flow_references()
176 # Condensate liquid always is removed first.
177 self.phase_separator.split_fraction[:, "vapor_outlet", "Vap"].fix(1.0)
178 self.phase_separator.split_fraction[:, "vapor_outlet", "Liq"].fix(0.0)
180 # Additional bounds and constraints
181 self._additional_constraints()
183 # Expand arcs
184 pyo.TransformationFactory("network.expand_arcs").apply_to(self)
186 self.cooler_to_separator_arc_expanded.flow_mol_equality.deactivate()
188 # add an expression for the degree of superheat
190 @self.Expression(self.flowsheet().time)
191 def degree_of_superheat(self,t):
192 return self.splitter.mixed_state[t].temperature - self.splitter.mixed_state[t].temperature_sat
194 def _create_split_flow_references(self):
195 self.outlet_idx = pyo.Set(initialize=self.outlet_list)
196 # Map each (t, o) to the outlet state's flow var
197 ref_map = {}
198 for o in self.outlet_list:
199 if o != "vent":
200 outlet_state_block = getattr(self.splitter, f"{o}_state")
201 for t in self.flowsheet().time:
202 ref_map[(t, o)] = outlet_state_block[t].flow_mol
204 return Reference(ref_map)
206 def _additional_constraints(self):
207 """
208 Additional constraints.
209 """
211 """
212 1) Set lower bounds on flow variables for all external ports.
213 """
214 [
215 state_block[t].flow_mol.setlb(0.0)
216 for state_block in (self.inlet_blocks + self.outlet_blocks)
217 for t in state_block
218 ]
220 """
221 2) Add overall material balance equation.
222 """
224 # Write phase-component balances
225 @self.Constraint(self.flowsheet().time, doc="Material balance equation")
226 def material_balance_equation(b, t):
227 pc_set = b.outlet_blocks[0].phase_component_set
228 comp_ls = b.outlet_blocks[0].component_list
229 phase_ls = b.outlet_blocks[0].phase_list
230 return (
231 0
232 == sum(
233 sum(
234 b.mixer.mixed_state[t].get_material_flow_terms(p, j)
235 - (b.phase_separator.mixed_state[t].get_material_flow_terms(p, j)
236 - b.splitter.vent_state[t].get_material_flow_terms(p, j))
237 for j in comp_ls
238 if (p, j) in pc_set
239 )
240 for p in phase_ls
241 )
242 - b.balance_flow_mol[t]
243 )
245 """
246 3) Add vent flow balance (scaled by the total known flows).
247 """
249 @self.Constraint(
250 self.flowsheet().time,
251 doc="Material balance of known inlet and outlet flows. "
252 "Only used for scaling the steam vent min flow, doesn't actually matter.",
253 )
254 def partial_material_balance_equation(b, t):
255 pc_set = b.inlet_blocks[0].phase_component_set
256 comp_ls = b.inlet_blocks[0].component_list
257 phase_ls = b.inlet_blocks[0].phase_list
258 return (
259 0
260 == sum(
261 sum(
262 sum(
263 o[t].get_material_flow_terms(p, j)
264 for o in b.inlet_blocks[:-1] + b.outlet_blocks[:-2]
265 )
266 for j in comp_ls
267 if (p, j) in pc_set
268 )
269 for p in phase_ls
270 )
271 - b._partial_total_flow_mol[t]
272 )
274 eps = 1e-5 # smoothing parameter; smaller = closer to exact max, larger = smoother
275 # calculate amount for vent flow: positive balance amount
276 @self.Constraint(self.flowsheet().time, doc="Steam vent flow balance.")
277 def vent_flow_balance(b, t):
278 return 0 == (
279 smooth_max(
280 b.balance_flow_mol[t] / (b._partial_total_flow_mol[t] + 1e-6),
281 0.0,
282 eps,
283 )
284 * (b._partial_total_flow_mol[t] + 1e-6)
285 - b.splitter.vent_state[t].flow_mol
286 )
288 # if balance is negative, vent will be zero. makeup will be required.
289 @self.Constraint(self.flowsheet().time, doc="Steam makeup balance")
290 def makeup_flow_balance(b, t):
291 return b.makeup_flow_mol[t] == (
292 b.splitter.vent_state[t].flow_mol - b.balance_flow_mol[t]
293 )
296 def calculate_scaling_factors(self):
297 super().calculate_scaling_factors()
298 [getattr(o, "calculate_scaling_factors")() for o in self.unit_ops]
300 def initialize_build(self, outlvl=idaeslog.NOTSET, **kwargs):
301 """
302 Initialize the Steam Header unit using Sequential Decomposition to determine optimal order
303 """
304 init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit")
306 init_log.info(
307 f"Starting {self.parent_block().__class__.__name__} initialization using Sequential Decomposition"
308 )
309 # Initialize the inverted values
310 initialise_inverted(self.cooler, "heat_duty")
311 initialise_inverted(self.cooler, "deltaP")
313 # create Sequential Decomposition object
314 seq = SequentialDecomposition()
315 seq.options.select_tear_method = "heuristic"
316 seq.options.tear_method = "Wegstein"
317 seq.options.iterLim = 1
319 # create computation graph
320 G = seq.create_graph(self)
321 heuristic_tear_set = seq.tear_set_arcs(G, method="heuristic")
322 # get calculation order
323 order = seq.calculation_order(G)
325 for o in heuristic_tear_set: 325 ↛ 326line 325 didn't jump to line 326 because the loop on line 325 never started
326 print(o.name)
328 print("Initialization order:")
329 for o in order:
330 print(o[0].name)
332 # define unit initialisation function
333 def init_unit(unit):
334 unit.initialize(outlvl=outlvl, **kwargs)
336 # run sequential decomposition
337 seq.run(self, init_unit)
339 def _get_stream_table_contents(self, time_point=0):
340 """
341 Create stream table showing all inlets, outlets, and liquid outlet
342 """
343 io_dict = {}
345 for inlet_name in self.inlet_list:
346 io_dict[inlet_name] = getattr(self, inlet_name)
348 for outlet_name in self.outlet_list:
349 io_dict[outlet_name] = getattr(self, outlet_name)
351 io_dict["condensate_outlet"] = self.condensate_outlet
353 return create_stream_table_dataframe(io_dict, time_point=time_point)