Coverage for backend/idaes_service/solver/custom/updated_pressure_changer.py: 43%
452 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +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
60__author__ = "Emmanuel Ogbe, Andrew Lee"
61_log = idaeslog.getLogger(__name__)
64class ThermodynamicAssumption(Enum):
65 """
66 Enum of supported thermodynamic assumptions.
67 """
69 isothermal = 1
70 isentropic = 2
71 pump = 3
72 adiabatic = 4
75class IsentropicPressureChangerInitializer(SingleControlVolumeUnitInitializer):
76 """
77 Initializer for isentropic pressure changer models (and derived types).
79 """
81 def initialization_routine(
82 self,
83 model: Block,
84 ):
85 """
86 Initialization routine for isentropic pressure changers.
88 This routine starts by initializing the inlet and outlet states as usual,
89 using the user provided operating conditions to estimate the outlet state.
90 The isentropic state is then initialized at the same conditions as the outlet.
91 Next, the pressure changer is solved with an isothermal assumption and fixed efficiency,
92 followed by a second solve with the isentropic constraints. Finally, if user-provided
93 performance constraints are present, these are activated and the model solved again.
95 Args:
96 model: model to be initialized
98 Returns:
99 Pyomo solver status object
101 """
102 init_log = idaeslog.getInitLogger(
103 model.name, self.get_output_level(), tag="unit"
104 )
105 solve_log = idaeslog.getSolveLogger(
106 model.name, self.get_output_level(), tag="unit"
107 )
109 # Create solver
110 solver = self._get_solver()
112 cv = model.control_volume
114 # check if performance curves exist and are active
115 activate_performance_curves = (
116 hasattr(model, "performance_curve")
117 and model.performance_curve.has_constraints()
118 and model.performance_curve.active
119 )
120 if activate_performance_curves:
121 model.performance_curve.deactivate()
122 # The performance curves will provide (maybe indirectly) efficiency
123 # and/or pressure ratio. To get through the standard isentropic
124 # pressure changer init, we'll see if the user provided a guess for
125 # pressure ratio or isentropic efficiency and fix them if needed. If
126 # not fixed and no guess provided, fill in something reasonable
127 # until the performance curves are turned on.
128 unfix_eff = {}
129 unfix_ratioP = {}
130 for t in model.flowsheet().time:
131 if not (
132 model.ratioP[t].fixed
133 or model.deltaP[t].fixed
134 or cv.properties_out[t].pressure.fixed
135 ):
136 if model.config.compressor:
137 if not (
138 value(model.ratioP[t]) >= 1.01
139 and value(model.ratioP[t]) <= 50
140 ):
141 model.ratioP[t] = 1.8
142 else:
143 if not (
144 value(model.ratioP[t]) >= 0.01
145 and value(model.ratioP[t]) <= 0.999
146 ):
147 model.ratioP[t] = 0.7
148 model.ratioP[t].fix()
149 unfix_ratioP[t] = True
150 if not model.efficiency_isentropic[t].fixed:
151 if not (
152 value(model.efficiency_isentropic[t]) >= 0.05
153 and value(model.efficiency_isentropic[t]) <= 1.0
154 ):
155 model.efficiency_isentropic[t] = 0.8
156 model.efficiency_isentropic[t].fix()
157 unfix_eff[t] = True
159 # Initialize state blocks
160 self.initialize_control_volume(cv, copy_inlet_state=False)
162 init_log.info_high("Initialization Step 1 Complete.")
163 # ---------------------------------------------------------------------
164 # Copy outlet state to isentropic state
166 # Map solution from inlet properties to outlet properties
167 state = to_json(
168 cv.properties_out,
169 wts=StoreSpec().value(),
170 return_dict=True,
171 )
172 from_json(
173 model.properties_isentropic,
174 sd=state,
175 wts=StoreSpec().value(only_not_fixed=True),
176 )
178 init_log.info_high("Initialization Step 2 Complete.")
180 # ---------------------------------------------------------------------
181 # Solve for isothermal conditions
182 if isinstance(
183 model.properties_isentropic[model.flowsheet().time.first()].temperature,
184 Var,
185 ):
186 model.properties_isentropic[:].temperature.fix()
187 elif isinstance(
188 model.properties_isentropic[model.flowsheet().time.first()].enth_mol,
189 Var,
190 ):
191 model.properties_isentropic[:].enth_mol.fix()
192 elif isinstance(
193 model.properties_isentropic[model.flowsheet().time.first()].temperature,
194 Expression,
195 ):
197 def tmp_rule(blk, t):
198 return (
199 blk.properties_isentropic[t].temperature
200 == blk.control_volume.properties_in[t].temperature
201 )
203 model.tmp_init_constraint = Constraint(
204 model.flowsheet().time, rule=tmp_rule
205 )
207 model.isentropic.deactivate()
209 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
210 res = solver.solve(model, tee=slc.tee)
211 init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res)))
213 # Revert changes for isothermal assumption
214 if isinstance(
215 model.properties_isentropic[model.flowsheet().time.first()].temperature,
216 Var,
217 ):
218 model.properties_isentropic[:].temperature.unfix()
219 elif isinstance(
220 model.properties_isentropic[model.flowsheet().time.first()].enth_mol,
221 Var,
222 ):
223 model.properties_isentropic[:].enth_mol.unfix()
224 elif isinstance(
225 model.properties_isentropic[model.flowsheet().time.first()].temperature,
226 Expression,
227 ):
228 model.del_component(model.tmp_init_constraint)
230 model.isentropic.activate()
232 # ---------------------------------------------------------------------
233 # Solve unit
234 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
235 res = solver.solve(model, tee=slc.tee)
236 init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res)))
238 if activate_performance_curves:
239 model.performance_curve.activate()
240 for t, v in unfix_eff.items():
241 if v:
242 model.efficiency_isentropic[t].unfix()
243 for t, v in unfix_ratioP.items():
244 if v:
245 model.ratioP[t].unfix()
246 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
247 res = solver.solve(model, tee=slc.tee)
248 init_log.info_high(f"Initialization Step 5 {idaeslog.condition(res)}.")
250 return res
253@declare_process_block_class("IsentropicPerformanceCurve")
254class IsentropicPerformanceCurveData(ProcessBlockData):
255 """Block that holds performance curves. Typically, these are in the form of
256 constraints that relate head, efficiency, or pressure ratio to volumetric
257 or mass flow. Additional variables can be included if needed, such as
258 speed. For convenience an option is provided to add head expressions to the
259 block. performance curves, and any additional variables, constraints, or
260 expressions can be added to this block either via callback provided to the
261 configuration, or after the model is constructed."""
263 CONFIG = ProcessBlockData.CONFIG(
264 doc="Configuration dictionary for the performance curve block."
265 )
266 CONFIG.declare(
267 "build_callback",
268 ConfigValue(
269 default=None, doc="Optional callback to add performance curve constraints"
270 ),
271 )
272 CONFIG.declare(
273 "build_head_expressions",
274 ConfigValue(
275 default=True,
276 domain=bool,
277 doc="If true add expressions for 'head' and 'head_isentropic'. "
278 "These expressions can be used in performance curve constraints.",
279 ),
280 )
282 def has_constraints(self):
283 for _ in self.component_data_objects(Constraint):
284 return True
285 return False
287 def build(self):
288 super().build()
289 if self.config.build_head_expressions:
290 try:
292 @self.Expression(self.flowsheet().time)
293 def head_isentropic(self, t): # units are energy/mass
294 b = self.parent_block()
295 if hasattr(b.control_volume.properties_in[t], "flow_mass"):
296 return (
297 b.work_isentropic[t]
298 / b.control_volume.properties_in[t].flow_mass
299 )
300 else:
301 return (
302 b.work_isentropic[t]
303 / b.control_volume.properties_in[t].flow_mol
304 / b.control_volume.properties_in[t].mw
305 )
307 @self.Expression(self.flowsheet().time)
308 def head(self, t): # units are energy/mass
309 b = self.parent_block()
310 if hasattr(b.control_volume.properties_in[t], "flow_mass"):
311 return (
312 b.work_mechanical[t]
313 / b.control_volume.properties_in[t].flow_mass
314 )
315 else:
316 return (
317 b.work_mechanical[t]
318 / b.control_volume.properties_in[t].flow_mol
319 / b.control_volume.properties_in[t].mw
320 )
322 except PropertyNotSupportedError:
323 _log.exception(
324 "flow_mass or flow_mol and mw are not supported by the "
325 "property package but are required for isentropic pressure"
326 " changer head calculation"
327 )
328 raise
330 if self.config.build_callback is not None:
331 self.config.build_callback(self)
334@declare_process_block_class("PressureChanger")
335class PressureChangerData(UnitModelBlockData):
336 """
337 Standard Compressor/Expander Unit Model Class
338 """
340 CONFIG = UnitModelBlockData.CONFIG()
342 CONFIG.declare(
343 "material_balance_type",
344 ConfigValue(
345 default=MaterialBalanceType.useDefault,
346 domain=In(MaterialBalanceType),
347 description="Material balance construction flag",
348 doc="""Indicates what type of mass balance should be constructed,
349**default** - MaterialBalanceType.useDefault.
350**Valid values:** {
351**MaterialBalanceType.useDefault - refer to property package for default
352balance type
353**MaterialBalanceType.none** - exclude material balances,
354**MaterialBalanceType.componentPhase** - use phase component balances,
355**MaterialBalanceType.componentTotal** - use total component balances,
356**MaterialBalanceType.elementTotal** - use total element balances,
357**MaterialBalanceType.total** - use total material balance.}""",
358 ),
359 )
360 CONFIG.declare(
361 "energy_balance_type",
362 ConfigValue(
363 default=EnergyBalanceType.useDefault,
364 domain=In(EnergyBalanceType),
365 description="Energy balance construction flag",
366 doc="""Indicates what type of energy balance should be constructed,
367**default** - EnergyBalanceType.useDefault.
368**Valid values:** {
369**EnergyBalanceType.useDefault - refer to property package for default
370balance type
371**EnergyBalanceType.none** - exclude energy balances,
372**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material,
373**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase,
374**EnergyBalanceType.energyTotal** - single energy balance for material,
375**EnergyBalanceType.energyPhase** - energy balances for each phase.}""",
376 ),
377 )
378 CONFIG.declare(
379 "momentum_balance_type",
380 ConfigValue(
381 default=MomentumBalanceType.pressureTotal,
382 domain=In(MomentumBalanceType),
383 description="Momentum balance construction flag",
384 doc="""Indicates what type of momentum balance should be
385constructed, **default** - MomentumBalanceType.pressureTotal.
386**Valid values:** {
387**MomentumBalanceType.none** - exclude momentum balances,
388**MomentumBalanceType.pressureTotal** - single pressure balance for material,
389**MomentumBalanceType.pressurePhase** - pressure balances for each phase,
390**MomentumBalanceType.momentumTotal** - single momentum balance for material,
391**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""",
392 ),
393 )
394 CONFIG.declare(
395 "has_phase_equilibrium",
396 ConfigValue(
397 default=False,
398 domain=Bool,
399 description="Phase equilibrium construction flag",
400 doc="""Indicates whether terms for phase equilibrium should be
401constructed, **default** = False.
402**Valid values:** {
403**True** - include phase equilibrium terms
404**False** - exclude phase equilibrium terms.}""",
405 ),
406 )
407 CONFIG.declare(
408 "compressor",
409 ConfigValue(
410 default=True,
411 domain=Bool,
412 description="Compressor flag",
413 doc="""Indicates whether this unit should be considered a
414 compressor (True (default), pressure increase) or an expander
415 (False, pressure decrease).""",
416 ),
417 )
418 CONFIG.declare(
419 "thermodynamic_assumption",
420 ConfigValue(
421 default=ThermodynamicAssumption.isothermal,
422 domain=In(ThermodynamicAssumption),
423 description="Thermodynamic assumption to use",
424 doc="""Flag to set the thermodynamic assumption to use for the unit.
425 - ThermodynamicAssumption.isothermal (default)
426 - ThermodynamicAssumption.isentropic
427 - ThermodynamicAssumption.pump
428 - ThermodynamicAssumption.adiabatic""",
429 ),
430 )
431 CONFIG.declare(
432 "property_package",
433 ConfigValue(
434 default=useDefault,
435 domain=is_physical_parameter_block,
436 description="Property package to use for control volume",
437 doc="""Property parameter object used to define property
438calculations, **default** - useDefault.
439**Valid values:** {
440**useDefault** - use default package from parent model or flowsheet,
441**PropertyParameterObject** - a PropertyParameterBlock object.}""",
442 ),
443 )
444 CONFIG.declare(
445 "property_package_args",
446 ConfigBlock(
447 implicit=True,
448 description="Arguments to use for constructing property packages",
449 doc="""A ConfigBlock with arguments to be passed to a property
450block(s) and used when constructing these,
451**default** - None.
452**Valid values:** {
453see property package for documentation.}""",
454 ),
455 )
456 CONFIG.declare(
457 "support_isentropic_performance_curves",
458 ConfigValue(
459 default=False,
460 domain=Bool,
461 doc="Include a block for performance curves, configure via"
462 " isentropic_performance_curves.",
463 ),
464 )
465 CONFIG.declare(
466 "isentropic_performance_curves",
467 IsentropicPerformanceCurveData.CONFIG(),
468 # doc included in IsentropicPerformanceCurveData
469 )
471 def build(self):
472 """
474 Args:
475 None
477 Returns:
478 None
479 """
480 # Call UnitModel.build
481 super().build()
483 # Add a control volume to the unit including setting up dynamics.
484 self.control_volume = ControlVolume0DBlock(
485 dynamic=self.config.dynamic,
486 has_holdup=self.config.has_holdup,
487 property_package=self.config.property_package,
488 property_package_args=self.config.property_package_args,
489 )
491 # Add geometry variables to control volume
492 if self.config.has_holdup: 492 ↛ 493line 492 didn't jump to line 493 because the condition on line 492 was never true
493 self.control_volume.add_geometry()
495 # Add inlet and outlet state blocks to control volume
496 self.control_volume.add_state_blocks(
497 has_phase_equilibrium=self.config.has_phase_equilibrium
498 )
500 # Add mass balance
501 # Set has_equilibrium is False for now
502 # TO DO; set has_equilibrium to True
503 self.control_volume.add_material_balances(
504 balance_type=self.config.material_balance_type,
505 has_phase_equilibrium=self.config.has_phase_equilibrium,
506 )
508 # Add energy balance
509 eb = self.control_volume.add_energy_balances(
510 balance_type=self.config.energy_balance_type, has_work_transfer=True
511 )
513 # add momentum balance
514 self.control_volume.add_momentum_balances(
515 balance_type=self.config.momentum_balance_type, has_pressure_change=True
516 )
518 # Add Ports
519 self.add_inlet_port()
520 self.add_outlet_port()
522 # Set Unit Geometry and holdup Volume
523 if self.config.has_holdup is True: 523 ↛ 524line 523 didn't jump to line 524 because the condition on line 523 was never true
524 self.volume = Reference(self.control_volume.volume[:])
526 # Construct performance equations
527 # Set references to balance terms at unit level
528 # Add Work transfer variable 'work'
529 # If the 'work' variable wasn't already built on the control volume but is needed, create it now.
530 if ( 530 ↛ 535line 530 didn't jump to line 535 because the condition on line 530 was never true
531 not hasattr(self.control_volume, "work")
532 and self.config.thermodynamic_assumption == ThermodynamicAssumption.pump
533 and eb is None
534 ):
535 units = (
536 self.control_volume.config.property_package.get_metadata().get_derived_units
537 )
538 self.control_volume.work = Var(
539 self.flowsheet().time,
540 domain=Reals,
541 initialize=0.0,
542 doc="Work transferred into control volume",
543 units=units("power"),
544 )
545 self.work_mechanical = Reference(self.control_volume.work[:])
547 # Add Momentum balance variable 'deltaP'
548 self.deltaP = Reference(self.control_volume.deltaP[:])
550 # Performance Variables
551 self.ratioP = Var(self.flowsheet().time, initialize=1.0, doc="Pressure Ratio")
553 # Pressure Ratio
554 @self.Constraint(self.flowsheet().time, doc="Pressure ratio constraint")
555 def ratioP_calculation(self, t):
556 return (
557 self.ratioP[t] * self.control_volume.properties_in[t].pressure
558 == self.control_volume.properties_out[t].pressure
559 )
561 # Construct equations for thermodynamic assumption
562 if self.config.thermodynamic_assumption == ThermodynamicAssumption.isothermal: 562 ↛ 563line 562 didn't jump to line 563 because the condition on line 562 was never true
563 self.add_isothermal()
564 elif self.config.thermodynamic_assumption == ThermodynamicAssumption.isentropic:
565 self.add_isentropic()
566 elif self.config.thermodynamic_assumption == ThermodynamicAssumption.pump:
567 self.add_pump()
568 elif self.config.thermodynamic_assumption == ThermodynamicAssumption.adiabatic: 568 ↛ exitline 568 didn't return from function 'build' because the condition on line 568 was always true
569 self.add_adiabatic()
571 def add_pump(self):
572 """
573 Add constraints for the incompressible fluid assumption
575 Args:
576 None
578 Returns:
579 None
580 """
581 units_meta = self.control_volume.config.property_package.get_metadata()
583 self.work_fluid = Var(
584 self.flowsheet().time,
585 initialize=1.0,
586 doc="Work required to increase the pressure of the liquid",
587 units=units_meta.get_derived_units("power"),
588 )
589 self.efficiency_pump = Var(
590 self.flowsheet().time, initialize=1.0, doc="Pump efficiency"
591 )
593 @self.Constraint(self.flowsheet().time, doc="Pump fluid work constraint")
594 def fluid_work_calculation(self, t):
595 return self.work_fluid[t] == (
596 (
597 self.control_volume.properties_out[t].pressure
598 - self.control_volume.properties_in[t].pressure
599 )
600 * self.control_volume.properties_out[t].flow_vol
601 )
603 # Actual work
604 @self.Constraint(
605 self.flowsheet().time, doc="Actual mechanical work calculation"
606 )
607 def actual_work(self, t):
608 if self.config.compressor: 608 ↛ 613line 608 didn't jump to line 613 because the condition on line 608 was always true
609 return self.work_fluid[t] == (
610 self.work_mechanical[t] * self.efficiency_pump[t]
611 )
612 else:
613 return self.work_mechanical[t] == (
614 self.work_fluid[t] * self.efficiency_pump[t]
615 )
617 def add_isothermal(self):
618 """
619 Add constraints for isothermal assumption.
621 Args:
622 None
624 Returns:
625 None
626 """
628 # Isothermal constraint
629 @self.Constraint(
630 self.flowsheet().time,
631 doc="For isothermal condition: Equate inlet and " "outlet temperature",
632 )
633 def isothermal(self, t):
634 return (
635 self.control_volume.properties_in[t].temperature
636 == self.control_volume.properties_out[t].temperature
637 )
639 def add_adiabatic(self):
640 """
641 Add constraints for adiabatic assumption.
643 Args:
644 None
646 Returns:
647 None
648 """
650 @self.Constraint(self.flowsheet().time)
651 def zero_work_equation(self, t):
652 return self.control_volume.work[t] == 0
654 def add_isentropic(self):
655 """
656 Add constraints for isentropic assumption.
658 Args:
659 None
661 Returns:
662 None
663 """
664 units_meta = self.control_volume.config.property_package.get_metadata()
666 # Get indexing sets from control volume
667 # Add isentropic variables
668 self.efficiency_isentropic = Var(
669 self.flowsheet().time,
670 initialize=0.8,
671 doc="Efficiency with respect to an isentropic process [-]",
672 )
673 self.work_isentropic = Var(
674 self.flowsheet().time,
675 initialize=0.0,
676 doc="Work input to unit if isentropic process",
677 units=units_meta.get_derived_units("power"),
678 )
680 # Build isentropic state block
681 tmp_dict = dict(**self.config.property_package_args)
682 tmp_dict["has_phase_equilibrium"] = self.config.has_phase_equilibrium
683 tmp_dict["defined_state"] = False
685 self.properties_isentropic = self.config.property_package.build_state_block(
686 self.flowsheet().time, doc="isentropic properties at outlet", **tmp_dict
687 )
689 # Connect isentropic state block properties
690 @self.Constraint(
691 self.flowsheet().time, doc="Pressure for isentropic calculations"
692 )
693 def isentropic_pressure(self, t):
694 return (
695 self.properties_isentropic[t].pressure
696 == self.control_volume.properties_out[t].pressure
697 )
699 # This assumes isentropic composition is the same as outlet
700 self.add_state_material_balances(
701 self.config.material_balance_type,
702 self.properties_isentropic,
703 self.control_volume.properties_out,
704 )
706 # This assumes isentropic entropy is the same as inlet
707 @self.Constraint(self.flowsheet().time, doc="Isentropic assumption")
708 def isentropic(self, t):
709 return (
710 self.properties_isentropic[t].entr_mol
711 == self.control_volume.properties_in[t].entr_mol
712 )
714 # Isentropic work
715 @self.Constraint(
716 self.flowsheet().time, doc="Calculate work of isentropic process"
717 )
718 def isentropic_energy_balance(self, t):
719 return self.work_isentropic[t] == (
720 sum(
721 self.properties_isentropic[t].get_enthalpy_flow_terms(p)
722 for p in self.properties_isentropic.phase_list
723 )
724 - sum(
725 self.control_volume.properties_in[t].get_enthalpy_flow_terms(p)
726 for p in self.control_volume.properties_in.phase_list
727 )
728 )
730 # Actual work
731 @self.Constraint(
732 self.flowsheet().time, doc="Actual mechanical work calculation"
733 )
734 def actual_work(self, t):
735 if self.config.compressor: 735 ↛ 740line 735 didn't jump to line 740 because the condition on line 735 was always true
736 return self.work_isentropic[t] == (
737 self.work_mechanical[t] * self.efficiency_isentropic[t]
738 )
739 else:
740 return self.work_mechanical[t] == (
741 self.work_isentropic[t] * self.efficiency_isentropic[t]
742 )
744 if self.config.support_isentropic_performance_curves: 744 ↛ 745line 744 didn't jump to line 745 because the condition on line 744 was never true
745 self.performance_curve = IsentropicPerformanceCurve(
746 **self.config.isentropic_performance_curves
747 )
749 def model_check(blk):
750 """
751 Check that pressure change matches with compressor argument (i.e. if
752 compressor = True, pressure should increase or work should be positive)
754 Args:
755 None
757 Returns:
758 None
759 """
760 if blk.config.compressor:
761 # Compressor
762 # Check that pressure does not decrease
763 if any(
764 blk.deltaP[t].fixed and (value(blk.deltaP[t]) < 0.0)
765 for t in blk.flowsheet().time
766 ):
767 _log.warning("{} Compressor set with negative deltaP.".format(blk.name))
768 if any(
769 blk.ratioP[t].fixed and (value(blk.ratioP[t]) < 1.0)
770 for t in blk.flowsheet().time
771 ):
772 _log.warning(
773 "{} Compressor set with ratioP less than 1.".format(blk.name)
774 )
775 if any(
776 blk.control_volume.properties_out[t].pressure.fixed
777 and (
778 value(blk.control_volume.properties_in[t].pressure)
779 > value(blk.control_volume.properties_out[t].pressure)
780 )
781 for t in blk.flowsheet().time
782 ):
783 _log.warning(
784 "{} Compressor set with pressure decrease.".format(blk.name)
785 )
786 # Check that work is not negative
787 if any(
788 blk.work_mechanical[t].fixed and (value(blk.work_mechanical[t]) < 0.0)
789 for t in blk.flowsheet().time
790 ):
791 _log.warning(
792 "{} Compressor maybe set with negative work.".format(blk.name)
793 )
794 else:
795 # Expander
796 # Check that pressure does not increase
797 if any(
798 blk.deltaP[t].fixed and (value(blk.deltaP[t]) > 0.0)
799 for t in blk.flowsheet().time
800 ):
801 _log.warning(
802 "{} Expander/turbine set with positive deltaP.".format(blk.name)
803 )
804 if any(
805 blk.ratioP[t].fixed and (value(blk.ratioP[t]) > 1.0)
806 for t in blk.flowsheet().time
807 ):
808 _log.warning(
809 "{} Expander/turbine set with ratioP greater "
810 "than 1.".format(blk.name)
811 )
812 if any(
813 blk.control_volume.properties_out[t].pressure.fixed
814 and (
815 value(blk.control_volume.properties_in[t].pressure)
816 < value(blk.control_volume.properties_out[t].pressure)
817 )
818 for t in blk.flowsheet().time
819 ):
820 _log.warning(
821 "{} Expander/turbine maybe set with pressure "
822 "increase.".format(blk.name),
823 )
824 # Check that work is not positive
825 if any(
826 blk.work_mechanical[t].fixed and (value(blk.work_mechanical[t]) > 0.0)
827 for t in blk.flowsheet().time
828 ):
829 _log.warning(
830 "{} Expander/turbine set with positive work.".format(blk.name)
831 )
833 # Run holdup block model checks
834 blk.control_volume.model_check()
836 # Run model checks on isentropic property block
837 try:
838 for t in blk.flowsheet().time:
839 blk.properties_in[t].model_check()
840 except AttributeError:
841 pass
843 def initialize_build(
844 blk,
845 state_args=None,
846 routine=None,
847 outlvl=idaeslog.NOTSET,
848 solver=None,
849 optarg=None,
850 ):
851 """
852 General wrapper for pressure changer initialization routines
854 Keyword Arguments:
855 routine : str stating which initialization routine to execute
856 * None - use routine matching thermodynamic_assumption
857 * 'isentropic' - use isentropic initialization routine
858 * 'isothermal' - use isothermal initialization routine
859 state_args : a dict of arguments to be passed to the property
860 package(s) to provide an initial state for
861 initialization (see documentation of the specific
862 property package) (default = {}).
863 outlvl : sets output level of initialization routine
864 optarg : solver options dictionary object (default=None, use
865 default solver options)
866 solver : str indicating which solver to use during
867 initialization (default = None, use default solver)
869 Returns:
870 None
871 """
872 if routine is None: 872 ↛ 877line 872 didn't jump to line 877 because the condition on line 872 was always true
873 # Use routine for specific type of unit
874 routine = blk.config.thermodynamic_assumption
876 # Call initialization routine
877 if routine is ThermodynamicAssumption.isentropic:
878 blk.init_isentropic(
879 state_args=state_args, outlvl=outlvl, solver=solver, optarg=optarg
880 )
881 elif routine is ThermodynamicAssumption.adiabatic:
882 blk.init_adiabatic(
883 state_args=state_args, outlvl=outlvl, solver=solver, optarg=optarg
884 )
885 else:
886 # Call the general initialization routine in UnitModelBlockData
887 super().initialize_build(
888 state_args=state_args, outlvl=outlvl, solver=solver, optarg=optarg
889 )
891 def init_adiabatic(blk, state_args, outlvl, solver, optarg):
892 """
893 Initialization routine for adiabatic pressure changers.
895 Keyword Arguments:
896 state_args : a dict of arguments to be passed to the property
897 package(s) to provide an initial state for
898 initialization (see documentation of the specific
899 property package) (default = {}).
900 outlvl : sets output level of initialization routine
901 optarg : solver options dictionary object (default={})
902 solver : str indicating which solver to use during
903 initialization (default = None)
905 Returns:
906 None
907 """
908 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
909 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
911 # Create solver
912 opt = get_solver(solver, optarg)
914 cv = blk.control_volume
915 t0 = blk.flowsheet().time.first()
917 if state_args is None: 917 ↛ 930line 917 didn't jump to line 930 because the condition on line 917 was always true
918 state_args = {}
919 state_dict = cv.properties_in[t0].define_port_members()
921 for k in state_dict.keys():
922 if state_dict[k].is_indexed(): 922 ↛ 923line 922 didn't jump to line 923 because the condition on line 922 was never true
923 state_args[k] = {}
924 for m in state_dict[k].keys():
925 state_args[k][m] = state_dict[k][m].value
926 else:
927 state_args[k] = state_dict[k].value
929 # Initialize state blocks
930 flags = cv.properties_in.initialize(
931 outlvl=outlvl,
932 optarg=optarg,
933 solver=solver,
934 hold_state=True,
935 state_args=state_args,
936 )
938 # Get initialisation guesses for outlet state
939 state_args_out = {}
940 # refresh state_dict (may have changed during initialization)
941 state_dict = cv.properties_in[t0].define_port_members()
942 for k in state_args:
943 if k == "pressure" and k not in state_args_out:
944 # Work out how to estimate outlet pressure
945 if cv.properties_out[t0].pressure.fixed:
946 # Fixed outlet pressure, use this value
947 state_args_out[k] = value(cv.properties_out[t0].pressure)
948 elif blk.deltaP[t0].fixed: 948 ↛ 949line 948 didn't jump to line 949 because the condition on line 948 was never true
949 state_args_out[k] = value(
950 cv.properties_in[t0].pressure + blk.deltaP[t0]
951 )
952 elif blk.ratioP[t0].fixed: 952 ↛ 953line 952 didn't jump to line 953 because the condition on line 952 was never true
953 state_args_out[k] = value(
954 cv.properties_in[t0].pressure * blk.ratioP[t0]
955 )
956 else:
957 # Not obvious what to do, use inlet state
958 state_args_out[k] = value(cv.properties_in[t0].pressure)
959 elif k not in state_args_out: 959 ↛ 942line 959 didn't jump to line 942 because the condition on line 959 was always true
960 # use the inlet state as a guess for the outlet state
961 if state_dict[k].is_indexed(): 961 ↛ 962line 961 didn't jump to line 962 because the condition on line 961 was never true
962 state_args_out[k] = {}
963 for m in state_dict[k].keys():
964 state_args_out[k][m] = state_dict[k][m].value
965 else:
966 state_args_out[k] = state_dict[k].value
968 cv.properties_out.initialize(
969 outlvl=outlvl,
970 optarg=optarg,
971 solver=solver,
972 hold_state=False,
973 state_args=state_args_out,
974 )
975 init_log.info_high("Initialization Step 1 Complete.")
977 # ---------------------------------------------------------------------
978 # Solve unit
979 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
980 res = opt.solve(blk, tee=slc.tee)
981 init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res)))
983 # ---------------------------------------------------------------------
984 # Release Inlet state
985 blk.control_volume.release_state(flags, outlvl)
987 if not check_optimal_termination(res): 987 ↛ 988line 987 didn't jump to line 988 because the condition on line 987 was never true
988 raise InitializationError(
989 f"{blk.name} failed to initialize successfully. Please check "
990 f"the output logs for more information."
991 )
993 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
995 def init_isentropic(blk, state_args, outlvl, solver, optarg):
996 """
997 Initialization routine for isentropic pressure changers.
999 Keyword Arguments:
1000 state_args : a dict of arguments to be passed to the property
1001 package(s) to provide an initial state for
1002 initialization (see documentation of the specific
1003 property package) (default = {}).
1004 outlvl : sets output level of initialization routine
1005 optarg : solver options dictionary object (default={})
1006 solver : str indicating which solver to use during
1007 initialization (default = None)
1009 Returns:
1010 None
1011 """
1012 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
1013 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
1015 # Create solver
1016 opt = get_solver(solver, optarg)
1018 cv = blk.control_volume
1019 t0 = blk.flowsheet().time.first()
1021 # performance curves exist and are active so initialize with them
1022 activate_performance_curves = (
1023 hasattr(blk, "performance_curve")
1024 and blk.performance_curve.has_constraints()
1025 and blk.performance_curve.active
1026 )
1027 if activate_performance_curves: 1027 ↛ 1028line 1027 didn't jump to line 1028 because the condition on line 1027 was never true
1028 blk.performance_curve.deactivate()
1029 # The performance curves will provide (maybe indirectly) efficiency
1030 # and/or pressure ratio. To get through the standard isentropic
1031 # pressure changer init, we'll see if the user provided a guess for
1032 # pressure ratio or isentropic efficiency and fix them if needed. If
1033 # not fixed and no guess provided, fill in something reasonable
1034 # until the performance curves are turned on.
1035 unfix_eff = {}
1036 unfix_ratioP = {}
1037 for t in blk.flowsheet().time:
1038 if not (
1039 blk.ratioP[t].fixed
1040 or blk.deltaP[t].fixed
1041 or cv.properties_out[t].pressure.fixed
1042 ):
1043 if blk.config.compressor:
1044 if not (
1045 value(blk.ratioP[t]) >= 1.01 and value(blk.ratioP[t]) <= 50
1046 ):
1047 blk.ratioP[t] = 1.8
1048 else:
1049 if not (
1050 value(blk.ratioP[t]) >= 0.01
1051 and value(blk.ratioP[t]) <= 0.999
1052 ):
1053 blk.ratioP[t] = 0.7
1054 blk.ratioP[t].fix()
1055 unfix_ratioP[t] = True
1056 if not blk.efficiency_isentropic[t].fixed:
1057 if not (
1058 value(blk.efficiency_isentropic[t]) >= 0.05
1059 and value(blk.efficiency_isentropic[t]) <= 1.0
1060 ):
1061 blk.efficiency_isentropic[t] = 0.8
1062 blk.efficiency_isentropic[t].fix()
1063 unfix_eff[t] = True
1065 if state_args is None: 1065 ↛ 1078line 1065 didn't jump to line 1078 because the condition on line 1065 was always true
1066 state_args = {}
1067 state_dict = cv.properties_in[t0].define_port_members()
1069 for k in state_dict.keys():
1070 if state_dict[k].is_indexed():
1071 state_args[k] = {}
1072 for m in state_dict[k].keys():
1073 state_args[k][m] = state_dict[k][m].value
1074 else:
1075 state_args[k] = state_dict[k].value
1077 # Initialize state blocks
1078 flags = cv.properties_in.initialize(
1079 outlvl=outlvl,
1080 optarg=optarg,
1081 solver=solver,
1082 hold_state=True,
1083 state_args=state_args,
1084 )
1086 # Get initialisation guesses for outlet and isentropic states
1087 state_args_out = {}
1088 # refresh state_dict (may have changed during initialization)
1089 state_dict = cv.properties_in[t0].define_port_members()
1090 for k in state_args:
1091 if k == "pressure" and k not in state_args_out:
1092 # Work out how to estimate outlet pressure
1093 if cv.properties_out[t0].pressure.fixed:
1094 # Fixed outlet pressure, use this value
1095 state_args_out[k] = value(cv.properties_out[t0].pressure)
1096 elif blk.deltaP[t0].fixed:
1097 state_args_out[k] = value(
1098 cv.properties_in[t0].pressure + blk.deltaP[t0]
1099 )
1100 elif blk.ratioP[t0].fixed:
1101 state_args_out[k] = value(
1102 cv.properties_in[t0].pressure * blk.ratioP[t0]
1103 )
1104 else:
1105 # Not obvious what to do, use inlet state
1106 state_args_out[k] = value(cv.properties_in[t0].pressure)
1107 elif k not in state_args_out: 1107 ↛ 1090line 1107 didn't jump to line 1090 because the condition on line 1107 was always true
1108 # use the inlet state as a guess for the outlet state
1109 if state_dict[k].is_indexed():
1110 state_args_out[k] = {}
1111 for m in state_dict[k].keys():
1112 state_args_out[k][m] = state_dict[k][m].value
1113 else:
1114 state_args_out[k] = state_dict[k].value
1116 cv.properties_out.initialize(
1117 outlvl=outlvl,
1118 optarg=optarg,
1119 solver=solver,
1120 hold_state=False,
1121 state_args=state_args_out,
1122 )
1124 init_log.info_high("Initialization Step 1 Complete.")
1125 # ---------------------------------------------------------------------
1126 # Initialize Isentropic block
1128 blk.properties_isentropic.initialize(
1129 outlvl=outlvl,
1130 optarg=optarg,
1131 solver=solver,
1132 state_args=state_args_out,
1133 )
1135 init_log.info_high("Initialization Step 2 Complete.")
1137 # ---------------------------------------------------------------------
1138 # Solve for isothermal conditions
1139 if isinstance(
1140 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1141 Var,
1142 ):
1143 blk.properties_isentropic[:].temperature.fix()
1144 elif isinstance( 1144 ↛ 1149line 1144 didn't jump to line 1149 because the condition on line 1144 was always true
1145 blk.properties_isentropic[blk.flowsheet().time.first()].enth_mol,
1146 Var,
1147 ):
1148 blk.properties_isentropic[:].enth_mol.fix()
1149 elif isinstance(
1150 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1151 Expression,
1152 ):
1154 def tmp_rule(self, t):
1155 return (
1156 blk.properties_isentropic[t].temperature
1157 == blk.control_volume.properties_in[t].temperature
1158 )
1160 blk.tmp_init_constraint = Constraint(blk.flowsheet().time, rule=tmp_rule)
1162 blk.isentropic.deactivate()
1164 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
1165 res = opt.solve(blk, tee=slc.tee)
1166 init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res)))
1168 if isinstance(
1169 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1170 Var,
1171 ):
1172 blk.properties_isentropic[:].temperature.unfix()
1173 elif isinstance( 1173 ↛ 1178line 1173 didn't jump to line 1178 because the condition on line 1173 was always true
1174 blk.properties_isentropic[blk.flowsheet().time.first()].enth_mol,
1175 Var,
1176 ):
1177 blk.properties_isentropic[:].enth_mol.unfix()
1178 elif isinstance(
1179 blk.properties_isentropic[blk.flowsheet().time.first()].temperature,
1180 Expression,
1181 ):
1182 blk.del_component(blk.tmp_init_constraint)
1184 blk.isentropic.activate()
1186 # ---------------------------------------------------------------------
1187 # Solve unit
1188 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
1189 res = opt.solve(blk, tee=slc.tee)
1190 init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res)))
1192 if activate_performance_curves: 1192 ↛ 1193line 1192 didn't jump to line 1193 because the condition on line 1192 was never true
1193 blk.performance_curve.activate()
1194 for t, v in unfix_eff.items():
1195 if v:
1196 blk.efficiency_isentropic[t].unfix()
1197 for t, v in unfix_ratioP.items():
1198 if v:
1199 blk.ratioP[t].unfix()
1200 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
1201 res = opt.solve(blk, tee=slc.tee)
1202 init_log.info_high(f"Initialization Step 5 {idaeslog.condition(res)}.")
1204 # ---------------------------------------------------------------------
1205 # Release Inlet state
1206 blk.control_volume.release_state(flags, outlvl)
1208 if not check_optimal_termination(res): 1208 ↛ 1209line 1208 didn't jump to line 1209 because the condition on line 1208 was never true
1209 raise InitializationError(
1210 f"{blk.name} failed to initialize successfully. Please check "
1211 f"the output logs for more information."
1212 )
1214 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
1216 def _get_performance_contents(self, time_point=0):
1217 var_dict = {}
1218 if hasattr(self, "deltaP"):
1219 var_dict["Mechanical Work"] = self.work_mechanical[time_point]
1220 if hasattr(self, "deltaP"):
1221 var_dict["Pressure Change"] = self.deltaP[time_point]
1222 if hasattr(self, "ratioP"):
1223 var_dict["Pressure Ratio"] = self.ratioP[time_point]
1224 if hasattr(self, "efficiency_pump"):
1225 var_dict["Efficiency"] = self.efficiency_pump[time_point]
1226 if hasattr(self, "efficiency_isentropic"):
1227 var_dict["Isentropic Efficiency"] = self.efficiency_isentropic[time_point]
1229 return {"vars": var_dict}
1231 def calculate_scaling_factors(self):
1232 super().calculate_scaling_factors()
1234 if hasattr(self, "work_fluid"):
1235 for t, v in self.work_fluid.items():
1236 iscale.set_scaling_factor(
1237 v,
1238 iscale.get_scaling_factor(
1239 self.control_volume.work[t], default=1, warning=True
1240 ),
1241 )
1243 if hasattr(self, "work_mechanical"):
1244 for t, v in self.work_mechanical.items():
1245 iscale.set_scaling_factor(
1246 v,
1247 iscale.get_scaling_factor(
1248 self.control_volume.work[t], default=1, warning=True
1249 ),
1250 )
1252 if hasattr(self, "work_isentropic"):
1253 for t, v in self.work_isentropic.items():
1254 iscale.set_scaling_factor(
1255 v,
1256 iscale.get_scaling_factor(
1257 self.control_volume.work[t], default=1, warning=True
1258 ),
1259 )
1261 if hasattr(self, "ratioP_calculation"):
1262 for t, c in self.ratioP_calculation.items():
1263 iscale.constraint_scaling_transform(
1264 c,
1265 iscale.get_scaling_factor(
1266 self.control_volume.properties_in[t].pressure,
1267 default=1,
1268 warning=True,
1269 ),
1270 overwrite=False,
1271 )
1273 if hasattr(self, "fluid_work_calculation"):
1274 for t, c in self.fluid_work_calculation.items():
1275 iscale.constraint_scaling_transform(
1276 c,
1277 iscale.get_scaling_factor(
1278 self.control_volume.deltaP[t], default=1, warning=True
1279 ),
1280 overwrite=False,
1281 )
1283 if hasattr(self, "actual_work"):
1284 for t, c in self.actual_work.items():
1285 iscale.constraint_scaling_transform(
1286 c,
1287 iscale.get_scaling_factor(
1288 self.control_volume.work[t], default=1, warning=True
1289 ),
1290 overwrite=False,
1291 )
1293 if hasattr(self, "isentropic_pressure"):
1294 for t, c in self.isentropic_pressure.items():
1295 iscale.constraint_scaling_transform(
1296 c,
1297 iscale.get_scaling_factor(
1298 self.control_volume.properties_in[t].pressure,
1299 default=1,
1300 warning=True,
1301 ),
1302 overwrite=False,
1303 )
1305 if hasattr(self, "isentropic"):
1306 for t, c in self.isentropic.items():
1307 iscale.constraint_scaling_transform(
1308 c,
1309 iscale.get_scaling_factor(
1310 self.control_volume.properties_in[t].entr_mol,
1311 default=1,
1312 warning=True,
1313 ),
1314 overwrite=False,
1315 )
1317 if hasattr(self, "isentropic_energy_balance"):
1318 for t, c in self.isentropic_energy_balance.items():
1319 iscale.constraint_scaling_transform(
1320 c,
1321 iscale.get_scaling_factor(
1322 self.control_volume.work[t], default=1, warning=True
1323 ),
1324 overwrite=False,
1325 )
1327 if hasattr(self, "zero_work_equation"):
1328 for t, c in self.zero_work_equation.items():
1329 iscale.constraint_scaling_transform(
1330 c,
1331 iscale.get_scaling_factor(
1332 self.control_volume.work[t], default=1, warning=True
1333 ),
1334 )
1336 if hasattr(self, "state_material_balances"):
1337 cvol = self.control_volume
1338 phase_list = cvol.properties_in.phase_list
1339 phase_component_set = cvol.properties_in.phase_component_set
1340 mb_type = cvol._constructed_material_balance_type
1341 if mb_type == MaterialBalanceType.componentPhase:
1342 for (t, p, j), c in self.state_material_balances.items():
1343 sf = iscale.get_scaling_factor(
1344 cvol.properties_in[t].get_material_flow_terms(p, j),
1345 default=1,
1346 warning=True,
1347 )
1348 iscale.constraint_scaling_transform(c, sf)
1349 elif mb_type == MaterialBalanceType.componentTotal:
1350 for (t, j), c in self.state_material_balances.items():
1351 sf = iscale.min_scaling_factor(
1352 [
1353 cvol.properties_in[t].get_material_flow_terms(p, j)
1354 for p in phase_list
1355 if (p, j) in phase_component_set
1356 ]
1357 )
1358 iscale.constraint_scaling_transform(c, sf)
1359 else:
1360 # There are some other material balance types but they create
1361 # constraints with different names.
1362 _log.warning(f"Unknown material balance type {mb_type}")
1365@declare_process_block_class("Turbine", doc="Isentropic turbine model")
1366class TurbineData(PressureChangerData):
1367 """
1368 Pressure changer with isentropic turbine options
1369 """
1371 default_initializer = IsentropicPressureChangerInitializer
1373 CONFIG = PressureChangerData.CONFIG()
1374 CONFIG.compressor = False
1375 CONFIG.get("compressor")._default = False
1376 CONFIG.get("compressor")._domain = In([False])
1377 CONFIG.thermodynamic_assumption = ThermodynamicAssumption.isentropic
1378 CONFIG.get("thermodynamic_assumption")._default = ThermodynamicAssumption.isentropic
1381@declare_process_block_class("Compressor", doc="Isentropic compressor model")
1382class CompressorData(PressureChangerData):
1383 """Pressure changer with isentropic compressor options"""
1385 default_initializer = IsentropicPressureChangerInitializer
1387 CONFIG = PressureChangerData.CONFIG()
1388 CONFIG.compressor = True
1389 CONFIG.get("compressor")._default = True
1390 CONFIG.get("compressor")._domain = In([True])
1391 CONFIG.thermodynamic_assumption = ThermodynamicAssumption.isentropic
1392 CONFIG.get("thermodynamic_assumption")._default = ThermodynamicAssumption.isentropic
1395@declare_process_block_class("Pump", doc="Pump model")
1396class PumpData(PressureChangerData):
1397 """Pressure changer with pump options"""
1399 CONFIG = PressureChangerData.CONFIG()
1400 CONFIG.compressor = True
1401 CONFIG.get("compressor")._default = True
1402 CONFIG.get("compressor")._domain = In([True])
1403 CONFIG.thermodynamic_assumption = ThermodynamicAssumption.pump
1404 CONFIG.get("thermodynamic_assumption")._default = ThermodynamicAssumption.pump