Coverage for backend/ahuora-builder/src/ahuora_builder/custom/thermal_utility_systems/desuperheater.py: 86%
168 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-05-13 02:47 +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 p_outlet = min(p_inlet_vap, p_inlet_liq) - pressure_loss
373 self.outlet_state[t0].pressure.set_value(p_outlet)
374 T_sat_outlet = value(self.outlet_state[t0].temperature_sat)
376 inlet_vap_enth_from_upstream = (
377 _value_or_none(self.inlet_vap_state[t0].enth_mol)
378 if inlet_vap_has_source
379 else None
380 )
381 inlet_vap_enth_from_fixed = (
382 _value_or_none(self.inlet_vap_state[t0].enth_mol)
383 if self.inlet_vap_state[t0].enth_mol.fixed
384 else None
385 )
386 h_inlet_vap = _pick_seed(
387 state_args_vap.get("enth_mol"),
388 inlet_vap_enth_from_upstream,
389 inlet_vap_enth_from_fixed,
390 _enthalpy_from_tp(T_sat_outlet + 20.0, p_inlet_vap),
391 )
393 inlet_liq_enth_from_upstream = (
394 _value_or_none(self.inlet_liq_state[t0].enth_mol)
395 if inlet_liq_has_source
396 else None
397 )
398 inlet_liq_enth_from_fixed = (
399 _value_or_none(self.inlet_liq_state[t0].enth_mol)
400 if self.inlet_liq_state[t0].enth_mol.fixed
401 else None
402 )
403 h_inlet_liq = _pick_seed(
404 state_args_liq.get("enth_mol"),
405 inlet_liq_enth_from_upstream,
406 inlet_liq_enth_from_fixed,
407 _enthalpy_from_tp(T_sat_outlet - 50.0, p_inlet_liq),
408 )
410 h_outlet = _pick_seed(
411 state_args_outlet.get("enth_mol"),
412 _enthalpy_from_tp(T_sat_outlet + deltaT_superheat, p_outlet) if deltaT_superheat is not None else None,
413 _enthalpy_from_tp(T_sat_outlet + 5, p_outlet),
414 )
416 #-----------------------------------------------------
417 # 2. Initialize state blocks using explicitly seeded values for all state variables
418 self.inlet_vap_state.initialize(
419 solver=solver,
420 optarg=optarg,
421 outlvl=outlvl,
422 state_args={
423 "flow_mol": f_inlet_vap,
424 "pressure": p_inlet_vap,
425 "enth_mol": h_inlet_vap,
426 },
427 )
428 init_log.info_high("Inlet vapour state initialization complete")
430 self.inlet_liq_state.initialize(
431 solver=solver,
432 optarg=optarg,
433 outlvl=outlvl,
434 state_args={
435 "flow_mol": f_inlet_liq,
436 "pressure": p_inlet_liq,
437 "enth_mol": h_inlet_liq,
438 },
439 )
440 init_log.info_high("Inlet liquid state initialization complete")
442 flags_outlet = self.outlet_state.initialize(
443 solver=solver,
444 optarg=optarg,
445 outlvl=outlvl,
446 state_args={
447 "flow_mol": f_inlet_vap + f_inlet_liq,
448 "pressure": p_outlet,
449 "enth_mol": h_outlet,
450 },
451 hold_state=True,
452 )
453 init_log.info_high("Outlet state initialization complete")
455 # #-----------------------------------------------------
456 # # 3. Deactivate constraints not specified in Step 2 then solve first pass
457 relaxed_eqns = [
458 self.eq_overall_material_balance,
459 self.eq_overall_energy_balance,
460 self.eq_overall_momentum_balance,
461 self.eq_desuperheating_temperature,
462 ]
464 for con in relaxed_eqns:
465 con.deactivate()
467 report_statistics(self)
469 dt = DiagnosticsToolbox(self)
470 dt.report_structural_issues()
471 dt.display_underconstrained_set()
473 active_constraints = sum(
474 1 for _ in self.component_data_objects(Constraint, active=True, descend_into=True)
475 )
477 if active_constraints > 0: 477 ↛ 478line 477 didn't jump to line 478 because the condition on line 477 was never true
478 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
479 res = opt.solve(self, tee=slc.tee)
480 if not check_optimal_termination(res):
481 dt.report_numerical_issues()
482 raise InitializationError(f"{self.name} failed relaxed initialization")
483 else:
484 init_log.info_high("Relaxed initialization pass skipped: no active constraints remained after seeding.")
486 # Restore full model
487 self.outlet_state.release_state(flags_outlet)
489 for con in relaxed_eqns:
490 con.activate()
492 dof = degrees_of_freedom(self)
493 if dof != 0: 493 ↛ 494line 493 didn't jump to line 494 because the condition on line 493 was never true
494 raise InitializationError(
495 f"{self.name} degrees of freedom were not 0 before final solve. DoF = {dof}"
496 )
498 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
499 res = opt.solve(self, tee=slc.tee)
500 if not check_optimal_termination(res): 500 ↛ 501line 500 didn't jump to line 501 because the condition on line 500 was never true
501 raise InitializationError(f"{self.name} failed final initialization")
503 init_log.info(f"Initialization complete: {idaeslog.condition(res)}")
506 def _get_performance_contents(self, time_point=0):
507 var_dict = {
508 "BFW temperature [K]": self.inlet_liq_state[time_point].temperature,
509 "Degree of superheat target [K]": self.deltaT_superheat[time_point],
510 "Heat loss [W]": self.heat_loss[time_point],
511 "Pressure loss [Pa]": self.pressure_loss[time_point],
512 }
513 expr_dict = {
514 "Liquid-to-vapour flow ratio": self.flow_ratio[time_point],
515 }
516 return {"vars": var_dict, "exprs": expr_dict}
518 # -----------------------------------------------------------------
519 # Common utilities
520 # -----------------------------------------------------------------
522 def calculate_scaling_factors(self):
523 super().calculate_scaling_factors()
526 def _build_state_blocks(
527 self,
528 stream_name_list: Iterable[str],
529 has_phase_equilibrium: bool,
530 is_defined_state: Optional[bool] = False,
531 is_build_port: Optional[bool] = False,
532 ) -> List[StateBlock]:
533 blocks: List[StateBlock] = []
535 base_args = dict(self.config.property_package_args)
536 base_args["has_phase_equilibrium"] = has_phase_equilibrium
537 base_args["defined_state"] = is_defined_state
539 for stream_name in stream_name_list:
540 args = dict(base_args)
541 args["doc"] = f"Thermophysical properties at {stream_name}"
542 sb = self.config.property_package.build_state_block(self.flowsheet().time, **args)
543 setattr(self, f"{stream_name}_state", sb)
544 blocks.append(sb)
546 if is_build_port: 546 ↛ 539line 546 didn't jump to line 539 because the condition on line 546 was always true
547 self.add_port(name=stream_name, block=sb)
549 return blocks
552 def _get_stream_table_contents(self, time_point=0):
553 io_dict = {name: getattr(self, name) for name in [*self.inlet_blocks, *self.outlet_blocks]}
554 return create_stream_table_dataframe(io_dict, time_point=time_point)