Coverage for backend/ahuora-builder/src/ahuora_builder/custom/thermal_utility_systems/simple_boiler.py: 91%
208 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-06-23 21:51 +0000
1"""Simple boiler unit model."""
3from typing import Iterable, List, Optional
4from pyomo.common.config import Bool, ConfigBlock, ConfigValue, In
5from pyomo.environ import (
6 Constraint,
7 Param,
8 Suffix,
9 Var,
10 check_optimal_termination,
11 value,
12 units as pyunits,
13)
15from idaes.core import StateBlock, UnitModelBlockData, declare_process_block_class, useDefault
16from idaes.core.scaling import CustomScalerBase
17from idaes.core.solvers import get_solver
18from idaes.core.util.config import is_physical_parameter_block
19from idaes.core.util.exceptions import InitializationError
20from idaes.core.util.model_diagnostics import DiagnosticsToolbox
21from idaes.core.util.model_statistics import degrees_of_freedom, report_statistics
22from idaes.core.util.tables import create_stream_table_dataframe
24import idaes.logger as idaeslog
26_log = idaeslog.getLogger(__name__)
28__author__ = "Ahuora Centre for Smart Energy Systems, University of Waikato, New Zealand"
29__Note__ = "When there is insufficient energy to raise feedwater to saturation, the model will fail in initialisation, if not enough energy is available to fully vaporise steam, saturated water will exit in the steam port."
31def _build_config(config: ConfigBlock) -> None:
32 """Declare config entries for SimpleBoiler."""
33 config.declare(
34 "dynamic",
35 ConfigValue(
36 domain=In([False]),
37 default=False,
38 description="Dynamic model flag - must be False",
39 ),
40 )
41 config.declare(
42 "has_holdup",
43 ConfigValue(
44 default=False,
45 domain=In([False]),
46 description="Holdup construction flag - must be False",
47 ),
48 )
49 config.declare(
50 "property_package",
51 ConfigValue(
52 default=useDefault,
53 domain=is_physical_parameter_block,
54 description="Property package used to build StateBlocks.",
55 ),
56 )
57 config.declare(
58 "property_package_args",
59 ConfigBlock(
60 implicit=True,
61 description="Arguments passed when constructing StateBlocks.",
62 ),
63 )
64 config.declare(
65 "has_phase_equilibrium",
66 ConfigValue(
67 default=False,
68 domain=Bool,
69 description="Default phase-equilibrium flag for inlet/outlet StateBlocks.",
70 ),
71 )
74class SimpleBoilerScaler(CustomScalerBase):
75 """Default scaler placeholder for the simple boiler model."""
77 DEFAULT_SCALING_FACTORS = {
78 "fuel_mass_flow": 1,
79 "fuel_calorific_value": 1e-6,
80 "heat_added": 1e-6,
81 }
84@declare_process_block_class("SimpleBoiler")
85class SimpleBoilerData(UnitModelBlockData):
86 """Simple boiler energy-side unit operation."""
88 default_scaler = SimpleBoilerScaler
90 CONFIG = ConfigBlock()
91 _build_config(CONFIG)
93 def build(self):
94 """Build the unit model structure and equations."""
95 super().build()
97 units_meta = self.config.property_package.get_metadata().get_derived_units
99 """
100 1. Build inlet and outlet state blocks and associate
101 with ports (where applicable)
102 """
103 self.inlet_blocks = self._build_state_blocks(
104 stream_name_list=["inlet"],
105 has_phase_equilibrium=self.config.has_phase_equilibrium,
106 is_defined_state=True,
107 is_build_port=True,
108 )
109 self.outlet_blocks = self._build_state_blocks(
110 stream_name_list=["steam", "blowdown"],
111 has_phase_equilibrium=self.config.has_phase_equilibrium,
112 is_defined_state=False,
113 is_build_port=True,
114 )
115 for sb in self.inlet_blocks + self.outlet_blocks:
116 for t in sb:
117 sb[t].flow_mol.setlb(0.0)
119 """
120 2. Create parameters, variables, references and expressions
121 """
122 self.blowdown_fraction = Var(
123 self.flowsheet().time,
124 initialize=0.05,
125 bounds=(0, None),
126 units=pyunits.dimensionless,
127 doc="Blowdown fraction relative to steam flow.",
128 )
130 self.fuel_mass_flow = Var(
131 self.flowsheet().time,
132 initialize=0.01,
133 bounds=(0, None),
134 units=pyunits.kg / pyunits.s,
135 doc="Fuel mass flow to the boiler.",
136 )
137 self.fuel_calorific_value = Var(
138 self.flowsheet().time,
139 initialize=40e6,
140 bounds=(0, None),
141 units=pyunits.J / pyunits.kg,
142 doc="Fuel calorific value.",
143 )
144 self.boiler_efficiency = Var(
145 self.flowsheet().time,
146 initialize=0.85,
147 bounds=(0, 1),
148 units=pyunits.dimensionless,
149 doc="Fraction of fuel energy transferred to the water/steam side.",
150 )
152 self.heat_added = Var(
153 self.flowsheet().time,
154 initialize=3.4e5,
155 bounds=(0, None),
156 units=pyunits.W,
157 doc="Useful heat added to the boiler water/steam side.",
158 )
160 self.feedwater_tds = Param(
161 self.flowsheet().time,
162 initialize=250,
163 mutable=True,
164 units=pyunits.ppm,
165 doc="Total dissolved solids in the boiler feedwater.",
166 )
167 self.max_boiler_tds = Param(
168 self.flowsheet().time,
169 initialize=3500,
170 mutable=True,
171 units=pyunits.ppm,
172 doc="Maximum allowable total dissolved solids in the boiler water.",
173 )
175 """
176 3. Declare constraints to define mass, energy, and momentum balances,
177 unit operation performance and other constraint
178 """
179 # Material balance: inlet splits into steam plus blowdown.
180 @self.Constraint(self.flowsheet().time, doc="Overall material balance")
181 def eq_overall_material_balance(b, t):
182 return (
183 b.inlet_state[t].flow_mol
184 == b.steam_state[t].flow_mol + b.blowdown_state[t].flow_mol
185 )
187 @self.Constraint(self.flowsheet().time, doc="Blowdown flow relation")
188 def eq_blowdown_flow(b, t):
189 return (
190 b.blowdown_state[t].flow_mol
191 == b.blowdown_fraction[t] * b.steam_state[t].flow_mol
192 )
194 @self.Constraint(self.flowsheet().time, doc="Blowdown fraction from TDS")
195 def eq_blowdown_fraction_tds(b, t):
196 feedwater_tds = pyunits.convert(
197 b.feedwater_tds[t],
198 to_units=pyunits.dimensionless,
199 )
200 max_boiler_tds = pyunits.convert(
201 b.max_boiler_tds[t],
202 to_units=pyunits.dimensionless,
203 )
204 return (
205 b.blowdown_fraction[t]
206 == feedwater_tds / (max_boiler_tds - feedwater_tds)
207 )
209 # Energy balance: fuel heat raises the total boiler-side enthalpy.
210 @self.Constraint(self.flowsheet().time, doc="Overall energy balance")
211 def eq_overall_energy_balance(b, t):
212 return (
213 b.inlet_state[t].flow_mol * b.inlet_state[t].enth_mol
214 + b.heat_added[t]
215 == b.steam_state[t].flow_mol * b.steam_state[t].enth_mol
216 + b.blowdown_state[t].flow_mol * b.blowdown_state[t].enth_mol
217 )
219 # Momentum: both outlet streams leave at inlet pressure in this simple model.
220 @self.Constraint(self.flowsheet().time, doc="Steam pressure equality")
221 def eq_steam_pressure(b, t):
222 return b.steam_state[t].pressure == b.inlet_state[t].pressure
224 @self.Constraint(self.flowsheet().time, doc="Blowdown pressure equality")
225 def eq_blowdown_pressure(b, t):
226 return b.blowdown_state[t].pressure == b.inlet_state[t].pressure
228 @self.Constraint(self.flowsheet().time, doc="Blowdown saturated liquid relation")
229 def eq_blowdown_saturated_liquid(b, t):
230 return b.blowdown_state[t].vapor_frac == 0
232 @self.Constraint(self.flowsheet().time, doc="Steam outlet must be at least saturated vapour")
233 def eq_steam_outlet_minimum_vapour_enthalpy(b, t):
234 return b.steam_state[t].enth_mol >= b.steam_state[t].enth_mol_sat_phase["Vap"]
236 # Performance relation tying the useful heat duty back to fuel input.
237 @self.Constraint(self.flowsheet().time, doc="Heat-added definition from fuel input")
238 def eq_heat_added_definition(b, t):
239 return b.heat_added[t] == (
240 b.fuel_mass_flow[t]
241 * b.fuel_calorific_value[t]
242 * b.boiler_efficiency[t]
243 )
244 self.scaling_factor = Suffix(direction=Suffix.EXPORT)
246 def initialize_build(self, state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None):
247 """
248 General wrapper for template initialization routines.
250 Keyword Arguments:
251 state_args : a dict of arguments to be passed to the property package(s)
252 outlvl : sets output level of initialization routine
253 optarg : solver options dictionary object
254 solver : str indicating which solver to use
256 Returns:
257 None
258 """
259 init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit")
260 solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit")
261 t0 = self.flowsheet().config.time.first()
263 opt = get_solver(solver, optarg)
265 pp = self.inlet_state[t0].params
266 state_args = {} if state_args is None else dict(state_args)
268 shared_state_args = {
269 key: state_args[key]
270 for key in ("flow_mol", "pressure", "enth_mol", "temperature")
271 if key in state_args
272 }
273 state_args_inlet = dict(shared_state_args)
274 state_args_inlet.update(state_args.get("inlet", {}))
275 state_args_steam = dict(state_args.get("steam", {}))
276 state_args_blowdown = dict(state_args.get("blowdown", {}))
278 # Step 1: prepare a few local helpers and normalise any caller-supplied state args.
279 # Keeping these utilities close to the initializer makes the seeding rules below
280 # easier to follow without having to jump around the file.
281 def _value_or_none(obj):
282 return value(obj, exception=False)
284 def _pick_seed(*candidates):
285 for candidate in candidates: 285 ↛ 288line 285 didn't jump to line 288 because the loop on line 285 didn't complete
286 if candidate is not None:
287 return candidate
288 return None
290 def _enthalpy_from_tp(temperature, pressure):
291 if temperature is None or pressure is None: 291 ↛ 292line 291 didn't jump to line 292 because the condition on line 291 was never true
292 return None
293 return value(pp.htpx(temperature * pyunits.K, pressure * pyunits.Pa))
295 def _safe_positive(value_, fallback):
296 if value_ is None or value_ <= 1e-9:
297 return fallback
298 return value_
300 def _saturated_phase_enthalpy(state, phase):
301 try:
302 return _value_or_none(state.enth_mol_sat_phase[phase])
303 except (AttributeError, KeyError):
304 return None
306 # Step 2: build starting guesses for the boiler state.
307 # The goal here is not to solve the model yet, only to give the property blocks
308 # physically reasonable values before the first relaxed solve.
309 inlet_has_source = len(list(self.inlet.sources())) > 0
311 # Start with the operating knobs that live on the unit itself.
312 feedwater_tds = value(
313 pyunits.convert(
314 self.feedwater_tds[t0],
315 to_units=pyunits.dimensionless,
316 )
317 )
318 max_boiler_tds = value(
319 pyunits.convert(
320 self.max_boiler_tds[t0],
321 to_units=pyunits.dimensionless,
322 )
323 )
324 blowdown_fraction_from_tds = (
325 feedwater_tds / (max_boiler_tds - feedwater_tds)
326 if max_boiler_tds > feedwater_tds
327 else None
328 )
329 blowdown_fraction = _pick_seed(
330 _value_or_none(self.blowdown_fraction[t0])
331 if self.blowdown_fraction[t0].fixed
332 else None,
333 blowdown_fraction_from_tds,
334 0.05,
335 )
336 self.blowdown_fraction[t0].set_value(blowdown_fraction)
337 fuel_mass_flow = _pick_seed(_value_or_none(self.fuel_mass_flow[t0]), 0.01)
338 fuel_calorific_value = _pick_seed(
339 _value_or_none(self.fuel_calorific_value[t0]),
340 40e6,
341 )
342 boiler_efficiency = _pick_seed(_value_or_none(self.boiler_efficiency[t0]), 0.85)
343 heat_added = fuel_mass_flow * fuel_calorific_value * boiler_efficiency
345 # Flow seeding:
346 # 1. Explicit state_args
347 # 2. Fixed values already present on the unit
348 # 3. Values propagated from connected upstream units
349 # 4. A simple internal estimate based on the blowdown-to-steam ratio
350 # 5. A conservative fallback if nothing else is available
351 inlet_flow_from_fixed = (
352 _value_or_none(self.inlet_state[t0].flow_mol)
353 if self.inlet_state[t0].flow_mol.fixed
354 else None
355 )
356 inlet_flow_from_upstream = (
357 _value_or_none(self.inlet_state[t0].flow_mol)
358 if inlet_has_source
359 else None
360 )
361 steam_flow_from_fixed = (
362 _value_or_none(self.steam_state[t0].flow_mol)
363 if self.steam_state[t0].flow_mol.fixed
364 else None
365 )
366 blowdown_flow_from_fixed = (
367 _value_or_none(self.blowdown_state[t0].flow_mol)
368 if self.blowdown_state[t0].flow_mol.fixed
369 else None
370 )
371 inlet_flow_from_steam = None
372 if steam_flow_from_fixed is not None:
373 inlet_flow_from_steam = steam_flow_from_fixed * (1 + blowdown_fraction)
374 inlet_flow_from_blowdown = None
375 if blowdown_flow_from_fixed is not None and abs(blowdown_fraction) >= 1e-8: 375 ↛ 376line 375 didn't jump to line 376 because the condition on line 375 was never true
376 inlet_flow_from_blowdown = (
377 blowdown_flow_from_fixed * (1 + blowdown_fraction) / blowdown_fraction
378 )
380 f_inlet = _pick_seed(
381 state_args_inlet.get("flow_mol"),
382 inlet_flow_from_fixed,
383 inlet_flow_from_upstream,
384 inlet_flow_from_steam,
385 inlet_flow_from_blowdown,
386 500.0,
387 )
389 # The outlet flows are split from the inlet flow using B = ratio * steam.
390 f_blowdown = _pick_seed(
391 state_args_blowdown.get("flow_mol"),
392 blowdown_flow_from_fixed,
393 f_inlet * blowdown_fraction / (1 + blowdown_fraction),
394 )
395 # Steam is the remaining flow after removing blowdown.
396 f_steam = _pick_seed(
397 state_args_steam.get("flow_mol"),
398 steam_flow_from_fixed,
399 f_inlet - f_blowdown,
400 0.95 * f_inlet,
401 )
403 # Pressure seeding:
404 # This simple boiler keeps all three streams at the same pressure.
405 inlet_pressure_from_fixed = (
406 _value_or_none(self.inlet_state[t0].pressure)
407 if self.inlet_state[t0].pressure.fixed
408 else None
409 )
410 inlet_pressure_from_upstream = (
411 _value_or_none(self.inlet_state[t0].pressure)
412 if inlet_has_source
413 else None
414 )
415 steam_pressure_from_fixed = (
416 _value_or_none(self.steam_state[t0].pressure)
417 if self.steam_state[t0].pressure.fixed
418 else None
419 )
420 blowdown_pressure_from_fixed = (
421 _value_or_none(self.blowdown_state[t0].pressure)
422 if self.blowdown_state[t0].pressure.fixed
423 else None
424 )
425 p_inlet = _pick_seed(
426 state_args_inlet.get("pressure"),
427 inlet_pressure_from_fixed,
428 inlet_pressure_from_upstream,
429 steam_pressure_from_fixed,
430 blowdown_pressure_from_fixed,
431 10e5,
432 )
433 p_steam = _pick_seed(state_args_steam.get("pressure"), steam_pressure_from_fixed, p_inlet)
434 p_blowdown = _pick_seed(
435 state_args_blowdown.get("pressure"),
436 blowdown_pressure_from_fixed,
437 p_inlet,
438 )
440 # Enthalpy seeding:
441 # Inlet enthalpy comes from explicit or propagated state data first, then from
442 # a T/P estimate. Blowdown is biased toward a saturated-liquid guess, and steam
443 # is back-calculated from the boiler-side energy balance.
444 inlet_enthalpy_from_fixed = (
445 _value_or_none(self.inlet_state[t0].enth_mol)
446 if self.inlet_state[t0].enth_mol.fixed
447 else None
448 )
449 inlet_enthalpy_from_upstream = (
450 _value_or_none(self.inlet_state[t0].enth_mol)
451 if inlet_has_source
452 else None
453 )
454 inlet_temperature_from_state = (
455 _value_or_none(self.inlet_state[t0].temperature)
456 if inlet_has_source or self.inlet_state[t0].enth_mol.fixed
457 else None
458 )
459 h_inlet = _pick_seed(
460 state_args_inlet.get("enth_mol"),
461 inlet_enthalpy_from_fixed,
462 inlet_enthalpy_from_upstream,
463 _enthalpy_from_tp(inlet_temperature_from_state, p_inlet),
464 _enthalpy_from_tp(298.15, p_inlet),
465 )
467 # Blowdown is treated as saturated liquid at the outlet pressure.
468 self.blowdown_state[t0].pressure.set_value(p_blowdown)
469 blowdown_sat_liq_enth = _saturated_phase_enthalpy(
470 self.blowdown_state[t0],
471 "Liq",
472 )
473 blowdown_sat_temp = _value_or_none(
474 getattr(self.blowdown_state[t0], "temperature_sat", None)
475 )
476 h_blowdown = _pick_seed(
477 state_args_blowdown.get("enth_mol"),
478 (
479 _value_or_none(self.blowdown_state[t0].enth_mol)
480 if self.blowdown_state[t0].enth_mol.fixed
481 else None
482 ),
483 blowdown_sat_liq_enth,
484 _enthalpy_from_tp(blowdown_sat_temp - 1e-4, p_blowdown)
485 if blowdown_sat_temp is not None
486 else None,
487 h_inlet,
488 )
490 # Steam is the energy-balance remainder after the blowdown enthalpy is fixed.
491 steam_enthalpy_from_fixed = (
492 _value_or_none(self.steam_state[t0].enth_mol)
493 if self.steam_state[t0].enth_mol.fixed
494 else None
495 )
496 f_steam_safe = _safe_positive(f_steam, 1.0)
497 h_steam_from_balance = (
498 (f_inlet * h_inlet + heat_added - f_blowdown * h_blowdown)
499 / f_steam_safe
500 )
501 self.steam_state[t0].pressure.set_value(p_steam)
502 steam_sat_temp = _value_or_none(getattr(self.steam_state[t0], "temperature_sat", None))
503 h_steam = _pick_seed(
504 state_args_steam.get("enth_mol"),
505 steam_enthalpy_from_fixed,
506 h_steam_from_balance,
507 _enthalpy_from_tp(steam_sat_temp + 10.0, p_steam)
508 if steam_sat_temp is not None
509 else None,
510 )
512 # Step 3: initialize the property blocks with the guesses above.
513 # The inlet is released immediately, while the outlet states are held so the
514 # relaxed unit solve can reconcile the model equations against the seeds.
515 self.inlet_state.initialize(
516 solver=solver,
517 optarg=optarg,
518 outlvl=outlvl,
519 state_args={
520 "flow_mol": f_inlet,
521 "pressure": p_inlet,
522 "enth_mol": h_inlet,
523 },
524 )
525 init_log.info_high("Inlet state initialization complete")
527 flags_steam = self.steam_state.initialize(
528 solver=solver,
529 optarg=optarg,
530 outlvl=outlvl,
531 state_args={
532 "flow_mol": f_steam,
533 "pressure": p_steam,
534 "enth_mol": h_steam,
535 },
536 hold_state=True,
537 )
538 init_log.info_high("Steam state initialization complete")
540 flags_blowdown = self.blowdown_state.initialize(
541 solver=solver,
542 optarg=optarg,
543 outlvl=outlvl,
544 state_args={
545 "flow_mol": f_blowdown,
546 "pressure": p_blowdown,
547 "enth_mol": h_blowdown,
548 },
549 hold_state=True,
550 )
551 init_log.info_high("Blowdown state initialization complete")
553 # Step 4: temporarily relax the unit equations so the first solve is easy.
554 # This mirrors the other thermal utility units: seed the states first, then let
555 # the solver clean up the structure before we restore the full model.
556 relaxed_eqns = [
557 self.eq_overall_material_balance,
558 self.eq_blowdown_flow,
559 self.eq_blowdown_fraction_tds,
560 self.eq_overall_energy_balance,
561 self.eq_steam_pressure,
562 self.eq_blowdown_pressure,
563 self.eq_blowdown_saturated_liquid,
564 self.eq_steam_outlet_minimum_vapour_enthalpy,
565 ]
567 for con in relaxed_eqns:
568 con.deactivate()
570 report_statistics(self)
572 dt = DiagnosticsToolbox(self)
573 dt.report_structural_issues()
574 dt.display_underconstrained_set()
576 active_constraints = sum(
577 1 for _ in self.component_data_objects(Constraint, active=True, descend_into=True)
578 )
580 if active_constraints > 0: 580 ↛ 587line 580 didn't jump to line 587 because the condition on line 580 was always true
581 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
582 res = opt.solve(self, tee=slc.tee)
583 if not check_optimal_termination(res): 583 ↛ 584line 583 didn't jump to line 584 because the condition on line 583 was never true
584 dt.report_numerical_issues()
585 raise InitializationError(f"{self.name} failed relaxed initialization")
586 else:
587 init_log.info_high(
588 "Relaxed initialization pass skipped: no active constraints remained after seeding."
589 )
591 # Step 5: restore the full model and solve the complete unit.
592 # At this point the state blocks have been seeded and the unit should be square.
593 self.steam_state.release_state(flags_steam)
594 self.blowdown_state.release_state(flags_blowdown)
596 for con in relaxed_eqns:
597 con.activate()
599 dof = degrees_of_freedom(self)
600 if dof != 0: 600 ↛ 601line 600 didn't jump to line 601 because the condition on line 600 was never true
601 raise InitializationError(
602 f"{self.name} degrees of freedom were not 0 before final solve. DoF = {dof}"
603 )
605 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
606 res = opt.solve(self, tee=slc.tee)
607 if not check_optimal_termination(res):
608 raise InitializationError(f"{self.name} failed final initialization")
610 init_log.info(f"Initialization complete: {idaeslog.condition(res)}")
612 def _get_performance_contents(self, time_point=0):
613 """Collect performance variables for reporting."""
614 var_dict = {
615 "Blowdown fraction": self.blowdown_fraction[time_point],
616 "Fuel mass flow": self.fuel_mass_flow[time_point],
617 "Fuel calorific value": self.fuel_calorific_value[time_point],
618 "Boiler efficiency": self.boiler_efficiency[time_point],
619 }
620 expr_dict = {
621 "Feedwater TDS": self.feedwater_tds[time_point],
622 "Maximum boiler TDS": self.max_boiler_tds[time_point],
623 "Heat added": self.heat_added[time_point],
624 }
625 return {"vars": var_dict, "exprs": expr_dict}
627 def calculate_scaling_factors(self):
628 super().calculate_scaling_factors()
630 def _build_state_blocks(
631 self,
632 stream_name_list: Iterable[str],
633 has_phase_equilibrium: bool,
634 is_defined_state: Optional[bool] = False,
635 is_build_port: Optional[bool] = False,
636 ) -> List[StateBlock]:
637 blocks: List[StateBlock] = []
639 base_args = dict(self.config.property_package_args)
640 base_args["has_phase_equilibrium"] = has_phase_equilibrium
641 base_args["defined_state"] = is_defined_state
643 for stream_name in stream_name_list:
644 args = dict(base_args)
645 args["doc"] = f"Thermophysical properties at {stream_name}"
646 sb = self.config.property_package.build_state_block(self.flowsheet().time, **args)
647 setattr(self, f"{stream_name}_state", sb)
648 blocks.append(sb)
650 if is_build_port: 650 ↛ 643line 650 didn't jump to line 643 because the condition on line 650 was always true
651 self.add_port(name=stream_name, block=sb)
653 return blocks
655 def _get_stream_table_contents(self, time_point=0):
656 io_dict = {
657 name: getattr(self, name)
658 for name in ["inlet", "steam", "blowdown"]
659 }
660 return create_stream_table_dataframe(io_dict, time_point=time_point)