Merge pull request #348 from kbrunik/h2new

PEM BOP Functionality

greenheart/simulation/greenheart_simulation.py

@@ -587,13 +587,14 @@ def energy_internals(
         )
         total_peripheral_energy = power_for_peripherals_kw_in * 365 * 24
         distributed_peripheral_power = total_peripheral_energy / high_count
-        for i in range(len(hopp_results["combined_hybrid_power_production_hopp"])):
-            r = (
-                hopp_results["combined_hybrid_power_production_hopp"][i]
-                - distributed_peripheral_power
-            )
-            if r > 0:
-                remaining_power_profile_in[i] = r
+
+        remaining_power_profile_in = np.where(
+            hopp_results["combined_hybrid_power_production_hopp"]
+                - distributed_peripheral_power > 0,
+            hopp_results["combined_hybrid_power_production_hopp"]
+                - distributed_peripheral_power,
+            0
+        )
 
         hopp_results_internal["combined_hybrid_power_production_hopp"] = tuple(
             remaining_power_profile_in
@@ -620,6 +621,12 @@ def energy_internals(
             verbose=verbose,
         )
 
+        # run electrolyzer bop model
+        electrolyzer_energy_consumption_bop_kw = he_elec.run_electrolyzer_bop(
+            greenheart_config,
+            electrolyzer_physics_results
+        )
+
         desal_results = he_elec.run_desal(
             hopp_config, electrolyzer_physics_results, design_scenario, verbose
         )
@@ -685,16 +692,18 @@ def energy_internals(
         h2_storage_energy_kwh = h2_storage_results["storage_energy"]
         h2_storage_power_kw = h2_storage_energy_kwh * (1.0 / (365 * 24))
 
+        total_accessory_power_renewable_kw = np.zeros(len(electrolyzer_energy_consumption_bop_kw))
+        total_accessory_power_renewable_kw += electrolyzer_energy_consumption_bop_kw
         # if transport is not HVDC and h2 storage is on shore, then power the storage from the grid
         if (design_scenario["transportation"] == "pipeline") and (
             design_scenario["h2_storage_location"] == "onshore"
         ):
-            total_accessory_power_renewable_kw = (
+            total_accessory_power_renewable_kw += (
                 desal_power_kw + h2_transport_compressor_power_kw
             )
             total_accessory_power_grid_kw = h2_storage_power_kw
         else:
-            total_accessory_power_renewable_kw = (
+            total_accessory_power_renewable_kw += (
                 desal_power_kw + h2_transport_compressor_power_kw + h2_storage_power_kw
             )
             total_accessory_power_grid_kw = 0.0
@@ -706,14 +715,13 @@ def energy_internals(
         grid_power_profile = np.zeros_like(
             hopp_results["combined_hybrid_power_production_hopp"]
         )
-        for i in range(len(hopp_results["combined_hybrid_power_production_hopp"])):
-            r = (
-                hopp_results["combined_hybrid_power_production_hopp"][i]
-                - total_accessory_power_renewable_kw
-            )
-            grid_power_profile[i] = total_accessory_power_grid_kw
-            if r > 0:
-                remaining_power_profile[i] = r
+        remaining_power_profile = np.where(
+            hopp_results["combined_hybrid_power_production_hopp"]
+                - total_accessory_power_renewable_kw > 0,
+            hopp_results["combined_hybrid_power_production_hopp"]
+                - total_accessory_power_renewable_kw,
+            0
+        )
 
         if verbose and not solver:
             print("\nEnergy/Power Results:")
@@ -755,6 +763,7 @@ def energy_internals(
                     desal_power_kw,
                     h2_transport_compressor_power_kw,
                     h2_storage_power_kw,
+                    electrolyzer_energy_consumption_bop_kw,
                     remaining_power_profile,
                 )
             else:
@@ -799,6 +808,7 @@ def simple_solver(initial_guess=0.0):
             desal_power_kw,
             h2_transport_compressor_power_kw,
             h2_storage_power_kw,
+            electrolyzer_bop_kw,
             remaining_power_profile,
         ) = energy_internals(
             power_for_peripherals_kw_in=initial_guess,
@@ -813,6 +823,7 @@ def simple_solver(initial_guess=0.0):
             desal_power_kw,
             h2_transport_compressor_power_kw,
             h2_storage_power_kw,
+            electrolyzer_bop_kw
         )
 
     #################### solving for energy needed for non-electrolyzer components ####################################

greenheart/simulation/technologies/hydrogen/electrolysis/PEM_BOP/BOP_efficiency_BOL.csv

@@ -0,0 +1,87 @@
+operating_ratio,efficiency
+0.937732223,4.285714286
+0.933250382,4.258241758
+0.929665402,4.230769231
+0.925182329,4.217032967
+0.920700488,4.18956044
+0.916217415,4.175824176
+0.911735574,4.148351648
+0.904564382,4.107142857
+0.896495097,4.07967033
+0.887531415,4.024725275
+0.878566501,3.983516484
+0.870498448,3.942307692
+0.862430395,3.901098901
+0.856154832,3.873626374
+0.847189918,3.832417582
+0.837328143,3.791208791
+0.827466368,3.75
+0.818501454,3.708791209
+0.805950328,3.653846154
+0.796985414,3.612637363
+0.786226778,3.571428571
+0.775468142,3.53021978
+0.757537082,3.461538462
+0.743192234,3.406593407
+0.730641108,3.351648352
+0.717193121,3.296703297
+0.703743902,3.255494505
+0.687605332,3.200549451
+0.674156113,3.159340659
+0.661603755,3.118131868
+0.646362046,3.063186813
+0.632017198,3.008241758
+0.614979303,2.980769231
+0.588977726,2.898351648
+0.570147341,2.857142857
+0.550420096,2.815934066
+0.530694082,2.760989011
+0.513656187,2.733516484
+0.490341497,2.692307692
+0.475096092,2.678571429
+0.449089588,2.651098901
+0.415906963,2.637362637
+0.388104272,2.637362637
+0.372856404,2.651098901
+0.348638693,2.678571429
+0.329802149,2.706043956
+0.309169418,2.760989011
+0.286741734,2.82967033
+0.266102843,2.953296703
+0.245462721,3.090659341
+0.228410043,3.228021978
+0.207766225,3.406593407
+0.194302222,3.53021978
+0.181735081,3.653846154
+0.173653477,3.763736264
+0.167365594,3.873626374
+0.16107648,3.997252747
+0.153889272,4.134615385
+0.146698369,4.313186813
+0.138608141,4.519230769
+0.131417237,4.697802198
+0.125120731,4.903846154
+0.120620411,5.082417582
+0.114320209,5.32967033
+0.110717982,5.494505495
+0.107116986,5.645604396
+0.103513527,5.824175824
+0.099907604,6.03021978
+0.097199773,6.222527473
+0.093595082,6.414835165
+0.089987927,6.634615385
+0.086377076,6.895604396
+0.084566107,7.087912088
+0.081858276,7.28021978
+0.079147982,7.5
+0.078235106,7.678571429
+0.076422904,7.884615385
+0.074610703,8.090659341
+0.070994924,8.406593407
+0.069177795,8.667582418
+0.067361898,8.914835165
+0.064645444,9.203296703
+0.06372764,9.436813187
+0.061910511,9.697802198
+0.061894496,9.876373626
+0.060088454,10.01373626

greenheart/simulation/technologies/hydrogen/electrolysis/PEM_BOP/PEM_BOP.py

@@ -0,0 +1,131 @@
+import os
+import numpy as np
+import pandas as pd
+import scipy.optimize
+from greenheart.simulation.technologies.hydrogen.electrolysis.PEM_tools import get_electrolyzer_BOL_efficiency
+
+file_path = os.path.dirname(os.path.abspath(__file__))
+
+
+def calc_efficiency_curve(operating_ratio, a, b, c, d):
+    """Calculates efficiency [kWh/kg] given operation ratio with flattened end curves.
+
+    Efficiency curve and general equation structure from Wang et. al (2023). See README.md
+    in PEM_BOP directory for more details.
+
+    Wang, X.; Star, A.G.; Ahluwalia, R.K. Performance of Polymer Electrolyte Membrane Water Electrolysis Systems: 
+    Configuration, Stack Materials, Turndown and Efficiency. Energies 2023, 16, 4964. 
+    https://doi.org/10.3390/en16134964
+
+    Args:
+        operating_ratio (list or np.array): Operation ratios.
+        a (float): Coefficient a.
+        b (float): Coefficient b.
+        c (float): Coefficient c.
+        d (float): Coefficient d.
+
+    Returns:
+        efficiency (list or np.array): Efficiency of electrolyzer BOP in kWh/kg.
+    """
+    efficiency = a + b * operating_ratio + c * operating_ratio**2 + d / operating_ratio
+
+    return efficiency
+
+
+def calc_efficiency(
+    operating_ratio, efficiency, min_ratio, max_ratio, min_efficiency, max_efficiency
+):
+    """Adjust efficiency list to not go above minimum or maximum operating ratios in BOP_efficiency_BOL.csv
+
+    Args:
+        operating_ratio (list or np.array): Operation ratios.
+        efficiency (list or np.array): Efficiencies calculated using curve fit.
+        min_ratio (float): Minimum operating ratio from the CSV.
+        max_ratio (float): Maximum operating ratio from the CSV.
+        min_efficiency (float): Efficiency at the minimum operating ratio.
+        max_efficiency (float): Efficiency at the maximum operating ratio.
+
+    Returns:
+        efficiency (list or np.array): Efficiencies limited with minimum and maximum values in kWh/kg.
+    """
+    efficiency = np.where(operating_ratio <= min_ratio, min_efficiency, efficiency)
+
+    efficiency = np.where(operating_ratio >= max_ratio, max_efficiency, efficiency)
+    return efficiency
+
+
+def calc_curve_coefficients():
+    """Calculates curve coefficients from BOP_efficiency_BOL.csv"""
+    df = pd.read_csv(os.path.join(file_path, "BOP_efficiency_BOL.csv"))
+    operating_ratios = df["operating_ratio"].values
+    efficiency = df["efficiency"].values
+
+    # Get min and max operating ratios
+    min_ratio_idx = df["operating_ratio"].idxmin()  # Index of minimum operating ratio
+    max_ratio_idx = df["operating_ratio"].idxmax()  # Index of maximum operating ratio
+
+    # Get the efficiency at the min and max operating ratios
+    min_efficiency = df["efficiency"].iloc[min_ratio_idx]
+    max_efficiency = df["efficiency"].iloc[max_ratio_idx]
+
+    # Get the actual min and max ratios
+    min_ratio = df["operating_ratio"].iloc[min_ratio_idx]
+    max_ratio = df["operating_ratio"].iloc[max_ratio_idx]
+
+    curve_coeff, curve_cov = scipy.optimize.curve_fit(
+        calc_efficiency_curve, operating_ratios, efficiency, p0=(1.0, 1.0, 1.0, 1.0)
+    )
+    return curve_coeff, min_ratio, max_ratio, min_efficiency, max_efficiency
+
+
+def pem_bop(
+    power_profile_to_electrolyzer_kw,
+    electrolyzer_rated_mw,
+    electrolyzer_turn_down_ratio,
+):
+    """
+    Calculate PEM balance of plant energy consumption at the beginning-of-life
+    based on power provided to the electrolyzer.
+
+    Args:
+        power_profile_to_electrolyzer_kw (list or np.array): Power profile to the electrolyzer in kW.
+        electrolyzer_rated_mw (float): The rating of the PEM electrolyzer in MW.
+        electrolyzer_turn_down_ratio (float): The electrolyzer turndown ratio.
+
+    Returns:
+        energy_consumption_bop_kwh (list or np.array): Energy consumed by electrolyzer BOP in kWh.
+    """
+    operating_ratios = power_profile_to_electrolyzer_kw / (electrolyzer_rated_mw * 1e3)
+
+    curve_coeff, min_ratio, max_ratio, min_efficiency, max_efficiency = (
+        calc_curve_coefficients()
+    )
+
+    efficiencies = calc_efficiency_curve(
+        operating_ratios,
+        *curve_coeff,
+    )  # kwh/kg
+
+    efficiencies = calc_efficiency(
+        operating_ratios,
+        efficiencies,
+        min_ratio,
+        max_ratio,
+        min_efficiency,
+        max_efficiency,
+    )
+
+    BOL_efficiency = get_electrolyzer_BOL_efficiency()  # kwh/kg
+
+    BOL_kg = (electrolyzer_rated_mw * 1000) / BOL_efficiency  # kg/hr
+
+    energy_consumption_bop_kwh = efficiencies * BOL_kg  # kwh
+
+    energy_consumption_bop_kwh = np.where(
+        power_profile_to_electrolyzer_kw
+        < electrolyzer_turn_down_ratio * electrolyzer_rated_mw * 1000,
+        0,
+        energy_consumption_bop_kwh,
+    )
+
+    return energy_consumption_bop_kwh

greenheart/simulation/technologies/hydrogen/electrolysis/PEM_BOP/README.md

@@ -0,0 +1,26 @@
+# Proton Exchange Membrane Water Electrolysis Balance-of-Plant
+
+This balance-of-plant (BOP) model is derived from Wang et. al (2023). It is represented as an overall BOP efficiency curve (kWh/kg) at different operating ratios at the beginning-of-life (BOL) of the electrolyzer. The operating ratios are the percentage of rated power provided to the electrolyzer.
+
+The BOP curve is largely dominated by electrical losses from the power electronics which includes a transformer and a rectifier to condition alternating current power, but there are additional mechanical and hydrogen losses.
+
+The model in GreenHEART calculates a curve fit based on the provided CSV it also limits the calculated efficiencies to the maximum and minimum operating ratio values in the CSV, making the overall function a piecewise implementation.
+
+**NOTE**: BOP assumes AC current as an input and assumes power electronics for AC to DC conversion. BOP efficiency curve has not been optimized for economies of scale or other electrical infrastructure connections.
+
+
+Citation for BOP model. Losses from BOP are shown in Figure 8. in Wang et. al (2023).
+```
+@Article{en16134964,
+AUTHOR = {Wang, Xiaohua and Star, Andrew G. and Ahluwalia, Rajesh K.},
+TITLE = {Performance of Polymer Electrolyte Membrane Water Electrolysis Systems: Configuration, Stack Materials, Turndown and Efficiency},
+JOURNAL = {Energies},
+VOLUME = {16},
+YEAR = {2023},
+NUMBER = {13},
+ARTICLE-NUMBER = {4964},
+URL = {https://www.mdpi.com/1996-1073/16/13/4964},
+ISSN = {1996-1073},
+DOI = {10.3390/en16134964}
+}
+```