1"""Functions to optimise hydrogen storage locations.
2
3References
4----------
5.. [#Musial19] Musial, W. D., Beiter, P. C., Nunemaker, J., Heimiller, D. M.,
6 Ahmann, J., and Busch, J. (2019). Oregon Offshore Wind Site Feasibility
7 and Cost Study. Technical Report NREL/TP-5000-74597. Golden, CO: National
8 Renewable Energy Laboratory. https://doi.org/10.2172/1570430.
9.. [#Dinh23] Dinh, Q. V., Dinh, V. N., Mosadeghi, H., Todesco Pereira, P. H.,
10 and Leahy, P. G. (2023). ‘A geospatial method for estimating the
11 levelised cost of hydrogen production from offshore wind’,
12 International Journal of Hydrogen Energy, 48(40), pp. 15000–15013.
13 https://doi.org/10.1016/j.ijhydene.2023.01.016.
14.. [#Pryor18] Pryor, S. C., Shepherd, T. J., and Barthelmie, R. J. (2018).
15 ‘Interannual variability of wind climates and wind turbine annual energy
16 production’, Wind Energy Science, 3(2), pp. 651–665.
17 https://doi.org/10.5194/wes-3-651-2018.
18.. [#Dinh21] Dinh, V. N., Leahy, P., McKeogh, E., Murphy, J., and Cummins, V.
19 (2021). ‘Development of a viability assessment model for hydrogen
20 production from dedicated offshore wind farms’, International Journal of
21 Hydrogen Energy, 46(48), pp. 24620–24631.
22 https://doi.org/10.1016/j.ijhydene.2020.04.232.
23.. [#HydrogenTools] Hydrogen Tools (no date) Lower and Higher Heating Values
24 of Fuels. Available at:
25 https://h2tools.org/hyarc/calculator-tools/lower-and-higher-heating-values-fuels
26 (Accessed: 25 April 2024).
27.. [#Dinh24a] Dinh, Q. V., Dinh, V. N., and Leahy, P. G. (2024). ‘A
28 differential evolution model for optimising the size and cost of
29 electrolysers coupled with offshore wind farms’.
30.. [#Dinh24b] Dinh, Q. V., Todesco Pereira, P. H., Dinh, V. N., Nagle, A. J.,
31 and Leahy, P. G. (2024). ‘Optimising the levelised cost of transmission
32 for green hydrogen and ammonia in new-build offshore energy
33 infrastructure: pipelines, tankers, and HVDC’.
34.. [#Baufume13] Baufumé, S., Grüger, F., Grube, T., Krieg, D., Linssen, J.,
35 Weber, M., Hake, J.-F., and Stolten, D. (2013). ‘GIS-based scenario
36 calculations for a nationwide German hydrogen pipeline infrastructure’,
37 International Journal of Hydrogen Energy, 38(10), pp. 3813–3829.
38 https://doi.org/10.1016/j.ijhydene.2012.12.147.
39.. [#IRENA22] International Renewable Energy Agency (2022). Global Hydrogen
40 Trade to Meet the 1.5°C Climate Goal: Technology Review of Hydrogen
41 Carriers. Abu Dhabi: International Renewable Energy Agency.
42 ISBN 978-92-9260-431-8. Available at:
43 https://www.irena.org/publications/2022/Apr/Global-hydrogen-trade-Part-II
44 (Accessed: 31 December 2023).
45.. [#Siemens] Siemens Energy (n.d.). Silyzer 300. Datasheet. Erlangen.
46 Available at:
47 https://assets.siemens-energy.com/siemens/assets/api/uuid:a193b68f-7ab4-4536-abe2-c23e01d0b526/datasheet-silyzer300.pdf
48 (Accessed: 31 December 2023).
49.. [#IEA19] International Energy Agency (2019). The Future of Hydrogen:
50 Seizing today’s opportunities. Paris: International Energy Agency.
51 Available at: https://www.iea.org/reports/the-future-of-hydrogen
52 (Accessed: 14 August 2023).
53"""
54
55from functools import partial
56
57import geopandas as gpd
58import numpy as np
59import pandas as pd
60from pyfluids import Fluid, FluidsList, Input
61from scipy import integrate
62from shapely.geometry import Point
63
64from h2ss import data as rd
65
66# NREL 15 MW reference turbine specifications
67REF_DIAMETER = 248 # m
68REF_HUB_HEIGHT = 149 # m
69REF_RATED_POWER = 15 # MW
70REF_V_CUT_IN = 3 # m s-1
71REF_V_RATED = 11 # m s-1
72REF_V_CUT_OUT = 25 # m s-1
73
74# power curve data
75pc = pd.DataFrame(
76 {
77 "wind_speed": range(REF_V_RATED + 2),
78 "power_curve": (
79 [0] * 4
80 + [0.499, 1.424, 2.732, 4.469, 6.643, 9.459, 12.975]
81 + [REF_RATED_POWER] * 2
82 ),
83 }
84)
85# for wind speeds between cut-in and rated
86pc["gradient"] = pc["power_curve"] - pc["power_curve"].shift(1)
87pc["intercept"] = pc["power_curve"] - pc["gradient"] * pc["wind_speed"]
88pc = pc[
89 (pc["wind_speed"] > REF_V_CUT_IN) & (pc["wind_speed"] < REF_V_RATED + 1)
90]
91
92
[docs]
93def ref_power_curve(v):
94 """Power curve for the reference wind turbine.
95
96 Parameters
97 ----------
98 v : float
99 Wind speed [m s⁻¹]
100
101 Returns
102 -------
103 float
104 Wind turbine power output at a given wind speed from the power curve
105 [MW]
106
107 Notes
108 -----
109 The NREL 15 MW reference wind turbine is used; see [#Musial19]_, Appendix
110 D, p. 84, for the data used to generate the power curve.
111 """
112 if v < REF_V_CUT_IN:
113 power_curve = 0
114 elif REF_V_CUT_IN <= v < REF_V_RATED:
115 pc_vals = pc[pc["wind_speed"] == np.trunc(v) + 1]
116 power_curve = (pc_vals["gradient"] * v + pc_vals["intercept"]).iloc[0]
117 elif REF_V_RATED <= v <= REF_V_CUT_OUT:
118 power_curve = REF_RATED_POWER
119 elif v > REF_V_CUT_OUT:
120 power_curve = 0
121 return power_curve
122
123
[docs]
124def weibull_probability_distribution(v, k, c):
125 """Weibull probability distribution function.
126
127 Parameters
128 ----------
129 v : float
130 Wind speed [m s⁻¹]
131 k : float
132 Shape (Weibull distribution parameter, k)
133 c : float
134 Scale (Weibull distribution parameter, C) [m s⁻¹]
135
136 Returns
137 -------
138 float
139 Weibull probability distribution function [s m⁻¹]
140
141 Notes
142 -----
143 The Weibull probability distribution function, :math:`f(v)` [s m⁻¹] is
144 based on Eqn. (1) of [#Dinh23]_, where :math:`k` and :math:`C`
145 [m s⁻¹] are the shape and scale Weibull distribution parameters,
146 respectively, and :math:`v` is the wind speed.
147
148 .. math::
149 f(v) = \\frac{k}{C} \\, \\left( \\frac{v}{C} \\right)^{k - 1}
150 \\, \\exp \\left( -\\left( \\frac{v}{C} \\right)^k \\right)
151
152 See also [#Pryor18]_, Eqn. (2).
153 """
154 return k / c * np.power((v / c), (k - 1)) * np.exp(-np.power((v / c), k))
155
156
[docs]
157def weibull_distribution(weibull_wf_data):
158 """Generate a power curve and Weibull distribution.
159
160 Parameters
161 ----------
162 weibull_wf_data : pandas.DataFrame
163 Dataframe of the Weibull distribution parameters for the wind farms
164
165 Returns
166 -------
167 pandas.DataFrame
168 Dataframe of the Weibull distribution for each wind farm for wind
169 speeds of 0 to 30 m s⁻¹ at an interval of 0.01
170 """
171 powercurve_weibull_data = {}
172 for n in weibull_wf_data["name"]:
173 powercurve_weibull_data[n] = {}
174 powercurve_weibull_data[n]["wind_speed"] = [
175 0 + 0.01 * n for n in range(3001)
176 ]
177 powercurve_weibull_data[n]["power_curve"] = []
178 powercurve_weibull_data[n][n] = []
179 for v in powercurve_weibull_data[n]["wind_speed"]:
180 powercurve_weibull_data[n]["power_curve"].append(
181 ref_power_curve(v=v)
182 )
183 powercurve_weibull_data[n][n].append(
184 weibull_probability_distribution(
185 v=v,
186 k=weibull_wf_data[weibull_wf_data["name"] == n][
187 ("k", "mean")
188 ].iloc[0],
189 c=weibull_wf_data[weibull_wf_data["name"] == n][
190 ("c", "mean")
191 ].iloc[0],
192 )
193 )
194 powercurve_weibull_data[n] = pd.DataFrame(powercurve_weibull_data[n])
195 powercurve_weibull_data = (
196 pd.concat(powercurve_weibull_data.values(), axis=1)
197 .T.drop_duplicates()
198 .T
199 )
200 return powercurve_weibull_data
201
202
[docs]
203def weibull_power_curve(v, k, c):
204 """Weibull probability distribution function multiplied by the power curve.
205
206 Parameters
207 ----------
208 v : float
209 Wind speed between cut-in and cut-out [m s⁻¹]
210 k : float
211 Shape (Weibull distribution parameter, k)
212 c : float
213 Scale (Weibull distribution parameter, C) [m s⁻¹]
214
215 Returns
216 -------
217 float
218 Weibull probability distribution multiplied by the power curve
219 [MW s m⁻¹]
220 """
221 if REF_V_CUT_IN <= v < REF_V_RATED:
222 pc_vals = pc[pc["wind_speed"] == np.trunc(v) + 1]
223 power_curve = (pc_vals["gradient"] * v + pc_vals["intercept"]).iloc[
224 0
225 ] * weibull_probability_distribution(v=v, k=k, c=c)
226 elif REF_V_RATED <= v <= REF_V_CUT_OUT:
227 power_curve = REF_RATED_POWER * weibull_probability_distribution(
228 v=v, k=k, c=c
229 )
230 return power_curve
231
232
[docs]
233def number_of_turbines(owf_cap, wt_power=REF_RATED_POWER):
234 """Number of reference wind turbines in the offshore wind farm.
235
236 Parameters
237 ----------
238 owf_cap : float
239 Maximum nameplate capacity of the proposed offshore wind farm [MW]
240 wt_power : float
241 Rated power of the reference wind turbine [MW]
242
243 Returns
244 -------
245 int
246 Number of wind turbines in the offshore wind farm comprising of
247 reference wind turbines
248
249 Notes
250 -----
251 The number of turbines, :math:`n` of an offshore wind farm was
252 determined using the floor division of the wind farm's maximum nameplate
253 capacity, :math:`P_{owf}` [MW] by the reference wind turbine's rated power,
254 :math:`P_{rated}` [MW].
255
256 .. math::
257 n = \\left\\lfloor \\frac{P_{owf}}{P_{rated}} \\right\\rfloor
258 """
259 return (owf_cap / wt_power).astype(int)
260
261
[docs]
262def annual_energy_production_function(
263 n_turbines, k, c, w_loss=0.1, acdc_loss=0.982
264):
265 """Annual energy production of the wind farm.
266
267 Parameters
268 ----------
269 n_turbines : int
270 Number of wind turbines in wind farm
271 k : float
272 Shape (Weibull distribution parameter)
273 c : float
274 Scale (Weibull distribution parameter) [m s⁻¹]
275 w_loss : float
276 Wake loss
277 acdc_loss : float
278 AC-DC conversion losses
279
280 Returns
281 -------
282 tuple[float, float, float]
283 Annual energy production of wind farm [MWh], integral [MW], and
284 absolute error [MW]
285
286 Notes
287 -----
288 The annual energy production, :math:`E_{annual}` [MWh], is based on Eqn.
289 (3) of [#Dinh23]_, where :math:`n` is the number of turbines in
290 the wind farm, :math:`w` is the wake loss, which is assumed to be a
291 constant value of 0.1, :math:`v_i` and :math:`v_o` [m s⁻¹] are the cut-in
292 and cut-out speeds of the wind turbine, respectively, :math:`P(v)` [MW] is
293 the wind turbine power output, :math:`f(v)` [s m⁻¹] is the Weibull
294 probability distribution function, and :math:`\\varepsilon` is the
295 AC-DC conversion loss.
296
297 .. math::
298 E_{annual} = 365 \\times 24 \\times n \\,
299 \\left( 1 - w \\right) \\, \\varepsilon \\,
300 \\int\\limits_{v_i}^{v_o} P(v) \\, f(v) \\,\\mathrm{d}v
301
302 In the function's implementation, both the limit and absolute error
303 tolerance for the integration have been increased.
304 """
305 integration = integrate.quad(
306 partial(weibull_power_curve, k=k, c=c),
307 a=REF_V_CUT_IN,
308 b=REF_V_CUT_OUT,
309 limit=100, # an upper bound on the number of subintervals used
310 epsabs=1.49e-6, # absolute error tolerance
311 )
312 aep = 365 * 24 * n_turbines * (1 - w_loss) * acdc_loss * integration[0]
313 return aep, integration[0], integration[1]
314
315
[docs]
316def annual_energy_production(weibull_wf_data):
317 """Annual energy production of the wind farms.
318
319 Parameters
320 ----------
321 weibull_wf_data : pandas.DataFrame
322 Dataframe of the Weibull distribution parameters for the wind farms
323
324 Returns
325 -------
326 pandas.DataFrame
327 Dataframe with the annual energy production for each wind farm
328 """
329 aep = []
330 for n in weibull_wf_data["name"]:
331 aep.append(
332 annual_energy_production_function(
333 n_turbines=weibull_wf_data[weibull_wf_data["name"] == n][
334 "n_turbines"
335 ].iloc[0],
336 k=weibull_wf_data[weibull_wf_data["name"] == n][
337 ("k", "mean")
338 ].iloc[0],
339 c=weibull_wf_data[weibull_wf_data["name"] == n][
340 ("c", "mean")
341 ].iloc[0],
342 )
343 )
344 aep = pd.DataFrame(aep)
345 aep.columns = ["AEP", "integral", "abserr"]
346 weibull_wf_data = pd.concat([weibull_wf_data, aep], axis=1)
347 return weibull_wf_data
348
349
[docs]
350def annual_hydrogen_production(aep, eta_conv=0.7, e_pcl=0.003053):
351 """Annual hydrogen production from the wind farm's energy generation.
352
353 Parameters
354 ----------
355 aep : float
356 Annual energy production of wind farm [MWh]
357 eta_conv : float
358 Conversion efficiency of the electrolyser
359 e_pcl : float
360 Electricity consumed by other parts of the hydrogen plant [MWh kg⁻¹]
361
362 Returns
363 -------
364 float
365 Annual hydrogen production [kg]
366
367 Notes
368 -----
369 Eqn. (4) of [#Dinh23]_, Eqn. (9) of [#Dinh21]_, and [#Dinh24a]_.
370 The value for the electricity consumed by other parts of the hydrogen
371 plant is from Table 3 of [#Dinh21]_ for proton exchange membrane (PEM)
372 electrolysers predicted for the year 2030.
373 This includes energy for water purification, hydrogen compression, and
374 losses.
375 The electrolyser conversion efficiency is based on [#Dinh24a]_.
376 See [#HydrogenTools]_ for the heating values.
377
378 .. math::
379 m_{annual} = \\dfrac{E_{annual}}{\\dfrac{LHV}{3,600 \\, \\eta} +
380 E_{plant}}
381
382 where :math:`m_{annual}` is the annual hydrogen production [kg],
383 :math:`E_{annual}` is the annual energy production of the wind farm [MWh],
384 :math:`LHV` is the lower heating value of hydrogen [MJ kg⁻¹], :math:`\\eta`
385 is the conversion efficiency of the electrolyser, and :math:`E_{plant}` is
386 the electricity consumed by other parts of the hydrogen plant [MWh kg⁻¹].
387 """
388 return aep / (119.96 / (3600 * eta_conv) + e_pcl)
389
390
[docs]
391def transmission_distance(
392 cavern_df, wf_data, injection_point_coords=(-6, -12, 53, 21)
393):
394 """Calculate the transmission distance to the injection point.
395
396 Parameters
397 ----------
398 cavern_df : geopandas.GeoDataFrame
399 Dataframe of potential caverns
400 wf_data : geopandas.GeoDataFrame
401 Geodataframe of the offshore wind farm data
402 injection_point_coords : tuple[float, float, float, float]
403 Injection point coordinates (lon-deg, lon-min, lat-deg, lat-min)
404
405 Returns
406 -------
407 tuple[geopandas.GeoDataFrame, geopandas.GeoSeries]
408 The cavern dataframe and injection point
409
410 Notes
411 -----
412 A hacky workaround was used to prevent
413 "DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is
414 deprecated, and will error in future. Ensure you extract a single
415 element from your array before performing this operation. (Deprecated
416 NumPy 1.25.)"
417 during the transformation of the injection point from a single-row series,
418 i.e. assigning a duplicate row and dropping it after reprojecting.
419 """
420 lond, lonm, latd, latm = injection_point_coords
421 injection_point = (
422 gpd.GeoSeries(
423 [
424 Point(lond + lonm / 60, latd + latm / 60),
425 Point(lond + lonm / 60, latd + latm / 60),
426 ],
427 crs=4326,
428 )
429 .to_crs(rd.CRS)
430 .drop_duplicates()
431 )
432 distance_ip = []
433 for j in range(len(cavern_df)):
434 distance_ip.append(
435 injection_point.distance(
436 cavern_df.iloc[[j]]["geometry"], align=False
437 ).values[0]
438 / 1000
439 )
440 cavern_df["distance_ip"] = distance_ip
441 distance_wf = {}
442 for i in range(len(wf_data)):
443 distance_wf[wf_data["name"][i]] = []
444 for j in range(len(cavern_df)):
445 distance_wf[wf_data["name"][i]].append(
446 (
447 wf_data.iloc[[i]]
448 .distance(cavern_df.iloc[[j]]["geometry"], align=False)
449 .values[0]
450 / 1000
451 + cavern_df.iloc[[j]]["distance_ip"].values[0]
452 )
453 )
454 cavern_df[f"dist_{wf_data['name'][i].replace(' ', '_')}"] = (
455 distance_wf[wf_data["name"][i]]
456 )
457 return cavern_df, injection_point
458
459
[docs]
460def electrolyser_capacity(
461 n_turbines, wt_power=REF_RATED_POWER, cap_ratio=0.83
462):
463 """
464 Calculate the electrolyser capacity for an offshore wind farm.
465
466 Parameters
467 ----------
468 n_turbines : int
469 Number of wind turbines in wind farm
470 wt_power : float
471 Rated power of the reference wind turbine [MW]
472 cap_ratio : float
473 Ratio of electrolyser capacity to the wind farm capacity
474
475 Returns
476 -------
477 int
478 Electrolyser capacity [MW]
479
480 Notes
481 -----
482 In the offshore electrolyser concept of [#Dinh24a]_, the electricity from
483 an offshore wind farm "is supplied to the electrolyser by medium voltage
484 alternating current (AC) cable. The electrolyser requires direct current
485 (DC) electricity input so this concept includes an AC-DC converter. The
486 hydrogen obtained is compressed and delivered onshore by hydrogen
487 pipelines".
488
489 The "energy from the offshore wind farm is only subjected to wake losses
490 and a small amount of electricity loss in the site network and the AC-DC
491 converter before being supplied to the electrolyser. The hydrogen obtained
492 will then be transported to onshore infrastructure by pipeline. The losses
493 incurred by pipeline transmission are very low" [#Dinh24a]_.
494
495 The electrolyser capacity, :math:`P_{electrolyser}` [MW] is rounded down
496 to an integer and is the product of the number of reference wind turbines
497 of the offshore wind farm, :math:`n`, the rated power of the
498 reference wind turbine, :math:`P_{rated}`, and the ratio of the
499 electrolyser capacity to the offshore wind farm capacity,
500 :math:`F_{electrolyser}`.
501
502 .. math::
503 P_{electrolyser} = \\left\\lfloor n \\, P_{rated} \\,
504 F_{electrolyser} \\right\\rfloor
505 """
506 return (n_turbines * wt_power * cap_ratio).astype(int)
507
508
[docs]
509def hydrogen_pipeline_density(pressure, temperature):
510 """Calculate the density of hydrogen in the pipelines.
511
512 Parameters
513 ----------
514 pressure : float
515 Pressure of hydrogen in the pipeline [Pa]
516 temperature : float
517 Temperature of hydrogen in the pipeline [°C]
518
519 Returns
520 -------
521 float
522 Pipeline hydrogen density [kg m⁻³]
523
524 Notes
525 -----
526 According to [#Baufume13]_, the envisaged future hydrogen pipeline system
527 is assumed to be operated at pressure levels up to 10 MPa.
528 They used an average operating pressure of 6.5 MPa for transmission
529 pipelines and calculated the density using the average soil temperature,
530 which was assumed to be the operating temperature of the hydrogen
531 [#Baufume13]_.
532 """
533 rho = (
534 Fluid(FluidsList.Hydrogen)
535 .with_state(
536 Input.pressure(pressure), Input.temperature(temperature + 273.15)
537 )
538 .density
539 )
540 # print(
541 # f"Density of hydrogen in the pipelines: {rho:.4f} "
542 # f"kg m\N{SUPERSCRIPT MINUS}\N{SUPERSCRIPT THREE}"
543 # )
544 return rho
545
546
[docs]
547def capex_pipeline(e_cap, p_rate=0.0055, pressure=100e5, temperature=10, u=15):
548 """Capital expenditure (CAPEX) for the pipeline.
549
550 Parameters
551 ----------
552 e_cap : float
553 Electrolyser capacity [MW]
554 p_rate : float
555 Electrolyser production rate [kg s⁻¹ MW⁻¹]
556 pressure : float
557 Pressure of hydrogen in the pipeline [Pa]
558 temperature : float
559 Temperature of hydrogen in the pipeline [°C]
560 u : float
561 Average fluid velocity [m s⁻¹]
562
563 Returns
564 -------
565 float
566 CAPEX of the pipeline per km of pipeline [€ km⁻¹]
567
568 Notes
569 -----
570 See Eqn. (18) of [#Dinh23]_ and Section 3.1 of [#Dinh24b]_, from which
571 the following text has been taken.
572
573 The estimation of offshore pipeline costs is based on the onshore pipeline
574 calculations of [#Baufume13]_, multiplied by a factor of two to reflect the
575 estimated cost scaling of onshore to expected offshore costs, as suggested
576 by [#IRENA22]_.
577
578 In the reference case, the resulting capital expenditure for transmission
579 pipelines and associated recompression stations was represented by a
580 second order polynomial function.
581
582 The electrolyser production rate was based on the Siemens - Silyzer 300
583 [#Siemens]_.
584
585 Because the electrolyser plant is assumed to be operating at full capacity
586 at all times, the CAPEX was calculated considering a 75% utilisation rate,
587 i.e. the pipeline capacity is 33% oversized [#IEA19]_.
588
589 .. math::
590 CAPEX = 2,000 \\, \\left( 16,000 \\, \\frac{P_{electrolyser}
591 \\times EPR}{\\rho_{H_2} \\, v_{H_2} \\, \\pi} + 1,197.2 \\,
592 \\sqrt{\\frac{P_{electrolyser} \\times EPR}{\\rho_{H_2} \\, v_{H_2}
593 \\, \\pi}} + 329 \\right)
594
595 where :math:`CAPEX` is the CAPEX of the pipeline per km of pipeline
596 [€ km⁻¹], :math:`P_{electrolyser}` is the electrolyser capacity [MW],
597 :math:`EPR` is the electrolyser production rate [kg s⁻¹ MW⁻¹],
598 :math:`\\rho_{H_2}` is the density of hydrogen [kg m⁻³], and
599 :math:`v_{H_2}` is the average fluid velocity [m s⁻¹].
600 """
601 rho = hydrogen_pipeline_density(pressure=pressure, temperature=temperature)
602 f = e_cap * p_rate / (rho * u * np.pi)
603 return 2000 * (16000 * f + 1197.2 * np.sqrt(f) + 329)
604
605
[docs]
606def lcot_pipeline_function(
607 capex,
608 d_transmission,
609 ahp,
610 opex_ratio=0.02,
611 discount_rate=0.08,
612 lifetime=30,
613):
614 """Levelised cost of transmission (LCOT) of hydrogen in pipelines.
615
616 Parameters
617 ----------
618 capex : float
619 Capital expenditure (CAPEX) of the pipeline [€ km⁻¹]
620 d_transmission : float
621 Pipeline transmission distance [km]
622 ahp : float
623 Annual hydrogen production [kg]
624 opex_ratio : float
625 Ratio of the operational expenditure (OPEX) to the CAPEX
626 discount_rate : float
627 Discount rate
628 lifetime : int
629 Lifetime of the pipeline [year]
630
631 Returns
632 -------
633 float
634 LCOT of hydrogen in the pipeline [€ kg⁻¹]
635
636 Notes
637 -----
638 The OPEX is calculated as a percentage of the CAPEX.
639 See the introduction of Section 3, Eqn. (1) and (2), and Section 3.5, Eqn.
640 (22) of [#Dinh24b]_; see Tables 2 and 3 for the assumptions and constants
641 used.
642
643 .. math::
644 LCOT = \\frac{\\mathrm{total\\ lifecycle\\ costs\\ of\\ all\\
645 components}}
646 {\\mathrm{lifetime\\ hydrogen\\ transported}}
647 .. math::
648 OPEX = CAPEX \\times F_{OPEX}
649 .. math::
650 LCOT = \\frac{\\left( CAPEX + \\displaystyle\\sum_{l=0}^{L}
651 \\frac{OPEX}{{(1 + DR)}^l} \\right) \\, d}
652 {\\displaystyle\\sum_{l=0}^{L} \\frac{AHP}{{(1 + DR)}^l}}
653
654 where :math:`LCOT` is the LCOT of hydrogen in pipelines [€ kg⁻¹],
655 :math:`CAPEX` is the CAPEX of the pipeline per km of pipeline [€ km⁻¹],
656 :math:`d` is the pipeline transmission distance [km], :math:`OPEX` is the
657 OPEX of the pipeline [€ km⁻¹], :math:`DR` is the discount rate,
658 :math:`AHP` is the annual hydrogen production [kg], :math:`L` is the
659 lifetime of the pipeline [year], and :math:`F_{OPEX}` is the ratio of the
660 OPEX to the CAPEX.
661 """
662 f = sum(
663 1 / np.power((1 + discount_rate), year) for year in range(lifetime + 1)
664 )
665 opex = capex * opex_ratio
666 return (capex + opex * f) * d_transmission / (ahp * f)
667
668
[docs]
669def lcot_pipeline(weibull_wf_data, cavern_df):
670 """Calculate the pipeline levelised cost of transmission.
671
672 Parameters
673 ----------
674 cavern_df : geopandas.GeoDataFrame
675 Dataframe of potential caverns
676 weibull_wf_data : pandas.DataFrame
677 Dataframe of the Weibull distribution parameters for the wind farms
678
679 Returns
680 -------
681 pandas.DataFrame
682 Dataframe of potential caverns
683 """
684 for wf in list(weibull_wf_data["name"]):
685 cavern_df[f"LCOT_{wf.replace(' ', '_')}"] = lcot_pipeline_function(
686 capex=weibull_wf_data[weibull_wf_data["name"].str.contains(wf)][
687 "CAPEX"
688 ].values[0],
689 d_transmission=cavern_df[f"dist_{wf.replace(' ', '_')}"],
690 ahp=weibull_wf_data[weibull_wf_data["name"].str.contains(wf)][
691 "AHP"
692 ].values[0],
693 )
694 return cavern_df
695
696
[docs]
697def rotor_area(diameter=REF_DIAMETER):
698 """Wind turbine rotor swept area.
699
700 Parameters
701 ----------
702 diameter : float
703 Wind turbine rotor diameter [m]
704
705 Returns
706 -------
707 float
708 Wind turbine rotor swept area [m²]
709
710 Notes
711 -----
712 .. math::
713 A = \\frac{\\pi \\, D^2}{4}
714
715 where :math:`A` is the wind turbine rotor swept area [m²] and :math:`D` is
716 the rotor diameter [m].
717 """
718 return np.pi * np.square(diameter) / 4
719
720
[docs]
721def power_wind_resource(v, rho=1.225, diameter=REF_DIAMETER):
722 """Total wind resource power passing through a wind turbine's rotor.
723
724 Parameters
725 ----------
726 v : float
727 Wind speed [m s⁻¹]
728 rho : float
729 Air density [kg m⁻³]
730 diameter : float
731 Wind turbine rotor diameter [m]
732
733 Returns
734 -------
735 float
736 Power contained in the wind resource [MW]
737
738 Notes
739 -----
740 .. math::
741 P_{wind} = \\frac{\\rho_{air} \\, A \\, v^3}{2 \\times 10^6}
742
743 where :math:`P_{wind}` is the power contained in the wind resource [MW],
744 :math:`\\rho_{air}` is the air density [kg m⁻³], :math:`A` is the rotor
745 swept area [m²], and :math:`v` is the wind speed [m s⁻¹].
746 """
747 return 0.5 * rho * rotor_area(diameter=diameter) * np.power(v, 3) / 1e6
748
749
[docs]
750def power_coefficient(v):
751 """Power coefficient of the reference wind turbine.
752
753 Parameters
754 ----------
755 v : float
756 Wind speed [m s⁻¹]
757
758 Returns
759 -------
760 float
761 Power coefficient
762
763 Notes
764 -----
765 This is specific to the power curve of the NREL 15 MW reference wind
766 turbine [#Musial19]_.
767
768 .. math::
769 C_p = \\frac{P}{P_{wind}} = \\frac{P \\times 2 \\times 10^6}
770 {\\rho_{air} \\, A \\, v^3}
771
772 where :math:`P` is the wind turbine power output [MW], :math:`P_{wind}` is
773 the power contained in the wind resource [MW], :math:`\\rho_{air}` is the
774 air density [kg m⁻³], :math:`A` is the rotor swept area [m²], :math:`v` is
775 the wind speed [m s⁻¹], and :math:`C_p` is the power coefficient of the
776 wind turbine.
777 """
778 if REF_V_CUT_IN <= v <= REF_V_CUT_OUT:
779 coeff = ref_power_curve(v=v) / power_wind_resource(v=v)
780 else:
781 coeff = 0
782 return coeff