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

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 

18 

19# Changing existing config block attributes 

20# pylint: disable=protected-access 

21 

22# Import Python libraries 

23 

24 

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 

38 

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 

60 

61 

62__author__ = "Emmanuel Ogbe, Andrew Lee" 

63_log = idaeslog.getLogger(__name__) 

64 

65 

66@declare_process_block_class("TurbineBase") 

67class TurbineBaseData(UnitModelBlockData): 

68 """ 

69 Standard Compressor/Expander Unit Model Class 

70 """ 

71 

72 CONFIG = UnitModelBlockData.CONFIG() 

73 

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 ) 

182 

183 def build(self): 

184 """ 

185 

186 Args: 

187 None 

188 

189 Returns: 

190 None 

191 """ 

192 # Call UnitModel.build 

193 super().build() 

194 

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 ) 

202 

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() 

206 

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 ) 

211 

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 ) 

219 

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 ) 

224 

225 # add momentum balance 

226 self.control_volume.add_momentum_balances( 

227 balance_type=self.config.momentum_balance_type, has_pressure_change=True 

228 ) 

229 

230 # Add Ports 

231 self.add_inlet_port() 

232 self.add_outlet_port() 

233 

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[:]) 

237 

238 self.work_mechanical = Reference(self.control_volume.work[:]) 

239 # self.work_mechanical.fix(-1000e3) 

240 # self.work_mechanical.unfix() 

241 

242 

243 # Add Momentum balance variable 'deltaP' 

244 self.deltaP = Reference(self.control_volume.deltaP[:]) 

245 

246 # Performance Variables 

247 self.ratioP = Var(self.flowsheet().time, initialize=1.0, doc="Pressure Ratio") 

248 

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 ) 

256 

257 units_meta = self.control_volume.config.property_package.get_metadata() 

258 

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 ) 

266 

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 # ) 

285 

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 ) 

292 

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 ) 

299 

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 ) 

308 

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 ) 

315 

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 ) 

322 

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 ) 

329 

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 ) 

336 

337 self.willans_efficiency = Var( 

338 self.flowsheet().time, 

339 initialize=1.0, 

340 doc="Willans efficiency", 

341 ) 

342 

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 

347 

348 self.properties_isentropic = self.config.property_package.build_state_block( 

349 self.flowsheet().time, doc="isentropic properties at outlet", **tmp_dict 

350 ) 

351 

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 ) 

361 

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 ) 

368 

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 ) 

376 

377 self.add_isentropic_work_definition() 

378 

379 if 'willans' in self.config.calculation_method: 

380 # Write isentropic efficiency eqn 

381 self.add_willans_line_relationship() 

382 

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() 

386 

387 elif self.config.calculation_method == "BPST_willans": 

388 self.calculate_BPST_willans_parameters() 

389 

390 elif self.config.calculation_method == "CT_willans": 

391 self.calculate_CT_willans_parameters() 

392 

393 # calculate slope and intercept 

394 self.calculate_willans_coefficients() 

395 

396 self.add_mechanical_and_isentropic_work_definition() 

397 self.add_electrical_work_definition() 

398 

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 

406 

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 

413 

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) 

420 

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 

428 

429 

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 

436 

437 

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) 

444 

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) 

452 

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 

459 

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 

466 

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]) 

474 

475 

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]) 

481 

482 def add_mechanical_and_isentropic_work_definition(self): 

483 

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 ) 

500 

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]) 

512 

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] 

520 

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 

540 

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 

551 

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) 

566 

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") 

572 

573 # Create solver 

574 opt = get_solver(solver, optarg) 

575 

576 cv = blk.control_volume 

577 t0 = blk.flowsheet().time.first() 

578 state_args_out = {} 

579 

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) 

582 

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] 

599 

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 ) 

615 

616 init_log.info_high("Initialization Step 1 Complete.") 

617 # --------------------------------------------------------------------- 

618 # Initialize Isentropic block 

619 

620 blk.properties_isentropic.initialize( 

621 outlvl=outlvl, 

622 optarg=optarg, 

623 solver=solver, 

624 state_args=state_args_out, 

625 ) 

626 

627 init_log.info_high("Initialization Step 2 Complete.") 

628 

629 # Skipping step 3 because Isothermal had problems. 

630 

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))) 

636 

637 

638 # --------------------------------------------------------------------- 

639 # Release Inlet state 

640 blk.control_volume.release_state(flags, outlvl) 

641 

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 ) 

647 

648 init_log.info(f"Initialization Complete: {idaeslog.condition(res)}") 

649 

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] 

662 

663 return {"vars": var_dict} 

664 

665 def calculate_scaling_factors(self): 

666 super().calculate_scaling_factors() 

667 

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 ) 

676 

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 ) 

685 

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 ) 

694 

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 ) 

706 

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 ) 

716 

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 ) 

726 

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 ) 

738 

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 ) 

750 

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 ) 

760 

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 ) 

769 

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}")