Coverage for backend/ahuora-compounds/ahuora_property_packages/saltwater/salt_water_property_package_VLE_V2.py: 11%

1""" 

2Salt water (water-NaCl) property package for IDAES unit models. 

3 

4This version is written specifically to be more robust for Heater-based 

5heating / evaporation cases that cross phase boundaries. The main numerical 

6changes relative to the earlier draft are: 

7 

81. Composition regularization uses a moderate flow epsilon rather than an 

9 extremely tiny value, which avoids huge Jacobian entries when liquid flow 

10 collapses toward zero. 

112. The outlet phase-split equations are written only when 

12 ``has_phase_equilibrium=True`` and ``defined_state=False`` in line with the 

13 IDAES custom property package pattern. 

143. The liquid-vapor and liquid-solid complementarity equations use smooth phase 

15 presence functions. This avoids forcing equilibrium when one of the phases 

16 is absent, which is essential for vapor-solid solutions such as 400 K, 

17 1 atm brine where no liquid can exist. 

184. A heuristic state-block initializer seeds the outlet state with a physically 

19 sensible phase regime (L, LV, VS, or LVS) based on the current T and P. 

20 

21The thermodynamic model remains deliberately simple: 

22* Modified Raoult law for water volatility. 

23* NaCl treated as nonvolatile. 

24* NaCl solubility expressed as a liquid-phase mass-fraction limit. 

25* Constant heat capacities and latent/solid offsets for enthalpy. 

26""" 

27 

28from math import log as ln 

29 

30import idaes.logger as idaeslog 

31import idaes.core.util.scaling as iscale 

32 

33from pyomo.common.config import ConfigValue, In 

34from pyomo.environ import ( 

35 Constraint, 

36 Expression, 

37 NonNegativeReals, 

38 Param, 

39 Set, 

40 Var, 

41 exp, 

42 value, 

43) 

44from pyomo.environ import units as pyunits 

45 

46from idaes.core import ( 

47 declare_process_block_class, 

48 EnergyBalanceType, 

49 MaterialBalanceType, 

50 MaterialFlowBasis, 

51 PhysicalParameterBlock, 

52 StateBlock, 

53 StateBlockData, 

54) 

55from idaes.core.base.components import Solute, Solvent 

56from idaes.core.base.phases import LiquidPhase, SolidPhase, VaporPhase, PhaseType as PT 

57from idaes.core.util.constants import Constants 

58from idaes.core.util.initialization import fix_state_vars, revert_state_vars 

59from idaes.core.util.misc import extract_data 

60 

61 

62_log = idaeslog.getLogger(__name__) 

63 

64 

65@declare_process_block_class("SaltWaterParameterBlock") 

66class SaltWaterParameterData(PhysicalParameterBlock): 

67 """Parameter block for an water-NaCl three-phase property package.""" 

68 

69 CONFIG = PhysicalParameterBlock.CONFIG() 

70 

71 CONFIG.declare( 

72 "activity_coefficient_model", 

73 ConfigValue( 

74 default="Ideal", 

75 domain=In(["Ideal", "Constant"]), 

76 description="Activity-coefficient model used in the modified Raoult law", 

77 doc=( 

78 "'Ideal' fixes gamma_i = 1. 'Constant' uses user-adjustable constant " 

79 "gamma_i parameters as a placeholder for a fuller model." 

80 ), 

81 ), 

82 ) 

83 

84 def build(self): 

85 super().build() 

86 self._state_block_class = SaltWaterStateBlock # noqa: F821 

87 

88 # Components 

89 self.water = Solvent(valid_phase_types=[PT.liquidPhase, PT.vaporPhase]) 

90 self.NaCl = Solute(valid_phase_types=[PT.liquidPhase, PT.solidPhase]) 

91 

92 # Phases 

93 self.Liq = LiquidPhase(component_list=["water", "NaCl"]) 

94 self.Vap = VaporPhase(component_list=["water"]) 

95 self.Sol = SolidPhase(component_list=["NaCl"]) 

96 

97 # Reference state 

98 self.temperature_ref = Param( 

99 initialize=298.15, 

100 units=pyunits.K, 

101 doc="Reference temperature for enthalpy calculations", 

102 ) 

103 self.pressure_ref = Param( 

104 initialize=101325.0, 

105 units=pyunits.Pa, 

106 doc="Reference pressure", 

107 ) 

108 

109 # Smoothing / scaling parameters. 

110 self.flow_epsilon = Param( 

111 initialize=1e-6, 

112 mutable=True, 

113 units=pyunits.kg / pyunits.s, 

114 doc=( 

115 "Regularization flow used in composition denominators and phase " 

116 "presence functions. Kept deliberately moderate to avoid extreme " 

117 "Jacobians near zero phase flow." 

118 ), 

119 ) 

120 self.flow_scale_ref = Param( 

121 initialize=1.0, 

122 mutable=True, 

123 units=pyunits.kg / pyunits.s, 

124 doc="Reference flow used to nondimensionalize complementarity constraints", 

125 ) 

126 self.pressure_scale_ref = Param( 

127 initialize=1e6, 

128 mutable=True, 

129 units=pyunits.Pa, 

130 doc="Reference pressure used to nondimensionalize complementarity constraints", 

131 ) 

132 

133 # Molecular weights 

134 self.mw_comp = Param( 

135 self.component_list, 

136 initialize=extract_data({"water": 18.01528e-3, "NaCl": 58.44e-3}), 

137 units=pyunits.kg / pyunits.mol, 

138 doc="Molecular weight", 

139 ) 

140 

141 # Placeholder activity coefficients for future model upgrades. 

142 self.gamma_constant = Param( 

143 self.component_list, 

144 mutable=True, 

145 initialize=extract_data({"water": 1.0, "NaCl": 1.0}), 

146 units=pyunits.dimensionless, 

147 doc="Constant activity coefficients used when model='Constant'", 

148 ) 

149 

150 self._phase_component_set_local = Set( 

151 initialize=[ 

152 ("Liq", "water"), 

153 ("Liq", "NaCl"), 

154 ("Vap", "water"), 

155 ("Sol", "NaCl"), 

156 ], 

157 dimen=2, 

158 ordered=True, 

159 doc="Valid phase-component pairs used for parameter indexing", 

160 ) 

161 

162 # Heat capacity data on a mass basis. 

163 cp_data = { 

164 ("Liq", "water"): 4182.0, 

165 ("Liq", "NaCl"): 850.0, 

166 ("Vap", "water"): 1864.0, 

167 ("Sol", "NaCl"): 850.0, 

168 } 

169 self.cp_mass_phase_comp_ref = Param( 

170 self._phase_component_set_local, 

171 initialize=extract_data(cp_data), 

172 units=pyunits.J / pyunits.kg / pyunits.K, 

173 doc="Reference mass heat capacity for each valid phase-component pair", 

174 ) 

175 

176 self.dh_vap_mass_comp = Param( 

177 ["water"], 

178 initialize=extract_data({"water": 2.257e6}), 

179 units=pyunits.J / pyunits.kg, 

180 doc="Latent heat offset used for vapor-phase water enthalpy", 

181 ) 

182 

183 self.dh_crystallization_mass_comp = Param( 

184 ["NaCl"], 

185 initialize=extract_data({"NaCl": -5.20e5}), 

186 units=pyunits.J / pyunits.kg, 

187 doc="Solid-phase enthalpy offset for NaCl relative to reference liquid state", 

188 ) 

189 

190 # Simple density data. 

191 self.dens_mass_phase_ref = Param( 

192 ["Liq", "Sol"], 

193 initialize=extract_data({"Liq": 1025.0, "Sol": 2160.0}), 

194 units=pyunits.kg / pyunits.m**3, 

195 doc="Representative phase densities for liquid brine and solid salt", 

196 ) 

197 

198 # NaCl solubility correlation coefficients (mass fraction). 

199 self.solubility_mass_frac_param_A1 = Param( 

200 initialize=0.2628, 

201 units=pyunits.dimensionless, 

202 ) 

203 self.solubility_mass_frac_param_A2 = Param( 

204 initialize=62.75e-6, 

205 units=pyunits.K**-1, 

206 ) 

207 self.solubility_mass_frac_param_A3 = Param( 

208 initialize=1.084e-6, 

209 units=pyunits.K**-2, 

210 ) 

211 

212 # Antoine parameters for water, valid over the range typically used in 

213 # low-to-moderate temperature brine unit operations. 

214 self.antoine_A = Param(initialize=8.07131, units=pyunits.dimensionless) 

215 self.antoine_B = Param(initialize=1730.63, units=pyunits.dimensionless) 

216 self.antoine_C = Param(initialize=233.426, units=pyunits.dimensionless) 

217 

218 # Default scaling 

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

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

221 self.set_default_scaling("flow_mass_phase_comp", 1.0, index=("Liq", "water")) 

222 self.set_default_scaling("flow_mass_phase_comp", 1e2, index=("Liq", "NaCl")) 

223 self.set_default_scaling("flow_mass_phase_comp", 1.0, index=("Vap", "water")) 

224 self.set_default_scaling("flow_mass_phase_comp", 1e2, index=("Sol", "NaCl")) 

225 

226 @classmethod 

227 def define_metadata(cls, obj): 

228 obj.add_default_units( 

229 { 

230 "time": pyunits.s, 

231 "length": pyunits.m, 

232 "mass": pyunits.kg, 

233 "amount": pyunits.mol, 

234 "temperature": pyunits.K, 

235 } 

236 ) 

237 

238 obj.add_properties( 

239 { 

240 "flow_mass_phase_comp": {"method": None}, 

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

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

243 "flow_mol_phase_comp": {"method": "_flow_mol_phase_comp"}, 

244 "mass_frac_phase_comp": {"method": "_mass_frac_phase_comp"}, 

245 "mole_frac_phase_comp": {"method": "_mole_frac_phase_comp"}, 

246 "enth_mass_phase_comp": {"method": "_enth_mass_phase_comp"}, 

247 "enth_mass_phase": {"method": "_enth_mass_phase"}, 

248 "pressure_sat_comp": {"method": "_pressure_sat_comp"}, 

249 "dens_mass_phase": {"method": "_dens_mass_phase"}, 

250 } 

251 ) 

252 

253 obj.define_custom_properties( 

254 { 

255 "flow_mass_phase": {"method": "_flow_mass_phase"}, 

256 "flow_mol_phase": {"method": "_flow_mol_phase"}, 

257 "enth_flow_phase": {"method": "_enth_flow_phase"}, 

258 "activity_coeff_phase_comp": {"method": "_activity_coeff_phase_comp"}, 

259 "modified_raoult_pressure_comp": {"method": "_modified_raoult_pressure_comp"}, 

260 "vle_residual_comp": {"method": "_vle_residual_comp"}, 

261 "solubility_mass_frac_phase_comp": {"method": "_solubility_mass_frac_phase_comp"}, 

262 "saturation_index_phase_comp": {"method": "_saturation_index_phase_comp"}, 

263 } 

264 ) 

265 

266 

267class _SaltWaterStateBlock(StateBlock): 

268 """Methods applied to indexed state blocks as a whole.""" 

269 

270 def fix_initialization_states(self): 

271 fix_state_vars(self) 

272 

273 def initialize( 

274 self, 

275 state_args=None, 

276 state_vars_fixed=False, 

277 hold_state=False, 

278 outlvl=idaeslog.NOTSET, 

279 solver=None, 

280 optarg=None, 

281 ): 

282 """ 

283 Lightweight initializer. 

284 

285 This package is mostly expression-based, but for outlet blocks with 

286 phase equilibrium enabled it helps to seed a physically sensible phase 

287 regime before the unit model solve. 

288 """ 

289 flags = None 

290 if not state_vars_fixed: 

291 flags = fix_state_vars(self, state_args) 

292 

293 # for k in self.values(): 

294 # if getattr(k.config, "has_phase_equilibrium", False) and not k.config.defined_state: 

295 # k._set_sensible_initial_guesses() 

296 

297 if hold_state: 

298 return flags 

299 

300 if not state_vars_fixed: 

301 self.release_state(flags, outlvl=outlvl) 

302 

303 def release_state(self, flags, outlvl=idaeslog.NOTSET): 

304 if flags is None: 

305 return 

306 revert_state_vars(self, flags) 

307 

308 

309@declare_process_block_class("SaltWaterStateBlock", block_class=_SaltWaterStateBlock) 

310class SaltWaterStateBlockData(StateBlockData): 

311 """State block for salt water with liquid/vapor/solid phases.""" 

312 

313 def build(self): 

314 super().build() 

315 self._make_state_vars() 

316 self._make_extra_vars() 

317 

318 

319 #if self.config.has_phase_equilibrium or self.config.defined_state: 

320 self._make_phase_split_equations() 

321 

322 def _make_state_vars(self): 

323 self.pressure = Var( 

324 domain=NonNegativeReals, 

325 initialize=101325.0, 

326 units=pyunits.Pa, 

327 doc="State pressure", 

328 ) 

329 

330 self.temperature = Var( 

331 domain=NonNegativeReals, 

332 initialize=298.15, 

333 bounds=(250.0, 650.0), 

334 units=pyunits.K, 

335 doc="State temperature", 

336 ) 

337 

338 self.mole_frac_comp = Var( 

339 self.component_list, 

340 initialize={ 

341 ("water"): 0.965, 

342 ("NaCl"): 0.035, 

343 }, 

344 bounds=(0, None), 

345 domain=NonNegativeReals, 

346 units=pyunits.dimensionless, 

347 doc="Mole fractions of components", 

348 ) 

349 

350 self.flow_mol = Var( 

351 domain=NonNegativeReals, 

352 initialize=100 / 18, 

353 units=pyunits.mol / pyunits.s, 

354 doc="Total molar flow rate",) 

355 

356 

357 

358 

359 def _make_extra_vars(self): 

360 # TODO: Refactor this into more build on demand properties. 

361 self.flow_mass_phase_comp = Var( 

362 self.phase_component_set, 

363 domain=NonNegativeReals, 

364 initialize={ 

365 ("Liq", "water"): 0.94, 

366 ("Liq", "NaCl"): 0.03, 

367 ("Vap", "water"): 0.02, 

368 ("Sol", "NaCl"): 0.005, 

369 }, 

370 units=pyunits.kg / pyunits.s, 

371 doc="Phase-component mass flow rate", 

372 ) 

373 @self.Constraint(self.component_list, doc="Mole fraction definition") 

374 def eq_mole_frac_comp(b, j): 

375 return b.mole_frac_comp[j] == sum( 

376 b.flow_mass_phase_comp[p, j] / b.params.mw_comp[j] for p in b.phase_list if (p, j) in b.phase_component_set 

377 ) / b.flow_mol 

378 

379 if self.config.defined_state is False: 

380 @self.Constraint(doc="Total molar flow rate definition") 

381 def total_molar_flow_definition(b): 

382 return b.flow_mol - sum( 

383 b.flow_mass_phase_comp[p, j] / b.params.mw_comp[j] for (p, j) in b.phase_component_set 

384 ) == 0 

385 

386 self.enth_mol = Var( 

387 initialize=1000, 

388 units=pyunits.J / pyunits.mol, 

389 doc="Molar enthalpy", 

390 ) 

391 self.enth_mol.fix(1) 

392 

393 self.entr_mol = Var( 

394 initialize=100, 

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

396 doc="Molar entropy", 

397 ) 

398 self.entr_mol.fix(1) 

399 

400 def rule_vapor_frac(b): 

401 return b.flow_mass_phase_comp["Vap", "water"] / ( 

402 b.flow_mass_phase_comp["Liq", "water"] + b.flow_mass_phase_comp["Vap", "water"] 

403 ) 

404 

405 self.vapor_frac = Expression(rule=rule_vapor_frac) 

406 

407 self.flow_mass = Var( 

408 domain=NonNegativeReals, 

409 initialize=100, 

410 units=pyunits.kg / pyunits.s, 

411 doc="Total mass flow rate", 

412 ) 

413 @self.Constraint(doc="Total mass flow rate definition") 

414 def total_mass_flow_definition(b): 

415 return b.flow_mass - sum( 

416 b.flow_mass_phase_comp[p, j] for (p, j) in b.phase_component_set 

417 ) == 0 

418 

419 def rule_flow_vol(b): 

420 return sum( 

421 b.flow_mass_phase_comp[p, j] / b.dens_mass_phase[p] 

422 for (p, j) in b.phase_component_set 

423 ) 

424 

425 self.flow_vol = Expression(rule=rule_flow_vol) 

426 

427 

428 

429 # ------------------------------------------------------------------------- 

430 # Build-on-demand properties 

431 def _flow_mass_phase(self): 

432 def rule_flow_mass_phase(b, p): 

433 return sum( 

434 b.flow_mass_phase_comp[p, j] 

435 for j in b.params.component_list 

436 if (p, j) in b.phase_component_set 

437 ) 

438 

439 self.flow_mass_phase = Expression(self.params.phase_list, rule=rule_flow_mass_phase) 

440 

441 def _flow_mol_phase_comp(self): 

442 def rule_flow_mol_phase_comp(b, p, j): 

443 return b.flow_mass_phase_comp[p, j] / b.params.mw_comp[j] 

444 

445 self.flow_mol_phase_comp = Expression( 

446 self.phase_component_set, rule=rule_flow_mol_phase_comp 

447 ) 

448 

449 def _flow_mol_phase(self): 

450 def rule_flow_mol_phase(b, p): 

451 return sum( 

452 b.flow_mol_phase_comp[p, j] 

453 for j in b.params.component_list 

454 if (p, j) in b.phase_component_set 

455 ) 

456 

457 self.flow_mol_phase = Expression(self.params.phase_list, rule=rule_flow_mol_phase) 

458 

459 def _mass_frac_phase_comp(self): 

460 def rule_mass_frac_phase_comp(b, p, j): 

461 denom = b.flow_mass_phase[p] + b.params.flow_epsilon 

462 return b.flow_mass_phase_comp[p, j] / denom 

463 

464 self.mass_frac_phase_comp = Expression( 

465 self.phase_component_set, rule=rule_mass_frac_phase_comp 

466 ) 

467 

468 def _mole_frac_phase_comp(self): 

469 def rule_mole_frac_phase_comp(b, p, j): 

470 denom = b.flow_mol_phase[p] + b.params.flow_epsilon / b.params.mw_comp[j] 

471 return b.flow_mol_phase_comp[p, j] / denom 

472 

473 self.mole_frac_phase_comp = Expression( 

474 self.phase_component_set, rule=rule_mole_frac_phase_comp 

475 ) 

476 

477 def _enth_mass_phase_comp(self): 

478 def rule_enth_mass_phase_comp(b, p, j): 

479 dT = b.temperature - b.params.temperature_ref 

480 if (p, j) == ("Liq", "water"): 

481 return b.params.cp_mass_phase_comp_ref[p, j] * dT 

482 if (p, j) == ("Liq", "NaCl"): 

483 return b.params.cp_mass_phase_comp_ref[p, j] * dT 

484 if (p, j) == ("Vap", "water"): 

485 return b.params.dh_vap_mass_comp["water"] + b.params.cp_mass_phase_comp_ref[p, j] * dT 

486 if (p, j) == ("Sol", "NaCl"): 

487 return b.params.dh_crystallization_mass_comp["NaCl"] + b.params.cp_mass_phase_comp_ref[p, j] * dT 

488 return 0 * pyunits.J / pyunits.kg 

489 

490 self.enth_mass_phase_comp = Expression( 

491 self.phase_component_set, rule=rule_enth_mass_phase_comp 

492 ) 

493 

494 def _enth_flow_phase(self): 

495 def rule_enth_flow_phase(b, p): 

496 return sum( 

497 b.flow_mass_phase_comp[p, j] * b.enth_mass_phase_comp[p, j] 

498 for j in b.params.component_list 

499 if (p, j) in b.phase_component_set 

500 ) 

501 

502 self.enth_flow_phase = Expression(self.params.phase_list, rule=rule_enth_flow_phase) 

503 

504 def _enth_mass_phase(self): 

505 def rule_enth_mass_phase(b, p): 

506 return b.enth_flow_phase[p] / (b.flow_mass_phase[p] + b.params.flow_epsilon) 

507 

508 self.enth_mass_phase = Expression(self.params.phase_list, rule=rule_enth_mass_phase) 

509 

510 def _pressure_sat_comp(self): 

511 def rule_pressure_sat_comp(b, j): 

512 if j == "water": 

513 t_c = (b.temperature - 273.15 * pyunits.K) / pyunits.K 

514 exponent = b.params.antoine_A - b.params.antoine_B / (t_c + b.params.antoine_C) 

515 return 133.322368 * pyunits.Pa * exp(ln(10.0) * exponent) 

516 return 1e-9 * pyunits.Pa 

517 

518 self.pressure_sat_comp = Expression( 

519 self.params.component_list, rule=rule_pressure_sat_comp 

520 ) 

521 

522 def _activity_coeff_phase_comp(self): 

523 def rule_activity_coeff_phase_comp(b, p, j): 

524 if p != "Liq": 

525 return 1.0 

526 if b.params.config.activity_coefficient_model == "Ideal": 

527 return 1.0 

528 return b.params.gamma_constant[j] 

529 

530 self.activity_coeff_phase_comp = Expression( 

531 self.phase_component_set, rule=rule_activity_coeff_phase_comp 

532 ) 

533 

534 def _modified_raoult_pressure_comp(self): 

535 def rule_modified_raoult_pressure_comp(b, j): 

536 if ("Liq", j) not in b.phase_component_set: 

537 return 0 * pyunits.Pa 

538 if j != "water": 

539 return 0 * pyunits.Pa 

540 return ( 

541 b.mole_frac_phase_comp["Liq", j] 

542 * b.activity_coeff_phase_comp["Liq", j] 

543 * b.pressure_sat_comp[j] 

544 ) 

545 

546 self.modified_raoult_pressure_comp = Expression( 

547 self.params.component_list, rule=rule_modified_raoult_pressure_comp 

548 ) 

549 

550 def _vle_residual_comp(self): 

551 def rule_vle_residual_comp(b, j): 

552 if j != "water": 

553 return 0 * pyunits.Pa 

554 vap_term = 0 * pyunits.Pa 

555 if ("Vap", j) in b.phase_component_set: 

556 vap_term = b.mole_frac_phase_comp["Vap", j] * b.pressure 

557 return vap_term - b.modified_raoult_pressure_comp[j] 

558 

559 self.vle_residual_comp = Expression( 

560 self.params.component_list, rule=rule_vle_residual_comp 

561 ) 

562 

563 def _solubility_mass_frac_phase_comp(self): 

564 def rule_solubility_mass_frac_phase_comp(b, p, j): 

565 if (p, j) != ("Liq", "NaCl"): 

566 return 0.0 

567 dt = b.temperature - 273.15 * pyunits.K 

568 return ( 

569 b.params.solubility_mass_frac_param_A1 

570 + b.params.solubility_mass_frac_param_A2 * dt 

571 + b.params.solubility_mass_frac_param_A3 * dt**2 

572 ) 

573 

574 self.solubility_mass_frac_phase_comp = Expression( 

575 [("Liq", "NaCl")], rule=rule_solubility_mass_frac_phase_comp 

576 ) 

577 

578 def _saturation_index_phase_comp(self): 

579 def rule_saturation_index_phase_comp(b, p, j): 

580 if (p, j) != ("Liq", "NaCl"): 

581 return 0.0 

582 return b.mass_frac_phase_comp[p, j] / ( 

583 b.solubility_mass_frac_phase_comp[p, j] + 1e-12 

584 ) 

585 

586 self.saturation_index_phase_comp = Expression( 

587 [("Liq", "NaCl")], rule=rule_saturation_index_phase_comp 

588 ) 

589 

590 def _dens_mass_phase(self): 

591 def rule_dens_mass_phase(b, p): 

592 if p == "Liq": 

593 return b.params.dens_mass_phase_ref["Liq"] 

594 if p == "Sol": 

595 return b.params.dens_mass_phase_ref["Sol"] 

596 return b.pressure / ((Constants.gas_constant / b.params.mw_comp["water"]) * b.temperature) 

597 

598 self.dens_mass_phase = Expression(self.params.phase_list, rule=rule_dens_mass_phase) 

599 

600 def _make_phase_split_equations(self): 

601 if not self.is_property_constructed("flow_mass_phase"): 

602 self._flow_mass_phase() 

603 if not self.is_property_constructed("flow_mol_phase_comp"): 

604 self._flow_mol_phase_comp() 

605 if not self.is_property_constructed("flow_mol_phase"): 

606 self._flow_mol_phase() 

607 if not self.is_property_constructed("mass_frac_phase_comp"): 

608 self._mass_frac_phase_comp() 

609 if not self.is_property_constructed("mole_frac_phase_comp"): 

610 self._mole_frac_phase_comp() 

611 if not self.is_property_constructed("pressure_sat_comp"): 

612 self._pressure_sat_comp() 

613 if not self.is_property_constructed("activity_coeff_phase_comp"): 

614 self._activity_coeff_phase_comp() 

615 if not self.is_property_constructed("modified_raoult_pressure_comp"): 

616 self._modified_raoult_pressure_comp() 

617 if not self.is_property_constructed("solubility_mass_frac_phase_comp"): 

618 self._solubility_mass_frac_phase_comp() 

619 

620 self.pressure_equil_slack = Var( 

621 domain=NonNegativeReals, 

622 initialize=1e3, 

623 units=pyunits.Pa, 

624 doc="Nonnegative slack for water modified-Raoult pressure inequality", 

625 ) 

626 self.pressure_equil_vap_slack = Var( 

627 domain=NonNegativeReals, 

628 initialize=1e3, 

629 units=pyunits.Pa, 

630 doc="Nonnegative slack for water modified-Raoult pressure inequality", 

631 ) 

632 self.solubility_slack = Var( 

633 domain=NonNegativeReals, 

634 initialize=0.05, 

635 units=pyunits.dimensionless, 

636 doc="Nonnegative slack for NaCl solubility inequality", 

637 ) 

638 

639 def rule_phase_presence(b, p): 

640 return b.flow_mass_phase[p] / (sum(b.flow_mass_phase[p] for p in b.params.phase_list) + b.params.flow_epsilon) 

641 

642 self.phase_presence = Expression(self.params.phase_list, rule=rule_phase_presence) 

643 

644 def rule_water_phase_presence(b, p): 

645 # amount of the total amount of water in the system that is in phase p 

646 if (p, "water") not in b.phase_component_set: 

647 return 0.0 

648 return b.flow_mass_phase_comp[p,"water"] / ( 

649 sum(b.flow_mass_phase_comp[phase,"water"] for phase in b.params.phase_list 

650 if (phase, "water") in b.phase_component_set) 

651 + b.params.flow_epsilon) 

652 

653 self.water_phase_presence = Expression(self.params.phase_list, rule=rule_water_phase_presence) 

654 

655 def rule_eq_pressure_slack(b): 

656 return (b.pressure_equil_vap_slack - b.pressure_equil_slack)/1e5 == (b.pressure - b.modified_raoult_pressure_comp["water"])/1e5 

657 

658 self.eq_pressure_equil_slack = Constraint(rule=rule_eq_pressure_slack) 

659 

660 def rule_eq_solubility_slack(b): 

661 return b.solubility_slack == ( 

662 b.solubility_mass_frac_phase_comp["Liq", "NaCl"] 

663 - b.mass_frac_phase_comp["Liq", "NaCl"] 

664 ) 

665 

666 self.eq_solubility_slack = Constraint(rule=rule_eq_solubility_slack) 

667 

668 def rule_vle_complementarity(b): 

669 return ( 

670 b.water_phase_presence["Liq"] 

671 * b.pressure_equil_slack 

672 / b.params.pressure_scale_ref 

673 == 0 

674 ) 

675 

676 self.eq_vle_complementarity = Constraint(rule=rule_vle_complementarity) 

677 

678 def rule_vle_vap_complementarity(b): 

679 return ( 

680 b.phase_presence["Vap"] 

681 * b.pressure_equil_vap_slack 

682 / b.params.pressure_scale_ref 

683 == 0 

684 ) 

685 

686 self.eq_vle_vap_complementarity = Constraint(rule=rule_vle_vap_complementarity) 

687 

688 def rule_sle_complementarity(b): 

689 return ( 

690 b.phase_presence["Liq"] 

691 * b.phase_presence["Sol"] 

692 * b.solubility_slack 

693 #* b.mole_frac_comp["NaCl"] # This is to say that the slack can be as large as you want if there is barely any Nacl 

694 == 0 

695 ) 

696 

697 self.eq_sle_complementarity = Constraint(rule=rule_sle_complementarity) 

698 

699 def _set_sensible_initial_guesses(self): 

700 """Heuristic initialization for phase-equilibrium outlet states.""" 

701 if not hasattr(self, "pressure_equil_slack"): 

702 return 

703 

704 # Ensure needed expressions exist. 

705 if not self.is_property_constructed("pressure_sat_comp"): 

706 self._pressure_sat_comp() 

707 if not self.is_property_constructed("solubility_mass_frac_phase_comp"): 

708 self._solubility_mass_frac_phase_comp() 

709 if not self.is_property_constructed("activity_coeff_phase_comp"): 

710 self._activity_coeff_phase_comp() 

711 

712 # Current totals are usually good proxies for the inlet copied by the 

713 # unit model initializer. 

714 total_water = value(self.flow_mass_phase_comp["Liq", "water"]) + value( 

715 self.flow_mass_phase_comp["Vap", "water"] 

716 ) 

717 total_nacl = value(self.flow_mass_phase_comp["Liq", "NaCl"]) + value( 

718 self.flow_mass_phase_comp["Sol", "NaCl"] 

719 ) 

720 

721 total_water = max(total_water, 1e-8) 

722 total_nacl = max(total_nacl, 0.0) 

723 

724 P = max(value(self.pressure), 1.0) 

725 psat = max(value(self.pressure_sat_comp["water"]), 1.0) 

726 gamma_w = value(self.activity_coeff_phase_comp["Liq", "water"]) 

727 wsat = min( 

728 0.999, 

729 max(1e-8, value(self.solubility_mass_frac_phase_comp["Liq", "NaCl"])), 

730 ) 

731 

732 mw_w = value(self.params.mw_comp["water"]) 

733 mw_s = value(self.params.mw_comp["NaCl"]) 

734 

735 # Water mole fraction in saturated liquid at the solubility limit. 

736 n_w_sat = (1.0 - wsat) / mw_w 

737 n_s_sat = wsat / mw_s 

738 xw_sat = n_w_sat / (n_w_sat + n_s_sat) 

739 

740 xw_eq = min(0.999, max(1e-6, P / max(gamma_w * psat, 1.0))) 

741 

742 # Case 1: No liquid can exist because even saturated brine has too high 

743 # a water partial pressure. Seed vapor + solid. 

744 if total_nacl > 1e-10 and xw_eq < xw_sat: 

745 liq_water = 1e-6 

746 liq_nacl = 1e-6 

747 vap_water = max(total_water - liq_water, 1e-6) 

748 sol_nacl = max(total_nacl - liq_nacl, 1e-8) 

749 

750 # Case 2: Potential LV or LVS state. 

751 elif psat > P: 

752 if total_nacl <= 1e-12: 

753 liq_water = 0.5 * total_water 

754 vap_water = total_water - liq_water 

755 liq_nacl = 0.0 

756 sol_nacl = 0.0 

757 else: 

758 n_salt_total = total_nacl / mw_s 

759 n_water_liq = xw_eq / max(1.0 - xw_eq, 1e-8) * n_salt_total 

760 liq_water_guess = n_water_liq * mw_w 

761 liq_water = min(max(liq_water_guess, 1e-6), max(total_water - 1e-6, 1e-6)) 

762 vap_water = max(total_water - liq_water, 1e-6) 

763 

764 # Try dissolving all salt first. 

765 liq_nacl_trial = total_nacl 

766 if liq_water + liq_nacl_trial > 0: 

767 w_nacl_trial = liq_nacl_trial / (liq_water + liq_nacl_trial) 

768 else: 

769 w_nacl_trial = 0.0 

770 

771 if w_nacl_trial <= wsat + 1e-8: 

772 liq_nacl = liq_nacl_trial 

773 sol_nacl = 0.0 

774 else: 

775 liq_nacl = wsat / max(1.0 - wsat, 1e-8) * liq_water 

776 liq_nacl = min(liq_nacl, total_nacl) 

777 sol_nacl = max(total_nacl - liq_nacl, 1e-8) 

778 

779 # Case 3: Liquid only or liquid + solid. 

780 else: 

781 vap_water = 1e-6 

782 liq_water = max(total_water - vap_water, 1e-6) 

783 max_dissolved_nacl = wsat / max(1.0 - wsat, 1e-8) * liq_water