Coverage for backend/idaes_service/solver/custom/thermal_utility_systems/willans_turbine.py: 66%
230 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +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
61__author__ = "Emmanuel Ogbe, Andrew Lee"
62_log = idaeslog.getLogger(__name__)
65@declare_process_block_class("TurbineBase")
66class TurbineBaseData(UnitModelBlockData):
67 """
68 Standard Compressor/Expander Unit Model Class
69 """
71 CONFIG = UnitModelBlockData.CONFIG()
73 CONFIG.declare(
74 "material_balance_type",
75 ConfigValue(
76 default=MaterialBalanceType.useDefault,
77 domain=In(MaterialBalanceType),
78 description="Material balance construction flag",
79 doc="""Indicates what type of mass balance should be constructed,
80**default** - MaterialBalanceType.useDefault.
81**Valid values:** {
82**MaterialBalanceType.useDefault - refer to property package for default
83balance type
84**MaterialBalanceType.none** - exclude material balances,
85**MaterialBalanceType.componentPhase** - use phase component balances,
86**MaterialBalanceType.componentTotal** - use total component balances,
87**MaterialBalanceType.elementTotal** - use total element balances,
88**MaterialBalanceType.total** - use total material balance.}""",
89 ),
90 )
91 CONFIG.declare(
92 "energy_balance_type",
93 ConfigValue(
94 default=EnergyBalanceType.useDefault,
95 domain=In(EnergyBalanceType),
96 description="Energy balance construction flag",
97 doc="""Indicates what type of energy balance should be constructed,
98**default** - EnergyBalanceType.useDefault.
99**Valid values:** {
100**EnergyBalanceType.useDefault - refer to property package for default
101balance type
102**EnergyBalanceType.none** - exclude energy balances,
103**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material,
104**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase,
105**EnergyBalanceType.energyTotal** - single energy balance for material,
106**EnergyBalanceType.energyPhase** - energy balances for each phase.}""",
107 ),
108 )
109 CONFIG.declare(
110 "momentum_balance_type",
111 ConfigValue(
112 default=MomentumBalanceType.pressureTotal,
113 domain=In(MomentumBalanceType),
114 description="Momentum balance construction flag",
115 doc="""Indicates what type of momentum balance should be
116constructed, **default** - MomentumBalanceType.pressureTotal.
117**Valid values:** {
118**MomentumBalanceType.none** - exclude momentum balances,
119**MomentumBalanceType.pressureTotal** - single pressure balance for material,
120**MomentumBalanceType.pressurePhase** - pressure balances for each phase,
121**MomentumBalanceType.momentumTotal** - single momentum balance for material,
122**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""",
123 ),
124 )
125 CONFIG.declare(
126 "has_phase_equilibrium",
127 ConfigValue(
128 default=False,
129 domain=Bool,
130 description="Phase equilibrium construction flag",
131 doc="""Indicates whether terms for phase equilibrium should be
132constructed, **default** = False.
133**Valid values:** {
134**True** - include phase equilibrium terms
135**False** - exclude phase equilibrium terms.}""",
136 ),
137 )
138 CONFIG.declare(
139 "property_package",
140 ConfigValue(
141 default=useDefault,
142 domain=is_physical_parameter_block,
143 description="Property package to use for control volume",
144 doc="""Property parameter object used to define property
145calculations, **default** - useDefault.
146**Valid values:** {
147**useDefault** - use default package from parent model or flowsheet,
148**PropertyParameterObject** - a PropertyParameterBlock object.}""",
149 ),
150 )
151 CONFIG.declare(
152 "property_package_args",
153 ConfigBlock(
154 implicit=True,
155 description="Arguments to use for constructing property packages",
156 doc="""A ConfigBlock with arguments to be passed to a property
157block(s) and used when constructing these,
158**default** - None.
159**Valid values:** {
160see property package for documentation.}""",
161 ),
162 )
163 CONFIG.declare(
164 "calculation_method",
165 ConfigValue(
166 default='isentropic',
167 domain=In(["isentropic", "simple_willans", "part_load_willans", "Tsat_willans", "BPST_willans", "CT_willans"]),
168 description="Calculation method used to model mechanical work",
169 doc="""Property parameter object used to define property
170calculations, **default** - 'isentropic'.
171**Valid values:** {
172**isentropic** - default method, uses isentropic efficiency to determine work
173**simple_willans** - simple willans line requring slope, intercept and max molar flow.
174**part_load_willans** - willans line with part load correction, a, b, c and max molar flow as parameters
175**Tsat_willans** - willans line with part load correction using saturation temperature. Requires only max molar flow
176**BPST_willans** - back pressure willans line with part load correction using pressure difference, requires max molar flow.
177**CT_willans** - condensing turbine willans line with part load correction using pressure difference, requires max molar flow.
178}""",
179 ),
180 )
182 def build(self):
183 """
185 Args:
186 None
188 Returns:
189 None
190 """
191 # Call UnitModel.build
192 super().build()
194 # Add a control volume to the unit including setting up dynamics.
195 self.control_volume = ControlVolume0DBlock(
196 dynamic=self.config.dynamic,
197 has_holdup=self.config.has_holdup,
198 property_package=self.config.property_package,
199 property_package_args=self.config.property_package_args,
200 )
202 # Add geometry variables to control volume
203 if self.config.has_holdup: 203 ↛ 204line 203 didn't jump to line 204 because the condition on line 203 was never true
204 self.control_volume.add_geometry()
206 # Add inlet and outlet state blocks to control volume
207 self.control_volume.add_state_blocks(
208 has_phase_equilibrium=self.config.has_phase_equilibrium
209 )
211 # Add mass balance
212 # Set has_equilibrium is False for now
213 # TO DO; set has_equilibrium to True
214 self.control_volume.add_material_balances(
215 balance_type=self.config.material_balance_type,
216 has_phase_equilibrium=self.config.has_phase_equilibrium,
217 )
219 # Add energy balance
220 eb = self.control_volume.add_energy_balances(
221 balance_type=self.config.energy_balance_type, has_work_transfer=True
222 )
224 # add momentum balance
225 self.control_volume.add_momentum_balances(
226 balance_type=self.config.momentum_balance_type, has_pressure_change=True
227 )
229 # Add Ports
230 self.add_inlet_port()
231 self.add_outlet_port()
233 # Set Unit Geometry and holdup Volume
234 if self.config.has_holdup is True: 234 ↛ 235line 234 didn't jump to line 235 because the condition on line 234 was never true
235 self.volume = Reference(self.control_volume.volume[:])
237 self.work_mechanical = Reference(self.control_volume.work[:])
238 # self.work_mechanical.fix(-1000e3)
239 # self.work_mechanical.unfix()
242 # Add Momentum balance variable 'deltaP'
243 self.deltaP = Reference(self.control_volume.deltaP[:])
245 # Performance Variables
246 self.ratioP = Var(self.flowsheet().time, initialize=1.0, doc="Pressure Ratio")
248 # Pressure Ratio
249 @self.Constraint(self.flowsheet().time, doc="Pressure ratio constraint")
250 def ratioP_calculation(self, t):
251 return (
252 self.ratioP[t] * self.control_volume.properties_in[t].pressure
253 == self.control_volume.properties_out[t].pressure
254 )
256 units_meta = self.control_volume.config.property_package.get_metadata()
258 # Get indexing sets from control volume
259 # Add isentropic variables
260 self.efficiency_isentropic = Var(
261 self.flowsheet().time,
262 initialize=0.5,
263 doc="Efficiency with respect to an isentropic process [-]",
264 )
266 # self.delta_h_is = Var(
267 # self.flowsheet().time,
268 # initialize=-100e3,
269 # doc="Enthalpy difference input to unit if isentropic process",
270 # units=units_meta.get_derived_units("energy") / units_meta.get_derived_units("amount"),
271 # )
272 # self.delta_h_act = Var(
273 # self.flowsheet().time,
274 # initialize=-100e3,
275 # doc="Enthalpy difference input to unit for actual process",
276 # units=units_meta.get_derived_units("energy") / units_meta.get_derived_units("amount"),
277 # )
278 # self.work_isentropic = Var(
279 # self.flowsheet().time,
280 # initialize=-100e3,
281 # doc="Work input to unit if isentropic process",
282 # units=units_meta.get_derived_units("power"),
283 # )
285 # Add motor/electrical work and efficiency variable
286 self.efficiency_motor = Var(
287 self.flowsheet().time,
288 initialize=1.0,
289 doc="Motor efficiency converting shaft work to electrical work [-]",
290 )
292 self.work_electrical = Var(
293 self.flowsheet().time,
294 initialize=1.0,
295 doc="Electrical work of a turbine [-]",
296 units=units_meta.get_derived_units("power")
297 )
299 # Add willans line parameters
300 if 'willans' in self.config.calculation_method:
301 self.willans_slope = Var(
302 self.flowsheet().time,
303 initialize=100,
304 doc="Slope of willans line",
305 units=units_meta.get_derived_units("energy") / units_meta.get_derived_units("amount"),
306 )
308 self.willans_intercept = Var(
309 self.flowsheet().time,
310 initialize=-100,
311 doc="Intercept of willans line",
312 units=units_meta.get_derived_units("power"),
313 )
315 self.willans_max_mol = Var(
316 self.flowsheet().time,
317 initialize=1.0,
318 doc="Max molar flow of willans line",
319 units=units_meta.get_derived_units("amount") / units_meta.get_derived_units("time"),
320 )
322 if self.config.calculation_method in ["part_load_willans", "Tsat_willans", "BPST_willans", "CT_willans"]:
323 self.willans_a = Var(
324 self.flowsheet().time,
325 initialize=1.0,
326 doc="Willans a coefficient",
327 )
329 self.willans_b = Var(
330 self.flowsheet().time,
331 initialize=1.0,
332 doc="Willans b coefficient",
333 units=units_meta.get_derived_units("power")
334 )
336 self.willans_efficiency = Var(
337 self.flowsheet().time,
338 initialize=1.0,
339 doc="Willans efficiency",
340 )
342 # Build isentropic state block
343 tmp_dict = dict(**self.config.property_package_args)
344 tmp_dict["has_phase_equilibrium"] = self.config.has_phase_equilibrium
345 tmp_dict["defined_state"] = False
347 self.properties_isentropic = self.config.property_package.build_state_block(
348 self.flowsheet().time, doc="isentropic properties at outlet", **tmp_dict
349 )
351 # Connect isentropic state block properties
352 @self.Constraint(
353 self.flowsheet().time, doc="Pressure for isentropic calculations"
354 )
355 def isentropic_pressure(self, t):
356 return (
357 self.properties_isentropic[t].pressure
358 == self.control_volume.properties_out[t].pressure
359 )
361 # This assumes isentropic composition is the same as outlet
362 self.add_state_material_balances(
363 self.config.material_balance_type,
364 self.properties_isentropic,
365 self.control_volume.properties_out,
366 )
368 # This assumes isentropic entropy is the same as inlet
369 @self.Constraint(self.flowsheet().time, doc="Isentropic assumption")
370 def isentropic(self, t):
371 return (
372 self.properties_isentropic[t].entr_mol
373 == self.control_volume.properties_in[t].entr_mol
374 )
376 self.add_isentropic_work_definition()
378 if 'willans' in self.config.calculation_method:
379 # Write isentropic efficiency eqn
380 self.add_willans_line_relationship()
382 if self.config.calculation_method in ["part_load_willans", "Tsat_willans", "BPST_willans", "CT_willans"]:
383 if self.config.calculation_method == "Tsat_willans": # use published values and dTsat to calculate willans a,b,c
384 self.calculate_Tsat_willans_parameters()
386 elif self.config.calculation_method == "BPST_willans":
387 self.calculate_BPST_willans_parameters()
389 elif self.config.calculation_method == "CT_willans":
390 self.calculate_CT_willans_parameters()
392 # calculate slope and intercept
393 self.calculate_willans_coefficients()
395 self.add_mechanical_and_isentropic_work_definition()
396 self.add_electrical_work_definition()
398 def calculate_CT_willans_parameters(self):
399 # a parameter
400 @self.Constraint(
401 self.flowsheet().time, doc="Willans CT a calculation"
402 )
403 def willans_CT_a_calculation(self, t):
404 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
406 # b parameter
407 @self.Constraint(
408 self.flowsheet().time, doc="Willans CT b calculation"
409 )
410 def willans_CT_b_calculation(self, t):
411 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
413 # c parameter
414 @self.Constraint(
415 self.flowsheet().time, doc="Willans CT efficiency calculation"
416 )
417 def willans_CT_efficiency_calculation(self, t):
418 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)
420 def calculate_BPST_willans_parameters(self):
421 # a parameter
422 @self.Constraint(
423 self.flowsheet().time, doc="Willans BPST a calculation"
424 )
425 def willans_BPST_a_calculation(self, t):
426 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
429 # b parameter
430 @self.Constraint(
431 self.flowsheet().time, doc="Willans BPST b calculation"
432 )
433 def willans_BPST_b_calculation(self, t):
434 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
437 # c parameter
438 @self.Constraint(
439 self.flowsheet().time, doc="Willans BPST c calculation"
440 )
441 def willans_BPST_efficiency_calculation(self, t):
442 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)
444 def calculate_Tsat_willans_parameters(self):
445 # a parameter
446 @self.Constraint(
447 self.flowsheet().time, doc="Willans Tsat a calculation"
448 )
449 def willans_Tsat_a_calculation(self, t):
450 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)
452 # b parameter
453 @self.Constraint(
454 self.flowsheet().time, doc="Willans Tsat b calculation"
455 )
456 def willans_Tsat_b_calculation(self, t):
457 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
459 # c parameter
460 @self.Constraint(
461 self.flowsheet().time, doc="Willans Tsat efficiency calculation"
462 )
463 def willans_Tsat_c_calculation(self, t):
464 return self.willans_efficiency[t] == 0.83333
466 def calculate_willans_coefficients(self):
467 # Calculate willans coefficients
468 @self.Constraint(
469 self.flowsheet().time, doc="Willans slope calculation"
470 )
471 def willans_slope_calculation(self, t):
472 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])
475 @self.Constraint(
476 self.flowsheet().time, doc="Willans intercept calculation"
477 )
478 def willans_intercept_calculation(self, t):
479 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])
481 def add_mechanical_and_isentropic_work_definition(self):
483 if 'willans' in self.config.calculation_method:
484 self.efficiency_isentropic = Expression(
485 self.flowsheet().time,
486 rule=lambda b, t: (
487 b.delta_h_act[t] / b.delta_h_is[t]
488 )
489 )
490 else:
491 # Mechanical work
492 @self.Constraint(
493 self.flowsheet().time, doc="Isentropic and mechanical work relationship"
494 )
495 def isentropic_and_mechanical_work_eq(self, t):
496 return self.work_mechanical[t] == (
497 self.delta_h_is[t] * self.control_volume.properties_in[t].flow_mol * self.efficiency_isentropic[t]
498 )
500 def add_willans_line_relationship(self):
501 @self.Constraint(
502 self.flowsheet().time, doc="Willans line and mechanical work relationship"
503 )
504 def willans_line_eq(self, t):
505 eps = 0.0001 # smoothing parameter; smaller = closer to exact max, larger = smoother
506 return self.work_mechanical[t] == smooth_min(
507 -(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]),
508 0.0,
509 eps
510 ) * (self.willans_slope[t] * self.willans_max_mol[t])
512 def add_electrical_work_definition(self):
513 # Electrical work
514 @self.Constraint(
515 self.flowsheet().time, doc="Calculate electrical work of turbine"
516 )
517 def electrical_energy_balance(self, t):
518 return self.work_electrical[t] == self.work_mechanical[t] * self.efficiency_motor[t]
520 def add_isentropic_work_definition(self):
521 self.delta_h_is = Expression(
522 self.flowsheet().time,
523 rule=lambda b, t: (
524 b.properties_isentropic[t].enth_mol - b.control_volume.properties_in[t].enth_mol
525 )
526 )
527 self.delta_h_act = Expression(
528 self.flowsheet().time,
529 rule=lambda b, t: (
530 b.control_volume.properties_out[t].enth_mol - b.control_volume.properties_in[t].enth_mol
531 )
532 )
533 # # Isentropic work
534 # @self.Constraint(
535 # self.flowsheet().time, doc="Calculate work of isentropic process"
536 # )
537 # def isentropic_energy_balance(self, t):
538 # 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
540 def initialize_build(
541 blk,
542 state_args=None,
543 routine=None,
544 outlvl=idaeslog.NOTSET,
545 solver=None,
546 optarg=None,
547 ):
548 """
549 General wrapper for pressure changer initialization routines
551 Keyword Arguments:
552 routine : str stating which initialization routine to execute
553 * None - use routine matching thermodynamic_assumption
554 * 'isentropic' - use isentropic initialization routine
555 * 'isothermal' - use isothermal initialization routine
556 state_args : a dict of arguments to be passed to the property
557 package(s) to provide an initial state for
558 initialization (see documentation of the specific
559 property package) (default = {}).
560 outlvl : sets output level of initialization routine
561 optarg : solver options dictionary object (default=None, use
562 default solver options)
563 solver : str indicating which solver to use during
564 initialization (default = None, use default solver)
566 Returns:
567 None
568 """
569 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
570 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
572 # Create solver
573 opt = get_solver(solver, optarg)
575 cv = blk.control_volume
576 t0 = blk.flowsheet().time.first()
577 state_args_out = {}
579 if state_args is None: 579 ↛ 592line 579 didn't jump to line 592 because the condition on line 579 was always true
580 state_args = {}
581 state_dict = cv.properties_in[t0].define_port_members()
583 for k in state_dict.keys():
584 if state_dict[k].is_indexed():
585 state_args[k] = {}
586 for m in state_dict[k].keys():
587 state_args[k][m] = state_dict[k][m].value
588 else:
589 state_args[k] = state_dict[k].value
591 # Get initialisation guesses for outlet and isentropic states
592 for k in state_args:
593 if k == "pressure" and k not in state_args_out:
594 # Work out how to estimate outlet pressure
595 if cv.properties_out[t0].pressure.fixed:
596 # Fixed outlet pressure, use this value
597 state_args_out[k] = value(cv.properties_out[t0].pressure)
598 elif blk.deltaP[t0].fixed: 598 ↛ 599line 598 didn't jump to line 599 because the condition on line 598 was never true
599 state_args_out[k] = value(state_args[k] + blk.deltaP[t0])
600 elif blk.ratioP[t0].fixed: 600 ↛ 601line 600 didn't jump to line 601 because the condition on line 600 was never true
601 state_args_out[k] = value(state_args[k] * blk.ratioP[t0])
602 else:
603 # Not obvious what to do, use inlet state
604 state_args_out[k] = state_args[k]
605 elif k not in state_args_out: 605 ↛ 592line 605 didn't jump to line 592 because the condition on line 605 was always true
606 state_args_out[k] = state_args[k]
608 # Initialize state blocks
609 flags = cv.properties_in.initialize(
610 outlvl=outlvl,
611 optarg=optarg,
612 solver=solver,
613 hold_state=True,
614 state_args=state_args,
615 )
616 cv.properties_out.initialize(
617 outlvl=outlvl,
618 optarg=optarg,
619 solver=solver,
620 hold_state=False,
621 state_args=state_args_out,
622 )
624 init_log.info_high("Initialization Step 1 Complete.")
625 # ---------------------------------------------------------------------
626 # Initialize Isentropic block
628 blk.properties_isentropic.initialize(
629 outlvl=outlvl,
630 optarg=optarg,
631 solver=solver,
632 state_args=state_args_out,
633 )
635 init_log.info_high("Initialization Step 2 Complete.")
637 # Skipping step 3 because Isothermal had problems.
639 # ---------------------------------------------------------------------
640 # Solve unit
641 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
642 res = opt.solve(blk, tee=slc.tee)
643 init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res)))
646 # ---------------------------------------------------------------------
647 # Release Inlet state
648 blk.control_volume.release_state(flags, outlvl)
650 if not check_optimal_termination(res): 650 ↛ 651line 650 didn't jump to line 651 because the condition on line 650 was never true
651 raise InitializationError(
652 f"{blk.name} failed to initialize successfully. Please check "
653 f"the output logs for more information."
654 )
656 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}")
658 def _get_performance_contents(self, time_point=0):
659 var_dict = {}
660 if hasattr(self, "deltaP"): 660 ↛ 662line 660 didn't jump to line 662 because the condition on line 660 was always true
661 var_dict["Mechanical Work"] = self.work_mechanical[time_point]
662 if hasattr(self, "deltaP"): 662 ↛ 664line 662 didn't jump to line 664 because the condition on line 662 was always true
663 var_dict["Electrical Work"] = self.work_electrical[time_point]
664 if hasattr(self, "deltaP"): 664 ↛ 666line 664 didn't jump to line 666 because the condition on line 664 was always true
665 var_dict["Pressure Change"] = self.deltaP[time_point]
666 if hasattr(self, "ratioP"): 666 ↛ 671line 666 didn't jump to line 671 because the condition on line 666 was always true
667 var_dict["Pressure Ratio"] = self.ratioP[time_point]
668 # if hasattr(self, "efficiency_isentropic"):
669 # var_dict["Isentropic Efficiency"] = self.efficiency_isentropic[time_point]
671 return {"vars": var_dict}
673 def calculate_scaling_factors(self):
674 super().calculate_scaling_factors()
676 if hasattr(self, "work_fluid"):
677 for t, v in self.work_fluid.items():
678 iscale.set_scaling_factor(
679 v,
680 iscale.get_scaling_factor(
681 self.control_volume.work[t], default=1, warning=True
682 ),
683 )
685 if hasattr(self, "work_mechanical"):
686 for t, v in self.work_mechanical.items():
687 iscale.set_scaling_factor(
688 v,
689 iscale.get_scaling_factor(
690 self.control_volume.work[t], default=1, warning=True
691 ),
692 )
694 if hasattr(self, "work_isentropic"):
695 for t, v in self.work_isentropic.items():
696 iscale.set_scaling_factor(
697 v,
698 iscale.get_scaling_factor(
699 self.control_volume.work[t], default=1, warning=True
700 ),
701 )
703 if hasattr(self, "ratioP_calculation"):
704 for t, c in self.ratioP_calculation.items():
705 iscale.constraint_scaling_transform(
706 c,
707 iscale.get_scaling_factor(
708 self.control_volume.properties_in[t].pressure,
709 default=1,
710 warning=True,
711 ),
712 overwrite=False,
713 )
715 if hasattr(self, "fluid_work_calculation"):
716 for t, c in self.fluid_work_calculation.items():
717 iscale.constraint_scaling_transform(
718 c,
719 iscale.get_scaling_factor(
720 self.control_volume.deltaP[t], default=1, warning=True
721 ),
722 overwrite=False,
723 )
725 if hasattr(self, "actual_work"):
726 for t, c in self.actual_work.items():
727 iscale.constraint_scaling_transform(
728 c,
729 iscale.get_scaling_factor(
730 self.control_volume.work[t], default=1, warning=True
731 ),
732 overwrite=False,
733 )
735 if hasattr(self, "isentropic_pressure"):
736 for t, c in self.isentropic_pressure.items():
737 iscale.constraint_scaling_transform(
738 c,
739 iscale.get_scaling_factor(
740 self.control_volume.properties_in[t].pressure,
741 default=1,
742 warning=True,
743 ),
744 overwrite=False,
745 )
747 if hasattr(self, "isentropic"):
748 for t, c in self.isentropic.items():
749 iscale.constraint_scaling_transform(
750 c,
751 iscale.get_scaling_factor(
752 self.control_volume.properties_in[t].entr_mol,
753 default=1,
754 warning=True,
755 ),
756 overwrite=False,
757 )
759 if hasattr(self, "isentropic_energy_balance"):
760 for t, c in self.isentropic_energy_balance.items():
761 iscale.constraint_scaling_transform(
762 c,
763 iscale.get_scaling_factor(
764 self.control_volume.work[t], default=1, warning=True
765 ),
766 overwrite=False,
767 )
769 if hasattr(self, "zero_work_equation"):
770 for t, c in self.zero_work_equation.items():
771 iscale.constraint_scaling_transform(
772 c,
773 iscale.get_scaling_factor(
774 self.control_volume.work[t], default=1, warning=True
775 ),
776 )
778 if hasattr(self, "state_material_balances"):
779 cvol = self.control_volume
780 phase_list = cvol.properties_in.phase_list
781 phase_component_set = cvol.properties_in.phase_component_set
782 mb_type = cvol._constructed_material_balance_type
783 if mb_type == MaterialBalanceType.componentPhase:
784 for (t, p, j), c in self.state_material_balances.items():
785 sf = iscale.get_scaling_factor(
786 cvol.properties_in[t].get_material_flow_terms(p, j),
787 default=1,
788 warning=True,
789 )
790 iscale.constraint_scaling_transform(c, sf)
791 elif mb_type == MaterialBalanceType.componentTotal:
792 for (t, j), c in self.state_material_balances.items():
793 sf = iscale.min_scaling_factor(
794 [
795 cvol.properties_in[t].get_material_flow_terms(p, j)
796 for p in phase_list
797 if (p, j) in phase_component_set
798 ]
799 )
800 iscale.constraint_scaling_transform(c, sf)
801 else:
802 # There are some other material balance types but they create
803 # constraints with different names.
804 _log.warning(f"Unknown material balance type {mb_type}")