Coverage for backend/ahuora-builder/src/ahuora_builder/custom/thermal_utility_systems/willans_turbine.py: 65%
224 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-05-13 02:47 +0000
« prev ^ index » next coverage.py v7.10.7, created at 2026-05-13 02:47 +0000
1#################################################################################
2# The Institute for the Design of Advanced Energy Systems Integrated Platform
3# Framework (IDAES IP) was produced under the DOE Institute for the
4# Design of Advanced Energy Systems (IDAES).
5#
6# Copyright (c) 2018-2024 by the software owners: The Regents of the
7# University of California, through Lawrence Berkeley National Laboratory,
8# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
9# University, West Virginia University Research Corporation, et al.
10# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
11# for full copyright and license information.
12#################################################################################
13"""
14Standard IDAES pressure changer model.
15"""
16# TODO: Missing docstrings
17# pylint: disable=missing-function-docstring
19# Changing existing config block attributes
20# pylint: disable=protected-access
22# Import Python libraries
25# Import Pyomo libraries
26from pyomo.environ import (
27 Block,
28 value,
29 Var,
30 Expression,
31 Constraint,
32 Reference,
33 check_optimal_termination,
34 Reals,
35 Param,
36)
37from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool
39# Import IDAES cores
40from idaes.core import (
41 ControlVolume0DBlock,
42 declare_process_block_class,
43 EnergyBalanceType,
44 MomentumBalanceType,
45 MaterialBalanceType,
46 ProcessBlockData,
47 UnitModelBlockData,
48 useDefault,
49)
50from idaes.core.util.exceptions import PropertyNotSupportedError, InitializationError
51from idaes.core.util.config import is_physical_parameter_block
52import idaes.logger as idaeslog
53from idaes.core.util import scaling as iscale
54from idaes.core.solvers import get_solver
55from idaes.core.initialization import SingleControlVolumeUnitInitializer
56from idaes.core.util import to_json, from_json, StoreSpec
57from idaes.core.util.math import smooth_max, safe_sqrt, sqrt, smooth_min
58from pyomo.environ import units as pyunits
59from ahuora_builder.state_args import extract_state_args
62__author__ = "Emmanuel Ogbe, Andrew Lee"
63_log = idaeslog.getLogger(__name__)
66@declare_process_block_class("TurbineBase")
67class TurbineBaseData(UnitModelBlockData):
68 """
69 Standard Compressor/Expander Unit Model Class
70 """
72 CONFIG = UnitModelBlockData.CONFIG()
74 CONFIG.declare(
75 "material_balance_type",
76 ConfigValue(
77 default=MaterialBalanceType.useDefault,
78 domain=In(MaterialBalanceType),
79 description="Material balance construction flag",
80 doc="""Indicates what type of mass balance should be constructed,
81**default** - MaterialBalanceType.useDefault.
82**Valid values:** {
83**MaterialBalanceType.useDefault - refer to property package for default
84balance type
85**MaterialBalanceType.none** - exclude material balances,
86**MaterialBalanceType.componentPhase** - use phase component balances,
87**MaterialBalanceType.componentTotal** - use total component balances,
88**MaterialBalanceType.elementTotal** - use total element balances,
89**MaterialBalanceType.total** - use total material balance.}""",
90 ),
91 )
92 CONFIG.declare(
93 "energy_balance_type",
94 ConfigValue(
95 default=EnergyBalanceType.useDefault,
96 domain=In(EnergyBalanceType),
97 description="Energy balance construction flag",
98 doc="""Indicates what type of energy balance should be constructed,
99**default** - EnergyBalanceType.useDefault.
100**Valid values:** {
101**EnergyBalanceType.useDefault - refer to property package for default
102balance type
103**EnergyBalanceType.none** - exclude energy balances,
104**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material,
105**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase,
106**EnergyBalanceType.energyTotal** - single energy balance for material,
107**EnergyBalanceType.energyPhase** - energy balances for each phase.}""",
108 ),
109 )
110 CONFIG.declare(
111 "momentum_balance_type",
112 ConfigValue(
113 default=MomentumBalanceType.pressureTotal,
114 domain=In(MomentumBalanceType),
115 description="Momentum balance construction flag",
116 doc="""Indicates what type of momentum balance should be
117constructed, **default** - MomentumBalanceType.pressureTotal.
118**Valid values:** {
119**MomentumBalanceType.none** - exclude momentum balances,
120**MomentumBalanceType.pressureTotal** - single pressure balance for material,
121**MomentumBalanceType.pressurePhase** - pressure balances for each phase,
122**MomentumBalanceType.momentumTotal** - single momentum balance for material,
123**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""",
124 ),
125 )
126 CONFIG.declare(
127 "has_phase_equilibrium",
128 ConfigValue(
129 default=False,
130 domain=Bool,
131 description="Phase equilibrium construction flag",
132 doc="""Indicates whether terms for phase equilibrium should be
133constructed, **default** = False.
134**Valid values:** {
135**True** - include phase equilibrium terms
136**False** - exclude phase equilibrium terms.}""",
137 ),
138 )
139 CONFIG.declare(
140 "property_package",
141 ConfigValue(
142 default=useDefault,
143 domain=is_physical_parameter_block,
144 description="Property package to use for control volume",
145 doc="""Property parameter object used to define property
146calculations, **default** - useDefault.
147**Valid values:** {
148**useDefault** - use default package from parent model or flowsheet,
149**PropertyParameterObject** - a PropertyParameterBlock object.}""",
150 ),
151 )
152 CONFIG.declare(
153 "property_package_args",
154 ConfigBlock(
155 implicit=True,
156 description="Arguments to use for constructing property packages",
157 doc="""A ConfigBlock with arguments to be passed to a property
158block(s) and used when constructing these,
159**default** - None.
160**Valid values:** {
161see property package for documentation.}""",
162 ),
163 )
164 CONFIG.declare(
165 "calculation_method",
166 ConfigValue(
167 default='isentropic',
168 domain=In(["isentropic", "simple_willans", "part_load_willans", "Tsat_willans", "BPST_willans", "CT_willans"]),
169 description="Calculation method used to model mechanical work",
170 doc="""Property parameter object used to define property
171calculations, **default** - 'isentropic'.
172**Valid values:** {
173**isentropic** - default method, uses isentropic efficiency to determine work
174**simple_willans** - simple willans line requring slope, intercept and max molar flow.
175**part_load_willans** - willans line with part load correction, a, b, c and max molar flow as parameters
176**Tsat_willans** - willans line with part load correction using saturation temperature. Requires only max molar flow
177**BPST_willans** - back pressure willans line with part load correction using pressure difference, requires max molar flow.
178**CT_willans** - condensing turbine willans line with part load correction using pressure difference, requires max molar flow.
179}""",
180 ),
181 )
183 def build(self):
184 """
186 Args:
187 None
189 Returns:
190 None
191 """
192 # Call UnitModel.build
193 super().build()
195 # Add a control volume to the unit including setting up dynamics.
196 self.control_volume = ControlVolume0DBlock(
197 dynamic=self.config.dynamic,
198 has_holdup=self.config.has_holdup,
199 property_package=self.config.property_package,
200 property_package_args=self.config.property_package_args,
201 )
203 # Add geometry variables to control volume
204 if self.config.has_holdup: 204 ↛ 205line 204 didn't jump to line 205 because the condition on line 204 was never true
205 self.control_volume.add_geometry()
207 # Add inlet and outlet state blocks to control volume
208 self.control_volume.add_state_blocks(
209 has_phase_equilibrium=self.config.has_phase_equilibrium
210 )
212 # Add mass balance
213 # Set has_equilibrium is False for now
214 # TO DO; set has_equilibrium to True
215 self.control_volume.add_material_balances(
216 balance_type=self.config.material_balance_type,
217 has_phase_equilibrium=self.config.has_phase_equilibrium,
218 )
220 # Add energy balance
221 eb = self.control_volume.add_energy_balances(
222 balance_type=self.config.energy_balance_type, has_work_transfer=True
223 )
225 # add momentum balance
226 self.control_volume.add_momentum_balances(
227 balance_type=self.config.momentum_balance_type, has_pressure_change=True
228 )
230 # Add Ports
231 self.add_inlet_port()
232 self.add_outlet_port()
234 # Set Unit Geometry and holdup Volume
235 if self.config.has_holdup is True: 235 ↛ 236line 235 didn't jump to line 236 because the condition on line 235 was never true
236 self.volume = Reference(self.control_volume.volume[:])
238 self.work_mechanical = Reference(self.control_volume.work[:])
239 # self.work_mechanical.fix(-1000e3)
240 # self.work_mechanical.unfix()
243 # Add Momentum balance variable 'deltaP'
244 self.deltaP = Reference(self.control_volume.deltaP[:])
246 # Performance Variables
247 self.ratioP = Var(self.flowsheet().time, initialize=1.0, doc="Pressure Ratio")
249 # Pressure Ratio
250 @self.Constraint(self.flowsheet().time, doc="Pressure ratio constraint")
251 def ratioP_calculation(self, t):
252 return (
253 self.ratioP[t] * self.control_volume.properties_in[t].pressure
254 == self.control_volume.properties_out[t].pressure
255 )
257 units_meta = self.control_volume.config.property_package.get_metadata()
259 # Get indexing sets from control volume
260 # Add isentropic variables
261 self.efficiency_isentropic = Var(
262 self.flowsheet().time,
263 initialize=0.5,
264 doc="Efficiency with respect to an isentropic process [-]",
265 )
267 # self.delta_h_is = Var(
268 # self.flowsheet().time,
269 # initialize=-100e3,
270 # doc="Enthalpy difference input to unit if isentropic process",
271 # units=units_meta.get_derived_units("energy") / units_meta.get_derived_units("amount"),
272 # )
273 # self.delta_h_act = Var(
274 # self.flowsheet().time,
275 # initialize=-100e3,
276 # doc="Enthalpy difference input to unit for actual process",
277 # units=units_meta.get_derived_units("energy") / units_meta.get_derived_units("amount"),
278 # )
279 # self.work_isentropic = Var(
280 # self.flowsheet().time,
281 # initialize=-100e3,
282 # doc="Work input to unit if isentropic process",
283 # units=units_meta.get_derived_units("power"),
284 # )
286 # Add motor/electrical work and efficiency variable
287 self.efficiency_motor = Var(
288 self.flowsheet().time,
289 initialize=1.0,
290 doc="Motor efficiency converting shaft work to electrical work [-]",
291 )
293 self.work_electrical = Var(
294 self.flowsheet().time,
295 initialize=1.0,
296 doc="Electrical work of a turbine [-]",
297 units=units_meta.get_derived_units("power")
298 )
300 # Add willans line parameters
301 if 'willans' in self.config.calculation_method:
302 self.willans_slope = Var(
303 self.flowsheet().time,
304 initialize=100,
305 doc="Slope of willans line",
306 units=units_meta.get_derived_units("energy") / units_meta.get_derived_units("amount"),
307 )
309 self.willans_intercept = Var(
310 self.flowsheet().time,
311 initialize=-100,
312 doc="Intercept of willans line",
313 units=units_meta.get_derived_units("power"),
314 )
316 self.willans_max_mol = Var(
317 self.flowsheet().time,
318 initialize=1.0,
319 doc="Max molar flow of willans line",
320 units=units_meta.get_derived_units("amount") / units_meta.get_derived_units("time"),
321 )
323 if self.config.calculation_method in ["part_load_willans", "Tsat_willans", "BPST_willans", "CT_willans"]:
324 self.willans_a = Var(
325 self.flowsheet().time,
326 initialize=1.0,
327 doc="Willans a coefficient",
328 )
330 self.willans_b = Var(
331 self.flowsheet().time,
332 initialize=1.0,
333 doc="Willans b coefficient",
334 units=units_meta.get_derived_units("power")
335 )
337 self.willans_efficiency = Var(
338 self.flowsheet().time,
339 initialize=1.0,
340 doc="Willans efficiency",
341 )
343 # Build isentropic state block
344 tmp_dict = dict(**self.config.property_package_args)
345 tmp_dict["has_phase_equilibrium"] = self.config.has_phase_equilibrium
346 tmp_dict["defined_state"] = False
348 self.properties_isentropic = self.config.property_package.build_state_block(
349 self.flowsheet().time, doc="isentropic properties at outlet", **tmp_dict
350 )
352 # Connect isentropic state block properties
353 @self.Constraint(
354 self.flowsheet().time, doc="Pressure for isentropic calculations"
355 )
356 def isentropic_pressure(self, t):
357 return (
358 self.properties_isentropic[t].pressure
359 == self.control_volume.properties_out[t].pressure
360 )
362 # This assumes isentropic composition is the same as outlet
363 self.add_state_material_balances(
364 self.config.material_balance_type,
365 self.properties_isentropic,
366 self.control_volume.properties_out,
367 )
369 # This assumes isentropic entropy is the same as inlet
370 @self.Constraint(self.flowsheet().time, doc="Isentropic assumption")
371 def isentropic(self, t):
372 return (
373 self.properties_isentropic[t].entr_mol
374 == self.control_volume.properties_in[t].entr_mol
375 )
377 self.add_isentropic_work_definition()
379 if 'willans' in self.config.calculation_method:
380 # Write isentropic efficiency eqn
381 self.add_willans_line_relationship()
383 if self.config.calculation_method in ["part_load_willans", "Tsat_willans", "BPST_willans", "CT_willans"]:
384 if self.config.calculation_method == "Tsat_willans": # use published values and dTsat to calculate willans a,b,c
385 self.calculate_Tsat_willans_parameters()
387 elif self.config.calculation_method == "BPST_willans":
388 self.calculate_BPST_willans_parameters()
390 elif self.config.calculation_method == "CT_willans":
391 self.calculate_CT_willans_parameters()
393 # calculate slope and intercept
394 self.calculate_willans_coefficients()
396 self.add_mechanical_and_isentropic_work_definition()
397 self.add_electrical_work_definition()
399 def calculate_CT_willans_parameters(self):
400 # a parameter
401 @self.Constraint(
402 self.flowsheet().time, doc="Willans CT a calculation"
403 )
404 def willans_CT_a_calculation(self, t):
405 return self.willans_a[t] == 1.288464 - 0.0015185 * (self.control_volume.properties_in[t].pressure / 1e5) / pyunits.Pa - 0.33415834 * (self.control_volume.properties_out[t].pressure / 1e5) / pyunits.Pa
407 # b parameter
408 @self.Constraint(
409 self.flowsheet().time, doc="Willans CT b calculation"
410 )
411 def willans_CT_b_calculation(self, t):
412 return self.willans_b[t] == (-437.7746025 + 29.00736723 * (self.control_volume.properties_in[t].pressure / 1e5) / pyunits.Pa + 10.35902331 * (self.control_volume.properties_out[t].pressure / 1e5) / pyunits.Pa) * 1000 * pyunits.W
414 # c parameter
415 @self.Constraint(
416 self.flowsheet().time, doc="Willans CT efficiency calculation"
417 )
418 def willans_CT_efficiency_calculation(self, t):
419 return self.willans_efficiency[t] == 1 / ((0.07886297 + 0.000528327 * (self.control_volume.properties_in[t].pressure / 1e5) / pyunits.Pa - 0.703153891 * (self.control_volume.properties_out[t].pressure / 1e5) / pyunits.Pa) + 1)
421 def calculate_BPST_willans_parameters(self):
422 # a parameter
423 @self.Constraint(
424 self.flowsheet().time, doc="Willans BPST a calculation"
425 )
426 def willans_BPST_a_calculation(self, t):
427 return self.willans_a[t] == 1.18795366 - 0.00029564 * (self.control_volume.properties_in[t].pressure / 1e5) / pyunits.Pa + 0.004647288 * (self.control_volume.properties_out[t].pressure / 1e5) / pyunits.Pa
430 # b parameter
431 @self.Constraint(
432 self.flowsheet().time, doc="Willans BPST b calculation"
433 )
434 def willans_BPST_b_calculation(self, t):
435 return self.willans_b[t] == (449.9767142 + 5.670176939 * (self.control_volume.properties_in[t].pressure / 1e5) / pyunits.Pa - 11.5045814 * (self.control_volume.properties_out[t].pressure / 1e5) / pyunits.Pa) * 1000 * pyunits.W
438 # c parameter
439 @self.Constraint(
440 self.flowsheet().time, doc="Willans BPST c calculation"
441 )
442 def willans_BPST_efficiency_calculation(self, t):
443 return self.willans_efficiency[t] == 1 / ((0.205149333 - 0.000695171 * (self.control_volume.properties_in[t].pressure / 1e5) / pyunits.Pa + 0.002844611 * (self.control_volume.properties_out[t].pressure / 1e5) / pyunits.Pa) + 1)
445 def calculate_Tsat_willans_parameters(self):
446 # a parameter
447 @self.Constraint(
448 self.flowsheet().time, doc="Willans Tsat a calculation"
449 )
450 def willans_Tsat_a_calculation(self, t):
451 return self.willans_a[t] == (1.155 + 0.000538 * (self.control_volume.properties_in[t].temperature_sat - self.control_volume.properties_out[t].temperature_sat) / pyunits.K)
453 # b parameter
454 @self.Constraint(
455 self.flowsheet().time, doc="Willans Tsat b calculation"
456 )
457 def willans_Tsat_b_calculation(self, t):
458 return self.willans_b[t] == (0 + 4.23 * (self.control_volume.properties_in[t].temperature_sat - self.control_volume.properties_out[t].temperature_sat) / pyunits.K)*1000 * pyunits.W
460 # c parameter
461 @self.Constraint(
462 self.flowsheet().time, doc="Willans Tsat efficiency calculation"
463 )
464 def willans_Tsat_c_calculation(self, t):
465 return self.willans_efficiency[t] == 0.83333
467 def calculate_willans_coefficients(self):
468 # Calculate willans coefficients
469 @self.Constraint(
470 self.flowsheet().time, doc="Willans slope calculation"
471 )
472 def willans_slope_calculation(self, t):
473 return self.willans_slope[t] == 1 / (self.willans_efficiency[t] * self.willans_a[t]) * ((self.control_volume.properties_in[t].enth_mol - self.properties_isentropic[t].enth_mol) - self.willans_b[t] / self.willans_max_mol[t])
476 @self.Constraint(
477 self.flowsheet().time, doc="Willans intercept calculation"
478 )
479 def willans_intercept_calculation(self, t):
480 return self.willans_intercept[t] == ((1 - self.willans_efficiency[t]) / (self.willans_efficiency[t] * self.willans_a[t])) * ((self.control_volume.properties_in[t].enth_mol - self.properties_isentropic[t].enth_mol) * self.willans_max_mol[t] - self.willans_b[t])
482 def add_mechanical_and_isentropic_work_definition(self):
484 if 'willans' in self.config.calculation_method:
485 self.efficiency_isentropic = Expression(
486 self.flowsheet().time,
487 rule=lambda b, t: (
488 b.delta_h_act[t] / b.delta_h_is[t]
489 )
490 )
491 else:
492 # Mechanical work
493 @self.Constraint(
494 self.flowsheet().time, doc="Isentropic and mechanical work relationship"
495 )
496 def isentropic_and_mechanical_work_eq(self, t):
497 return self.work_mechanical[t] == (
498 self.delta_h_is[t] * self.control_volume.properties_in[t].flow_mol * self.efficiency_isentropic[t]
499 )
501 def add_willans_line_relationship(self):
502 @self.Constraint(
503 self.flowsheet().time, doc="Willans line and mechanical work relationship"
504 )
505 def willans_line_eq(self, t):
506 eps = 0.0001 # smoothing parameter; smaller = closer to exact max, larger = smoother
507 return self.work_mechanical[t] == smooth_min(
508 -(self.willans_slope[t] * self.control_volume.properties_in[t].flow_mol - self.willans_intercept[t]) / (self.willans_slope[t] * self.willans_max_mol[t]),
509 0.0,
510 eps
511 ) * (self.willans_slope[t] * self.willans_max_mol[t])
513 def add_electrical_work_definition(self):
514 # Electrical work
515 @self.Constraint(
516 self.flowsheet().time, doc="Calculate electrical work of turbine"
517 )
518 def electrical_energy_balance(self, t):
519 return self.work_electrical[t] == self.work_mechanical[t] * self.efficiency_motor[t]
521 def add_isentropic_work_definition(self):
522 self.delta_h_is = Expression(
523 self.flowsheet().time,
524 rule=lambda b, t: (
525 b.properties_isentropic[t].enth_mol - b.control_volume.properties_in[t].enth_mol
526 )
527 )
528 self.delta_h_act = Expression(
529 self.flowsheet().time,
530 rule=lambda b, t: (
531 b.control_volume.properties_out[t].enth_mol - b.control_volume.properties_in[t].enth_mol
532 )
533 )
534 # # Isentropic work
535 # @self.Constraint(
536 # self.flowsheet().time, doc="Calculate work of isentropic process"
537 # )
538 # def isentropic_energy_balance(self, t):
539 # return self.work_isentropic[t] == ( self.properties_isentropic[t].enth_mol - self.control_volume.properties_in[t].enth_mol ) * self.control_volume.properties_in[t].flow_mol
541 def initialize_build(
542 blk,
543 state_args=None,
544 routine=None,
545 outlvl=idaeslog.NOTSET,
546 solver=None,
547 optarg=None,
548 ):
549 """
550 General wrapper for pressure changer initialization routines
552 Keyword Arguments:
553 routine : str stating which initialization routine to execute
554 * None - use routine matching thermodynamic_assumption
555 * 'isentropic' - use isentropic initialization routine
556 * 'isothermal' - use isothermal initialization routine
557 state_args : a dict of arguments to be passed to the property
558 package(s) to provide an initial state for
559 initialization (see documentation of the specific
560 property package) (default = {}).
561 outlvl : sets output level of initialization routine
562 optarg : solver options dictionary object (default=None, use
563 default solver options)
564 solver : str indicating which solver to use during
565 initialization (default = None, use default solver)
567 Returns:
568 None
569 """
570 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
571 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
573 # Create solver
574 opt = get_solver(solver, optarg)
576 cv = blk.control_volume
577 t0 = blk.flowsheet().time.first()
578 state_args_out = {}
580 if state_args is None: 580 ↛ 584line 580 didn't jump to line 584 because the condition on line 580 was always true
581 state_args = extract_state_args(cv.properties_in[t0], use_state_vars=False)
583 # Get initialisation guesses for outlet and isentropic states
584 for k in state_args:
585 if k == "pressure" and k not in state_args_out:
586 # Work out how to estimate outlet pressure
587 if cv.properties_out[t0].pressure.fixed:
588 # Fixed outlet pressure, use this value
589 state_args_out[k] = value(cv.properties_out[t0].pressure)
590 elif blk.deltaP[t0].fixed: 590 ↛ 591line 590 didn't jump to line 591 because the condition on line 590 was never true
591 state_args_out[k] = value(state_args[k] + blk.deltaP[t0])
592 elif blk.ratioP[t0].fixed: 592 ↛ 593line 592 didn't jump to line 593 because the condition on line 592 was never true
593 state_args_out[k] = value(state_args[k] * blk.ratioP[t0])
594 else:
595 # Not obvious what to do, use inlet state
596 state_args_out[k] = state_args[k]
597 elif k not in state_args_out: 597 ↛ 584line 597 didn't jump to line 584 because the condition on line 597 was always true
598 state_args_out[k] = state_args[k]
600 # Initialize state blocks
601 flags = cv.properties_in.initialize(
602 outlvl=outlvl,
603 optarg=optarg,
604 solver=solver,
605 hold_state=True,
606 state_args=state_args,
607 )
608 cv.properties_out.initialize(
609 outlvl=outlvl,
610 optarg=optarg,
611 solver=solver,
612 hold_state=False,
613 state_args=state_args_out,
614 )
616 init_log.info_high("Initialization Step 1 Complete.")
617 # ---------------------------------------------------------------------
618 # Initialize Isentropic block
620 blk.properties_isentropic.initialize(
621 outlvl=outlvl,
622 optarg=optarg,
623 solver=solver,
624 state_args=state_args_out,
625 )
627 init_log.info_high("Initialization Step 2 Complete.")
629 # Skipping step 3 because Isothermal had problems.
631 # ---------------------------------------------------------------------
632 # Solve unit
633 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
634 res = opt.solve(blk, tee=slc.tee)
635 init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res)))
638 # ---------------------------------------------------------------------
639 # Release Inlet state
640 blk.control_volume.release_state(flags, outlvl)
642 if not check_optimal_termination(res): 642 ↛ 643line 642 didn't jump to line 643 because the condition on line 642 was never true
643 raise InitializationError(
644 f"{blk.name} failed to initialize successfully. Please check "
645 f"the output logs for more information."
646 )
648 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
650 def _get_performance_contents(self, time_point=0):
651 var_dict = {}
652 if hasattr(self, "deltaP"): 652 ↛ 654line 652 didn't jump to line 654 because the condition on line 652 was always true
653 var_dict["Mechanical Work"] = self.work_mechanical[time_point]
654 if hasattr(self, "deltaP"): 654 ↛ 656line 654 didn't jump to line 656 because the condition on line 654 was always true
655 var_dict["Electrical Work"] = self.work_electrical[time_point]
656 if hasattr(self, "deltaP"): 656 ↛ 658line 656 didn't jump to line 658 because the condition on line 656 was always true
657 var_dict["Pressure Change"] = self.deltaP[time_point]
658 if hasattr(self, "ratioP"): 658 ↛ 663line 658 didn't jump to line 663 because the condition on line 658 was always true
659 var_dict["Pressure Ratio"] = self.ratioP[time_point]
660 # if hasattr(self, "efficiency_isentropic"):
661 # var_dict["Isentropic Efficiency"] = self.efficiency_isentropic[time_point]
663 return {"vars": var_dict}
665 def calculate_scaling_factors(self):
666 super().calculate_scaling_factors()
668 if hasattr(self, "work_fluid"):
669 for t, v in self.work_fluid.items():
670 iscale.set_scaling_factor(
671 v,
672 iscale.get_scaling_factor(
673 self.control_volume.work[t], default=1, warning=True
674 ),
675 )
677 if hasattr(self, "work_mechanical"):
678 for t, v in self.work_mechanical.items():
679 iscale.set_scaling_factor(
680 v,
681 iscale.get_scaling_factor(
682 self.control_volume.work[t], default=1, warning=True
683 ),
684 )
686 if hasattr(self, "work_isentropic"):
687 for t, v in self.work_isentropic.items():
688 iscale.set_scaling_factor(
689 v,
690 iscale.get_scaling_factor(
691 self.control_volume.work[t], default=1, warning=True
692 ),
693 )
695 if hasattr(self, "ratioP_calculation"):
696 for t, c in self.ratioP_calculation.items():
697 iscale.constraint_scaling_transform(
698 c,
699 iscale.get_scaling_factor(
700 self.control_volume.properties_in[t].pressure,
701 default=1,
702 warning=True,
703 ),
704 overwrite=False,
705 )
707 if hasattr(self, "fluid_work_calculation"):
708 for t, c in self.fluid_work_calculation.items():
709 iscale.constraint_scaling_transform(
710 c,
711 iscale.get_scaling_factor(
712 self.control_volume.deltaP[t], default=1, warning=True
713 ),
714 overwrite=False,
715 )
717 if hasattr(self, "actual_work"):
718 for t, c in self.actual_work.items():
719 iscale.constraint_scaling_transform(
720 c,
721 iscale.get_scaling_factor(
722 self.control_volume.work[t], default=1, warning=True
723 ),
724 overwrite=False,
725 )
727 if hasattr(self, "isentropic_pressure"):
728 for t, c in self.isentropic_pressure.items():
729 iscale.constraint_scaling_transform(
730 c,
731 iscale.get_scaling_factor(
732 self.control_volume.properties_in[t].pressure,
733 default=1,
734 warning=True,
735 ),
736 overwrite=False,
737 )
739 if hasattr(self, "isentropic"):
740 for t, c in self.isentropic.items():
741 iscale.constraint_scaling_transform(
742 c,
743 iscale.get_scaling_factor(
744 self.control_volume.properties_in[t].entr_mol,
745 default=1,
746 warning=True,
747 ),
748 overwrite=False,
749 )
751 if hasattr(self, "isentropic_energy_balance"):
752 for t, c in self.isentropic_energy_balance.items():
753 iscale.constraint_scaling_transform(
754 c,
755 iscale.get_scaling_factor(
756 self.control_volume.work[t], default=1, warning=True
757 ),
758 overwrite=False,
759 )
761 if hasattr(self, "zero_work_equation"):
762 for t, c in self.zero_work_equation.items():
763 iscale.constraint_scaling_transform(
764 c,
765 iscale.get_scaling_factor(
766 self.control_volume.work[t], default=1, warning=True
767 ),
768 )
770 if hasattr(self, "state_material_balances"):
771 cvol = self.control_volume
772 phase_list = cvol.properties_in.phase_list
773 phase_component_set = cvol.properties_in.phase_component_set
774 mb_type = cvol._constructed_material_balance_type
775 if mb_type == MaterialBalanceType.componentPhase:
776 for (t, p, j), c in self.state_material_balances.items():
777 sf = iscale.get_scaling_factor(
778 cvol.properties_in[t].get_material_flow_terms(p, j),
779 default=1,
780 warning=True,
781 )
782 iscale.constraint_scaling_transform(c, sf)
783 elif mb_type == MaterialBalanceType.componentTotal:
784 for (t, j), c in self.state_material_balances.items():
785 sf = iscale.min_scaling_factor(
786 [
787 cvol.properties_in[t].get_material_flow_terms(p, j)
788 for p in phase_list
789 if (p, j) in phase_component_set
790 ]
791 )
792 iscale.constraint_scaling_transform(c, sf)
793 else:
794 # There are some other material balance types but they create
795 # constraints with different names.
796 _log.warning(f"Unknown material balance type {mb_type}")