Coverage for backend/ahuora-builder/src/ahuora_builder/custom/updated_pressure_changer.py: 41%
431 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-05-13 02:47 +0000
1# THis is unncecessary when https://github.com/IDAES/idaes-pse/pull/1556/ is merged
2# and in the production of the idaes release
3#################################################################################
4# The Institute for the Design of Advanced Energy Systems Integrated Platform
5# Framework (IDAES IP) was produced under the DOE Institute for the
6# Design of Advanced Energy Systems (IDAES).
7#
8# Copyright (c) 2018-2024 by the software owners: The Regents of the
9# University of California, through Lawrence Berkeley National Laboratory,
10# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
11# University, West Virginia University Research Corporation, et al.
12# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
13# for full copyright and license information.
14#################################################################################
15"""
16Standard IDAES pressure changer model.
17"""
18# TODO: Missing docstrings
19# pylint: disable=missing-function-docstring
21# Changing existing config block attributes
22# pylint: disable=protected-access
24# Import Python libraries
25from enum import Enum
27# Import Pyomo libraries
28from pyomo.environ import (
29 Block,
30 value,
31 Var,
32 Expression,
33 Constraint,
34 Reference,
35 check_optimal_termination,
36 Reals,
37)
38from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool
40# Import IDAES cores
41from idaes.core import (
42 ControlVolume0DBlock,
43 declare_process_block_class,
44 EnergyBalanceType,
45 MomentumBalanceType,
46 MaterialBalanceType,
47 ProcessBlockData,
48 UnitModelBlockData,
49 useDefault,
50)
51from idaes.core.util.exceptions import PropertyNotSupportedError, InitializationError
52from idaes.core.util.config import is_physical_parameter_block
53import idaes.logger as idaeslog
54from idaes.core.util import scaling as iscale
55from idaes.core.solvers import get_solver
56from idaes.core.initialization import SingleControlVolumeUnitInitializer
57from idaes.core.util import to_json, from_json, StoreSpec
58from ahuora_builder.state_args import extract_state_args
61__author__ = "Emmanuel Ogbe, Andrew Lee"
62_log = idaeslog.getLogger(__name__)
65class ThermodynamicAssumption(Enum):
66 """
67 Enum of supported thermodynamic assumptions.
68 """
70 isothermal = 1
71 isentropic = 2
72 pump = 3
73 adiabatic = 4
76class IsentropicPressureChangerInitializer(SingleControlVolumeUnitInitializer):
77 """
78 Initializer for isentropic pressure changer models (and derived types).
80 """
82 def initialization_routine(
83 self,
84 model: Block,
85 ):
86 """
87 Initialization routine for isentropic pressure changers.
89 This routine starts by initializing the inlet and outlet states as usual,
90 using the user provided operating conditions to estimate the outlet state.
91 The isentropic state is then initialized at the same conditions as the outlet.
92 Next, the pressure changer is solved with an isothermal assumption and fixed efficiency,
93 followed by a second solve with the isentropic constraints. Finally, if user-provided
94 performance constraints are present, these are activated and the model solved again.
96 Args:
97 model: model to be initialized
99 Returns:
100 Pyomo solver status object
102 """
103 init_log = idaeslog.getInitLogger(
104 model.name, self.get_output_level(), tag="unit"
105 )
106 solve_log = idaeslog.getSolveLogger(
107 model.name, self.get_output_level(), tag="unit"
108 )
110 # Create solver
111 solver = self._get_solver()
113 cv = model.control_volume
115 # check if performance curves exist and are active
116 activate_performance_curves = (
117 hasattr(model, "performance_curve")
118 and model.performance_curve.has_constraints()
119 and model.performance_curve.active
120 )
121 if activate_performance_curves:
122 model.performance_curve.deactivate()
123 # The performance curves will provide (maybe indirectly) efficiency
124 # and/or pressure ratio. To get through the standard isentropic
125 # pressure changer init, we'll see if the user provided a guess for
126 # pressure ratio or isentropic efficiency and fix them if needed. If
127 # not fixed and no guess provided, fill in something reasonable
128 # until the performance curves are turned on.
129 unfix_eff = {}
130 unfix_ratioP = {}
131 for t in model.flowsheet().time:
132 if not (
133 model.ratioP[t].fixed
134 or model.deltaP[t].fixed
135 or cv.properties_out[t].pressure.fixed
136 ):
137 if model.config.compressor:
138 if not (
139 value(model.ratioP[t]) >= 1.01
140 and value(model.ratioP[t]) <= 50
141 ):
142 model.ratioP[t] = 1.8
143 else:
144 if not (
145 value(model.ratioP[t]) >= 0.01
146 and value(model.ratioP[t]) <= 0.999
147 ):
148 model.ratioP[t] = 0.7
149 model.ratioP[t].fix()
150 unfix_ratioP[t] = True
151 if not model.efficiency_isentropic[t].fixed:
152 if not (
153 value(model.efficiency_isentropic[t]) >= 0.05
154 and value(model.efficiency_isentropic[t]) <= 1.0
155 ):
156 model.efficiency_isentropic[t] = 0.8
157 model.efficiency_isentropic[t].fix()
158 unfix_eff[t] = True
160 # Initialize state blocks
161 self.initialize_control_volume(cv, copy_inlet_state=False)
163 init_log.info_high("Initialization Step 1 Complete.")
164 # ---------------------------------------------------------------------
165 # Copy outlet state to isentropic state
167 # Map solution from inlet properties to outlet properties
168 state = to_json(
169 cv.properties_out,
170 wts=StoreSpec().value(),
171 return_dict=True,
172 )
173 from_json(
174 model.properties_isentropic,
175 sd=state,
176 wts=StoreSpec().value(only_not_fixed=True),
177 )
179 init_log.info_high("Initialization Step 2 Complete.")
181 # ---------------------------------------------------------------------
182 # Solve for isothermal conditions
183 if isinstance(
184 model.properties_isentropic[model.flowsheet().time.first()].temperature,
185 Var,
186 ):
187 model.properties_isentropic[:].temperature.fix()
188 elif isinstance(
189 model.properties_isentropic[model.flowsheet().time.first()].enth_mol,
190 Var,
191 ):
192 model.properties_isentropic[:].enth_mol.fix()
193 elif isinstance(
194 model.properties_isentropic[model.flowsheet().time.first()].temperature,
195 Expression,
196 ):
198 def tmp_rule(blk, t):
199 return (
200 blk.properties_isentropic[t].temperature
201 == blk.control_volume.properties_in[t].temperature
202 )
204 model.tmp_init_constraint = Constraint(
205 model.flowsheet().time, rule=tmp_rule
206 )
208 model.isentropic.deactivate()
210 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
211 res = solver.solve(model, tee=slc.tee)
212 init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res)))
214 # Revert changes for isothermal assumption
215 if isinstance(
216 model.properties_isentropic[model.flowsheet().time.first()].temperature,
217 Var,
218 ):
219 model.properties_isentropic[:].temperature.unfix()
220 elif isinstance(
221 model.properties_isentropic[model.flowsheet().time.first()].enth_mol,
222 Var,
223 ):
224 model.properties_isentropic[:].enth_mol.unfix()
225 elif isinstance(
226 model.properties_isentropic[model.flowsheet().time.first()].temperature,
227 Expression,
228 ):
229 model.del_component(model.tmp_init_constraint)
231 model.isentropic.activate()
233 # ---------------------------------------------------------------------
234 # Solve unit
235 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
236 res = solver.solve(model, tee=slc.tee)
237 init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res)))
239 if activate_performance_curves:
240 model.performance_curve.activate()
241 for t, v in unfix_eff.items():
242 if v:
243 model.efficiency_isentropic[t].unfix()
244 for t, v in unfix_ratioP.items():
245 if v:
246 model.ratioP[t].unfix()
247 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
248 res = solver.solve(model, tee=slc.tee)
249 init_log.info_high(f"Initialization Step 5 {idaeslog.condition(res)}.")
251 return res
254@declare_process_block_class("IsentropicPerformanceCurve")
255class IsentropicPerformanceCurveData(ProcessBlockData):
256 """Block that holds performance curves. Typically, these are in the form of
257 constraints that relate head, efficiency, or pressure ratio to volumetric
258 or mass flow. Additional variables can be included if needed, such as
259 speed. For convenience an option is provided to add head expressions to the
260 block. performance curves, and any additional variables, constraints, or
261 expressions can be added to this block either via callback provided to the
262 configuration, or after the model is constructed."""
264 CONFIG = ProcessBlockData.CONFIG(
265 doc="Configuration dictionary for the performance curve block."
266 )
267 CONFIG.declare(
268 "build_callback",
269 ConfigValue(
270 default=None, doc="Optional callback to add performance curve constraints"
271 ),
272 )
273 CONFIG.declare(
274 "build_head_expressions",
275 ConfigValue(
276 default=True,
277 domain=bool,
278 doc="If true add expressions for 'head' and 'head_isentropic'. "
279 "These expressions can be used in performance curve constraints.",
280 ),
281 )
283 def has_constraints(self):
284 for _ in self.component_data_objects(Constraint):
285 return True
286 return False
288 def build(self):
289 super().build()
290 if self.config.build_head_expressions:
291 try:
293 @self.Expression(self.flowsheet().time)
294 def head_isentropic(self, t): # units are energy/mass
295 b = self.parent_block()
296 if hasattr(b.control_volume.properties_in[t], "flow_mass"):
297 return (
298 b.work_isentropic[t]
299 / b.control_volume.properties_in[t].flow_mass
300 )
301 else:
302 return (
303 b.work_isentropic[t]
304 / b.control_volume.properties_in[t].flow_mol
305 / b.control_volume.properties_in[t].mw
306 )
308 @self.Expression(self.flowsheet().time)
309 def head(self, t): # units are energy/mass
310 b = self.parent_block()
311 if hasattr(b.control_volume.properties_in[t], "flow_mass"):
312 return (
313 b.work_mechanical[t]
314 / b.control_volume.properties_in[t].flow_mass
315 )
316 else:
317 return (
318 b.work_mechanical[t]
319 / b.control_volume.properties_in[t].flow_mol
320 / b.control_volume.properties_in[t].mw
321 )
323 except PropertyNotSupportedError:
324 _log.exception(
325 "flow_mass or flow_mol and mw are not supported by the "
326 "property package but are required for isentropic pressure"
327 " changer head calculation"
328 )
329 raise
331 if self.config.build_callback is not None:
332 self.config.build_callback(self)
335@declare_process_block_class("PressureChanger")
336class PressureChangerData(UnitModelBlockData):
337 """
338 Standard Compressor/Expander Unit Model Class
339 """
341 CONFIG = UnitModelBlockData.CONFIG()
343 CONFIG.declare(
344 "material_balance_type",
345 ConfigValue(
346 default=MaterialBalanceType.useDefault,
347 domain=In(MaterialBalanceType),
348 description="Material balance construction flag",
349 doc="""Indicates what type of mass balance should be constructed,
350**default** - MaterialBalanceType.useDefault.
351**Valid values:** {
352**MaterialBalanceType.useDefault - refer to property package for default
353balance type
354**MaterialBalanceType.none** - exclude material balances,
355**MaterialBalanceType.componentPhase** - use phase component balances,
356**MaterialBalanceType.componentTotal** - use total component balances,
357**MaterialBalanceType.elementTotal** - use total element balances,
358**MaterialBalanceType.total** - use total material balance.}""",
359 ),
360 )
361 CONFIG.declare(
362 "energy_balance_type",
363 ConfigValue(
364 default=EnergyBalanceType.useDefault,
365 domain=In(EnergyBalanceType),
366 description="Energy balance construction flag",
367 doc="""Indicates what type of energy balance should be constructed,
368**default** - EnergyBalanceType.useDefault.
369**Valid values:** {
370**EnergyBalanceType.useDefault - refer to property package for default
371balance type
372**EnergyBalanceType.none** - exclude energy balances,
373**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material,
374**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase,
375**EnergyBalanceType.energyTotal** - single energy balance for material,
376**EnergyBalanceType.energyPhase** - energy balances for each phase.}""",
377 ),
378 )
379 CONFIG.declare(
380 "momentum_balance_type",
381 ConfigValue(
382 default=MomentumBalanceType.pressureTotal,
383 domain=In(MomentumBalanceType),
384 description="Momentum balance construction flag",
385 doc="""Indicates what type of momentum balance should be
386constructed, **default** - MomentumBalanceType.pressureTotal.
387**Valid values:** {
388**MomentumBalanceType.none** - exclude momentum balances,
389**MomentumBalanceType.pressureTotal** - single pressure balance for material,
390**MomentumBalanceType.pressurePhase** - pressure balances for each phase,
391**MomentumBalanceType.momentumTotal** - single momentum balance for material,
392**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""",
393 ),
394 )
395 CONFIG.declare(
396 "has_phase_equilibrium",
397 ConfigValue(
398 default=False,
399 domain=Bool,
400 description="Phase equilibrium construction flag",
401 doc="""Indicates whether terms for phase equilibrium should be
402constructed, **default** = False.
403**Valid values:** {
404**True** - include phase equilibrium terms
405**False** - exclude phase equilibrium terms.}""",
406 ),
407 )
408 CONFIG.declare(
409 "compressor",
410 ConfigValue(
411 default=True,
412 domain=Bool,
413 description="Compressor flag",
414 doc="""Indicates whether this unit should be considered a
415 compressor (True (default), pressure increase) or an expander
416 (False, pressure decrease).""",
417 ),
418 )
419 CONFIG.declare(
420 "thermodynamic_assumption",
421 ConfigValue(
422 default=ThermodynamicAssumption.isothermal,
423 domain=In(ThermodynamicAssumption),
424 description="Thermodynamic assumption to use",
425 doc="""Flag to set the thermodynamic assumption to use for the unit.
426 - ThermodynamicAssumption.isothermal (default)
427 - ThermodynamicAssumption.isentropic
428 - ThermodynamicAssumption.pump
429 - ThermodynamicAssumption.adiabatic""",
430 ),
431 )
432 CONFIG.declare(
433 "property_package",
434 ConfigValue(
435 default=useDefault,
436 domain=is_physical_parameter_block,
437 description="Property package to use for control volume",
438 doc="""Property parameter object used to define property
439calculations, **default** - useDefault.
440**Valid values:** {
441**useDefault** - use default package from parent model or flowsheet,
442**PropertyParameterObject** - a PropertyParameterBlock object.}""",
443 ),
444 )
445 CONFIG.declare(
446 "property_package_args",
447 ConfigBlock(
448 implicit=True,
449 description="Arguments to use for constructing property packages",
450 doc="""A ConfigBlock with arguments to be passed to a property
451block(s) and used when constructing these,
452**default** - None.
453**Valid values:** {
454see property package for documentation.}""",
455 ),
456 )
457 CONFIG.declare(
458 "support_isentropic_performance_curves",
459 ConfigValue(
460 default=False,
461 domain=Bool,
462 doc="Include a block for performance curves, configure via"
463 " isentropic_performance_curves.",
464 ),
465 )
466 CONFIG.declare(
467 "isentropic_performance_curves",
468 IsentropicPerformanceCurveData.CONFIG(),
469 # doc included in IsentropicPerformanceCurveData
470 )
472 def build(self):
473 """
475 Args:
476 None
478 Returns:
479 None
480 """
481 # Call UnitModel.build
482 super().build()
484 # Add a control volume to the unit including setting up dynamics.
485 self.control_volume = ControlVolume0DBlock(
486 dynamic=self.config.dynamic,
487 has_holdup=self.config.has_holdup,
488 property_package=self.config.property_package,
489 property_package_args=self.config.property_package_args,
490 )
492 # Add geometry variables to control volume
493 if self.config.has_holdup: 493 ↛ 494line 493 didn't jump to line 494 because the condition on line 493 was never true
494 self.control_volume.add_geometry()
496 # Add inlet and outlet state blocks to control volume
497 self.control_volume.add_state_blocks(
498 has_phase_equilibrium=self.config.has_phase_equilibrium
499 )
501 # Add mass balance
502 # Set has_equilibrium is False for now
503 # TO DO; set has_equilibrium to True
504 self.control_volume.add_material_balances(
505 balance_type=self.config.material_balance_type,
506 has_phase_equilibrium=self.config.has_phase_equilibrium,
507 )
509 # Add energy balance
510 eb = self.control_volume.add_energy_balances(
511 balance_type=self.config.energy_balance_type, has_work_transfer=True
512 )
514 # add momentum balance
515 self.control_volume.add_momentum_balances(
516 balance_type=self.config.momentum_balance_type, has_pressure_change=True
517 )
519 # Add Ports
520 self.add_inlet_port()
521 self.add_outlet_port()
523 # Set Unit Geometry and holdup Volume
524 if self.config.has_holdup is True: 524 ↛ 525line 524 didn't jump to line 525 because the condition on line 524 was never true
525 self.volume = Reference(self.control_volume.volume[:])
527 # Construct performance equations
528 # Set references to balance terms at unit level
529 # Add Work transfer variable 'work'
530 # If the 'work' variable wasn't already built on the control volume but is needed, create it now.
531 if ( 531 ↛ 536line 531 didn't jump to line 536 because the condition on line 531 was never true
532 not hasattr(self.control_volume, "work")
533 and self.config.thermodynamic_assumption == ThermodynamicAssumption.pump
534 and eb is None
535 ):
536 units = (
537 self.control_volume.config.property_package.get_metadata().get_derived_units
538 )
539 self.control_volume.work = Var(
540 self.flowsheet().time,
541 domain=Reals,
542 initialize=0.0,
543 doc="Work transferred into control volume",
544 units=units("power"),
545 )
546 self.work_mechanical = Reference(self.control_volume.work[:])
548 # Add Momentum balance variable 'deltaP'
549 self.deltaP = Reference(self.control_volume.deltaP[:])
551 # Performance Variables
552 self.ratioP = Var(self.flowsheet().time, initialize=1.0, doc="Pressure Ratio")
554 # Pressure Ratio
555 @self.Constraint(self.flowsheet().time, doc="Pressure ratio constraint")
556 def ratioP_calculation(self, t):
557 return (
558 self.ratioP[t] * self.control_volume.properties_in[t].pressure
559 == self.control_volume.properties_out[t].pressure
560 )
562 # Construct equations for thermodynamic assumption
563 if self.config.thermodynamic_assumption == ThermodynamicAssumption.isothermal: 563 ↛ 564line 563 didn't jump to line 564 because the condition on line 563 was never true
564 self.add_isothermal()
565 elif self.config.thermodynamic_assumption == ThermodynamicAssumption.isentropic:
566 self.add_isentropic()
567 elif self.config.thermodynamic_assumption == ThermodynamicAssumption.pump:
568 self.add_pump()
569 elif self.config.thermodynamic_assumption == ThermodynamicAssumption.adiabatic: 569 ↛ exitline 569 didn't return from function 'build' because the condition on line 569 was always true
570 self.add_adiabatic()
572 def add_pump(self):
573 """
574 Add constraints for the incompressible fluid assumption
576 Args:
577 None
579 Returns:
580 None
581 """
582 units_meta = self.control_volume.config.property_package.get_metadata()
584 self.work_fluid = Var(
585 self.flowsheet().time,
586 initialize=1.0,
587 doc="Work required to increase the pressure of the liquid",
588 units=units_meta.get_derived_units("power"),
589 )
590 self.efficiency_pump = Var(
591 self.flowsheet().time, initialize=1.0, doc="Pump efficiency"
592 )
594 @self.Constraint(self.flowsheet().time, doc="Pump fluid work constraint")
595 def fluid_work_calculation(self, t):
596 return self.work_fluid[t] == (
597 (
598 self.control_volume.properties_out[t].pressure
599 - self.control_volume.properties_in[t].pressure
600 )
601 * self.control_volume.properties_out[t].flow_vol
602 )
604 # Actual work
605 @self.Constraint(
606 self.flowsheet().time, doc="Actual mechanical work calculation"
607 )
608 def actual_work(self, t):
609 if self.config.compressor: 609 ↛ 614line 609 didn't jump to line 614 because the condition on line 609 was always true
610 return self.work_fluid[t] == (
611 self.work_mechanical[t] * self.efficiency_pump[t]
612 )
613 else:
614 return self.work_mechanical[t] == (
615 self.work_fluid[t] * self.efficiency_pump[t]
616 )
618 def add_isothermal(self):
619 """
620 Add constraints for isothermal assumption.
622 Args:
623 None
625 Returns:
626 None
627 """
629 # Isothermal constraint
630 @self.Constraint(
631 self.flowsheet().time,
632 doc="For isothermal condition: Equate inlet and " "outlet temperature",
633 )
634 def isothermal(self, t):
635 return (
636 self.control_volume.properties_in[t].temperature
637 == self.control_volume.properties_out[t].temperature
638 )
640 def add_adiabatic(self):
641 """
642 Add constraints for adiabatic assumption.
644 Args:
645 None
647 Returns:
648 None
649 """
651 @self.Constraint(self.flowsheet().time)
652 def zero_work_equation(self, t):
653 return self.control_volume.work[t] == 0
655 def add_isentropic(self):
656 """
657 Add constraints for isentropic assumption.
659 Args:
660 None
662 Returns:
663 None
664 """
665 units_meta = self.control_volume.config.property_package.get_metadata()
667 # Get indexing sets from control volume
668 # Add isentropic variables
669 self.efficiency_isentropic = Var(
670 self.flowsheet().time,
671 initialize=0.8,
672 doc="Efficiency with respect to an isentropic process [-]",
673 )
674 self.work_isentropic = Var(
675 self.flowsheet().time,
676 initialize=0.0,
677 doc="Work input to unit if isentropic process",
678 units=units_meta.get_derived_units("power"),
679 )
681 # Build isentropic state block
682 tmp_dict = dict(**self.config.property_package_args)
683 tmp_dict["has_phase_equilibrium"] = self.config.has_phase_equilibrium
684 tmp_dict["defined_state"] = False
686 self.properties_isentropic = self.config.property_package.build_state_block(
687 self.flowsheet().time, doc="isentropic properties at outlet", **tmp_dict
688 )
690 # Connect isentropic state block properties
691 @self.Constraint(
692 self.flowsheet().time, doc="Pressure for isentropic calculations"
693 )
694 def isentropic_pressure(self, t):
695 return (
696 self.properties_isentropic[t].pressure
697 == self.control_volume.properties_out[t].pressure
698 )
700 # This assumes isentropic composition is the same as outlet
701 self.add_state_material_balances(
702 self.config.material_balance_type,
703 self.properties_isentropic,
704 self.control_volume.properties_out,
705 )
707 # This assumes isentropic entropy is the same as inlet
708 @self.Constraint(self.flowsheet().time, doc="Isentropic assumption")
709 def isentropic(self, t):
710 return (
711 self.properties_isentropic[t].entr_mol
712 == self.control_volume.properties_in[t].entr_mol
713 )
715 # Isentropic work
716 @self.Constraint(
717 self.flowsheet().time, doc="Calculate work of isentropic process"
718 )
719 def isentropic_energy_balance(self, t):
720 return self.work_isentropic[t] == (
721 sum(
722 self.properties_isentropic[t].get_enthalpy_flow_terms(p)
723 for p in self.properties_isentropic.phase_list
724 )
725 - sum(
726 self.control_volume.properties_in[t].get_enthalpy_flow_terms(p)
727 for p in self.control_volume.properties_in.phase_list
728 )
729 )
731 # Actual work
732 @self.Constraint(
733 self.flowsheet().time, doc="Actual mechanical work calculation"
734 )
735 def actual_work(self, t):
736 if self.config.compressor: 736 ↛ 741line 736 didn't jump to line 741 because the condition on line 736 was always true
737 return self.work_isentropic[t] == (
738 self.work_mechanical[t] * self.efficiency_isentropic[t]
739 )
740 else:
741 return self.work_mechanical[t] == (
742 self.work_isentropic[t] * self.efficiency_isentropic[t]
743 )
745 if self.config.support_isentropic_performance_curves: 745 ↛ 746line 745 didn't jump to line 746 because the condition on line 745 was never true
746 self.performance_curve = IsentropicPerformanceCurve(
747 **self.config.isentropic_performance_curves
748 )
750 def model_check(blk):
751 """
752 Check that pressure change matches with compressor argument (i.e. if
753 compressor = True, pressure should increase or work should be positive)
755 Args:
756 None
758 Returns:
759 None
760 """
761 if blk.config.compressor:
762 # Compressor
763 # Check that pressure does not decrease
764 if any(
765 blk.deltaP[t].fixed and (value(blk.deltaP[t]) < 0.0)
766 for t in blk.flowsheet().time
767 ):
768 _log.warning("{} Compressor set with negative deltaP.".format(blk.name))
769 if any(
770 blk.ratioP[t].fixed and (value(blk.ratioP[t]) < 1.0)
771 for t in blk.flowsheet().time
772 ):
773 _log.warning(
774 "{} Compressor set with ratioP less than 1.".format(blk.name)
775 )
776 if any(
777 blk.control_volume.properties_out[t].pressure.fixed
778 and (
779 value(blk.control_volume.properties_in[t].pressure)
780 > value(blk.control_volume.properties_out[t].pressure)
781 )
782 for t in blk.flowsheet().time
783 ):
784 _log.warning(
785 "{} Compressor set with pressure decrease.".format(blk.name)
786 )
787 # Check that work is not negative
788 if any(
789 blk.work_mechanical[t].fixed and (value(blk.work_mechanical[t]) < 0.0)
790 for t in blk.flowsheet().time
791 ):
792 _log.warning(
793 "{} Compressor maybe set with negative work.".format(blk.name)
794 )
795 else:
796 # Expander
797 # Check that pressure does not increase
798 if any(
799 blk.deltaP[t].fixed and (value(blk.deltaP[t]) > 0.0)
800 for t in blk.flowsheet().time
801 ):
802 _log.warning(
803 "{} Expander/turbine set with positive deltaP.".format(blk.name)
804 )
805 if any(
806 blk.ratioP[t].fixed and (value(blk.ratioP[t]) > 1.0)
807 for t in blk.flowsheet().time
808 ):
809 _log.warning(
810 "{} Expander/turbine set with ratioP greater "
811 "than 1.".format(blk.name)
812 )
813 if any(
814 blk.control_volume.properties_out[t].pressure.fixed
815 and (
816 value(blk.control_volume.properties_in[t].pressure)
817 < value(blk.control_volume.properties_out[t].pressure)
818 )
819 for t in blk.flowsheet().time
820 ):
821 _log.warning(
822 "{} Expander/turbine maybe set with pressure "
823 "increase.".format(blk.name),
824 )
825 # Check that work is not positive
826 if any(
827 blk.work_mechanical[t].fixed and (value(blk.work_mechanical[t]) > 0.0)
828 for t in blk.flowsheet().time
829 ):
830 _log.warning(
831 "{} Expander/turbine set with positive work.".format(blk.name)
832 )
834 # Run holdup block model checks
835 blk.control_volume.model_check()
837 # Run model checks on isentropic property block
838 try:
839 for t in blk.flowsheet().time:
840 blk.properties_in[t].model_check()
841 except AttributeError:
842 pass
844 def initialize_build(
845 blk,
846 state_args=None,
847 routine=None,
848 outlvl=idaeslog.NOTSET,
849 solver=None,
850 optarg=None,
851 ):
852 """
853 General wrapper for pressure changer initialization routines
855 Keyword Arguments:
856 routine : str stating which initialization routine to execute
857 * None - use routine matching thermodynamic_assumption
858 * 'isentropic' - use isentropic initialization routine
859 * 'isothermal' - use isothermal initialization routine
860 state_args : a dict of arguments to be passed to the property
861 package(s) to provide an initial state for
862 initialization (see documentation of the specific
863 property package) (default = {}).
864 outlvl : sets output level of initialization routine
865 optarg : solver options dictionary object (default=None, use
866 default solver options)
867 solver : str indicating which solver to use during
868 initialization (default = None, use default solver)
870 Returns:
871 None
872 """
873 if routine is None: 873 ↛ 878line 873 didn't jump to line 878 because the condition on line 873 was always true
874 # Use routine for specific type of unit
875 routine = blk.config.thermodynamic_assumption
877 # Call initialization routine
878 if routine is ThermodynamicAssumption.isentropic:
879 blk.init_isentropic(
880 state_args=state_args, outlvl=outlvl, solver=solver, optarg=optarg
881 )
882 elif routine is ThermodynamicAssumption.adiabatic:
883 blk.init_adiabatic(
884 state_args=state_args, outlvl=outlvl, solver=solver, optarg=optarg
885 )
886 else:
887 # Call the general initialization routine in UnitModelBlockData
888 super().initialize_build(
889 state_args=state_args, outlvl=outlvl, solver=solver, optarg=optarg
890 )
892 def init_adiabatic(blk, state_args, outlvl, solver, optarg):
893 """
894 Initialization routine for adiabatic pressure changers.
896 Keyword Arguments:
897 state_args : a dict of arguments to be passed to the property
898 package(s) to provide an initial state for
899 initialization (see documentation of the specific
900 property package) (default = {}).
901 outlvl : sets output level of initialization routine
902 optarg : solver options dictionary object (default={})
903 solver : str indicating which solver to use during
904 initialization (default = None)
906 Returns:
907 None
908 """
909 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
910 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
912 # Create solver
913 opt = get_solver(solver, optarg)
915 cv = blk.control_volume
916 t0 = blk.flowsheet().time.first()
918 if state_args is None: 918 ↛ 922line 918 didn't jump to line 922 because the condition on line 918 was always true
919 state_args = extract_state_args(cv.properties_in[t0], use_state_vars=False)
921 # Initialize state blocks
922 flags = cv.properties_in.initialize(
923 outlvl=outlvl,
924 optarg=optarg,
925 solver=solver,
926 hold_state=True,
927 state_args=state_args,
928 )
930 # Get initialisation guesses for outlet state
931 state_args_out = {}
932 # refresh inlet values (may have changed during initialization)
933 inlet_state_args = extract_state_args(cv.properties_in[t0], use_state_vars=False)
934 for k in state_args:
935 if k == "pressure" and k not in state_args_out:
936 # Work out how to estimate outlet pressure
937 if cv.properties_out[t0].pressure.fixed:
938 # Fixed outlet pressure, use this value
939 state_args_out[k] = value(cv.properties_out[t0].pressure)
940 elif blk.deltaP[t0].fixed:
941 state_args_out[k] = value(
942 cv.properties_in[t0].pressure + blk.deltaP[t0]
943 )
944 elif blk.ratioP[t0].fixed: 944 ↛ 945line 944 didn't jump to line 945 because the condition on line 944 was never true
945 state_args_out[k] = value(
946 cv.properties_in[t0].pressure * blk.ratioP[t0]
947 )
948 else:
949 # Not obvious what to do, use inlet state
950 state_args_out[k] = value(cv.properties_in[t0].pressure)
951 elif k not in state_args_out: 951 ↛ 934line 951 didn't jump to line 934 because the condition on line 951 was always true
952 # use the inlet state as a guess for the outlet state
953 state_args_out[k] = inlet_state_args[k]
955 cv.properties_out.initialize(
956 outlvl=outlvl,
957 optarg=optarg,
958 solver=solver,
959 hold_state=False,
960 state_args=state_args_out,
961 )
962 init_log.info_high("Initialization Step 1 Complete.")
964 # ---------------------------------------------------------------------
965 # Solve unit
966 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
967 res = opt.solve(blk, tee=slc.tee)
968 init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res)))
970 # ---------------------------------------------------------------------
971 # Release Inlet state
972 blk.control_volume.release_state(flags, outlvl)
974 if not check_optimal_termination(res): 974 ↛ 975line 974 didn't jump to line 975 because the condition on line 974 was never true
975 raise InitializationError(
976 f"{blk.name} failed to initialize successfully. Please check "
977 f"the output logs for more information."
978 )
980 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
982 def init_isentropic(blk, state_args, outlvl, solver, optarg):
983 """
984 Initialization routine for isentropic pressure changers.
986 Keyword Arguments:
987 state_args : a dict of arguments to be passed to the property
988 package(s) to provide an initial state for
989 initialization (see documentation of the specific
990 property package) (default = {}).
991 outlvl : sets output level of initialization routine
992 optarg : solver options dictionary object (default={})
993 solver : str indicating which solver to use during
994 initialization (default = None)
996 Returns:
997 None
998 """
999 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
1000 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
1002 # Create solver
1003 opt = get_solver(solver, optarg)
1005 cv = blk.control_volume
1006 t0 = blk.flowsheet().time.first()
1008 # performance curves exist and are active so initialize with them
1009 activate_performance_curves = (
1010 hasattr(blk, "performance_curve")
1011 and blk.performance_curve.has_constraints()
1012 and blk.performance_curve.active
1013 )
1014 if activate_performance_curves: 1014 ↛ 1015line 1014 didn't jump to line 1015 because the condition on line 1014 was never true
1015 blk.performance_curve.deactivate()
1016 # The performance curves will provide (maybe indirectly) efficiency
1017 # and/or pressure ratio. To get through the standard isentropic
1018 # pressure changer init, we'll see if the user provided a guess for
1019 # pressure ratio or isentropic efficiency and fix them if needed. If
1020 # not fixed and no guess provided, fill in something reasonable
1021 # until the performance curves are turned on.
1022 unfix_eff = {}
1023 unfix_ratioP = {}
1024 for t in blk.flowsheet().time:
1025 if not (
1026 blk.ratioP[t].fixed
1027 or blk.deltaP[t].fixed
1028 or cv.properties_out[t].pressure.fixed
1029 ):
1030 if blk.config.compressor:
1031 if not (
1032 value(blk.ratioP[t]) >= 1.01 and value(blk.ratioP[t]) <= 50
1033 ):
1034 blk.ratioP[t] = 1.8
1035 else:
1036 if not (
1037 value(blk.ratioP[t]) >= 0.01
1038 and value(blk.ratioP[t]) <= 0.999
1039 ):
1040 blk.ratioP[t] = 0.7
1041 blk.ratioP[t].fix()
1042 unfix_ratioP[t] = True
1043 if not blk.efficiency_isentropic[t].fixed:
1044 if not (
1045 value(blk.efficiency_isentropic[t]) >= 0.05
1046 and value(blk.efficiency_isentropic[t]) <= 1.0
1047 ):
1048 blk.efficiency_isentropic[t] = 0.8
1049 blk.efficiency_isentropic[t].fix()
1050 unfix_eff[t] = True
1052 if state_args is None: 1052 ↛ 1056line 1052 didn't jump to line 1056 because the condition on line 1052 was always true
1053 state_args = extract_state_args(cv.properties_in[t0], use_state_vars=False)
1055 # Initialize state blocks
1056 flags = cv.properties_in.initialize(
1057 outlvl=outlvl,
1058 optarg=optarg,
1059 solver=solver,
1060 hold_state=True,
1061 state_args=state_args,
1062 )
1064 # Get initialisation guesses for outlet and isentropic states
1065 state_args_out = {}
1066 # refresh inlet values (may have changed during initialization)
1067 inlet_state_args = extract_state_args(cv.properties_in[t0], use_state_vars=False)
1068 for k in state_args:
1069 if k == "pressure" and k not in state_args_out:
1070 # Work out how to estimate outlet pressure
1071 if cv.properties_out[t0].pressure.fixed:
1072 # Fixed outlet pressure, use this value
1073 state_args_out[k] = value(cv.properties_out[t0].pressure)
1074 elif blk.deltaP[t0].fixed:
1075 state_args_out[k] = value(
1076 cv.properties_in[t0].pressure + blk.deltaP[t0]
1077 )
1078 elif blk.ratioP[t0].fixed:
1079 state_args_out[k] = value(
1080 cv.properties_in[t0].pressure * blk.ratioP[t0]
1081 )
1082 else:
1083 # Not obvious what to do, use inlet state
1084 state_args_out[k] = value(cv.properties_in[t0].pressure)
1085 elif k not in state_args_out: 1085 ↛ 1068line 1085 didn't jump to line 1068 because the condition on line 1085 was always true
1086 # use the inlet state as a guess for the outlet state
1087 state_args_out[k] = inlet_state_args[k]
1089 cv.properties_out.initialize(
1090 outlvl=outlvl,
1091 optarg=optarg,
1092 solver=solver,
1093 hold_state=False,
1094 state_args=state_args_out,
1095 )
1097 init_log.info_high("Initialization Step 1 Complete.")
1098 # ---------------------------------------------------------------------
1099 # Initialize Isentropic block
1101 blk.properties_isentropic.initialize(
1102 outlvl=outlvl,
1103 optarg=optarg,
1104 solver=solver,
1105 state_args=state_args_out,
1106 )
1108 init_log.info_high("Initialization Step 2 Complete.")
1110 # ---------------------------------------------------------------------
1111 # Solve for isothermal conditions
1112 if isinstance(
1113 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1114 Var,
1115 ):
1116 blk.properties_isentropic[:].temperature.fix()
1117 elif isinstance( 1117 ↛ 1122line 1117 didn't jump to line 1122 because the condition on line 1117 was always true
1118 blk.properties_isentropic[blk.flowsheet().time.first()].enth_mol,
1119 Var,
1120 ):
1121 blk.properties_isentropic[:].enth_mol.fix()
1122 elif isinstance(
1123 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1124 Expression,
1125 ):
1127 def tmp_rule(self, t):
1128 return (
1129 blk.properties_isentropic[t].temperature
1130 == blk.control_volume.properties_in[t].temperature
1131 )
1133 blk.tmp_init_constraint = Constraint(blk.flowsheet().time, rule=tmp_rule)
1135 blk.isentropic.deactivate()
1137 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
1138 res = opt.solve(blk, tee=slc.tee)
1139 init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res)))
1141 if isinstance(
1142 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1143 Var,
1144 ):
1145 blk.properties_isentropic[:].temperature.unfix()
1146 elif isinstance( 1146 ↛ 1151line 1146 didn't jump to line 1151 because the condition on line 1146 was always true
1147 blk.properties_isentropic[blk.flowsheet().time.first()].enth_mol,
1148 Var,
1149 ):
1150 blk.properties_isentropic[:].enth_mol.unfix()
1151 elif isinstance(
1152 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1153 Expression,
1154 ):
1155 blk.del_component(blk.tmp_init_constraint)
1157 blk.isentropic.activate()
1159 # ---------------------------------------------------------------------
1160 # Solve unit
1161 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
1162 res = opt.solve(blk, tee=slc.tee)
1163 init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res)))
1165 if activate_performance_curves: 1165 ↛ 1166line 1165 didn't jump to line 1166 because the condition on line 1165 was never true
1166 blk.performance_curve.activate()
1167 for t, v in unfix_eff.items():
1168 if v:
1169 blk.efficiency_isentropic[t].unfix()
1170 for t, v in unfix_ratioP.items():
1171 if v:
1172 blk.ratioP[t].unfix()
1173 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
1174 res = opt.solve(blk, tee=slc.tee)
1175 init_log.info_high(f"Initialization Step 5 {idaeslog.condition(res)}.")
1177 # ---------------------------------------------------------------------
1178 # Release Inlet state
1179 blk.control_volume.release_state(flags, outlvl)
1181 if not check_optimal_termination(res): 1181 ↛ 1182line 1181 didn't jump to line 1182 because the condition on line 1181 was never true
1182 raise InitializationError(
1183 f"{blk.name} failed to initialize successfully. Please check "
1184 f"the output logs for more information."
1185 )
1187 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
1189 def _get_performance_contents(self, time_point=0):
1190 var_dict = {}
1191 if hasattr(self, "deltaP"):
1192 var_dict["Mechanical Work"] = self.work_mechanical[time_point]
1193 if hasattr(self, "deltaP"):
1194 var_dict["Pressure Change"] = self.deltaP[time_point]
1195 if hasattr(self, "ratioP"):
1196 var_dict["Pressure Ratio"] = self.ratioP[time_point]
1197 if hasattr(self, "efficiency_pump"):
1198 var_dict["Efficiency"] = self.efficiency_pump[time_point]
1199 if hasattr(self, "efficiency_isentropic"):
1200 var_dict["Isentropic Efficiency"] = self.efficiency_isentropic[time_point]
1202 return {"vars": var_dict}
1204 def calculate_scaling_factors(self):
1205 super().calculate_scaling_factors()
1207 if hasattr(self, "work_fluid"):
1208 for t, v in self.work_fluid.items():
1209 iscale.set_scaling_factor(
1210 v,
1211 iscale.get_scaling_factor(
1212 self.control_volume.work[t], default=1, warning=True
1213 ),
1214 )
1216 if hasattr(self, "work_mechanical"):
1217 for t, v in self.work_mechanical.items():
1218 iscale.set_scaling_factor(
1219 v,
1220 iscale.get_scaling_factor(
1221 self.control_volume.work[t], default=1, warning=True
1222 ),
1223 )
1225 if hasattr(self, "work_isentropic"):
1226 for t, v in self.work_isentropic.items():
1227 iscale.set_scaling_factor(
1228 v,
1229 iscale.get_scaling_factor(
1230 self.control_volume.work[t], default=1, warning=True
1231 ),
1232 )
1234 if hasattr(self, "ratioP_calculation"):
1235 for t, c in self.ratioP_calculation.items():
1236 iscale.constraint_scaling_transform(
1237 c,
1238 iscale.get_scaling_factor(
1239 self.control_volume.properties_in[t].pressure,
1240 default=1,
1241 warning=True,
1242 ),
1243 overwrite=False,
1244 )
1246 if hasattr(self, "fluid_work_calculation"):
1247 for t, c in self.fluid_work_calculation.items():
1248 iscale.constraint_scaling_transform(
1249 c,
1250 iscale.get_scaling_factor(
1251 self.control_volume.deltaP[t], default=1, warning=True
1252 ),
1253 overwrite=False,
1254 )
1256 if hasattr(self, "actual_work"):
1257 for t, c in self.actual_work.items():
1258 iscale.constraint_scaling_transform(
1259 c,
1260 iscale.get_scaling_factor(
1261 self.control_volume.work[t], default=1, warning=True
1262 ),
1263 overwrite=False,
1264 )
1266 if hasattr(self, "isentropic_pressure"):
1267 for t, c in self.isentropic_pressure.items():
1268 iscale.constraint_scaling_transform(
1269 c,
1270 iscale.get_scaling_factor(
1271 self.control_volume.properties_in[t].pressure,
1272 default=1,
1273 warning=True,
1274 ),
1275 overwrite=False,
1276 )
1278 if hasattr(self, "isentropic"):
1279 for t, c in self.isentropic.items():
1280 iscale.constraint_scaling_transform(
1281 c,
1282 iscale.get_scaling_factor(
1283 self.control_volume.properties_in[t].entr_mol,
1284 default=1,
1285 warning=True,
1286 ),
1287 overwrite=False,
1288 )
1290 if hasattr(self, "isentropic_energy_balance"):
1291 for t, c in self.isentropic_energy_balance.items():
1292 iscale.constraint_scaling_transform(
1293 c,
1294 iscale.get_scaling_factor(
1295 self.control_volume.work[t], default=1, warning=True
1296 ),
1297 overwrite=False,
1298 )
1300 if hasattr(self, "zero_work_equation"):
1301 for t, c in self.zero_work_equation.items():
1302 iscale.constraint_scaling_transform(
1303 c,
1304 iscale.get_scaling_factor(
1305 self.control_volume.work[t], default=1, warning=True
1306 ),
1307 )
1309 if hasattr(self, "state_material_balances"):
1310 cvol = self.control_volume
1311 phase_list = cvol.properties_in.phase_list
1312 phase_component_set = cvol.properties_in.phase_component_set
1313 mb_type = cvol._constructed_material_balance_type
1314 if mb_type == MaterialBalanceType.componentPhase:
1315 for (t, p, j), c in self.state_material_balances.items():
1316 sf = iscale.get_scaling_factor(
1317 cvol.properties_in[t].get_material_flow_terms(p, j),
1318 default=1,
1319 warning=True,
1320 )
1321 iscale.constraint_scaling_transform(c, sf)
1322 elif mb_type == MaterialBalanceType.componentTotal:
1323 for (t, j), c in self.state_material_balances.items():
1324 sf = iscale.min_scaling_factor(
1325 [
1326 cvol.properties_in[t].get_material_flow_terms(p, j)
1327 for p in phase_list
1328 if (p, j) in phase_component_set
1329 ]
1330 )
1331 iscale.constraint_scaling_transform(c, sf)
1332 else:
1333 # There are some other material balance types but they create
1334 # constraints with different names.
1335 _log.warning(f"Unknown material balance type {mb_type}")
1338@declare_process_block_class("Turbine", doc="Isentropic turbine model")
1339class TurbineData(PressureChangerData):
1340 """
1341 Pressure changer with isentropic turbine options
1342 """
1344 default_initializer = IsentropicPressureChangerInitializer
1346 CONFIG = PressureChangerData.CONFIG()
1347 CONFIG.compressor = False
1348 CONFIG.get("compressor")._default = False
1349 CONFIG.get("compressor")._domain = In([False])
1350 CONFIG.thermodynamic_assumption = ThermodynamicAssumption.isentropic
1351 CONFIG.get("thermodynamic_assumption")._default = ThermodynamicAssumption.isentropic
1354@declare_process_block_class("Compressor", doc="Isentropic compressor model")
1355class CompressorData(PressureChangerData):
1356 """Pressure changer with isentropic compressor options"""
1358 default_initializer = IsentropicPressureChangerInitializer
1360 CONFIG = PressureChangerData.CONFIG()
1361 CONFIG.compressor = True
1362 CONFIG.get("compressor")._default = True
1363 CONFIG.get("compressor")._domain = In([True])
1364 CONFIG.thermodynamic_assumption = ThermodynamicAssumption.isentropic
1365 CONFIG.get("thermodynamic_assumption")._default = ThermodynamicAssumption.isentropic
1368@declare_process_block_class("Pump", doc="Pump model")
1369class PumpData(PressureChangerData):
1370 """Pressure changer with pump options"""
1372 CONFIG = PressureChangerData.CONFIG()
1373 CONFIG.compressor = True
1374 CONFIG.get("compressor")._default = True
1375 CONFIG.get("compressor")._domain = In([True])
1376 CONFIG.thermodynamic_assumption = ThermodynamicAssumption.pump
1377 CONFIG.get("thermodynamic_assumption")._default = ThermodynamicAssumption.pump