Coverage for backend/ahuora-builder/src/ahuora_builder/custom/thermal_utility_systems/desuperheater.py: 86%
169 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-06-23 21:51 +0000
1"""Desuperheater unit model."""
3from typing import Dict, Iterable, List, Optional
5from pyomo.common.config import Bool, ConfigBlock, ConfigValue, In
6import pyomo.environ as pyo
7from pyomo.environ import (
8 Constraint,
9 Expression,
10 NonNegativeReals,
11 Suffix,
12 Var,
13 Param,
14 value,
15 units as pyunits,
16 check_optimal_termination,
17)
19from idaes.core import StateBlock, UnitModelBlockData, declare_process_block_class, useDefault
20from idaes.core.initialization import ModularInitializerBase
21from idaes.core.solvers import get_solver
22from idaes.core.util.config import is_physical_parameter_block
23from idaes.core.util.exceptions import InitializationError
24from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme
25from idaes.core.util.tables import create_stream_table_dataframe
26from idaes.core.util.model_statistics import degrees_of_freedom, report_statistics
27from idaes.core.util.model_diagnostics import DiagnosticsToolbox
28from idaes.core.util.math import smooth_min
30import idaes.logger as idaeslog
32_log = idaeslog.getLogger(__name__)
34__author__ = "Ahuora Centre for Smart Energy Systems, University of Waikato, New Zealand"
37def _build_config(config: ConfigBlock) -> None:
38 """Declare config entries for Desuperheater."""
39 config.declare(
40 "dynamic",
41 ConfigValue(
42 domain=In([False]),
43 default=False,
44 description="Dynamic model flag - must be False",
45 ),
46 )
47 config.declare(
48 "has_holdup",
49 ConfigValue(
50 default=False,
51 domain=In([False]),
52 description="Holdup construction flag - must be False",
53 ),
54 )
55 config.declare(
56 "property_package",
57 ConfigValue(
58 default=useDefault,
59 domain=is_physical_parameter_block,
60 description="Property package used to build StateBlocks.",
61 ),
62 )
63 config.declare(
64 "property_package_args",
65 ConfigBlock(
66 implicit=True,
67 description="Arguments passed when constructing StateBlocks.",
68 ),
69 )
70 config.declare(
71 "has_phase_equilibrium",
72 ConfigValue(
73 default=False,
74 domain=Bool,
75 description="Default phase-equilibrium flag for inlet/outlet StateBlocks.",
76 ),
77 )
80class DesuperheaterScaler(CustomScalerBase):
81 """
82 Default modular scaler for the generic unit model.
83 This Scaler relies on the associated property and reaction packages,
84 either through user provided options (submodel_scalers argument) or by default
85 Scalers assigned to the packages.
86 """
88 DEFAULT_SCALING_FACTORS = {
89 "var1": 1e-3,
90 "var2": 1e-3,
91 }
94@declare_process_block_class("Desuperheater")
95class DesuperheaterData(UnitModelBlockData):
96 """Desuperheater unit operation."""
98 default_scaler = DesuperheaterScaler
100 CONFIG = ConfigBlock()
101 _build_config(CONFIG)
103 def build(self):
104 """Build the unit model structure and equations."""
105 super().build()
107 units_meta = self.config.property_package.get_metadata().get_derived_units
109 """
110 1. Build inlet, outlet and internal state blocks and associate
111 with ports (where applicable)
112 """
113 self.inlet_blocks = self._build_state_blocks(
114 stream_name_list=["inlet_vap", "inlet_liq"],
115 has_phase_equilibrium=self.config.has_phase_equilibrium,
116 is_defined_state=True,
117 is_build_port=True,
118 )
119 self.outlet_blocks = self._build_state_blocks(
120 stream_name_list=["outlet"],
121 has_phase_equilibrium=self.config.has_phase_equilibrium,
122 is_defined_state=False,
123 is_build_port=True,
124 )
126 """
127 2. Create parameters, variables, references and expressions
128 """
129 # State variables
130 self.heat_loss = Var(
131 self.flowsheet().time,
132 initialize=0,
133 bounds=(0, None),
134 units=units_meta("power"),
135 doc="Heat loss.",
136 )
137 self.pressure_loss = Var(
138 self.flowsheet().time,
139 initialize=0,
140 bounds=(0, None),
141 units=units_meta("pressure"),
142 doc="Pressure loss.",
143 )
144 self.deltaT_superheat = Var(
145 self.flowsheet().time,
146 initialize=1e-4,
147 bounds=(0, None),
148 units=units_meta("temperature"),
149 doc="Target outlet superheat.",
150 )
151 self.eps = Param(
152 initialize=1e-6,
153 mutable=True,
154 units=units_meta("pressure"),
155 doc="Small number to prevent division by zero in momentum balance.",
156 )
158 # Non-normal state variables, other parameters and expression
159 self.flow_ratio = Expression(
160 self.flowsheet().time,
161 rule=lambda b, t: (
162 b.inlet_liq_state[t].flow_mol / (b.inlet_vap_state[t].flow_mol + 1e-9)
163 ),
164 doc="Ratio of injected water flow to inlet steam flow.",
165 )
167 """
168 3. Declare constraints to define mass, energy, and momentum balances,
169 unit operation performance and other constraints
170 """
171 # a) Material balance equations
172 @self.Constraint(self.flowsheet().time, doc="Overall material balance")
173 def eq_overall_material_balance(b, t):
174 return b.outlet_state[t].flow_mol == b.inlet_vap_state[t].flow_mol + b.inlet_liq_state[t].flow_mol
176 # b) Energy balance equations
177 @self.Constraint(self.flowsheet().time, doc="Overall energy balance")
178 def eq_overall_energy_balance(b, t):
179 return (
180 b.inlet_vap_state[t].flow_mol * b.inlet_vap_state[t].enth_mol + b.inlet_liq_state[t].flow_mol * b.inlet_liq_state[t].enth_mol
181 ==
182 b.outlet_state[t].flow_mol * b.outlet_state[t].enth_mol + b.heat_loss[t]
183 )
185 # c) Momentum balance equations
186 @self.Constraint(self.flowsheet().time, doc="Overall momentum balance")
187 def eq_overall_momentum_balance(b, t):
188 return (
189 smooth_min(b.inlet_vap_state[t].pressure, b.inlet_liq_state[t].pressure, eps=b.eps)
190 ==
191 b.outlet_state[t].pressure + b.pressure_loss[t]
192 )
194 # d) Performance equations and other constraints
195 @self.Constraint(self.flowsheet().time, doc="Outlet superheat relation")
196 def eq_desuperheating_temperature(b, t):
197 return b.outlet_state[t].temperature == b.outlet_state[t].temperature_sat + b.deltaT_superheat[t]
199 self.scaling_factor = Suffix(direction=Suffix.EXPORT)
201 def initialize_build(self, state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None):
202 """
203 General wrapper for template initialization routines
205 Keyword Arguments:
206 state_args : a dict of arguments to be passed to the property
207 package(s) to provide an initial state for
208 initialization (see documentation of the specific
209 property package) (default = {}).
210 outlvl : sets output level of initialization routine
211 optarg : solver options dictionary object (default=None)
212 solver : str indicating which solver to use during
213 initialization (default = None)
215 Returns: None
216 """
217 init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit")
218 solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit")
219 t0 = self.flowsheet().config.time.first() # use first time point for initialisation
221 # Set solver options
222 opt = get_solver(solver, optarg)
224 pp = self.inlet_vap_state[t0].params # proporty package calculation block
225 state_args = {} if state_args is None else dict(state_args)
227 shared_state_args = {
228 key: state_args[key]
229 for key in ("flow_mol", "pressure", "enth_mol", "temperature")
230 if key in state_args
231 }
232 state_args_vap = dict(shared_state_args)
233 state_args_vap.update(state_args.get("inlet_vap", {}))
234 state_args_liq = dict(state_args.get("inlet_liq", {}))
235 state_args_outlet = dict(state_args.get("outlet", {}))
237 # Helper functions
238 def _value_or_none(obj): # returns the value of a Var or Param if it is fixed, otherwise returns None
239 return value(obj, exception=False)
241 def _pick_seed(*candidates): # loops through different options of seeding and picks the first (in order of preference)
242 for candidate in candidates: 242 ↛ 245line 242 didn't jump to line 245 because the loop on line 242 didn't complete
243 if candidate is not None:
244 return candidate
245 return None
247 def _enthalpy_from_tp(temperature, pressure):
248 if temperature is None or pressure is None: 248 ↛ 249line 248 didn't jump to line 249 because the condition on line 248 was never true
249 return None
250 return value(pp.htpx(temperature * pyunits.K, pressure * pyunits.Pa))
252 inlet_vap_has_source = len(list(self.inlet_vap.sources())) > 0
253 inlet_liq_has_source = len(list(self.inlet_liq.sources())) > 0
255 pressure_loss = _value_or_none(self.pressure_loss[t0])
256 deltaT_superheat = _value_or_none(self.deltaT_superheat[t0])
258 #-----------------------------------------------------
259 # 1. Build state seeds. Prefer explicit state_args, then connected/fixed
260 # values already present on the unit, then infer from local specs, then
261 # fall back to generic guesses.
263 # Inlets flow ranking
264 # 1. Explicit state_args
265 # 2. A fixed inlet flow on this unit
266 # 3. A propagated upstream inlet flow
267 # 4. Back-calculate from fixed outlet flow(s) plus either liquid or vap inlet flow
268 # 5. Nominal fallback (could be replaced with enthalpy-based estimate in future)
269 inlet_vap_flow_from_fixed = (
270 _value_or_none(self.inlet_vap_state[t0].flow_mol)
271 if self.inlet_vap_state[t0].flow_mol.fixed
272 else None
273 )
274 inlet_vap_flow_from_upstream = (
275 _value_or_none(self.inlet_vap_state[t0].flow_mol)
276 if inlet_vap_has_source
277 else None
278 )
279 inlet_liq_flow_from_fixed = (
280 _value_or_none(self.inlet_liq_state[t0].flow_mol)
281 if self.inlet_liq_state[t0].flow_mol.fixed
282 else None
283 )
284 inlet_liq_flow_from_upstream = (
285 _value_or_none(self.inlet_liq_state[t0].flow_mol)
286 if inlet_liq_has_source
287 else None
288 )
289 outlet_flow_from_fixed = (
290 _value_or_none(self.outlet_state[t0].flow_mol)
291 if self.outlet_state[t0].flow_mol.fixed
292 else None
293 )
294 inlet_vap_flow_from_outlet = None
295 if outlet_flow_from_fixed is not None and inlet_liq_flow_from_fixed is not None: 295 ↛ 296line 295 didn't jump to line 296 because the condition on line 295 was never true
296 inlet_vap_flow_from_outlet = outlet_flow_from_fixed - inlet_liq_flow_from_fixed
298 f_inlet_vap = _pick_seed(
299 state_args_vap.get("flow_mol"),
300 inlet_vap_flow_from_fixed,
301 inlet_vap_flow_from_upstream,
302 inlet_vap_flow_from_outlet,
303 500.0,
304 )
306 inlet_liq_flow_from_outlet = None
307 if outlet_flow_from_fixed is not None and inlet_vap_flow_from_fixed is not None: 307 ↛ 308line 307 didn't jump to line 308 because the condition on line 307 was never true
308 inlet_liq_flow_from_outlet = outlet_flow_from_fixed - inlet_vap_flow_from_fixed
309 f_inlet_liq = _pick_seed(
310 state_args_liq.get("flow_mol"),
311 inlet_liq_flow_from_fixed,
312 inlet_liq_flow_from_upstream,
313 inlet_liq_flow_from_outlet,
314 0.1 * f_inlet_vap, # TODO: make this an enthalpy based calculation if they are known, may need to be done at the end
315 50.0,
316 )
318 # Pressure ranking
319 # 1. Explicit state_args
320 # 2. A propagated upstream inlet pressure
321 # 3. A fixed inlet pressure on this unit
322 # 4. Assume liq and vapor same as fixed outlet pressure + pressure loss
323 # 5. Nominal fallback
324 inlet_vap_pressure_from_upstream = (
325 _value_or_none(self.inlet_vap_state[t0].pressure)
326 if inlet_vap_has_source
327 else None
328 )
329 inlet_vap_pressure_from_fixed = (
330 _value_or_none(self.inlet_vap_state[t0].pressure)
331 if self.inlet_vap_state[t0].pressure.fixed
332 else None
333 )
334 inlet_vap_pressure_from_outlet = (
335 _value_or_none(self.outlet_state[t0].pressure) + pressure_loss
336 if self.outlet_state[t0].pressure.fixed
337 else None
338 )
339 p_inlet_vap = _pick_seed(
340 state_args_vap.get("pressure"),
341 inlet_vap_pressure_from_upstream,
342 inlet_vap_pressure_from_fixed,
343 inlet_vap_pressure_from_outlet,
344 10e5,
345 )
347 inlet_liq_pressure_from_fixed = (
348 _value_or_none(self.inlet_liq_state[t0].pressure)
349 if self.inlet_liq_state[t0].pressure.fixed
350 else None
351 )
352 inlet_liq_pressure_from_upstream = (
353 _value_or_none(self.inlet_liq_state[t0].pressure)
354 if inlet_liq_has_source
355 else None
356 )
358 p_inlet_liq = _pick_seed(
359 state_args_liq.get("pressure"),
360 inlet_liq_pressure_from_upstream,
361 inlet_liq_pressure_from_fixed,
362 p_inlet_vap,
363 10e5,
364 )
366 # Enthalpy ranking
367 # 1. Explicit state_args
368 # 2. A propagated upstream inlet enthalpy
369 # 3. A fixed inlet enthalpy on this unit
370 # 4. Assume superheat/subcooling on inlet pressure
371 # Set the outlet pressure first so the property block can report the corresponding saturation temperature.
372 seeded_p_outlet = (
373 _value_or_none(self.outlet_state[t0].pressure)
374 if self.outlet_state[t0].pressure.fixed
375 else None
376 )
377 default_p_outlet =min(p_inlet_vap, p_inlet_liq) - pressure_loss
378 p_outlet = _pick_seed(
379 seeded_p_outlet,
380 default_p_outlet
381 )
382 T_sat_outlet = value(self.outlet_state[t0].temperature_sat)
384 inlet_vap_enth_from_upstream = (
385 _value_or_none(self.inlet_vap_state[t0].enth_mol)
386 if inlet_vap_has_source
387 else None
388 )
389 inlet_vap_enth_from_fixed = (
390 _value_or_none(self.inlet_vap_state[t0].enth_mol)
391 if self.inlet_vap_state[t0].enth_mol.fixed
392 else None
393 )
394 h_inlet_vap = _pick_seed(
395 state_args_vap.get("enth_mol"),
396 inlet_vap_enth_from_upstream,
397 inlet_vap_enth_from_fixed,
398 _enthalpy_from_tp(T_sat_outlet + 20.0, p_inlet_vap),
399 )
401 inlet_liq_enth_from_upstream = (
402 _value_or_none(self.inlet_liq_state[t0].enth_mol)
403 if inlet_liq_has_source
404 else None
405 )
406 inlet_liq_enth_from_fixed = (
407 _value_or_none(self.inlet_liq_state[t0].enth_mol)
408 if self.inlet_liq_state[t0].enth_mol.fixed
409 else None
410 )
411 h_inlet_liq = _pick_seed(
412 state_args_liq.get("enth_mol"),
413 inlet_liq_enth_from_upstream,
414 inlet_liq_enth_from_fixed,
415 _enthalpy_from_tp(T_sat_outlet - 50.0, p_inlet_liq),
416 )
418 h_outlet = _pick_seed(
419 state_args_outlet.get("enth_mol"),
420 _enthalpy_from_tp(T_sat_outlet + deltaT_superheat, p_outlet) if deltaT_superheat is not None else None,
421 _enthalpy_from_tp(T_sat_outlet + 5, p_outlet),
422 )
424 #-----------------------------------------------------
425 # 2. Initialize state blocks using explicitly seeded values for all state variables
426 self.inlet_vap_state.initialize(
427 solver=solver,
428 optarg=optarg,
429 outlvl=outlvl,
430 state_args={
431 "flow_mol": f_inlet_vap,
432 "pressure": p_inlet_vap,
433 "enth_mol": h_inlet_vap,
434 },
435 )
436 init_log.info_high("Inlet vapour state initialization complete")
438 self.inlet_liq_state.initialize(
439 solver=solver,
440 optarg=optarg,
441 outlvl=outlvl,
442 state_args={
443 "flow_mol": f_inlet_liq,
444 "pressure": p_inlet_liq,
445 "enth_mol": h_inlet_liq,
446 },
447 )
448 init_log.info_high("Inlet liquid state initialization complete")
450 flags_outlet = self.outlet_state.initialize(
451 solver=solver,
452 optarg=optarg,
453 outlvl=outlvl,
454 state_args={
455 "flow_mol": f_inlet_vap + f_inlet_liq,
456 "pressure": p_outlet,
457 "enth_mol": h_outlet,
458 },
459 hold_state=True,
460 )
461 init_log.info_high("Outlet state initialization complete")
463 # #-----------------------------------------------------
464 # # 3. Deactivate constraints not specified in Step 2 then solve first pass
465 relaxed_eqns = [
466 self.eq_overall_material_balance,
467 self.eq_overall_energy_balance,
468 self.eq_overall_momentum_balance,
469 self.eq_desuperheating_temperature,
470 ]
472 for con in relaxed_eqns:
473 con.deactivate()
475 report_statistics(self)
477 dt = DiagnosticsToolbox(self)
478 dt.report_structural_issues()
479 dt.display_underconstrained_set()
481 active_constraints = sum(
482 1 for _ in self.component_data_objects(Constraint, active=True, descend_into=True)
483 )
485 if active_constraints > 0: 485 ↛ 486line 485 didn't jump to line 486 because the condition on line 485 was never true
486 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
487 res = opt.solve(self, tee=slc.tee)
488 if not check_optimal_termination(res):
489 dt.report_numerical_issues()
490 raise InitializationError(f"{self.name} failed relaxed initialization")
491 else:
492 init_log.info_high("Relaxed initialization pass skipped: no active constraints remained after seeding.")
494 # Restore full model
495 self.outlet_state.release_state(flags_outlet)
497 for con in relaxed_eqns:
498 con.activate()
500 dof = degrees_of_freedom(self)
501 if dof != 0: 501 ↛ 502line 501 didn't jump to line 502 because the condition on line 501 was never true
502 raise InitializationError(
503 f"{self.name} degrees of freedom were not 0 before final solve. DoF = {dof}"
504 )
506 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
507 res = opt.solve(self, tee=slc.tee)
508 if not check_optimal_termination(res): 508 ↛ 509line 508 didn't jump to line 509 because the condition on line 508 was never true
509 raise InitializationError(f"{self.name} failed final initialization")
511 init_log.info(f"Initialization complete: {idaeslog.condition(res)}")
514 def _get_performance_contents(self, time_point=0):
515 var_dict = {
516 "BFW temperature [K]": self.inlet_liq_state[time_point].temperature,
517 "Degree of superheat target [K]": self.deltaT_superheat[time_point],
518 "Heat loss [W]": self.heat_loss[time_point],
519 "Pressure loss [Pa]": self.pressure_loss[time_point],
520 }
521 expr_dict = {
522 "Liquid-to-vapour flow ratio": self.flow_ratio[time_point],
523 }
524 return {"vars": var_dict, "exprs": expr_dict}
526 # -----------------------------------------------------------------
527 # Common utilities
528 # -----------------------------------------------------------------
530 def calculate_scaling_factors(self):
531 super().calculate_scaling_factors()
534 def _build_state_blocks(
535 self,
536 stream_name_list: Iterable[str],
537 has_phase_equilibrium: bool,
538 is_defined_state: Optional[bool] = False,
539 is_build_port: Optional[bool] = False,
540 ) -> List[StateBlock]:
541 blocks: List[StateBlock] = []
543 base_args = dict(self.config.property_package_args)
544 base_args["has_phase_equilibrium"] = has_phase_equilibrium
545 base_args["defined_state"] = is_defined_state
547 for stream_name in stream_name_list:
548 args = dict(base_args)
549 args["doc"] = f"Thermophysical properties at {stream_name}"
550 sb = self.config.property_package.build_state_block(self.flowsheet().time, **args)
551 setattr(self, f"{stream_name}_state", sb)
552 blocks.append(sb)
554 if is_build_port: 554 ↛ 547line 554 didn't jump to line 547 because the condition on line 554 was always true
555 self.add_port(name=stream_name, block=sb)
557 return blocks
560 def _get_stream_table_contents(self, time_point=0):
561 io_dict = {name: getattr(self, name) for name in [*self.inlet_blocks, *self.outlet_blocks]}
562 return create_stream_table_dataframe(io_dict, time_point=time_point)