add pem bop tests

greenheart/simulation/greenheart_simulation.py

@@ -625,8 +625,6 @@ def energy_internals(
             greenheart_config,
             electrolyzer_physics_results
         )
-        print("ElectrolyzEr bop",electrolyzer_energy_consumption_bop_kw[1])
-        print("Electrolyzer Energy Consumption BOP",np.mean(electrolyzer_energy_consumption_bop_kw))
 
         desal_results = he_elec.run_desal(
             hopp_config, electrolyzer_physics_results, design_scenario, verbose

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

@@ -0,0 +1,120 @@
+import os
+import numpy as np
+import pandas as pd
+import scipy.optimize
+import matplotlib.pyplot as plt
+
+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.
+
+    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.
+    """
+    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.
+    """
+    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,
+):
+    from greenheart.tools.eco.electrolysis import get_electrolyzer_BOL_efficiency
+
+    """Calculate PEM balance of plant energy consumption 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/tools/eco/electrolysis.py

@@ -32,7 +32,7 @@
     PEM_H2_Clusters as PEMClusters,
 )
 
-from greenheart.simulation.technologies.hydrogen.electrolysis.PEM_BOP import pem_bop
+from greenheart.simulation.technologies.hydrogen.electrolysis.PEM_BOP.PEM_BOP import pem_bop
 
 # from electrolyzer import run_electrolyzer
 
@@ -453,7 +453,8 @@ def run_electrolyzer_bop(
     if "include_bop_power" in plant_config["electrolyzer"]:
         if plant_config['electrolyzer']['include_bop_power']:
             energy_consumption_bop = pem_bop(electrolyzer_physics_results["power_to_electrolyzer_kw"],
-                                            plant_config['electrolyzer']["rating"])
+                                            plant_config["electrolyzer"]["rating"],
+                                            plant_config["electrolyzer"]["turndown_ratio"])
         else:
             energy_consumption_bop = np.zeros(len(electrolyzer_physics_results["power_to_electrolyzer_kw"]))
     else:

tests/greenheart/test_hydrogen/test_PEM_bop.py

@@ -0,0 +1,45 @@
+from pytest import approx, fixture
+import numpy as np
+
+from greenheart.simulation.technologies.hydrogen.electrolysis.PEM_BOP.PEM_BOP import (
+    pem_bop,
+)
+
+
+@fixture
+def bop_energy():
+    power_profile_kw = np.array(
+        [
+            0,
+            9999,  # just below turndown ratio
+            10000,  # exactly at turndown ratio
+            82746,
+            93774,  # max power in csv. largest operating ratio
+            100000,  # full power
+        ]
+    )
+
+    electrolyzer_rating_mw = 100  # MW
+    turndown_ratio = 0.1
+
+    bop_energy = pem_bop(
+        power_profile_to_electrolyzer_kw=power_profile_kw,
+        electrolyzer_rated_mw=electrolyzer_rating_mw,
+        electrolyzer_turn_down_ratio=turndown_ratio,
+    )
+    return bop_energy
+
+
+def test_bop_energy(subtests, bop_energy):
+    with subtests.test("No power"):
+        assert bop_energy[0] == 0
+    with subtests.test("below turndown"):
+        assert bop_energy[1] == 0
+    with subtests.test("at turndown"):
+        assert bop_energy[2] == approx(11032.668, 1e-2)
+    with subtests.test("mid-range power"):
+        assert bop_energy[3] == approx(6866.87, 1e-2)
+    with subtests.test("max power in curve"):
+        assert bop_energy[4] == approx(7847.85)
+    with subtests.test("full power"):
+        assert bop_energy[5] == approx(7847.85)