Coverage for backend/pinch_service/OpenPinch/src/analysis/power_cogeneration_analysis.py: 6%
175 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +0000
1from __future__ import annotations
3from typing import TYPE_CHECKING
4from ..utils import *
5from ..classes import *
6from .support_methods import *
8if TYPE_CHECKING:
9 from ..classes import *
11__all__ = ["get_power_cogeneration_above_pinch"]
13#######################################################################################################
14# Public API --- TODO
15#######################################################################################################
17def get_power_cogeneration_above_pinch(z: Zone):
18 """Calculate the power cogeneration potential above pinch for a given process zone."""
20 # === Step 1: Prepare turbine and model parameters ===
21 turbine_params = _prepare_turbine_parameters(z.config)
22 if turbine_params is None:
23 return z
25 # === Step 2: Preprocess utilities ===
26 utility_data = _preprocess_utilities(z, turbine_params)
27 if utility_data is None:
28 return z
30 # === Step 3: Solve turbine work and efficiency ===
31 w_total, Wmax = _solve_turbine_work(turbine_params, utility_data)
33 # === Step 4: Assign back to zone ===
34 z.work_target = w_total
35 z.turbine_efficiency_target = w_total / Wmax if Wmax > ZERO else 0.0
37 return z
40#######################################################################################################
41# Helper functions
42#######################################################################################################
44def _prepare_turbine_parameters(config: Configuration):
45 """Load and sanitize turbine parameters from config."""
47 P_in = float(config.P_TURBINE_BOX) # bar
48 T_in = float(config.T_TURBINE_BOX) # °C
49 min_eff = float(config.MIN_EFF) # minimum isentropic efficiency (decimal)
50 model = config.COMBOBOX
52 load_frac = min(max(config.LOAD, 0), 1) # clamp between 0 and 1
53 mech_eff = min(max(config.MECH_EFF, 0), 1) # clamp between 0 and 1
55 return {
56 "P_in": P_in,
57 "T_in": T_in,
58 "min_eff": min_eff,
59 "model": model,
60 "load_frac": load_frac,
61 "mech_eff": mech_eff,
62 "CONDESATE_FLASH_CORRECTION": getattr(config, "CONDESATE_FLASH_CORRECTION", False), # optional
63 }
66def _preprocess_utilities(z: Zone, turbine_params: dict):
67 """Identify eligible hot utilities above pinch and prepare initial turbine stream data."""
69 P_in = turbine_params["P_in"]
71 Q_HU = 0.0
72 HU_num = 0
73 u: Stream
74 for u in z.hot_utilities:
75 if u.t_supply < T_CRIT and u.heat_flow > ZERO:
76 HU_num += 1
77 Q_HU += u.heat_flow
79 if Q_HU < ZERO:
80 return None # No turbine opportunity
82 # Initialize arrays
83 P_out = [P_in]
84 Q_users = [0]
85 w_k = [0]
86 w_isen_k = [0]
87 m_k = [0]
88 eff_k = [0]
89 dh_is_k = [0]
90 h_out = [h_pT(P_in, turbine_params["T_in"])]
91 h_tar = [hL_p(P_in)]
92 h_sat = [hV_p(P_in)]
93 turbine = [0]
95 m_in_est = 0.0 # estimated steam inlet mass flow
96 Q_boiler = 0.0 # total boiler heat required
98 s = 0
99 for u in z.hot_utilities:
100 if u.t_supply < T_CRIT and u.heat_flow > ZERO:
101 P_out_n = psat_T(u.t_target) if abs(u.t_supply - u.t_target) < 1 + ZERO else psat_T(u.t_target + u.dt_cont * 2)
103 if P_in >= P_out_n:
104 s += 1
105 for lst in [w_k, w_isen_k, eff_k, turbine]:
106 lst.append(0)
108 P_out.append(P_out_n)
109 Q_users.append(u.heat_flow)
110 Q_boiler += u.heat_flow
112 h_out.append(hV_p(P_out_n))
113 h_tar.append(hL_p(P_out_n))
114 h_sat.append(hV_p(P_out_n))
116 dh_is = h_out[0] - h_ps(P_out_n, s_ph(P_out[0], h_out[0]))
117 dh_is_k.append(dh_is)
118 mflow = u.heat_flow / (h_out[0] - dh_is - h_tar[-1])
119 m_k.append(mflow)
120 m_in_est += mflow
122 return {
123 "P_out": P_out,
124 "Q_users": Q_users,
125 "w_k": w_k,
126 "w_isen_k": w_isen_k,
127 "m_k": m_k,
128 "eff_k": eff_k,
129 "dh_is_k": dh_is_k,
130 "h_out": h_out,
131 "h_tar": h_tar,
132 "h_sat": h_sat,
133 "turbine": turbine,
134 "m_in_est": m_in_est,
135 "Q_boiler": Q_boiler,
136 "s": s + 1,
137 }
140def _solve_turbine_work(turbine_params: dict, utility_data: dict):
141 """Solve turbine mass flow, work production, and efficiency based on utilities and turbine parameters."""
143 P_out = utility_data["P_out"]
144 Q_users = utility_data["Q_users"]
145 w_k = utility_data["w_k"]
146 w_isen_k = utility_data["w_isen_k"]
147 m_k = utility_data["m_k"]
148 eff_k = utility_data["eff_k"]
149 dh_is_k = utility_data["dh_is_k"]
150 h_out = utility_data["h_out"]
151 h_tar = utility_data["h_tar"]
152 h_sat = utility_data["h_sat"]
153 s = utility_data["s"]
155 model = turbine_params["model"]
156 load_frac = turbine_params["load_frac"]
157 n_mech = turbine_params["mech_eff"]
158 min_eff = turbine_params["min_eff"]
159 flash_correction = turbine_params.get("CONDESATE_FLASH_CORRECTION", False)
161 m_in_est = utility_data["m_in_est"]
162 i = 0
164 while True:
165 m_in = m_in_est
166 m_in_k = m_in_est
167 m_in_est = 0
169 for j in range(1, s):
170 m_in_k -= m_k[j - 1]
171 dh_is_k[j] = h_out[j - 1] - h_ps(P_out[j], s_ph(P_out[j - 1], h_out[j - 1]))
172 w_isen_k[j] = m_in_k * dh_is_k[j]
173 m_max = m_in_k / load_frac if load_frac > 0 else 0
175 if m_in_k > 0:
176 if model == 'Sun & Smith (2015)':
177 w_k[j] = Work_SunModel(P_out[j-1], h_out[j-1], P_out[j], h_sat[j], m_in_k, m_max, dh_is_k[j], n_mech)
178 elif model == 'Medina-Flores et al. (2010)':
179 w_k[j] = Work_MedinaModel(P_out[j-1], m_in_k, dh_is_k[j])
180 elif model == 'Varbanov et al. (2004)':
181 w_k[j] = Work_THM(P_out[j-1], h_out[j-1], P_out[j], h_sat[j], m_in_k, dh_is_k[j], n_mech)
182 elif model == 'Fixed Isentropic Turbine':
183 w_k[j] = w_isen_k[j] * min_eff
184 else:
185 w_k[j] = 0
187 if w_isen_k[j] > 0:
188 eff_k[j] = w_k[j] / w_isen_k[j]
189 else:
190 eff_k[j] = min_eff
192 if eff_k[j] <= min_eff:
193 w_k[j] = min_eff * w_isen_k[j]
195 if m_in_k > 0:
196 h_out[j] = h_out[j-1] - w_k[j] / (m_in_k * n_mech)
197 else:
198 h_out[j] = h_out[j-1]
200 if m_in_k > 0:
201 if flash_correction:
202 Q_flash = m_k[j-1] * (h_tar[j-1] - h_tar[j])
203 m_k[j] = (Q_users[j] - Q_flash) / (h_out[j] - h_tar[j])
204 else:
205 m_k[j] = Q_users[j] / (h_out[j] - h_tar[j])
206 else:
207 m_k[j] = 0
209 m_in_est += m_k[j]
211 i += 1
212 if abs(m_in - m_in_est) < ZERO or i >= 3:
213 break
215 w_total = sum(w_k)
216 Wmax = sum(w_isen_k)
218 return w_total, Wmax
223# def get_power_cogeneration_above_pinch(z: Zone): # type: ignore
224# config: Configuration = z.config
225# P_in = float(config.P_TURBINE_BOX) # bar
226# T_in = float(config.T_TURBINE_BOX) # C
227# MinEff = float(config.MIN_EFF) # %
228# SelectedModel = config.COMBOBOX
230# load_frac = config.LOAD
231# n_mech = config.MECH_EFF
233# if load_frac > 1:
234# load_frac = 1
235# elif load_frac < 0:
236# load_frac = 0
238# if n_mech > 1:
239# n_mech = 1
240# elif n_mech < 0:
241# n_mech = 0
243# Q_HU = 0
244# for Utility_k in z.hot_utilities:
245# Q_HU += Utility_k.heat_flow
246# HU_num = len(z.hot_utilities)
247# if Q_HU < ZERO:
248# return z
250# if SelectedModel == 'Sun & Smith (2015)':
251# SunCoef = [ [ [None for k in range(3)] for j in range(3)] for i in range(2)]
252# Set_Coeff(SunCoef, None)
253# elif SelectedModel == 'Varbanov et al. (2004)':
254# VarCoef = [ [ [None for k in range(4)] for j in range(2)] for i in range(2)]
255# Set_Coeff(None, VarCoef)
257# P_out = [P_in] # bar
258# Q_users = [0]
259# w_k = [0]
260# w_isen_k = [0]
261# m_k = [0] # kg/s
262# eff_k = [0]
263# dh_is_k = [0]
264# h_out = [h_pT(P_in, T_in)] # kJ/kg for turbine
265# h_tar = [hL_p(P_out[0])] # kJ/kg for after condensing
266# h_sat = [hV_p(P_out[0])] # kJ/kg
267# turbine = [0]
269# m_in_est = 0 # kg/s
270# Q_boiler = 0 # kW
272# # Preparation of data for turbine calculation
273# s = 0
274# for i in range(HU_num):
275# Utility_i = z.hot_utilities[i]
276# if Utility_i.t_supply < T_CRIT:
277# P_out_n = psat_T(Utility_i.t_target) if abs(Utility_i.t_supply - Utility_i.t_target) < 1 + ZERO \
278# else psat_T(Utility_i.t_target + Utility_i.dt_cont * 2) # bar
279# if P_in >= P_out_n and Utility_i.heat_flow > ZERO:
280# s += 1
281# for l in [w_k, w_isen_k, eff_k, turbine]:
282# l.append(0)
284# P_out.append(P_out_n) # bar
285# Q_users.append(Utility_i.heat_flow) # kW
286# Q_boiler += Q_users[s] # kW
288# h_out.append(hV_p(P_out[s])) # kJ/kg
289# h_tar.append(hL_p(P_out[s])) # kJ/kg
290# h_sat.append(hV_p(P_out[s])) # kJ/kg
292# dh_is_k.append(h_out[0] - h_ps(P_out[s], s_ph(P_out[0], h_out[0]))) # kJ/kg
293# m_k.append(Q_users[s] / (h_out[0] - dh_is_k[s] - h_tar[s])) # kg/s
294# m_in_est += m_k[s] # kg/s
295# s += 1
297# # Turbine calculation
298# i = 0
299# while True:
300# m_in = m_in_est
301# m_in_k = m_in_est
302# m_in_est = 0
304# for j in range(1, s):
305# m_in_k -= m_k[j - 1] # kg/s
306# dh_is_k[j] = h_out[j - 1] - h_ps(P_out[j], s_ph(P_out[j - 1], h_out[j - 1])) # kJ/kg
307# w_isen_k[j] = m_in_k * dh_is_k[j] # kW
308# m_max = m_in_k / load_frac # kg/s
310# # Determine the work production based on various Turbine Models
311# if m_in_k > 0:
312# if SelectedModel == 'Sun & Smith (2015)':
313# w_k[j] = Work_SunModel(P_out[j - 1], h_out[j - 1], P_out[j], h_sat[j], m_in_k, m_max, dh_is_k[j], n_mech, SunCoef) # kW
314# elif SelectedModel == 'Medina-Flores et al. (2010)':
315# w_k[j] = Work_MedinaModel(P_out[j - 1], m_in_k, dh_is_k[j]) # kW
316# elif SelectedModel == 'Varbanov et al. (2004)':
317# w_k[j] = Work_THM(P_out[j - 1], h_out[j - 1], P_out[j], h_sat[j], m_in_k, dh_is_k[j], n_mech, VarCoef) # kW
318# elif SelectedModel == 'Fixed Isentropic Turbine':
319# w_k[j] = w_isen_k[j] * MinEff # kW
320# else:
321# w_k[j] = 0
323# if w_isen_k[j] > 0:
324# eff_k[j] = w_k[j] / w_isen_k[j]
325# else:
326# eff_k[j] = MinEff
327# if eff_k[j] <= MinEff:
328# w_k[j] = MinEff * w_isen_k[j] # kW
329# if m_in_k > 0:
330# h_out[j] = h_out[j - 1] - w_k[j] / (m_in_k * n_mech)
331# else:
332# h_out[j] = h_out[j - 1] # kJ/kg
334# # Correct for high pressure condensate flash, if selected
335# if m_in_k > 0:
336# if config.CONDESATE_FLASH_CORRECTION:
337# Q_flash = m_k[j - 1] * (h_tar[j - 1] - h_tar[j])
338# m_k[j] = (Q_users[j] - Q_flash) / (h_out[j] - h_tar[j])
339# else:
340# m_k[j] = Q_users[j] / (h_out[j] - h_tar[j])
341# else:
342# m_k[j] = 0 # kg/s
344# m_in_est += m_k[j] # kg/s
346# i += 1
347# if (abs(m_in - m_in_est) < ZERO and i >= 100) or i >= 3:
348# break
350# w = 0
351# Wmax = 0
352# for j in range(s):
353# w += w_k[j]
354# Wmax += w_isen_k[j]
356# z.work_target = w # kW
357# if Wmax > 0:
358# z.turbine_efficiency_target = w / Wmax
359# else:
360# z.turbine_efficiency_target = 0 # %
362# return z
364# #######################################################################################################
365# # Helper Functions
366# #######################################################################################################
368def Work_MedinaModel(P_in, m, dh_is):
369 """'Determines power generation based on the thermodynamic model of Medina-Flores & Picón-Núñez (2010).
370 """
371 # Reference for Turbine Model 1: Medina-Flores & Picón-Núñez (2010)
372 # Modelling the power production of single and multiple extraction steam turbines.
373 # Chemical Engineering Science 65, 2811-2820
374 # Part load not used
375 A0 = 185.4 + 43.3 * (P_in * 0.1) # a0 in kW, P in bar
376 b0 = 1.2057 + 0.0075 * (P_in * 0.1) # b0 in dimensionless, P in bar
377 return (m * dh_is - A0) / b0 # kW
379def Work_SunModel(P_in, h_in, P_out, h_sat, m, m_max, dh_is, n_mech, t_type=1):
380 """Determines power generation based on the correlation model of Sun and Smith (2015).
381 """
382 # 'Reference for Turbine Model 2: Sun & Smith (2015)
383 # ' Performance Modeling of New and Existing Steam Turbines
384 # ' I&CE 54, 1908-1915
386 coeff = {
387 "BPST": {
388 "a": [1.18795366, -0.00029564, 0.004647288],
389 "b": [449.9767142, 5.670176939, -11.5045814],
390 "c": [0.205149333, -0.000695171, 0.002844611],
391 },
392 "CT": {
393 "a": [1.314991261, -0.001634725, -0.367975103],
394 "b": [-437.7746025, 29.00736723, 10.35902331],
395 "c": [0.07886297, 0.000528327, -0.703153891],
396 }
397 }
399 # Model coefficients where P in bar
400 A0 = coeff[t_type]["a"][0] + coeff[t_type]["a"][1] * (P_in) + coeff[t_type]["a"][2] * (P_out)
401 b0 = coeff[t_type]["b"][0] + coeff[t_type]["b"][1] * (P_in) + coeff[t_type]["b"][2] * (P_out)
402 c0 = coeff[t_type]["c"][0] + coeff[t_type]["c"][1] * (P_in) + coeff[t_type]["c"][2] * (P_out)
404 # Willan's line coefficients and predicted work (after isentropic and mechanical efficiency loss)
405 W_int = c0 / A0 * (m_max * dh_is - b0)
406 n = (1 + c0) / A0 * (dh_is - b0 / m_max)
407 w_act = n * m - W_int
408 h_out = h_in - w_act / (n_mech * m)
410 if h_out <= h_sat + ZERO and t_type == 1:
411 w_act = Work_SunModel(P_in, h_in, P_out, h_sat, m, m_max, dh_is, n_mech, 2)
412 return w_act
414def Work_THM(P_in, h_in, P_out, h_sat, m, dh_is, n_mech, t_size=1, t_type=1):
415 """Determines power generation based on the Turbine Hardware Model of Varbanov et al. (2004).
416 """
417 # 'Reference for Turbine Model 3: Varbanov et al. (2004)
418 # ' Modelling and Optimization of Utility Systems
419 # ' Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 561–578
420 # Part load not used
421 coeff = {
422 "BPST": {
423 "<2MW": [0, 0.00108, 1.097, 0.00172],
424 ">2MW": [0, 0.00423, 1.155, 0.000538],
425 },
426 "CT": {
427 "<2MW": [0, 0.000662, 1.191, 0.000759],
428 ">2MW": [-0.463, 0.00353, 1.22, 0.000148],
429 }
430 }
432 dT_sat = Tsat_p(P_in) - Tsat_p(P_out)
433 a = (coeff[t_type][t_size][0] + coeff[t_type][t_size][1] * dT_sat) * 1000
434 b = coeff[t_type][t_size][2] + coeff[t_type][t_size][3] * dT_sat
435 w_max = (dh_is * m - a) / b
437 if w_max > 2000 and t_size == 1:
438 t_size = 2
439 w_max = Work_THM(P_in, h_in, P_out, h_sat, m, dh_is, n_mech, t_size)
441 h_out = h_in - w_max / (n_mech * m)
442 if h_out <= h_sat + ZERO and t_type == 1:
443 t_type = 2
444 w_max = Work_THM(P_in, h_in, P_out, h_sat, m, dh_is, n_mech, t_size, t_type)
446 return w_max
450def Set_Coeff(SunCoef=None, VarCoef=None):
451 """Sets the model coefficients."""
454 if VarCoef != None:
455 # Model coefficients for the Varbanov et al. model
456 # BPST
457 # < 2MW
458 VarCoef[0][0][0] = 0
459 VarCoef[0][0][1] = 0.00108
460 VarCoef[0][0][2] = 1.097
461 VarCoef[0][0][3] = 0.00172
463 # > 2MW
464 VarCoef[0][1][0] = 0
465 VarCoef[0][1][1] = 0.00423
466 VarCoef[0][1][2] = 1.155
467 VarCoef[0][1][3] = 0.000538
469 # CT
470 # < 2MW
471 VarCoef[1][0][0] = 0
472 VarCoef[1][0][1] = 0.000662
473 VarCoef[1][0][2] = 1.191
474 VarCoef[1][0][3] = 0.000759
476 # > 2MW
477 VarCoef[1][1][0] = -0.463
478 VarCoef[1][1][1] = 0.00353
479 VarCoef[1][1][2] = 1.22
480 VarCoef[1][1][3] = 0.000148