Coverage for backend/ahuora-compounds/ahuora_property_packages/modular/builder/common_parsers.py: 15%

1from math import floor 

2from typing import Any, Dict, List 

3from ahuora_compounds.packages.chemsep import ChemsepCompound as Compound 

4from pyomo.environ import units as pyunits 

5from idaes.models.properties.modular_properties.state_definitions import FTPx 

6from idaes.models.properties.modular_properties.phase_equil.bubble_dew import (LogBubbleDew) 

7from idaes.core import LiquidPhase, VaporPhase, Component, PhaseType as PT 

8from idaes.models.properties.modular_properties.phase_equil import (SmoothVLE) 

9from idaes.models.properties.modular_properties.phase_equil.forms import log_fugacity 

10from idaes.models.properties.modular_properties.eos.ceos import Cubic, CubicType 

11from idaes.models.properties.modular_properties.pure import RPP4, Perrys 

12from ahuora_property_packages.modular.builder.data.chem_sep import ChemSep 

13from pyomo.common.fileutils import this_file_dir 

14from ahuora_property_packages.types import States 

15import csv 

16 

17def build_base_units(compounds: List[Compound], valid_states: List[States]) -> Dict[str, Any]: 

18 return { 

19 'time': pyunits.s, 

20 'length': pyunits.m, 

21 'mass': pyunits.kg, 

22 'amount': pyunits.mol, 

23 'temperature': pyunits.K, 

24 } 

25 

26def build_bubble_dew_method(compounds: List[Compound], valid_states: List[States]) -> Any: 

27 return LogBubbleDew 

28 

29def build_components(compounds: List[Compound], valid_states: List[States]) -> Dict[str, Any]: 

30 

31 def _serialise_component(compound: Compound) -> Dict[str, Any]: 

32 

33 # configuration default to all components 

34 config = { 

35 "type": Component, 

36 "parameter_data": { 

37 "mw": (compound.MolecularWeight.value, pyunits.kg/pyunits.kilomol), 

38 "pressure_crit": (compound.CriticalPressure.value, pyunits.Pa), 

39 "temperature_crit": (compound.CriticalTemperature.value, pyunits.K), 

40 "omega": compound.AcentricityFactor.value, 

41 } 

42 } 

43 

44 valid_phase = _valid_phases(compound) 

45 if valid_phase != PT.vaporPhase: 

46 config["phase_equilibrium_form"] = {("Vap", "Liq"): log_fugacity} 

47 else: 

48 config["valid_phase_types"] = valid_phase 

49 

50 # Energies of Formation 

51 if compound.HeatOfFormation is not None: # this does not work when passed 

52 config["parameter_data"].update({ 

53 "enth_mol_form_vap_comp_ref": (compound.HeatOfFormation.value, pyunits.J/pyunits.kilomol) 

54 }) 

55 else: 

56 raise ValueError("No Heat of Formation Data") 

57 

58 if compound.AbsEntropy is not None: 

59 config["parameter_data"].update({ 

60 "entr_mol_form_vap_comp_ref": (-1 * compound.AbsEntropy.value, pyunits.J/pyunits.kilomol/pyunits.K) 

61 }) 

62 else: 

63 raise ValueError("No Absolute Entropy Data") 

64 

65 # Ideal Gas Molar Calculations 

66 # All three properties intrinsically linked together 

67 if compound.RPPHeatCapacityCp is not None: 

68 if compound.RPPHeatCapacityCp["eqno"] == 4: 

69 config["enth_mol_ig_comp"] = RPP4 

70 config["entr_mol_ig_comp"] = RPP4 

71 config["parameter_data"].update({"cp_mol_ig_comp_coeff": { 

72 "A": (float(compound.RPPHeatCapacityCp["A"]), pyunits.J / pyunits.kilomol / pyunits.K), 

73 "B": (float(compound.RPPHeatCapacityCp["B"]), pyunits.J / pyunits.kilomol / pyunits.K**2), 

74 "C": (float(compound.RPPHeatCapacityCp["C"]), pyunits.J / pyunits.kilomol / pyunits.K**3), 

75 "D": (float(compound.RPPHeatCapacityCp["D"]), pyunits.J / pyunits.kilomol / pyunits.K**4), 

76 }}) 

77 elif compound.RPPHeatCapacityCp["eqno"] == 100 or compound.RPPHeatCapacityCp["eqno"] == 5: 

78 config["enth_mol_ig_comp"] = ChemSep 

79 config["entr_mol_ig_comp"] = ChemSep 

80 config["parameter_data"].update({"cp_mol_ig_comp_coeff": { 

81 "A": (compound.RPPHeatCapacityCp["A"], pyunits.J / pyunits.kilomol / pyunits.K), 

82 "B": (compound.RPPHeatCapacityCp["B"], pyunits.J / pyunits.kilomol / pyunits.K**2), 

83 "C": (compound.RPPHeatCapacityCp["C"], pyunits.J / pyunits.kilomol / pyunits.K**3), 

84 "D": (compound.RPPHeatCapacityCp["D"], pyunits.J / pyunits.kilomol / pyunits.K**4), 

85 "E": (compound.RPPHeatCapacityCp["E"], pyunits.J / pyunits.kilomol / pyunits.K**5), 

86 }}) 

87 else: 

88 raise ValueError(f"Invalid equation number for heat capacity {compound.RPPHeatCapacityCp['eqno']}") 

89 else: 

90 raise ValueError("No Heat Capacity Data") 

91 

92 # Saturation Pressure (Vapor) 

93 if compound.AntoineVaporPressure is not None: 

94 if compound.AntoineVaporPressure["eqno"] == 10: 

95 config["pressure_sat_comp"] = ChemSep 

96 config["parameter_data"].update({"pressure_sat_comp_coeff": { 

97 "A": (compound.AntoineVaporPressure["A"], None), 

98 "B": (compound.AntoineVaporPressure["B"], pyunits.K), 

99 "C": (compound.AntoineVaporPressure["C"], pyunits.K), 

100 }}) 

101 else: 

102 raise ValueError("No Antoine Vapor Pressure Equation Data") 

103 else: 

104 raise ValueError("No Antoine Vapor Pressure Data") 

105 

106 # Liquid Density 

107 if compound.LiquidDensity is not None: 

108 if compound.LiquidDensity["eqno"] == 105: 

109 config["dens_mol_liq_comp"] = Perrys 

110 config["parameter_data"].update({"dens_mol_liq_comp_coeff": { 

111 "eqn_type": 1, 

112 "1": (compound.LiquidDensity["A"], pyunits.kmol / pyunits.m**3), 

113 "2": (compound.LiquidDensity["B"], None), 

114 "3": (compound.LiquidDensity["C"], pyunits.K), 

115 "4": (compound.LiquidDensity["D"], None), 

116 }}) 

117 elif compound.LiquidDensity["eqno"] == 106: 

118 config["dens_mol_liq_comp"] = ChemSep 

119 config["parameter_data"].update({"dens_mol_liq_comp_coeff": { 

120 "A": (compound.LiquidDensity["A"], pyunits.kmol / pyunits.m**3), 

121 "B": (compound.LiquidDensity["B"], None), 

122 "C": (compound.LiquidDensity["C"], None), 

123 "D": (compound.LiquidDensity["D"], None), 

124 "E": (compound.LiquidDensity["E"], None), 

125 }}) 

126 else: 

127 raise ValueError("No Liquid Density Equation Data") 

128 else: 

129 raise ValueError("No Liquid Density Data") 

130 

131 # Liquid Heat Capacity & Entropy / Enthalpy 

132 if compound.LiquidHeatCapacityCp is not None: 

133 if compound.LiquidHeatCapacityCp["eqno"] == 100: 

134 # Uses correct equations to calculate 

135 config["enth_mol_liq_comp"] = Perrys 

136 config["entr_mol_liq_comp"] = Perrys 

137 config["parameter_data"].update({"cp_mol_liq_comp_coeff": { 

138 "1": (compound.LiquidHeatCapacityCp["A"], pyunits.J / pyunits.kilomol / pyunits.K**1), 

139 "2": (compound.LiquidHeatCapacityCp["B"], pyunits.J / pyunits.kilomol / pyunits.K**2), 

140 "3": (compound.LiquidHeatCapacityCp["C"], pyunits.J / pyunits.kilomol / pyunits.K**3), 

141 "4": (compound.LiquidHeatCapacityCp["D"], pyunits.J / pyunits.kilomol / pyunits.K**4), 

142 "5": (compound.LiquidHeatCapacityCp["E"], pyunits.J / pyunits.kilomol / pyunits.K**5), 

143 }}) 

144 

145 # ASSUMPTION: Molar heat of formation, liq is zero - given the semi-okay by Ben 

146 config["parameter_data"].update({ 

147 "enth_mol_form_liq_comp_ref": (0, pyunits.J / pyunits.kilomol) 

148 }) 

149 

150 config["parameter_data"].update({ 

151 "entr_mol_form_liq_comp_ref": (0, pyunits.J / pyunits.kilomol / pyunits.K) 

152 }) 

153 elif compound.LiquidHeatCapacityCp["eqno"] == 16: 

154 # Uses correct equations to calculate 

155 config["enth_mol_liq_comp"] = ChemSep 

156 config["entr_mol_liq_comp"] = ChemSep 

157 config["parameter_data"].update({"cp_mol_liq_comp_coeff": { 

158 "A": (compound.LiquidHeatCapacityCp["A"], pyunits.J / pyunits.kilomol / pyunits.K**1), 

159 "B": (compound.LiquidHeatCapacityCp["B"], pyunits.J / pyunits.kilomol / pyunits.K**2), 

160 "C": (compound.LiquidHeatCapacityCp["C"], pyunits.J / pyunits.kilomol / pyunits.K**3), 

161 "D": (compound.LiquidHeatCapacityCp["D"], pyunits.J / pyunits.kilomol / pyunits.K**4), 

162 "E": (compound.LiquidHeatCapacityCp["E"], pyunits.J / pyunits.kilomol / pyunits.K**5), 

163 }}) 

164 

165 # ASSUMPTION: Molar heat of formation, liq is zero - given the semi-okay by Ben 

166 config["parameter_data"].update({ 

167 "enth_mol_form_liq_comp_ref": (0, pyunits.J / pyunits.kilomol) 

168 }) 

169 

170 config["parameter_data"].update({ 

171 "entr_mol_form_liq_comp_ref": (0, pyunits.J / pyunits.kilomol / pyunits.K) 

172 }) 

173 else: 

174 raise ValueError(f"No Liquid Heat Capacity Equation Data {compound.LiquidHeatCapacityCp['eqno']}") 

175 else: 

176 # Compound only exists in vapor phase 

177 config["parameter_data"].update({"valid_phase_types": PT.vaporPhase}) 

178 

179 return config 

180 

181 def _valid_phases(compound: Compound) -> PT: 

182 if compound.NormalMeltingPointTemperature.value > compound.NormalBoilingPointTemperature.value: 

183 # no liquid phase exists (sublimation) 

184 return PT.vaporPhase 

185 # Assumption: Anything above hydrogen can exist as both liquid and vapor 

186 elif compound.NormalBoilingPointTemperature.value >= 21: 

187 return [PT.liquidPhase, PT.vaporPhase] 

188 else: 

189 return PT.vaporPhase 

190 

191 components_output = {} 

192 for compound in compounds: 

193 components_output[compound.CompoundID.value] = _serialise_component(compound) 

194 return components_output 

195 

196def build_phase_equilibrium_state(compounds: List[Compound], valid_states: List[States]) -> Dict[str, Any]: 

197 return {("Vap", "Liq"): SmoothVLE} 

198 

199def build_phases(compounds: List[Compound], valid_states: List[States]) -> Dict[str, Any]: 

200 phases = {} 

201 for state in valid_states: 

202 if state == "Liq": 

203 phases["Liq"] = { 

204 "type": LiquidPhase, 

205 "equation_of_state": Cubic, 

206 "equation_of_state_options": {"type": CubicType.PR}, 

207 } 

208 elif state == "Vap": 

209 phases["Vap"] = { 

210 "type": VaporPhase, 

211 "equation_of_state": Cubic, 

212 "equation_of_state_options": {"type": CubicType.PR}, 

213 } 

214 return phases 

215 

216def build_phases_in_equilibrium(compounds: List[Compound], valid_states: List[States]) -> List: 

217 return [("Vap", "Liq")] 

218 

219def build_pressure_ref(compounds: List[Compound], valid_states: List[States]) -> tuple: 

220 return (101325, pyunits.Pa) 

221 

222def build_state_bounds(compounds: List[Compound], valid_states: List[States]) -> Dict[str, Any]: 

223 """ 

224 State bounds are used to find optimal solution 

225 TODO: need to find a way to dynamically determine these 

226 """ 

227 min_melting_point = min([compound.NormalMeltingPointTemperature.value for compound in compounds]) 

228 min([compound.CriticalTemperature.value for compound in compounds]) 

229 

230 # TODO: Refactor this logic, need a more versatile approach 

231 return { 

232 "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), 

233 "temperature": (max(min_melting_point-50,1), 300, 3000, pyunits.K), 

234 "pressure": (5e4, 1e5, 1e6, pyunits.Pa), 

235 } 

236 

237def build_state_definition(compounds: List[Compound], valid_states: List[States]) -> Any: 

238 return FTPx 

239 

240def build_temperature_ref(compounds: List[Compound], valid_states: List[States]) -> tuple: 

241 return (298.15, pyunits.K) 

242 

243def build_pr_kappa(compounds: List[Compound], valid_states: List[States]) -> Dict[str, Any]: 

244 kappa_parameters = {} 

245 compound_id_map = {} 

246 

247 for i, compound1 in enumerate(compounds): 

248 compound_id_map[str(floor(compound1.LibraryIndex.value))] = compound1 

249 for j, compound2 in enumerate(compounds): 

250 kappa_parameters[(compound1.CompoundID.value, compound2.CompoundID.value)] = 0.000 

251 # Setting all interactions initially to zero 

252 

253 # TODO: Adjust method so the multiple kappa values for single pair are supported 

254 # Open and read the interaction data file 

255 file = open(this_file_dir() + "/data/pr.dat", 'r') 

256 reader = csv.reader(file, delimiter=';') 

257 for row in reader: 

258 

259 # Skip invalid rows 

260 if len(row) < 4: 

261 continue 

262 

263 # Extract ID1, ID2, kappa (k12) and comment 

264 id1, id2, kappa, _ = row 

265 

266 try: 

267 kappa_value = float(kappa) 

268 except ValueError: 

269 continue 

270 

271 # Check if the compound IDs exist in the compound list 

272 if id1 in compound_id_map and id2 in compound_id_map: 

273 compound1 = compound_id_map[id1] 

274 compound2 = compound_id_map[id2] 

275 kappa_parameters[(compound1.CompoundID.value, compound2.CompoundID.value)] = kappa_value 

276 kappa_parameters[(compound2.CompoundID.value, compound1.CompoundID.value)] = kappa_value 

277 

278 return {"PR_kappa": kappa_parameters}