|
| 1 | +import numpy as np |
| 2 | +from CoolProp.CoolProp import PropsSI as PSI |
| 3 | +from tespy.components import ( |
| 4 | + Compressor, |
| 5 | + CycleCloser, |
| 6 | + HeatExchanger, |
| 7 | + Merge, |
| 8 | + Motor, |
| 9 | + MovingBoundaryHeatExchanger, |
| 10 | + PowerSource, |
| 11 | + Sink, |
| 12 | + Source, |
| 13 | + Splitter, |
| 14 | + Valve, |
| 15 | +) |
| 16 | +from tespy.connections import Connection, PowerConnection |
| 17 | +from tespy.networks import Network |
| 18 | + |
| 19 | +# "Bus" Parameters |
| 20 | +# Heat source |
| 21 | +T_ev_in = 25 |
| 22 | +T_ev_out = 15 |
| 23 | +p_ev = 1 |
| 24 | +# Water/steam side |
| 25 | +T_co_in = 45 |
| 26 | +T_co_out = 105 |
| 27 | +p_co = 1 |
| 28 | +# T_co_out = 65 |
| 29 | +# p_co = 0.2 |
| 30 | + |
| 31 | +# Default component parameters |
| 32 | +td_pinch = 5 |
| 33 | +td_dew = 5 |
| 34 | +eta_s = 0.8 |
| 35 | + |
| 36 | +wf = "Butane" |
| 37 | + |
| 38 | +# Preliminary calculations, actually optimization variables... |
| 39 | +p_ev_wf = PSI("P", "Q", 1, "T", T_ev_out - td_pinch + 273.15, wf) * 1e-5 |
| 40 | +# p_co_wf = PSI('P', 'Q', 0, 'T', PSI('T', 'Q', 1, 'P', p_co * 1e5, 'water') + td_pinch, wf) * 1e-5 |
| 41 | +p_co_wf = ( |
| 42 | + PSI("P", "Q", 0, "T", PSI("T", "Q", 1, "P", p_co * 1e5, "water") + td_pinch + td_dew, wf) * 1e-5 |
| 43 | +) # hier müsste der Druck wirklich sukzessive gesenkt werden... |
| 44 | +pr_tot = p_co_wf / p_ev_wf |
| 45 | +pr_per_stage = np.power(pr_tot, 1 / 2) |
| 46 | +p_mid_wf = p_ev_wf * pr_per_stage |
| 47 | + |
| 48 | +# Create the network |
| 49 | +nw = Network(T_unit="C", p_unit="bar", h_unit="kJ / kg", m_unit="t / h", iterinfo=False) |
| 50 | + |
| 51 | +# Create the components |
| 52 | +# Ammonia cycle |
| 53 | +compressor_1 = Compressor("Compressor 1") |
| 54 | +compressor_2 = Compressor("Compressor 2") |
| 55 | +evaporator = HeatExchanger("Evaporator") |
| 56 | +expansion_valve_1 = Valve("Expansion Valve 1") |
| 57 | +expansion_valve_2 = Valve("Expansion Valve 2") |
| 58 | +cycle_closer = CycleCloser("Cycle Closer") |
| 59 | +condenser = MovingBoundaryHeatExchanger("Condenser") |
| 60 | +splitter = Splitter("Splitter") |
| 61 | +merge = Merge("Merge") |
| 62 | +economizer = MovingBoundaryHeatExchanger("Economizer") |
| 63 | + |
| 64 | +# Water/steam side |
| 65 | +water_in = Source("Water In") |
| 66 | +steam_out = Sink("Steam Out") |
| 67 | + |
| 68 | +# Heat source |
| 69 | +heat_source_in = Source("Heat Source In") |
| 70 | +heat_source_out = Sink("Heat Source Out") |
| 71 | + |
| 72 | +motor_1 = Motor("Motor 1") |
| 73 | +motor_2 = Motor("Motor 2") |
| 74 | +grid_1 = PowerSource("Grid 1") |
| 75 | +grid_2 = PowerSource("Grid 2") |
| 76 | + |
| 77 | +# Create connections |
| 78 | +# Ammonia cycle |
| 79 | +c1 = Connection(cycle_closer, "out1", compressor_1, "in1", label="01") |
| 80 | +c2 = Connection(compressor_1, "out1", merge, "in1", label="02") |
| 81 | +c3 = Connection(merge, "out1", compressor_2, "in1", label="03") |
| 82 | +c4 = Connection(compressor_2, "out1", condenser, "in1", label="04") |
| 83 | +c5 = Connection(condenser, "out1", splitter, "in1", label="05") |
| 84 | +c5a = Connection(splitter, "out1", economizer, "in1", label="05a") |
| 85 | +c6 = Connection(economizer, "out1", expansion_valve_1, "in1", label="06") |
| 86 | +c7 = Connection(expansion_valve_1, "out1", evaporator, "in2", label="07") |
| 87 | +c0cc = Connection(evaporator, "out2", cycle_closer, "in1", label="07cc") |
| 88 | +c5b = Connection(splitter, "out2", expansion_valve_2, "in1", label="05b") |
| 89 | +c8 = Connection(expansion_valve_2, "out1", economizer, "in2", label="08") |
| 90 | +c9 = Connection(economizer, "out2", merge, "in2", label="09") |
| 91 | + |
| 92 | +# Water/steam side |
| 93 | +c10 = Connection(water_in, "out1", condenser, "in2", label="10") |
| 94 | +c11 = Connection(condenser, "out2", steam_out, "in1", label="11") |
| 95 | + |
| 96 | +# Heat source |
| 97 | +c20 = Connection(heat_source_in, "out1", evaporator, "in1", label="20") |
| 98 | +c21 = Connection(evaporator, "out1", heat_source_out, "in1", label="21") |
| 99 | + |
| 100 | +# Power connections |
| 101 | +e01 = PowerConnection(grid_1, "power", motor_1, "power_in", label="E1") |
| 102 | +w01 = PowerConnection(motor_1, "power_out", compressor_1, "power", label="W1") |
| 103 | +e02 = PowerConnection(grid_2, "power", motor_2, "power_in", label="E2") |
| 104 | +w02 = PowerConnection(motor_2, "power_out", compressor_2, "power", label="W2") |
| 105 | + |
| 106 | +# Add connections to network |
| 107 | +nw.add_conns(c1, c2, c3, c4, c5, c5a, c5b, c6, c7, c8, c9, c0cc, c10, c11, c20, c21, e01, e02, w01, w02) |
| 108 | + |
| 109 | +# Set component parameters |
| 110 | +compressor_1.set_attr(eta_s=eta_s) |
| 111 | +compressor_2.set_attr(eta_s=eta_s) |
| 112 | +condenser.set_attr(pr1=1, pr2=1, td_pinch=td_pinch) |
| 113 | +evaporator.set_attr(pr1=1, pr2=1, ttd_min=td_pinch) |
| 114 | +economizer.set_attr(pr1=1, pr2=1, td_pinch=td_pinch) |
| 115 | +motor_1.set_attr(eta=0.985) |
| 116 | +motor_2.set_attr(eta=0.985) |
| 117 | + |
| 118 | +# Set connection parameters |
| 119 | +# Water/steam side |
| 120 | +c10.set_attr(fluid={"water": 1}, p=p_co, T=T_co_in, m=1) # Water inlet |
| 121 | +c11.set_attr(T=T_co_out) # Steam outlet |
| 122 | + |
| 123 | +# Heat source |
| 124 | +c20.set_attr(fluid={"air": 1}, p=p_ev, T=T_ev_in) # Heat source inlet |
| 125 | +c21.set_attr(T=T_ev_out) # Heat source outlet |
| 126 | + |
| 127 | +# Butane cycle |
| 128 | +c7.set_attr(fluid={wf: 1}) # Pure butane |
| 129 | +c1.set_attr(p=p_ev_wf) |
| 130 | +c2.set_attr(p=p_mid_wf) |
| 131 | +c4.set_attr(p=p_co_wf) |
| 132 | + |
| 133 | +# Generate starting values |
| 134 | +condenser.set_attr(td_pinch=None) |
| 135 | +evaporator.set_attr(ttd_min=None) |
| 136 | +economizer.set_attr(td_pinch=None) |
| 137 | + |
| 138 | +c1.set_attr(T=T_ev_in - td_pinch) |
| 139 | +c4.set_attr(td_dew=td_dew) |
| 140 | +c5.set_attr(T=T_co_in + td_pinch) |
| 141 | +c6.set_attr(T=T_co_in) |
| 142 | +nw.solve("design") |
| 143 | + |
| 144 | +# Solve real design |
| 145 | +condenser.set_attr(td_pinch=td_pinch) # For condenser |
| 146 | +c5.set_attr(T=None) |
| 147 | +evaporator.set_attr(ttd_min=td_pinch) |
| 148 | +c1.set_attr(T=None) |
| 149 | +economizer.set_attr(td_pinch=td_pinch) |
| 150 | +c6.set_attr(T=None) |
| 151 | + |
| 152 | +nw.solve("design") |
| 153 | + |
| 154 | +# Exergy analysis via ExerPy |
| 155 | +from exerpy import ExergyAnalysis |
| 156 | + |
| 157 | +p0 = 1.000e5 |
| 158 | +T0 = 25 + 273.15 |
| 159 | + |
| 160 | +ean = ExergyAnalysis.from_tespy(nw, Tamb=T0, pamb=p0, split_physical_exergy=True) |
| 161 | +fuel = {"inputs": ["E1", "E2"], "outputs": []} |
| 162 | +product = {"inputs": ["11"], "outputs": ["10"]} |
| 163 | +loss = {"inputs": ["21"], "outputs": ["20"]} |
| 164 | +ean.analyse(E_F=fuel, E_P=product, E_L=loss) |
| 165 | +df_component_results, df_material_connection_results, df_non_material_connection_results = ean.exergy_results( |
| 166 | + print_results=False |
| 167 | +) |
| 168 | + |
| 169 | +from exerpy import EconomicAnalysis, ExergoeconomicAnalysis |
| 170 | + |
| 171 | +# Define the CEPCI values for cost correction. |
| 172 | +CEPCI_2013 = 567.3 |
| 173 | +CEPCI_2023 = 797.9 |
| 174 | +CEPCI_factor = CEPCI_2023 / CEPCI_2013 |
| 175 | + |
| 176 | +# Define values for electricity price and full load hours. |
| 177 | +elec_price_cent_kWh = 40.0 # cent/kWh |
| 178 | +tau = 5500 # hours/year |
| 179 | + |
| 180 | +# Define economic parameters. |
| 181 | +r_n = 0.02 # Cost elevation rate |
| 182 | +i_eff = 0.08 # Interest rate |
| 183 | +n = 20 # Number of years |
| 184 | +omc_relative = 0.03 # Relative operation and maintenance costs (compared to PEC) |
| 185 | + |
| 186 | +import pandas as pd |
| 187 | +from tabulate import tabulate |
| 188 | + |
| 189 | +k_lookup = {"Evaporator": 100, "Economizer": 750, "Condenser": 1000} |
| 190 | + |
| 191 | +PEC_computed = {} |
| 192 | +for comp in nw.comps["object"]: |
| 193 | + name = comp.label |
| 194 | + PEC = 0.0 # Default PEC |
| 195 | + # --- Heat Exchangers --- |
| 196 | + if isinstance(comp, (HeatExchanger, MovingBoundaryHeatExchanger)): |
| 197 | + kA = getattr(comp.kA, "val", 0.0) |
| 198 | + if kA: |
| 199 | + k = k_lookup.get(name) |
| 200 | + if k is None: |
| 201 | + raise KeyError(f"No k-value defined for heat exchanger '{name}'") |
| 202 | + A = kA / k |
| 203 | + PEC = 15526 * (A / 42) ** 0.8 * CEPCI_factor |
| 204 | + else: |
| 205 | + PEC = 0.0 |
| 206 | + PEC_computed[name] = PEC |
| 207 | + |
| 208 | + # --- Compressors (and Fans) --- |
| 209 | + elif isinstance(comp, Compressor): |
| 210 | + VM = getattr(comp.inl[0].v, "val", 0.0) |
| 211 | + PEC = 19850 * ((VM * 3600) / 279.8) ** 0.73 |
| 212 | + PEC *= CEPCI_factor # Adjust PEC cost. |
| 213 | + PEC_computed[name] = PEC |
| 214 | + |
| 215 | + # --- Other components --- |
| 216 | + else: |
| 217 | + PEC_computed[name] = 0.0 |
| 218 | + |
| 219 | + |
| 220 | +for comp in ean.components.values(): |
| 221 | + if comp.__class__.__name__ == "Motor": |
| 222 | + name = comp.name |
| 223 | + # Retrieve the electrical input power X from the attribute "energy_flow_1" |
| 224 | + X = getattr(comp, "E_F", None) |
| 225 | + if X is not None: |
| 226 | + PEC = 10710 * (X / 250000) ** 0.65 |
| 227 | + PEC *= CEPCI_factor # Adjust PEC cost. |
| 228 | + else: |
| 229 | + PEC = 0.0 |
| 230 | + PEC_computed[name] = PEC |
| 231 | + |
| 232 | +# ------------------------------ |
| 233 | +# Economic Analysis |
| 234 | +# ------------------------------ |
| 235 | +# Convert electricity price from cent/kWh to €/GJ. |
| 236 | +# 1 kWh = 3.6 MJ and 1 GJ = 277.78 kWh. |
| 237 | +elec_cost_eur_per_kWh = elec_price_cent_kWh / 100.0 |
| 238 | +elec_cost_eur_per_GJ = elec_cost_eur_per_kWh * 277.78 |
| 239 | + |
| 240 | +econ_pars = {"tau": tau, "i_eff": i_eff, "n": n, "r_n": r_n} |
| 241 | +components_order = list(PEC_computed.keys()) |
| 242 | +PEC_list = [PEC_computed[comp] for comp in components_order] |
| 243 | +# Multiply each PEC by 6.32 to obtain TCI. |
| 244 | +TCI_list = [pec * 6.32 for pec in PEC_list] |
| 245 | +OMC_relative = [omc_relative if pec > 0 else 0.0 for pec in TCI_list] |
| 246 | + |
| 247 | +econ_analysis = EconomicAnalysis(econ_pars) |
| 248 | +Z_CC, Z_OMC, Z_total = econ_analysis.compute_component_costs(TCI_list, OMC_relative) |
| 249 | + |
| 250 | +# Create a DataFrame to display PEC, TCI, annual OMC, and Z for each component |
| 251 | +component_costs_df = pd.DataFrame( |
| 252 | + { |
| 253 | + "Component": components_order, |
| 254 | + "PEC [EUR]": [round(pec, 2) for pec in PEC_list], |
| 255 | + "CC [EUR]": [round(tci, 2) for tci in TCI_list], |
| 256 | + "Z_CC [EUR/h]": [round(z_cc, 2) for z_cc in Z_CC], |
| 257 | + "Z_OMC [EUR/h]": [round(omc, 2) for omc in Z_OMC], |
| 258 | + "Z [EUR/h]": [round(z, 2) for z in Z_total], |
| 259 | + } |
| 260 | +) |
| 261 | + |
| 262 | +# Calculate totals |
| 263 | +total_pec = sum(PEC_list) |
| 264 | +total_tci = sum(TCI_list) |
| 265 | +total_z_cc = sum(Z_CC) |
| 266 | +total_z_omc = sum(Z_OMC) |
| 267 | +total_z = sum(Z_total) |
| 268 | + |
| 269 | +# Add a total row |
| 270 | +component_costs_df.loc[len(component_costs_df)] = [ |
| 271 | + "TOTAL", |
| 272 | + round(total_pec, 2), |
| 273 | + round(total_tci, 2), |
| 274 | + round(total_z_cc, 2), |
| 275 | + round(total_z_omc, 2), |
| 276 | + round(total_z, 2), |
| 277 | +] |
| 278 | + |
| 279 | +# Print the component costs table without separators |
| 280 | +print("\nComponent Investment Costs (Year 2023):") |
| 281 | +print(tabulate(component_costs_df, headers="keys", tablefmt="psql", floatfmt=".2f")) |
| 282 | + |
| 283 | +# Build the exergoeconomic cost dictionary. |
| 284 | +Exe_Eco_Costs = {} |
| 285 | +for comp, z in zip(components_order, Z_total, strict=False): |
| 286 | + Exe_Eco_Costs[f"{comp}_Z"] = z |
| 287 | +Exe_Eco_Costs["10_c"] = 0.0 |
| 288 | +Exe_Eco_Costs["20_c"] = 0.0 |
| 289 | +Exe_Eco_Costs["E2_c"] = elec_cost_eur_per_GJ |
| 290 | + |
| 291 | +# ------------------------------ |
| 292 | +# Exergoeconomic Analysis |
| 293 | +# ------------------------------ |
| 294 | +exergoeco_analysis = ExergoeconomicAnalysis(ean) |
| 295 | +exergoeco_analysis.run(Exe_Eco_Costs=Exe_Eco_Costs, Tamb=ean.Tamb) |
| 296 | +# Unpack four DataFrames; we only use the component results. |
| 297 | +df_comp, df_mat1, df_mat2, df_non_mat = exergoeco_analysis.exergoeconomic_results() |
0 commit comments