Coverage for backend/idaes_service/solver/custom/hda_ideal

1##############################################################################

2# Institute for the Design of Advanced Energy Systems Process Systems

4# software owners: The Regents of the University of California, through

5# Lawrence Berkeley National Laboratory, National Technology & Engineering

6# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia

9# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and

10# license information, respectively. Both files are also available online

11# at the URL "https://github.com/IDAES/idaes-pse".

12##############################################################################

13"""

14Example ideal parameter block for the VLE calucations for a

15Benzene-Toluene-o-Xylene system.

16"""

19# Import Pyomo libraries

20from pyomo.environ import (

21 Constraint,

22 Expression,

23 log,

24 NonNegativeReals,

25 Var,

26 Set,

27 Param,

28 sqrt,

29 log10,

30 units as pyunits,

31)

32from pyomo.util.calc_var_value import calculate_variable_from_constraint

33from pyomo.common.config import ConfigValue

35# Import IDAES cores

36from idaes.core import (

37 declare_process_block_class,

38 MaterialFlowBasis,

39 PhysicalParameterBlock,

40 StateBlockData,

41 StateBlock,

42 MaterialBalanceType,

43 EnergyBalanceType,

44 Component,

45 LiquidPhase,

46 VaporPhase,

47)

48from idaes.core.util.constants import Constants as const

49from idaes.core.util.initialization import fix_state_vars, solve_indexed_blocks

50from idaes.core.initialization import InitializerBase

51from idaes.core.util.misc import add_object_reference

52from idaes.core.util.model_statistics import number_unfixed_variables

53from idaes.core.util.misc import extract_data

54from idaes.core.solvers import get_solver

55import idaes.core.util.scaling as iscale

56import idaes.logger as idaeslog

58# Set up logger

59_log = idaeslog.getLogger(__name__)

62class HDAInitializer(InitializerBase):

63 """

64 Initializer for HDA Property package.

66 """

68 CONFIG = InitializerBase.CONFIG()

69 CONFIG.declare(

70 "solver",

71 ConfigValue(default=None, domain=str, description="Initialization solver"),

72 )

73 CONFIG.declare(

74 "solver_options",

75 ConfigValue(default=None, description="Initialization solver options"),

76 )

78 def initialization_routine(self, blk):

79 init_log = idaeslog.getInitLogger(

80 blk.name, self.config.output_level, tag="properties"

81 )

82 solve_log = idaeslog.getSolveLogger(

83 blk.name, self.config.output_level, tag="properties"

84 )

86 # Set solver

87 solver = get_solver(self.config.solver, self.config.solver_options)

89 # ---------------------------------------------------------------------

90 # If present, initialize bubble and dew point calculations

91 for k in blk.keys():

92 if hasattr(blk[k], "eq_temperature_dew"):

93 calculate_variable_from_constraint(

94 blk[k].temperature_dew, blk[k].eq_temperature_dew

95 )

97 if hasattr(blk[k], "eq_pressure_dew"):

98 calculate_variable_from_constraint(

99 blk[k].pressure_dew, blk[k].eq_pressure_dew

100 )

101

102 init_log.info_high(

103 "Initialization Step 1 - Dew and bubble points " "calculation completed."

104 )

105

106 # ---------------------------------------------------------------------

107 # If flash, initialize T1 and Teq

108 for k in blk.keys():

109 if blk[k].config.has_phase_equilibrium and not blk[k].config.defined_state:

110 blk[k]._t1.value = max(

111 blk[k].temperature.value, blk[k].temperature_bubble.value

112 )

113 blk[k]._teq.value = min(blk[k]._t1.value, blk[k].temperature_dew.value)

114

115 init_log.info_high(

116 "Initialization Step 2 - Equilibrium temperature " " calculation completed."

117 )

118

119 # ---------------------------------------------------------------------

120 # Initialize flow rates and compositions

121 free_vars = 0

122 for k in blk.keys():

123 free_vars += number_unfixed_variables(blk[k])

124 if free_vars > 0:

125 try:

126 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:

127 res = solve_indexed_blocks(solver, [blk], tee=slc.tee)

128 except:

129 res = None

130 else:

131 res = None

132

133 init_log.info("Initialization Complete")

134

135 return res

136

137

138@declare_process_block_class("HDAParameterBlock")

139class HDAParameterData(PhysicalParameterBlock):

140 CONFIG = PhysicalParameterBlock.CONFIG()

141

142 def build(self):

143 """

144 Callable method for Block construction.

145 """

146 super(HDAParameterData, self).build()

147

148 self._state_block_class = IdealStateBlock

149

150 self.benzene = Component()

151 self.toluene = Component()

152 self.methane = Component()

153 self.hydrogen = Component()

154

155 self.Liq = LiquidPhase()

156 self.Vap = VaporPhase()

157

158 # List of components in each phase (optional)

159 self.phase_comp = {"Liq": self.component_list, "Vap": self.component_list}

160

161 # List of phase equilibrium index

162 self.phase_equilibrium_idx = Set(initialize=[1, 2, 3, 4])

163

164 self.phase_equilibrium_list = {

165 1: ["benzene", ("Vap", "Liq")],

166 2: ["toluene", ("Vap", "Liq")],

167 3: ["hydrogen", ("Vap", "Liq")],

168 4: ["methane", ("Vap", "Liq")],

169 }

170

171 # Thermodynamic reference state

172 self.pressure_ref = Param(

173 mutable=True, default=101325, units=pyunits.Pa, doc="Reference pressure"

174 )

175 self.temperature_ref = Param(

176 mutable=True, default=298.15, units=pyunits.K, doc="Reference temperature"

177 )

178

179 # Source: The Properties of Gases and Liquids (1987)

180 # 4th edition, Chemical Engineering Series - Robert C. Reid

181 pressure_crit_data = {

182 "benzene": 48.9e5,

183 "toluene": 41e5,

184 "hydrogen": 12.9e5,

185 "methane": 46e5,

186 }

187

188 self.pressure_crit = Param(

189 self.component_list,

190 within=NonNegativeReals,

191 mutable=True,

192 units=pyunits.Pa,

193 initialize=extract_data(pressure_crit_data),

194 doc="Critical pressure",

195 )

196

197 # Source: The Properties of Gases and Liquids (1987)

198 # 4th edition, Chemical Engineering Series - Robert C. Reid

199 temperature_crit_data = {

200 "benzene": 562.2,

201 "toluene": 591.8,

202 "hydrogen": 33.0,

203 "methane": 190.4,

204 }

205

206 self.temperature_crit = Param(

207 self.component_list,

208 within=NonNegativeReals,

209 mutable=True,

210 units=pyunits.K,

211 initialize=extract_data(temperature_crit_data),

212 doc="Critical temperature",

213 )

214

215 # Source: The Properties of Gases and Liquids (1987)

216 # 4th edition, Chemical Engineering Series - Robert C. Reid

217 mw_comp_data = {

218 "benzene": 78.1136e-3,

219 "toluene": 92.1405e-3,

220 "hydrogen": 2.016e-3,

221 "methane": 16.043e-3,

222 }

223

224 self.mw_comp = Param(

225 self.component_list,

226 mutable=True,

227 units=pyunits.kg / pyunits.mol,

228 initialize=extract_data(mw_comp_data),

229 doc="molecular weight",

230 )

231

232 # Constants for liquid densities

233 # Source: Perry's Chemical Engineers Handbook

234 # - Robert H. Perry (Cp_liq)

235 dens_liq_data = {

236 ("benzene", "1"): 1.0162,

237 ("benzene", "2"): 0.2655,

238 ("benzene", "3"): 562.16,

239 ("benzene", "4"): 0.28212,

240 ("toluene", "1"): 0.8488,

241 ("toluene", "2"): 0.26655,

242 ("toluene", "3"): 591.8,

243 ("toluene", "4"): 0.2878,

244 ("hydrogen", "1"): 5.414,

245 ("hydrogen", "2"): 0.34893,

246 ("hydrogen", "3"): 33.19,

247 ("hydrogen", "4"): 0.2706,

248 ("methane", "1"): 2.9214,

249 ("methane", "2"): 0.28976,

250 ("methane", "3"): 190.56,

251 ("methane", "4"): 0.28881,

252 }

253

254 self.dens_liq_param_1 = Param(

255 self.component_list,

256 mutable=True,

257 initialize={c: v for (c, j), v in dens_liq_data.items() if j == "1"},

258 doc="Parameter 1 to compute liquid densities",

259 units=pyunits.kmol * pyunits.m**-3,

260 )

261

262 self.dens_liq_param_2 = Param(

263 self.component_list,

264 mutable=True,

265 initialize={c: v for (c, j), v in dens_liq_data.items() if j == "2"},

266 doc="Parameter 2 to compute liquid densities",

267 units=pyunits.dimensionless,

268 )

269

270 self.dens_liq_param_3 = Param(

271 self.component_list,

272 mutable=True,

273 initialize={c: v for (c, j), v in dens_liq_data.items() if j == "3"},

274 doc="Parameter 3 to compute liquid densities",

275 units=pyunits.K,

276 )

277

278 self.dens_liq_param_4 = Param(

279 self.component_list,

280 mutable=True,

281 initialize={c: v for (c, j), v in dens_liq_data.items() if j == "4"},

282 doc="Parameter 4 to compute liquid densities",

283 units=pyunits.dimensionless,

284 )

285

286 # Boiling point at standard pressure

287 # Source: Perry's Chemical Engineers Handbook

288 # - Robert H. Perry (Cp_liq)

289 bp_data = {

290 ("benzene"): 353.25,

291 ("toluene"): 383.95,

292 ("hydrogen"): 20.45,

293 ("methane"): 111.75,

294 }

295

296 self.temperature_boil = Param(

297 self.component_list,

298 mutable=True,

299 units=pyunits.K,

300 initialize=extract_data(bp_data),

301 doc="Pure component boiling points at standard pressure",

302 )

303

304 # Constants for specific heat capacity, enthalpy

305 # Sources: The Properties of Gases and Liquids (1987)

306 # 4th edition, Chemical Engineering Series - Robert C. Reid

307 # Perry's Chemical Engineers Handbook

308 # - Robert H. Perry (Cp_liq)

309 cp_ig_data = {

310 ("Liq", "benzene", "1"): 1.29e5,

311 ("Liq", "benzene", "2"): -1.7e2,

312 ("Liq", "benzene", "3"): 6.48e-1,

313 ("Liq", "benzene", "4"): 0,

314 ("Liq", "benzene", "5"): 0,

315 ("Vap", "benzene", "1"): -3.392e1,

316 ("Vap", "benzene", "2"): 4.739e-1,

317 ("Vap", "benzene", "3"): -3.017e-4,

318 ("Vap", "benzene", "4"): 7.130e-8,

319 ("Vap", "benzene", "5"): 0,

320 ("Liq", "toluene", "1"): 1.40e5,

321 ("Liq", "toluene", "2"): -1.52e2,

322 ("Liq", "toluene", "3"): 6.95e-1,

323 ("Liq", "toluene", "4"): 0,

324 ("Liq", "toluene", "5"): 0,

325 ("Vap", "toluene", "1"): -2.435e1,

326 ("Vap", "toluene", "2"): 5.125e-1,

327 ("Vap", "toluene", "3"): -2.765e-4,

328 ("Vap", "toluene", "4"): 4.911e-8,

329 ("Vap", "toluene", "5"): 0,

330 ("Liq", "hydrogen", "1"): 0, # 6.6653e1,

331 ("Liq", "hydrogen", "2"): 0, # 6.7659e3,

332 ("Liq", "hydrogen", "3"): 0, # -1.2363e2,

333 ("Liq", "hydrogen", "4"): 0, # 4.7827e2, # Eqn 2

334 ("Liq", "hydrogen", "5"): 0,

335 ("Vap", "hydrogen", "1"): 2.714e1,

336 ("Vap", "hydrogen", "2"): 9.274e-3,

337 ("Vap", "hydrogen", "3"): -1.381e-5,

338 ("Vap", "hydrogen", "4"): 7.645e-9,

339 ("Vap", "hydrogen", "5"): 0,

340 ("Liq", "methane", "1"): 0, # 6.5708e1,

341 ("Liq", "methane", "2"): 0, # 3.8883e4,

342 ("Liq", "methane", "3"): 0, # -2.5795e2,

343 ("Liq", "methane", "4"): 0, # 6.1407e2, # Eqn 2

344 ("Liq", "methane", "5"): 0,

345 ("Vap", "methane", "1"): 1.925e1,

346 ("Vap", "methane", "2"): 5.213e-2,

347 ("Vap", "methane", "3"): 1.197e-5,

348 ("Vap", "methane", "4"): -1.132e-8,

349 ("Vap", "methane", "5"): 0,

350 }

351

352 self.cp_ig_1 = Param(

353 self.phase_list,

354 self.component_list,

355 mutable=True,

356 initialize={(p, c): v for (p, c, j), v in cp_ig_data.items() if j == "1"},

357 doc="Parameter 1 to compute Cp_comp",

358 units=pyunits.J / pyunits.mol / pyunits.K,

359 )

360

361 self.cp_ig_2 = Param(

362 self.phase_list,

363 self.component_list,

364 mutable=True,

365 initialize={(p, c): v for (p, c, j), v in cp_ig_data.items() if j == "2"},

366 doc="Parameter 2 to compute Cp_comp",

367 units=pyunits.J / pyunits.mol / pyunits.K**2,

368 )

369

370 self.cp_ig_3 = Param(

371 self.phase_list,

372 self.component_list,

373 mutable=True,

374 initialize={(p, c): v for (p, c, j), v in cp_ig_data.items() if j == "3"},

375 doc="Parameter 3 to compute Cp_comp",

376 units=pyunits.J / pyunits.mol / pyunits.K**3,

377 )

378

379 self.cp_ig_4 = Param(

380 self.phase_list,

381 self.component_list,

382 mutable=True,

383 initialize={(p, c): v for (p, c, j), v in cp_ig_data.items() if j == "4"},

384 doc="Parameter 4 to compute Cp_comp",

385 units=pyunits.J / pyunits.mol / pyunits.K**4,

386 )

387

388 self.cp_ig_5 = Param(

389 self.phase_list,

390 self.component_list,

391 mutable=True,

392 initialize={(p, c): v for (p, c, j), v in cp_ig_data.items() if j == "5"},

393 doc="Parameter 5 to compute Cp_comp",

394 units=pyunits.J / pyunits.mol / pyunits.K**5,

395 )

396

397 # Source: The Properties of Gases and Liquids (1987)

398 # 4th edition, Chemical Engineering Series - Robert C. Reid

399 # fitted to Antoine form

400 # H2, Methane from NIST webbook

401 pressure_sat_coeff_data = {

402 ("benzene", "A"): 4.202,

403 ("benzene", "B"): 1322,

404 ("benzene", "C"): -38.56,

405 ("toluene", "A"): 4.216,

406 ("toluene", "B"): 1435,

407 ("toluene", "C"): -43.33,

408 ("hydrogen", "A"): 3.543,

409 ("hydrogen", "B"): 99.40,

410 ("hydrogen", "C"): 7.726,

411 ("methane", "A"): 3.990,

412 ("methane", "B"): 443.0,

413 ("methane", "C"): -0.49,

414 }

415

416 self.pressure_sat_coeff_A = Param(

417 self.component_list,

418 mutable=True,

419 initialize={

420 c: v for (c, j), v in pressure_sat_coeff_data.items() if j == "A"

421 },

422 doc="Parameter A to compute saturated pressure",

423 units=pyunits.dimensionless,

424 )

425

426 self.pressure_sat_coeff_B = Param(

427 self.component_list,

428 mutable=True,

429 initialize={

430 c: v for (c, j), v in pressure_sat_coeff_data.items() if j == "B"

431 },

432 doc="Parameter B to compute saturated pressure",

433 units=pyunits.K,

434 )

435

436 self.pressure_sat_coeff_C = Param(

437 self.component_list,

438 mutable=True,

439 initialize={

440 c: v for (c, j), v in pressure_sat_coeff_data.items() if j == "C"

441 },

442 doc="Parameter C to compute saturated pressure",

443 units=pyunits.K,

444 )

445

446 # Source: The Properties of Gases and Liquids (1987)

447 # 4th edition, Chemical Engineering Series - Robert C. Reid

448 dh_vap = {"benzene": 3.387e4, "toluene": 3.8262e4, "hydrogen": 0, "methane": 0}

449

450 self.dh_vap = Param(

451 self.component_list,

452 mutable=True,

453 units=pyunits.J / pyunits.mol,

454 initialize=extract_data(dh_vap),

455 doc="heat of vaporization",

456 )

457

458 # Set default scaling factors

459 self.set_default_scaling("flow_mol", 1e3)

460 self.set_default_scaling("flow_mol_phase_comp", 1e3)

461 self.set_default_scaling("flow_mol_phase", 1e3)

462 self.set_default_scaling("material_flow_terms", 1e3)

463 self.set_default_scaling("enthalpy_flow_terms", 1e-2)

464 self.set_default_scaling("mole_frac_comp", 1e1)

465 self.set_default_scaling("temperature", 1e-2)

466 self.set_default_scaling("temperature_dew", 1e-2)

467 self.set_default_scaling("temperature_bubble", 1e-2)

468 self.set_default_scaling("pressure", 1e-5)

469 self.set_default_scaling("pressure_sat", 1e-5)

470 self.set_default_scaling("pressure_dew", 1e-5)

471 self.set_default_scaling("pressure_bubble", 1e-5)

472 self.set_default_scaling("mole_frac_phase_comp", 1e1)

473 self.set_default_scaling("enth_mol_phase", 1e-3, index="Liq")

474 self.set_default_scaling("enth_mol_phase", 1e-4, index="Vap")

475 self.set_default_scaling("enth_mol", 1e-3)

476 self.set_default_scaling("entr_mol_phase", 1e-2)

477 self.set_default_scaling("entr_mol", 1e-2)

478

479 @classmethod

480 def define_metadata(cls, obj):

481 """Define properties supported and units."""

482 obj.add_properties(

483 {

484 "flow_mol": {"method": None},

485 "flow_mol_phase_comp": {"method": None},

486 "mole_frac_comp": {"method": None},

487 "temperature": {"method": None},

488 "pressure": {"method": None},

489 "flow_mol_phase": {"method": None},

490 "dens_mol_phase": {"method": "_dens_mol_phase"},

491 "pressure_sat": {"method": "_pressure_sat"},

492 "mole_frac_phase_comp": {"method": "_mole_frac_phase"},

493 "energy_internal_mol_phase_comp": {

494 "method": "_energy_internal_mol_phase_comp"

495 },

496 "energy_internal_mol_phase": {"method": "_energy_internal_mol_phase"},

497 "enth_mol_phase_comp": {"method": "_enth_mol_phase_comp"},

498 "enth_mol_phase": {"method": "_enth_mol_phase"},

499 "entr_mol_phase_comp": {"method": "_entr_mol_phase_comp"},

500 "entr_mol_phase": {"method": "_entr_mol_phase"},

501 "temperature_bubble": {"method": "_temperature_bubble"},

502 "temperature_dew": {"method": "_temperature_dew"},

503 "pressure_bubble": {"method": "_pressure_bubble"},

504 "pressure_dew": {"method": "_pressure_dew"},

505 "fug_phase_comp": {"method": "_fug_phase_comp"},

506 }

507 )

508

509 obj.define_custom_properties(

510 {

511 # Enthalpy of vaporization

512 "dh_vap": {"method": "_dh_vap", "units": obj.derived_units.ENERGY_MOLE},

513 # Entropy of vaporization

514 "ds_vap": {

515 "method": "_ds_vap",

516 "units": obj.derived_units.ENTROPY_MOLE,

517 },

518 }

519 )

520

521 obj.add_default_units(

522 {

523 "time": pyunits.s,

524 "length": pyunits.m,

525 "mass": pyunits.kg,

526 "amount": pyunits.mol,

527 "temperature": pyunits.K,

528 }

529 )

530

531

532class _IdealStateBlock(StateBlock):

533 """

534 This Class contains methods which should be applied to Property Blocks as a

535 whole, rather than individual elements of indexed Property Blocks.

536 """

537

538 default_initializer = HDAInitializer

539

540 def fix_initialization_states(blk):

541 """

542 Fixes state variables for state blocks.

543

544 Returns:

545 None

546 """

547

548 # Fix state variables

549 fix_state_vars(blk)

550

551 # Also need to deactivate sum of mole fraction constraint

552 for k in blk.values():

553 if not k.config.defined_state:

554 k.equilibrium_constraint.deactivate()

555

556

557@declare_process_block_class("IdealStateBlock", block_class=_IdealStateBlock)

558class IdealStateBlockData(StateBlockData):

559 """An example property package for ideal VLE."""

560

561 def build(self):

562 """Callable method for Block construction."""

563 super().build()

564

565 # Add state variables

566 self.flow_mol_phase_comp = Var(

567 self._params.phase_list,

568 self._params.component_list,

569 initialize=0.5,

570 units=pyunits.mol / pyunits.s,

571 bounds=(1e-12, 100),

572 doc="Phase-component molar flow rates",

573 )

574

575 self.pressure = Var(

576 initialize=101325,

577 bounds=(100000, 1000000),

578 units=pyunits.Pa,

579 domain=NonNegativeReals,

580 doc="State pressure",

581 )

582 self.temperature = Var(

583 initialize=298.15,

584 units=pyunits.K,

585 bounds=(298, 1000),

586 domain=NonNegativeReals,

587 doc="State temperature",

588 )

589

590 # Add supporting variables

591 def flow_mol_phase(b, p):

592 return sum(b.flow_mol_phase_comp[p, j] for j in b._params.component_list)

593

594 self.flow_mol_phase = Expression(

595 self._params.phase_list, rule=flow_mol_phase, doc="Phase molar flow rates"

596 )

597

598 def flow_mol(b):

599 return sum(

600 b.flow_mol_phase_comp[p, j]

601 for j in b._params.component_list

602 for p in b._params.phase_list

603 )

604

605 self.flow_mol = Expression(rule=flow_mol, doc="Total molar flowrate")

606

607 def mole_frac_phase_comp(b, p, j):

608 return b.flow_mol_phase_comp[p, j] / b.flow_mol_phase[p]

609

610 self.mole_frac_phase_comp = Expression(

611 self._params.phase_list,

612 self._params.component_list,

613 rule=mole_frac_phase_comp,

614 doc="Phase mole fractions",

615 )

616

617 def mole_frac_comp(b, j):

618 return (

619 sum(b.flow_mol_phase_comp[p, j] for p in b._params.phase_list)

620 / b.flow_mol

621 )

622

623 self.mole_frac_comp = Expression(

624 self._params.component_list,

625 rule=mole_frac_comp,

626 doc="Mixture mole fractions",

627 )

628

629 # Reaction Stoichiometry

630 add_object_reference(

631 self, "phase_equilibrium_list_ref", self._params.phase_equilibrium_list

632 )

633

634 if self.config.has_phase_equilibrium and self.config.defined_state is False:

635 # Definition of equilibrium temperature for smooth VLE

636 self._teq = Var(

637 initialize=self.temperature.value,

638 units=pyunits.K,

639 doc="Temperature for calculating phase equilibrium",

640 )

641 self._t1 = Var(

642 initialize=self.temperature.value,

643 units=pyunits.K,

644 doc="Intermediate temperature for calculating Teq",

645 )

646

647 self.eps_1 = Param(

648 default=0.01,

649 units=pyunits.K,

650 mutable=True,

651 doc="Smoothing parameter for Teq",

652 )

653 self.eps_2 = Param(

654 default=0.0005,

655 units=pyunits.K,

656 mutable=True,

657 doc="Smoothing parameter for Teq",

658 )

659

660 # PSE paper Eqn 13

661 def rule_t1(b):

662 return b._t1 == 0.5 * (

663 b.temperature

664 + b.temperature_bubble

665 + sqrt((b.temperature - b.temperature_bubble) ** 2 + b.eps_1**2)

666 )

667

668 self._t1_constraint = Constraint(rule=rule_t1)

669

670 # PSE paper Eqn 14

671 # TODO : Add option for supercritical extension

672 def rule_teq(b):

673 return b._teq == 0.5 * (

674 b._t1

675 + b.temperature_dew

676 - sqrt((b._t1 - b.temperature_dew) ** 2 + b.eps_2**2)

677 )

678

679 self._teq_constraint = Constraint(rule=rule_teq)

680

681 def rule_tr_eq(b, i):

682 return b._teq / b._params.temperature_crit[i]

683

684 self._tr_eq = Expression(

685 self._params.component_list,

686 rule=rule_tr_eq,

687 doc="Component reduced temperatures",

688 )

689

690 def rule_equilibrium(b, i):

691 return b.fug_phase_comp["Liq", i] == b.fug_phase_comp["Vap", i]

692

693 self.equilibrium_constraint = Constraint(

694 self._params.component_list, rule=rule_equilibrium

695 )

696

697 # -----------------------------------------------------------------------------

698 # Property Methods

699 def _dens_mol_phase(self):

700 self.dens_mol_phase = Var(

701 self._params.phase_list,

702 initialize=1.0,

703 units=pyunits.mol * pyunits.m**-3,

704 doc="Molar density",

705 )

706

707 def rule_dens_mol_phase(b, p):

708 if p == "Vap":

709 return b._dens_mol_vap()

710 else:

711 return b._dens_mol_liq()

712

713 self.eq_dens_mol_phase = Constraint(

714 self._params.phase_list, rule=rule_dens_mol_phase

715 )

716

717 def _energy_internal_mol_phase_comp(self):

718 self.energy_internal_mol_phase_comp = Var(

719 self._params.phase_list,

720 self._params.component_list,

721 units=pyunits.J / pyunits.mol,

722 doc="Phase-component molar specific internal energies",

723 )

724

725 def rule_energy_internal_mol_phase_comp(b, p, j):

726 if p == "Vap":

727 return b.energy_internal_mol_phase_comp[p, j] == b.enth_mol_phase_comp[

728 p, j

729 ] - const.gas_constant * (b.temperature - b._params.temperature_ref)

730 else:

731 return (

732 b.energy_internal_mol_phase_comp[p, j]

733 == b.enth_mol_phase_comp[p, j]

734 )

735

736 self.eq_energy_internal_mol_phase_comp = Constraint(

737 self._params.phase_list,

738 self._params.component_list,

739 rule=rule_energy_internal_mol_phase_comp,

740 )

741

742 def _energy_internal_mol_phase(self):

743 self.energy_internal_mol_phase = Var(

744 self._params.phase_list,

745 units=pyunits.J / pyunits.mol,

746 doc="Phase molar specific internal energies",

747 )

748

749 def rule_energy_internal_mol_phase(b, p):

750 return b.energy_internal_mol_phase[p] == sum(

751 b.energy_internal_mol_phase_comp[p, i] * b.mole_frac_phase_comp[p, i]

752 for i in b._params.component_list

753 )

754

755 self.eq_energy_internal_mol_phase = Constraint(

756 self._params.phase_list, rule=rule_energy_internal_mol_phase

757 )

758

759 def _enth_mol_phase_comp(self):

760 self.enth_mol_phase_comp = Var(

761 self._params.phase_list,

762 self._params.component_list,

763 initialize=7e5,

764 units=pyunits.J / pyunits.mol,

765 doc="Phase-component molar specific enthalpies",

766 )

767

768 def rule_enth_mol_phase_comp(b, p, j):

769 if p == "Vap":

770 return b._enth_mol_comp_vap(j)

771 else:

772 return b._enth_mol_comp_liq(j)

773

774 self.eq_enth_mol_phase_comp = Constraint(

775 self._params.phase_list,

776 self._params.component_list,

777 rule=rule_enth_mol_phase_comp,

778 )

779

780 def _enth_mol_phase(self):

781 self.enth_mol_phase = Var(

782 self._params.phase_list,

783 initialize=7e5,

784 units=pyunits.J / pyunits.mol,

785 doc="Phase molar specific enthalpies",

786 )

787

788 def rule_enth_mol_phase(b, p):

789 return b.enth_mol_phase[p] == sum(

790 b.enth_mol_phase_comp[p, i] * b.mole_frac_phase_comp[p, i]

791 for i in b._params.component_list

792 )

793

794 self.eq_enth_mol_phase = Constraint(

795 self._params.phase_list, rule=rule_enth_mol_phase

796 )

797

798 def _entr_mol_phase_comp(self):

799 self.entr_mol_phase_comp = Var(

800 self._params.phase_list,

801 self._params.component_list,

802 units=pyunits.J / pyunits.mol / pyunits.K,

803 doc="Phase-component molar specific entropies",

804 )

805

806 def rule_entr_mol_phase_comp(b, p, j):

807 if p == "Vap":

808 return b._entr_mol_comp_vap(j)

809 else:

810 return b._entr_mol_comp_liq(j)

811

812 self.eq_entr_mol_phase_comp = Constraint(

813 self._params.phase_list,

814 self._params.component_list,

815 rule=rule_entr_mol_phase_comp,

816 )

817

818 def _entr_mol_phase(self):

819 self.entr_mol_phase = Var(

820 self._params.phase_list,

821 units=pyunits.J / pyunits.mol / pyunits.K,

822 doc="Phase molar specific enthropies",

823 )

824

825 def rule_entr_mol_phase(b, p):

826 return b.entr_mol_phase[p] == sum(

827 b.entr_mol_phase_comp[p, i] * b.mole_frac_phase_comp[p, i]

828 for i in b._params.component_list

829 )

830

831 self.eq_entr_mol_phase = Constraint(

832 self._params.phase_list, rule=rule_entr_mol_phase

833 )

834

835 # -----------------------------------------------------------------------------

836 # General Methods

837 def get_material_flow_terms(self, p, j):

838 """Create material flow terms for control volume."""

839 if not self.is_property_constructed("material_flow_terms"):

840 try:

841

842 def rule_material_flow_terms(blk, p, j):

843 return blk.flow_mol_phase_comp[p, j]

844

845 self.material_flow_terms = Expression(

846 self.params.phase_list,

847 self.params.component_list,

848 rule=rule_material_flow_terms,

849 )

850 except AttributeError:

851 self.del_component(self.material_flow_terms)

852

853 if j in self.params.component_list:

854 return self.material_flow_terms[p, j]

855 else:

856 return 0

857

858 def get_enthalpy_flow_terms(self, p):

859 """Create enthalpy flow terms."""

860 if not self.is_property_constructed("enthalpy_flow_terms"):

861 try:

862

863 def rule_enthalpy_flow_terms(blk, p):

864 return blk.flow_mol_phase[p] * blk.enth_mol_phase[p]

865

866 self.enthalpy_flow_terms = Expression(

867 self.params.phase_list, rule=rule_enthalpy_flow_terms

868 )

869 except AttributeError:

870 self.del_component(self.enthalpy_flow_terms)

871 return self.enthalpy_flow_terms[p]

872

873 def get_material_density_terms(self, p, j):

874 """Create material density terms."""

875 if not self.is_property_constructed("material_density_terms"):

876 try:

877

878 def rule_material_density_terms(b, p, j):

879 return self.dens_mol_phase[p] * self.mole_frac_phase_comp[p, j]

880

881 self.material_density_terms = Expression(

882 self.params.phase_list,

883 self.params.component_list,

884 rule=rule_material_density_terms,

885 )

886 except AttributeError:

887 self.del_component(self.material_density_terms)

888

889 if j in self.params.component_list:

890 return self.material_density_terms[p, j]

891 else:

892 return 0

893

894 def get_enthalpy_density_terms(self, p):

895 """Create energy density terms."""

896 if not self.is_property_constructed("enthalpy_density_terms"):

897 try:

898

899 def rule_energy_density_terms(b, p):

900 return self.dens_mol_phase[p] * self.energy_internal_mol_phase[p]

901

902 self.energy_density_terms = Expression(

903 self.params.phase_list, rule=rule_energy_density_terms

904 )

905 except AttributeError:

906 self.del_component(self.energy_density_terms)

907 return self.enthalpy_density_terms[p]

908

909 def default_material_balance_type(self):

910 return MaterialBalanceType.componentPhase

911

912 def default_energy_balance_type(self):

913 return EnergyBalanceType.enthalpyTotal

914

915 def get_material_flow_basis(b):

916 return MaterialFlowBasis.molar

917

918 def define_state_vars(self):

919 """Define state vars."""

920 return {

921 "flow_mol_phase_comp": self.flow_mol_phase_comp,

922 "temperature": self.temperature,

923 "pressure": self.pressure,

924 }

925

926 # Property package utility functions

927 def calculate_bubble_point_temperature(self, clear_components=True):

928 """ "To compute the bubble point temperature of the mixture."""

929

930 if hasattr(self, "eq_temperature_bubble"):

931 # Do not delete components if the block already has the components

932 clear_components = False

933

934 calculate_variable_from_constraint(

935 self.temperature_bubble, self.eq_temperature_bubble

936 )

937

938 return self.temperature_bubble.value

939

940 if clear_components is True:

941 self.del_component(self.eq_temperature_bubble)

942 self.del_component(self._p_sat_bubbleT)

943 self.del_component(self.temperature_bubble)

944

945 def calculate_dew_point_temperature(self, clear_components=True):

946 """ "To compute the dew point temperature of the mixture."""

947

948 if hasattr(self, "eq_temperature_dew"):

949 # Do not delete components if the block already has the components

950 clear_components = False

951

952 calculate_variable_from_constraint(

953 self.temperature_dew, self.eq_temperature_dew

954 )

955

956 return self.temperature_dew.value

957

958 # Delete the var/constraint created in this method that are part of the

959 # IdealStateBlock if the user desires

960 if clear_components is True:

961 self.del_component(self.eq_temperature_dew)

962 self.del_component(self._p_sat_dewT)

963 self.del_component(self.temperature_dew)

964

965 def calculate_bubble_point_pressure(self, clear_components=True):

966 """ "To compute the bubble point pressure of the mixture."""

967

968 if hasattr(self, "eq_pressure_bubble"):

969 # Do not delete components if the block already has the components

970 clear_components = False

971

972 calculate_variable_from_constraint(

973 self.pressure_bubble, self.eq_pressure_bubble

974 )

975

976 return self.pressure_bubble.value

977

978 # Delete the var/constraint created in this method that are part of the

979 # IdealStateBlock if the user desires

980 if clear_components is True:

981 self.del_component(self.eq_pressure_bubble)

982 self.del_component(self._p_sat_bubbleP)

983 self.del_component(self.pressure_bubble)

984

985 def calculate_dew_point_pressure(self, clear_components=True):

986 """ "To compute the dew point pressure of the mixture."""

987

988 if hasattr(self, "eq_pressure_dew"):

989 # Do not delete components if the block already has the components

990 clear_components = False

991

992 calculate_variable_from_constraint(self.pressure_dew, self.eq_pressure_dew)

993

994 return self.pressure_dew.value

995

996 # Delete the var/constraint created in this method that are part of the

997 # IdealStateBlock if the user desires

998 if clear_components is True:

999 self.del_component(self.eq_pressure_dew)

1000 self.del_component(self._p_sat_dewP)

1001 self.del_component(self.pressure_dew)

1002

1003 # -----------------------------------------------------------------------------

1004 # Bubble and Dew Points

1005 # Ideal-Ideal properties allow for the simplifications below

1006 # Other methods require more complex equations with shadow compositions

1007

1008 # For future work, propose the following:

1009 # Core class writes a set of constraints Phi_L_i == Phi_V_i

1010 # Phi_L_i and Phi_V_i make calls to submethods which add shadow compositions

1011 # as needed

1012 def _temperature_bubble(self):

1013 self.temperature_bubble = Param(

1014 initialize=33.0, units=pyunits.K, doc="Bubble point temperature"

1015 )

1016

1017 def _temperature_dew(self):

1018

1019 self.temperature_dew = Var(

1020 initialize=298.15, units=pyunits.K, doc="Dew point temperature"

1021 )

1022

1023 def rule_psat_dew(b, j):

1024 return (

1025 1e5

1026 * pyunits.Pa

1027 * 10

1028 ** (

1029 b._params.pressure_sat_coeff_A[j]

1030 - b._params.pressure_sat_coeff_B[j]

1031 / (b.temperature_dew + b._params.pressure_sat_coeff_C[j])

1032 )

1033 )

1034

1035 try:

1036 # Try to build expression

1037 self._p_sat_dewT = Expression(

1038 self._params.component_list, rule=rule_psat_dew

1039 )

1040

1041 def rule_temp_dew(b):

1042 return (

1043 b.pressure

1044 * sum(

1045 b.mole_frac_comp[i] / b._p_sat_dewT[i]

1046 for i in ["toluene", "benzene"]

1047 )

1048 - 1

1049 == 0

1050 )

1051

1052 self.eq_temperature_dew = Constraint(rule=rule_temp_dew)

1053 except AttributeError:

1054 # If expression fails, clean up so that DAE can try again later

1055 # Deleting only var/expression as expression construction will fail

1056 # first; if it passes then constraint construction will not fail.

1057 self.del_component(self.temperature_dew)

1058 self.del_component(self._p_sat_dewT)

1059

1060 def _pressure_bubble(self):

1061 self.pressure_bubble = Param(

1062 initialize=1e8, units=pyunits.Pa, doc="Bubble point pressure"

1063 )

1064

1065 def _pressure_dew(self):

1066 self.pressure_dew = Var(

1067 initialize=298.15, units=pyunits.Pa, doc="Dew point pressure"

1068 )

1069

1070 def rule_psat_dew(b, j):

1071 return (

1072 1e5

1073 * pyunits.Pa

1074 * 10

1075 ** (

1076 b._params.pressure_sat_coeff_A[j]

1077 - b._params.pressure_sat_coeff_B[j]

1078 / (b.temperature + b._params.pressure_sat_coeff_C[j])

1079 )

1080 )

1081

1082 try:

1083 # Try to build expression

1084 self._p_sat_dewP = Expression(

1085 self._params.component_list, rule=rule_psat_dew

1086 )

1087

1088 def rule_pressure_dew(b):

1089 return (

1090 b.pressure_dew

1091 * sum(

1092 b.mole_frac_comp[i] / b._p_sat_dewP[i]

1093 for i in ["toluene", "benzene"]

1094 )

1095 - 1

1096 == 0

1097 )

1098

1099 self.eq_pressure_dew = Constraint(rule=rule_pressure_dew)

1100 except AttributeError:

1101 # If expression fails, clean up so that DAE can try again later

1102 # Deleting only var/expression as expression construction will fail

1103 # first; if it passes then constraint construction will not fail.

1104 self.del_component(self.pressure_dew)

1105 self.del_component(self._p_sat_dewP)

1106

1107 # -----------------------------------------------------------------------------

1108 # Liquid phase properties

1109 def _dens_mol_liq(b):

1110 return b.dens_mol_phase["Liq"] == 1e3 * sum(

1111 b.mole_frac_phase_comp["Liq", j]

1112 * b._params.dens_liq_param_1[j]

1113 / b._params.dens_liq_param_2[j]

1114 ** (

1115 1

1116 + (1 - b.temperature / b._params.dens_liq_param_3[j])

1117 ** b._params.dens_liq_param_4[j]

1118 )

1119 for j in ["benzene", "toluene"]

1120 )

1121

1122 def _fug_phase_comp(self):

1123 def fug_phase_comp_rule(b, p, i):

1124 if p == "Liq":

1125 if i in ["hydrogen", "methane"]:

1126 return b.mole_frac_phase_comp["Liq", i]

1127 else:

1128 return b.pressure_sat[i] * b.mole_frac_phase_comp["Liq", i]

1129 else:

1130 if i in ["hydrogen", "methane"]:

1131 return 1e-6

1132 else:

1133 return b.mole_frac_phase_comp["Vap", i] * b.pressure

1134

1135 self.fug_phase_comp = Expression(

1136 self._params.phase_list,

1137 self._params.component_list,

1138 rule=fug_phase_comp_rule,

1139 )

1140

1141 def _pressure_sat(self):

1142 self.pressure_sat = Var(

1143 self._params.component_list,

1144 initialize=101325,

1145 units=pyunits.Pa,

1146 doc="Vapor pressure",

1147 )

1148

1149 def rule_P_sat(b, j):

1150 return (

1151 (

1152 log10(b.pressure_sat[j] / pyunits.Pa * 1e-5)

1153 - b._params.pressure_sat_coeff_A[j]

1154 )

1155 * (b._teq + b._params.pressure_sat_coeff_C[j])

1156 ) == -b._params.pressure_sat_coeff_B[j]

1157

1158 self.eq_pressure_sat = Constraint(self._params.component_list, rule=rule_P_sat)

1159

1160 def _enth_mol_comp_liq(b, j):

1161 return b.enth_mol_phase_comp["Liq", j] * 1e3 == (

1162 (b._params.cp_ig_5["Liq", j] / 5)

1163 * (b.temperature**5 - b._params.temperature_ref**5)

1164 + (b._params.cp_ig_4["Liq", j] / 4)

1165 * (b.temperature**4 - b._params.temperature_ref**4)

1166 + (b._params.cp_ig_3["Liq", j] / 3)

1167 * (b.temperature**3 - b._params.temperature_ref**3)

1168 + (b._params.cp_ig_2["Liq", j] / 2)

1169 * (b.temperature**2 - b._params.temperature_ref**2)

1170 + b._params.cp_ig_1["Liq", j] * (b.temperature - b._params.temperature_ref)

1171 )

1172

1173 def _entr_mol_comp_liq(b, j):

1174 return b.entr_mol_phase_comp["Liq", j] * 1e3 == (

1175 (

1176 (b._params.cp_ig_5["Liq", j] / 4)

1177 * (b.temperature**4 - b._params.temperature_ref**4)

1178 + (b._params.cp_ig_4["Liq", j] / 3)

1179 * (b.temperature**3 - b._params.temperature_ref**3)

1180 + (b._params.cp_ig_3["Liq", j] / 2)

1181 * (b.temperature**2 - b._params.temperature_ref**2)

1182 + b._params.cp_ig_2["Liq", j]

1183 * (b.temperature - b._params.temperature_ref)

1184 + b._params.cp_ig_1["Liq", j]

1185 * log(b.temperature / b._params.temperature_ref)

1186 )

1187 - const.gas_constant

1188 * log(

1189 b.mole_frac_phase_comp["Liq", j] * b.pressure / b._params.pressure_ref

1190 )

1191 )

1192

1193 # -----------------------------------------------------------------------------

1194 # Vapour phase properties

1195 def _dens_mol_vap(b):

1196 return b.pressure == (

1197 b.dens_mol_phase["Vap"] * const.gas_constant * b.temperature

1198 )

1199

1200 def _dh_vap(self):

1201 # heat of vaporization

1202 add_object_reference(self, "dh_vap", self._params.dh_vap)

1203

1204 def _ds_vap(self):

1205 # entropy of vaporization = dh_Vap/T_boil

1206 # TODO : something more rigorous would be nice

1207 self.ds_vap = Var(

1208 self._params.component_list,

1209 initialize=86,

1210 units=pyunits.J / pyunits.mol / pyunits.K,

1211 doc="Entropy of vaporization",

1212 )

1213

1214 def rule_ds_vap(b, j):

1215 return b.dh_vap[j] == (b.ds_vap[j] * b._params.temperature_boil[j])

1216

1217 self.eq_ds_vap = Constraint(self._params.component_list, rule=rule_ds_vap)

1218

1219 def _enth_mol_comp_vap(b, j):

1220 return b.enth_mol_phase_comp["Vap", j] == b.dh_vap[j] + (

1221 (b._params.cp_ig_5["Vap", j] / 5)

1222 * (b.temperature**5 - b._params.temperature_ref**5)

1223 + (b._params.cp_ig_4["Vap", j] / 4)

1224 * (b.temperature**4 - b._params.temperature_ref**4)

1225 + (b._params.cp_ig_3["Vap", j] / 3)

1226 * (b.temperature**3 - b._params.temperature_ref**3)

1227 + (b._params.cp_ig_2["Vap", j] / 2)

1228 * (b.temperature**2 - b._params.temperature_ref**2)

1229 + b._params.cp_ig_1["Vap", j] * (b.temperature - b._params.temperature_ref)

1230 )

1231

1232 def _entr_mol_comp_vap(b, j):

1233 return b.entr_mol_phase_comp["Vap", j] == (

1234 b.ds_vap[j]

1235 + (

1236 (b._params.cp_ig_5["Vap", j] / 4)

1237 * (b.temperature**4 - b._params.temperature_ref**4)

1238 + (b._params.cp_ig_4["Vap", j] / 3)

1239 * (b.temperature**3 - b._params.temperature_ref**3)

1240 + (b._params.cp_ig_3["Vap", j] / 2)

1241 * (b.temperature**2 - b._params.temperature_ref**2)

1242 + b._params.cp_ig_2["Vap", j]

1243 * (b.temperature - b._params.temperature_ref)

1244 + b._params.cp_ig_1["Vap", j]

1245 * log(b.temperature / b._params.temperature_ref)

1246 )

1247 - const.gas_constant

1248 * log(

1249 b.mole_frac_phase_comp["Vap", j] * b.pressure / b._params.pressure_ref

1250 )

1251 )

1252

1253 def calculate_scaling_factors(self):

1254 # Get default scale factors

1255 super().calculate_scaling_factors()

1256

1257 is_two_phase = len(self._params.phase_list) == 2

1258 sf_flow = iscale.get_scaling_factor(self.flow_mol, default=1, warning=True)

1259 sf_T = iscale.get_scaling_factor(self.temperature, default=1, warning=True)

1260 sf_P = iscale.get_scaling_factor(self.pressure, default=1, warning=True)

1261

1262 if self.is_property_constructed("_teq"):

1263 iscale.set_scaling_factor(self._teq, sf_T)

1264 if self.is_property_constructed("_teq_constraint"):

1265 iscale.constraint_scaling_transform(

1266 self._teq_constraint, sf_T, overwrite=False

1267 )

1268

1269 if self.is_property_constructed("_t1"):

1270 iscale.set_scaling_factor(self._t1, sf_T)

1271 if self.is_property_constructed("_t1_constraint"):

1272 iscale.constraint_scaling_transform(

1273 self._t1_constraint, sf_T, overwrite=False

1274 )

1275

1276 if self.is_property_constructed("_mole_frac_pdew"):

1277 iscale.set_scaling_factor(self._mole_frac_pdew, 1e3)

1278 iscale.constraint_scaling_transform(

1279 self._sum_mole_frac_pdew, 1e3, overwrite=False

1280 )

1281

1282 if self.is_property_constructed("total_flow_balance"):

1283 iscale.constraint_scaling_transform(

1284 self.total_flow_balance, sf_flow, overwrite=False

1285 )

1286

1287 if self.is_property_constructed("component_flow_balances"):

1288 for i, c in self.component_flow_balances.items():

1289 if is_two_phase:

1290 s = iscale.get_scaling_factor(

1291 self.mole_frac_comp[i], default=1, warning=True

1292 )

1293 s *= sf_flow

1294 iscale.constraint_scaling_transform(c, s, overwrite=False)

1295 else:

1296 s = iscale.get_scaling_factor(

1297 self.mole_frac_comp[i], default=1, warning=True

1298 )

1299 iscale.constraint_scaling_transform(c, s, overwrite=False)

1300

1301 if self.is_property_constructed("dens_mol_phase"):

1302 for c in self.eq_dens_mol_phase.values():

1303 iscale.constraint_scaling_transform(c, sf_P, overwrite=False)

1304

1305 if self.is_property_constructed("dens_mass_phase"):

1306 for p, c in self.eq_dens_mass_phase.items():

1307 sf = iscale.get_scaling_factor(

1308 self.dens_mass_phase[p], default=1, warning=True

1309 )

1310 iscale.constraint_scaling_transform(c, sf, overwrite=False)

1311

1312 if self.is_property_constructed("enth_mol_phase"):

1313 for p, c in self.eq_enth_mol_phase.items():

1314 sf = iscale.get_scaling_factor(

1315 self.enth_mol_phase[p], default=1, warning=True

1316 )

1317 iscale.constraint_scaling_transform(c, sf, overwrite=False)

1318

1319 if self.is_property_constructed("enth_mol"):

1320 sf = iscale.get_scaling_factor(self.enth_mol, default=1, warning=True)

1321 sf *= sf_flow

1322 iscale.constraint_scaling_transform(self.eq_enth_mol, sf, overwrite=False)

1323

1324 if self.is_property_constructed("entr_mol_phase"):

1325 for p, c in self.eq_entr_mol_phase.items():

1326 sf = iscale.get_scaling_factor(

1327 self.entr_mol_phase[p], default=1, warning=True

1328 )

1329 iscale.constraint_scaling_transform(c, sf, overwrite=False)

1330

1331 if self.is_property_constructed("entr_mol"):

1332 sf = iscale.get_scaling_factor(self.entr_mol, default=1, warning=True)

1333 sf *= sf_flow

1334 iscale.constraint_scaling_transform(self.eq_entr_mol, sf, overwrite=False)

1335

1336 if self.is_property_constructed("gibbs_mol_phase"):

1337 for p, c in self.eq_gibbs_mol_phase.items():

1338 sf = iscale.get_scaling_factor(

1339 self.gibbs_mol_phase[p], default=1, warning=True

1340 )

1341 iscale.constraint_scaling_transform(c, sf, overwrite=False)

Coverage for backend/idaes_service/solver/custom/hda_ideal_VLE.py: 26%

382 statements