diff --git a/src/models/pump_detailed.py b/src/models/pump_detailed.py new file mode 100644 index 00000000..45e09454 --- /dev/null +++ b/src/models/pump_detailed.py @@ -0,0 +1,500 @@ +import pandas as pd +from enum import Enum, auto +from numpy import polyfit + +from pyomo.environ import ( + Var, + Param, + Set, + check_optimal_termination, + value, + units as pyunits, + Reals, +) +from pyomo.common.config import ConfigValue, In + +from idaes.models.unit_models.pressure_changer import PumpData +from idaes.core import declare_process_block_class +from idaes.core.util.constants import Constants +from idaes.core.util.exceptions import InitializationError +from idaes.core.util.exceptions import ConfigurationError +import idaes.logger as idaeslog + +from watertap.core import InitializationMixin +from watertap.costing.unit_models.pump import cost_pump +from watertap.core.solvers import get_solver + +_log = idaeslog.getLogger(__name__) + + +# ------- Validation functions ---------------------- +def _validate_curve_data(filepath): + if isinstance(filepath, str): + try: + df = pd.read_csv(filepath) + except Exception as e: + raise ValueError(f"Failed to read CSV file '{filepath}': {e}") + required_cols = ["flow (m3/s)", "head (m)", "efficiency (-)"] + if all(col in df.columns for col in required_cols): + return df + raise ValueError( + f"CSV file must contain columns: {required_cols}, but found: {list(df.columns)}" + ) + raise ValueError("Must be a string filepath to a CSV file") + + +def _validate_surrogate_coeffs(coeffs): + if not isinstance(coeffs, dict): + raise ValueError("Surrogate coefficients must be provided as a dictionary.") + + expected_keys = {0, 1, 2, 3} + if set(coeffs.keys()) != expected_keys: + raise ValueError("Surrogate coefficient keys must be exactly {0, 1, 2, 3}.") + + return coeffs + + +# -------------------------------------------------- + + +# ---------- Enums for configuration options ---------- +class Efficiency(Enum): + Fixed = auto() # default is constant efficiency + Flow = auto() # flow-only correlation + + +class PumpCurveDataType(Enum): + DataSet = auto() + # user provides flow, head, and efficiency data to fit a curve and surrogate coefficients are calculated via polyfit + SurrogateCoefficent = auto() + # curve fit correlation based on flow and head is provided by the user as surrogate coefficients. + + +# -------------------------------------------------- + + +@declare_process_block_class("PumpDetailed") +class PumpIsothermalData(InitializationMixin, PumpData): + """ + Detailed Isothermal Pump Unit Model Class + + * User can specify whether to use a variable efficiency based on flow and head, or to use a constant efficiency at the best efficiency point (BEP). + * This function builds on the standard pump model by adding variables and constraints to represent pump performance curves. + * The model includes options for defining the pump curve based on flow and head at the BEP, as well as scaling the efficiency curve accordingly. + * This allows for more accurate representation of pump performance under varying operating conditions. + """ + + CONFIG = PumpData.CONFIG() + + CONFIG.declare( + "variable_efficiency", + ConfigValue( + default=Efficiency.Fixed, + domain=In(Efficiency), + description="Variable pump efficiency flag", + doc="""Indicates whether pump efficiency should be variable or fixed. + ** Fixed: The default is constant efficiency at the best efficiency point (BEP) + ** Flow: Efficiency is a function of flow only at reference speed based on pump performance curve data provided. + """, + ), + ) + + # Config defining if surrogate coefficients should be calculated based or hardcoded. + CONFIG.declare( + "pump_curve_data_type", + ConfigValue( + default=PumpCurveDataType.SurrogateCoefficent, + domain=In(PumpCurveDataType), + description="Type of pump curve data", + doc="""Indicates whether the pump curve data is provided as a dataset or as surrogate coefficients.""", + ), + ) + + CONFIG.declare( + "head_surrogate_coeffs", + ConfigValue( + default=None, + domain=_validate_surrogate_coeffs, + description="Surrogate coefficients for the pump head curve based on flow only", + doc="""Coefficients for the pump head curve surrogate based on flow only. + The head is calculated as a function of flow using a cubic polynomial with these coefficients.""", + ), + ) + + CONFIG.declare( + "efficiency_surrogate_coeffs", + ConfigValue( + default=None, + domain=_validate_surrogate_coeffs, + description="Surrogate coefficients for the pump efficiency curve based on flow only", + doc="""Coefficients for the pump efficiency curve surrogate based on flow only. + The efficiency is calculated as a function of flow using a cubic polynomial with these coefficients.""", + ), + ) + + CONFIG.declare( + "pump_curves", + ConfigValue( + default=None, + domain=_validate_curve_data, + description="Head and efficiency curve at reference speed", + doc="Data from digitization of head vs. flow and efficiency vs. flow curves at a reference speed from a pump datasheet", + ), + ) + + def _validate_config(self): + if self.config.variable_efficiency is Efficiency.Flow: + if self.config.pump_curve_data_type == PumpCurveDataType.DataSet: + if self.config.pump_curves is None: + raise ConfigurationError( + "pump_curves must be provided as a CSV filepath when pump_curve_data_type is DataSet." + ) + elif ( + self.config.pump_curve_data_type + == PumpCurveDataType.SurrogateCoefficent + ): + if ( + not self.config.head_surrogate_coeffs + or not self.config.efficiency_surrogate_coeffs + ): + raise ConfigurationError( + "surrogate_coeffs must be provided for the pump head curve and efficiency curve when pump_curve_data_type is set to SurrogateCoefficent." + ) + + # ----------------------------------------------- + + def build(self): + super().build() + self._validate_config() + + # Isothermal pump set up + if hasattr(self.control_volume, "enthalpy_balances"): + self.control_volume.del_component(self.control_volume.enthalpy_balances) + + @self.control_volume.Constraint( + self.flowsheet().config.time, doc="Isothermal constraint" + ) + def isothermal_balance(b, t): + return b.properties_in[t].temperature == b.properties_out[t].temperature + + if self.config.variable_efficiency is not Efficiency.Fixed: + # Variable efficiency pump set-up + #### Design point variables #### + self.design_flow = Var( + initialize=1.0, + bounds=(0, 10000), + doc="""Design flowrate of the centrifugal pump. + This could be the flowrate at the best efficiency point (BEP) or user selected operating point. + Used to build the system curve.""", + units=pyunits.m**3 / pyunits.s, + ) + + self.design_head = Var( + initialize=1.0, + bounds=(0, 10000), + doc="""Design head of the centrifugal pump. + This could be the head at the best efficiency point (BEP) or user selected operating point. + Used to build the system curve.""", + units=pyunits.m, + ) + + self.design_efficiency = Var( + initialize=0.8, + bounds=(0, 1), + doc="""Design efficiency of the centrifugal pump. + This could be the efficiency at the best efficiency point (BEP) or user selected operating point.""", + units=pyunits.dimensionless, + ) + + self.design_speed_fraction = Var( + initialize=0.8, + bounds=(0, 1), + doc="""Design speed fraction of the centrifugal pump. + This could be the speed fraction at the best efficiency point (BEP) or user selected operating point.""", + units=pyunits.dimensionless, + ) + + ### System curve variables ### + # design_head = system_curve_geometric_head + system_curve_flow_constant * (design_flow)**2 + + self.system_curve_geometric_head = Var( + initialize=0.0, + bounds=(0, 10000), + doc="""Geometric head constant for the pump, that represents the static head component used to define the system curve.""", + units=pyunits.m, + ) + + self.system_curve_flow_constant = Var( + initialize=1.0, + bounds=(0, 10000), + doc="""Geometric flow constant for the pump, represents the major and minor losses in pump. Used to define the system curve.""", + units=pyunits.m * (pyunits.m**3 / pyunits.s) ** (-2), + ) + + # Constraints connecting inlet and outlet conditions to the design point head and flow, used to solve for the system curve constants + @self.Constraint( + doc="Design head is the pressure difference across the pump at the design point" + ) + def design_head_constraint(b): + return b.design_head == b.control_volume.deltaP[0] / ( + b.control_volume.properties_out[0].dens_mass_phase["Liq"] + * Constants.acceleration_gravity + ) + + @self.Constraint( + doc="Design flow is the flow through the pump at the design point" + ) + def design_flow_constraint(b): + return b.design_flow == pyunits.convert( + b.control_volume.properties_in[0].flow_vol_phase["Liq"], + to_units=pyunits.m**3 / pyunits.s, + ) + + # Constraints to calculate system curve constants based on design point + @self.Constraint(doc="System curve head constant calculation") + def system_curve_head_constant_calculation(b): + return ( + b.design_head + == b.system_curve_geometric_head + + b.system_curve_flow_constant * (b.design_flow**2) + ) + + #### Pump curve variables #### + self.ref_flow = Var( + initialize=1.0, + bounds=(0, 10000), + doc="Reference flowrate for the pump on the pump curve from specification sheet", + units=pyunits.m**3 / pyunits.s, + ) + + self.ref_head = Var( + initialize=10.0, + bounds=(0, 10000), + doc="Reference head for the pump on the pump curve from specification sheet.", + units=pyunits.m, + ) + + self.ref_efficiency = Var( + initialize=0.8, + bounds=(0, 1), + doc="Reference efficiency for the pump on the pump curve from specification sheet.", + units=pyunits.dimensionless, + ) + + self.ref_speed_fraction = Var( + initialize=1, + bounds=(0, 1), + doc="Reference speed fraction for the pump on the pump curve from specification sheet.", + units=pyunits.dimensionless, + ) + + # Add efficiency variables for more and VFD efficiency. Only defined for variable efficiency + self.vfd_efficiency = Param( + initialize=0.97, + mutable=True, + doc="Variable frequency drive (VFD) efficiency", + units=pyunits.dimensionless, + ) + + self.motor_efficiency = Param( + initialize=0.95, + mutable=True, + doc="Motor efficiency", + units=pyunits.dimensionless, + ) + + if self.config.pump_curve_data_type == PumpCurveDataType.DataSet: + self.surrogate_index = Set( + initialize=range(4), + doc="Index for surrogate coefficients for a cubic polynomial fit", + ) + + # pump_curves is converted from a filepath name to DataFrame from the validator + curves_df = self.config.pump_curves + + p_head = polyfit( + curves_df["flow (m3/s)"], + curves_df["head (m)"], + 3, + ) + p_eff = polyfit( + curves_df["flow (m3/s)"], + curves_df["efficiency (-)"], + 3, + ) + head_surrogate_coeffs = { + i: float(p_head[3 - i]) for i in self.surrogate_index + } + efficiency_surrogate_coeffs = { + i: float(p_eff[3 - i]) for i in self.surrogate_index + } + + elif ( + self.config.pump_curve_data_type + == PumpCurveDataType.SurrogateCoefficent + ): + self.surrogate_index = Set( + initialize=[0, 1, 2, 3], doc="Index for surrogate coefficients" + ) + + head_surrogate_coeffs = self.config.head_surrogate_coeffs + efficiency_surrogate_coeffs = self.config.efficiency_surrogate_coeffs + + # Also would be nice for fit to be in terms of ft and gpm, or m and m3/hr or m3/s + + # Note: Coefficients are dimensionless because they're applied to flow in m³/s (numeric value) + # The polynomial produces head in meters + self.head_surrogate_coefficients = Param( + self.surrogate_index, + initialize=head_surrogate_coeffs, + doc="Coefficients for the head surrogate based on flow only (flow in m³/s, output in m)", + units=pyunits.dimensionless, + ) + + self.efficiency_surrogate_coefficients = Param( + self.surrogate_index, + initialize=efficiency_surrogate_coeffs, + doc="Coefficients for the efficiency surrogate based on flow only", + units=pyunits.dimensionless, + ) + + # Constraints to calculate reference point variables by solving the system curve and pump curve simultaneously at the reference point. + @self.Constraint(doc="Reference head equality") + def ref_head_constraint(b): + return ( + b.ref_head + == b.system_curve_geometric_head + + b.system_curve_flow_constant * (b.ref_flow**2) + ) + + @self.Constraint( + doc="Reference head calculation calculated using the pump curve surrogate" + ) + def ref_head_surrogate(b): + # Convert flow to dimensionless value in m³/s for polynomial + # flow_val = pyunits.convert(b.ref_flow, to_units=pyunits.m**3/pyunits.s) + return ( + b.ref_head + == ( + b.head_surrogate_coefficients[0] + + b.head_surrogate_coefficients[1] + * (b.ref_flow / (pyunits.m**3 / pyunits.s)) + + b.head_surrogate_coefficients[2] + * ((b.ref_flow / (pyunits.m**3 / pyunits.s)) ** 2) + + b.head_surrogate_coefficients[3] + * ((b.ref_flow / (pyunits.m**3 / pyunits.s)) ** 3) + ) + * pyunits.m + ) + + # Calculate the reference efficiency based on the surrogate coefficients and reference flow + @self.Constraint( + doc="Reference efficiency calculated calculated using the pump curve surrogate" + ) + def ref_efficiency_constraint(b): + # Convert flow to dimensionless value in m³/s for polynomial + return b.ref_efficiency == ( + b.efficiency_surrogate_coefficients[0] + + b.efficiency_surrogate_coefficients[1] + * (b.ref_flow / (pyunits.m**3 / pyunits.s)) + + b.efficiency_surrogate_coefficients[2] + * ((b.ref_flow / (pyunits.m**3 / pyunits.s)) ** 2) + + b.efficiency_surrogate_coefficients[3] + * ((b.ref_flow / (pyunits.m**3 / pyunits.s)) ** 3) + ) + + @self.Constraint( + doc="Design speed fraction calculation using affinity laws" + ) + def design_speed_fraction_constraint(b): + return b.design_flow == b.ref_flow * ( + b.design_speed_fraction / b.ref_speed_fraction + ) + + # Expression to calculate the efficiency at the design point based as a function of speed fraction + @self.Constraint(doc="Design efficiency calculation") + def design_efficiency_constraint(b): + return b.design_efficiency == b.ref_efficiency + + @self.Constraint( + doc="Overall efficiency calculation including motor efficiency and VFD efficiency" + ) + def overall_efficiency_constraint(b): + return ( + b.efficiency_pump[0] + == b.design_efficiency * b.motor_efficiency * b.vfd_efficiency + ) + + elif self.config.variable_efficiency is Efficiency.Fixed: + # Fixed efficiency pump set-up (efficiency is constant at the best efficiency point) + # User should directly fix efficiency_pump[0] and deltaP + pass + + def initialize_build( + self, + state_args=None, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + ): + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") + opt = get_solver(solver, optarg) + + if state_args is None: + state_args = {} + state_dict = self.control_volume.properties_in[ + self.flowsheet().config.time.first() + ].define_port_members() + + for k in state_dict.keys(): + if state_dict[k].is_indexed(): + state_args[k] = {} + for m in state_dict[k].keys(): + state_args[k][m] = state_dict[k][m].value + else: + state_args[k] = state_dict[k].value + + flags = self.control_volume.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args, + ) + + init_log.info_high("Initialization Step 1 Complete.") + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + + init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res))) + + self.control_volume.release_state(flags, outlvl + 1) + init_log.info("Initialization Complete: {}".format(idaeslog.condition(res))) + + if not check_optimal_termination(res): + # It's possible initialization fails if the initial flowrate + # if m.fs.unit.design_speed_fraction.fixed and m.fs.unit.design_head.fixed: + # design_speed = m.fs.unit.design_speed_fraction.value + # m.fs.unit.design_speed_fraction.unfix() + # m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].value * .75) + # m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix(m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].value * .75) + # m.fs.unit.control_volume.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].unfix() # This should remain fixed the whole time, but unfixing b/c it might cause property model to fail. + # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + # res = opt.solve(self, tee=slc.tee) + # # Then refix speed, try again + # m.fs.unit.design_speed_fraction.fix(design_speed) + # m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].unfix() + # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + # res = opt.solve(self, tee=slc.tee) + # if not check_optimal_termination(res): + # raise InitializationError(f"Unit model {self.name} failed to initialize") + # else: + raise InitializationError(f"Unit model {self.name} failed to initialize") + + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + + @property + def default_costing_method(self): + return cost_pump diff --git a/src/models/tests/test_pump_curves_data.csv b/src/models/tests/test_pump_curves_data.csv new file mode 100644 index 00000000..f813bd15 --- /dev/null +++ b/src/models/tests/test_pump_curves_data.csv @@ -0,0 +1,10 @@ +flow (m3/s),efficiency (-),head (m) +0.095,0.59,92.9929878 +0.118,0.68,90.40823171 +0.144,0.75,87.62012195 +0.169,0.8,84.26768293 +0.175,0.825,82.07073171 +0.18,0.83,81.71707317 +0.202,0.83,75.71371951 +0.215,0.8,72.25060976 +0.233,0.755,64.12865854 diff --git a/src/models/tests/test_pump_curves_data_uf.csv b/src/models/tests/test_pump_curves_data_uf.csv new file mode 100644 index 00000000..87dfe030 --- /dev/null +++ b/src/models/tests/test_pump_curves_data_uf.csv @@ -0,0 +1,11 @@ +flow (m3/s),efficiency (-),head (m) +0.063083523,0.38,92.98780488 +0.094625284,0.52,88.41463415 +0.126167045,0.64,85.36585366 +0.157708806,0.72,83.84146341 +0.189250568,0.79,79.26829268 +0.220792329,0.81,73.17073171 +0.242997729,0.83,66.46341463 +0.25233409,0.82,64.02439024 +0.315417613,0.71,44.20731707 +0.331188494,0.63,36.58536585 diff --git a/src/models/tests/test_pump_detailed.py b/src/models/tests/test_pump_detailed.py new file mode 100644 index 00000000..3105a73c --- /dev/null +++ b/src/models/tests/test_pump_detailed.py @@ -0,0 +1,492 @@ +import pytest +import os +from pyomo.environ import ( + ConcreteModel, + assert_optimal_termination, + units as pyunits, + value, + Reals, +) +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.scaling import calculate_scaling_factors +from idaes.core.util.exceptions import ConfigurationError +from idaes.core.util.constants import Constants +from pyomo.util.check_units import assert_units_consistent +from watertap.property_models.seawater_prop_pack import SeawaterParameterBlock +from watertap.core.solvers import get_solver +from models.pump_detailed import PumpDetailed, Efficiency, PumpCurveDataType + + +solver = get_solver() + + +# Build function with design flow and head as inputs +def build_pump_w_flow_head(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.SurrogateCoefficent, + head_surrogate_coeffs={0: 114.22, 1: -410.6, 2: 2729.2, 3: -8089.1}, + efficiency_surrogate_coeffs={0: 0.389, 1: -0.535, 2: 41.373, 3: -138.82}, + ) + + # Input flow and head + feed_flow_vol = 0.126 * pyunits.m**3 / pyunits.s + pump_head = 60.96 * pyunits.m + density = 1000 * pyunits.kg / pyunits.m**3 + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_TDS = 0.035 + feed_pressure_in = 101325 * pyunits.Pa + feed_pressure_out = ( + feed_pressure_in + pump_head * density * Constants.acceleration_gravity + ) + feed_temperature = 273.15 + 25 + + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + m.fs.unit.outlet.pressure[0].fix(feed_pressure_out) + + m.fs.unit.system_curve_geometric_head.fix(4.57) + m.fs.unit.ref_speed_fraction.fix(1.0) + + assert_units_consistent(m) + + return m + + +# Build function with design flow and speed as inputs +def build_pump_w_flow_speed(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.SurrogateCoefficent, + head_surrogate_coeffs={0: 114.22, 1: -410.6, 2: 2729.2, 3: -8089.1}, + efficiency_surrogate_coeffs={0: 0.389, 1: -0.535, 2: 41.373, 3: -138.82}, + ) + # Input flow and speed + feed_flow_vol = 0.126 * pyunits.m**3 / pyunits.s + design_speed_fraction = 0.829 + density = 1000 * pyunits.kg / pyunits.m**3 + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_TDS = 0.035 + + feed_pressure_in = 101325 * pyunits.Pa + feed_temperature = 273.15 + 25 + + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + + m.fs.unit.design_speed_fraction.fix(design_speed_fraction) + m.fs.unit.system_curve_geometric_head.fix(4.57) + m.fs.unit.ref_speed_fraction.fix(1.0) + + calculate_scaling_factors(m) + assert_units_consistent(m) + + return m + + +@pytest.mark.component +def test_fixed_eff_pump(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Fixed, + ) + + # Input flow and speed + feed_flow_vol = 0.126 * pyunits.m**3 / pyunits.s + density = 1000 * pyunits.kg / pyunits.m**3 + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_TDS = 0.035 + + feed_pressure_in = 101325 * pyunits.Pa + feed_temperature = 273.15 + 25 + + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + + m.fs.unit.efficiency_pump.fix(0.85) + m.fs.unit.deltaP.fix(500000) + + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + results = solver.solve(m) + assert_optimal_termination(results) + + assert pytest.approx(m.fs.unit.work_mechanical[0].value, rel=1e-3) == 72411 + + +# Three tests for different combinations of inputs for variable efficiency +@pytest.mark.component +def test_pump_w_flow_head(): + m = build_pump_w_flow_head() + + assert hasattr(m.fs.unit, "inlet") + assert hasattr(m.fs.unit, "outlet") + assert hasattr(m.fs.unit, "deltaP") # this is just a reference + assert hasattr(m.fs.unit.control_volume, "deltaP") + + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + results = solver.solve(m) + assert_optimal_termination(results) + + +@pytest.mark.component +def test_pump_w_flow_speed(): + m = build_pump_w_flow_speed() + + assert hasattr(m.fs.unit, "inlet") + assert hasattr(m.fs.unit, "outlet") + assert hasattr(m.fs.unit, "deltaP") # this is just a reference + assert hasattr(m.fs.unit.control_volume, "deltaP") + + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + results = solver.solve(m) + assert_optimal_termination(results) + + +@pytest.mark.component +def test_pump_w_head_speed(): + m = build_pump_w_flow_head() + + assert hasattr(m.fs.unit, "inlet") + assert hasattr(m.fs.unit, "outlet") + assert hasattr(m.fs.unit, "deltaP") # this is just a reference + assert hasattr(m.fs.unit.control_volume, "deltaP") + + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + m.fs.unit.design_speed_fraction.fix(0.8184) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].unfix() + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].unfix() + m.fs.unit.control_volume.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].fix() + + results = solver.solve(m) + assert_optimal_termination(results) + + +# Test for passing a dataset for the pump curves +@pytest.mark.component +def test_data_points(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + pump_curves_filepath = os.path.join( + os.path.dirname(__file__), "test_pump_curves_data.csv" + ) + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + # flow in m3/s and head in m + pump_curves=pump_curves_filepath, + ) + + # Input flow and head + feed_flow_vol = 0.1435 * pyunits.m**3 / pyunits.s + pump_head = 87 * pyunits.m + density = 1000 * pyunits.kg / pyunits.m**3 + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_TDS = 0.035 + feed_pressure_in = 101325 * pyunits.Pa + feed_pressure_out = ( + feed_pressure_in + pump_head * density * Constants.acceleration_gravity + ) + feed_temperature = 273.15 + 25 + + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + m.fs.unit.outlet.pressure[0].fix(feed_pressure_out) + + m.fs.unit.system_curve_geometric_head.fix(0) + m.fs.unit.ref_speed_fraction.fix(1.0) + + assert hasattr(m.fs.unit, "inlet") + assert hasattr(m.fs.unit, "surrogate_index") + + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + results = solver.solve(m) + assert_optimal_termination(results) + + +@pytest.mark.component +def test_low_speed(): + # Doubles as low speed (50%) test + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + pump_curves=os.path.join( + os.path.dirname(__file__), "test_pump_curves_data_uf.csv" + ), + ) + # Input flow and head for initial solve + feed_flow_vol = 0.12 * pyunits.m**3 / pyunits.s + pump_head = 21.3 * pyunits.m + density = 1000 * pyunits.kg / pyunits.m**3 + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_TDS = 0.035 + + feed_pressure_in = 101325 * pyunits.Pa + feed_pressure_out = ( + feed_pressure_in + pump_head * density * Constants.acceleration_gravity + ) + feed_temperature = 273.15 + 25 + + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + m.fs.unit.outlet.pressure[0].fix(feed_pressure_out) + + m.fs.unit.system_curve_geometric_head.fix(0) + m.fs.unit.ref_speed_fraction.fix(1.0) + + calculate_scaling_factors(m) + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + # Now apply the pump speed + test_pump_speed = 0.5 + m.fs.unit.design_speed_fraction.fix(test_pump_speed) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].unfix() + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].unfix() + m.fs.unit.control_volume.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].fix() + calculate_scaling_factors(m) + + results = solver.solve(m) + assert_optimal_termination(results) + + assert m.fs.unit.efficiency_pump[0].value == pytest.approx(0.64, abs=0.02) + + +@pytest.mark.component +def test_negative_inlet_pressure(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + pump_curves=os.path.join( + os.path.dirname(__file__), + "test_pump_curves_data_uf.csv", + ), + ) + + # Input flow and head for initial solve + feed_flow_vol = ( + 0.21 * pyunits.m**3 / pyunits.s + ) # Actual flow vol will be calculated in second solve + density = 1000 * pyunits.kg / pyunits.m**3 + feed_pressure_in = -12 * pyunits.psi # 101325 * pyunits.Pa + + feed_pressure_out = 50 * pyunits.psi + feed_mass_frac_TDS = 0.0005 + feed_temperature = 298.15 + geometric_head = pyunits.convert( + 12 * pyunits.psi / (density * Constants.acceleration_gravity), + to_units=pyunits.m, + ) + + # Fix pump characteristics + m.fs.unit.system_curve_geometric_head.fix(geometric_head) + m.fs.unit.ref_speed_fraction.fix(1.0) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + m.fs.unit.outlet.pressure[0].fix(feed_pressure_out) + + # Change the bounds for the inlet pressure + m.fs.unit.control_volume.properties_in[0].pressure.setlb(None) + m.fs.unit.control_volume.properties_in[0].pressure.domain = Reals + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + calculate_scaling_factors(m) + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + solver = get_solver() + test_pump_speed = 0.91 + m.fs.unit.design_speed_fraction.fix(test_pump_speed) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].unfix() + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].unfix() + m.fs.unit.control_volume.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].fix() + calculate_scaling_factors(m) + results = solver.solve(m) + assert_optimal_termination(results) + + assert value( + pyunits.convert(m.fs.unit.work_mechanical[0], to_units=pyunits.kW) + ) == pytest.approx(166.54, rel=1e-3) + + +# Test an invalid surrogate coefficient case +@pytest.mark.component +def test_invalid_surrogate_coefficients(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + with pytest.raises( + ValueError, + match=r"Surrogate coefficient keys must be exactly \{0, 1, 2, 3\}\.", + ): + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.SurrogateCoefficent, + head_surrogate_coeffs={0: 114.22, 1: -410.6, 2: 2729.2}, + efficiency_surrogate_coeffs={0: 0.389, 1: -0.535}, + ) + + +# Test an invalid filepath name +@pytest.mark.component +def test_invalid_filepath(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + with pytest.raises( + ValueError, + match=r"Failed to read CSV file '.*': .*", + ): + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + pump_curves=os.path.join(os.path.dirname(__file__), "DNE.csv"), + ) + + +# Test missing filepath for dataset mode +@pytest.mark.component +def test_missing_filepath_for_dataset_mode(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + with pytest.raises( + ConfigurationError, + match=r"pump_curves must be provided as a CSV filepath when pump_curve_data_type is DataSet\.", + ): + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + pump_curves=None, + ) + + +# Test with variable_efficiency not assigned using proper class +@pytest.mark.component +def test_invalid_efficiency_type(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + with pytest.raises( + ValueError, + match="'InvalidType' is not a valid Efficiency", + ): + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency="InvalidType", + pump_curve_data_type=PumpCurveDataType.SurrogateCoefficent, + ) + + +# Test surrogate mode selected but no surrogate coefficients provided +@pytest.mark.component +def test_missing_surrogate_values(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + with pytest.raises( + ConfigurationError, + match="surrogate_coeffs must be provided for the pump head curve and efficiency curve when pump_curve_data_type is set to SurrogateCoefficent.", + ): + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.SurrogateCoefficent, + ) diff --git a/src/models/tests/test_uf_pump.py b/src/models/tests/test_uf_pump.py new file mode 100644 index 00000000..632a35c1 --- /dev/null +++ b/src/models/tests/test_uf_pump.py @@ -0,0 +1,170 @@ +# file is meant to be moved into the case study days test files in the Testing PR +# Comparison of the uf pumps on the two days. + +import pytest +import os +from pyomo.environ import ( + ConcreteModel, + assert_optimal_termination, + units as pyunits, + value, + Reals, +) +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.scaling import calculate_scaling_factors +from pyomo.util.check_units import assert_units_consistent +from watertap.property_models.seawater_prop_pack import SeawaterParameterBlock +from watertap.core.solvers import get_solver +from models.pump_detailed import PumpDetailed, Efficiency, PumpCurveDataType + +solver = get_solver() + + +@pytest.mark.unit +def test_uf_pump_8_19_full_speed(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + pump_curves=os.path.join( + os.path.dirname(__file__), + "test_pump_curves_data_uf.csv", + ), + ) + + # Input flow and head for initial solve + feed_flow_vol = ( + 0.21 * pyunits.m**3 / pyunits.s + ) # Actual flow vol will be calculated in second solve + density = 1000 * pyunits.kg / pyunits.m**3 + feed_pressure_in = -12 * pyunits.psi # 101325 * pyunits.Pa + + m.fs.unit.control_volume.properties_in[0].pressure.setlb(None) + m.fs.unit.control_volume.properties_in[0].pressure.domain = Reals + + feed_pressure_out = 50 * pyunits.psi + feed_mass_frac_TDS = 0.0005 + feed_temperature = 298.15 + + geometric_head = pyunits.convert( + 12 * pyunits.psi / (density * 9.81 * pyunits.m / pyunits.s**2), + to_units=pyunits.m, + ) + + # Fix pump characteristics + m.fs.unit.system_curve_geometric_head.fix(geometric_head) + m.fs.unit.ref_speed_fraction.fix(1.0) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + m.fs.unit.outlet.pressure[0].fix(feed_pressure_out) + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + + calculate_scaling_factors(m) + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + solver = get_solver() + test_pump_speed = 0.91 + m.fs.unit.design_speed_fraction.fix(test_pump_speed) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].unfix() + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].unfix() + m.fs.unit.control_volume.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].fix() + + calculate_scaling_factors(m) + results = solver.solve(m) + assert_optimal_termination(results) + expected_power = 175 # kW + + print( + "Calculated power: ", + value(pyunits.convert(m.fs.unit.work_mechanical[0], to_units=pyunits.kW)), + ) + assert value( + pyunits.convert(m.fs.unit.work_mechanical[0], to_units=pyunits.kW) + ) == pytest.approx(expected_power, rel=0.15) + + +@pytest.mark.unit +def test_uf_pump_8_19_low_speed(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.unit = PumpDetailed( + property_package=m.fs.properties, + variable_efficiency=Efficiency.Flow, + pump_curve_data_type=PumpCurveDataType.DataSet, + pump_curves=os.path.join( + os.path.dirname(__file__), + "test_pump_curves_data_uf.csv", + ), + ) + + # Input flow and head for initial solve + feed_flow_vol = ( + 0.21 * pyunits.m**3 / pyunits.s + ) # Actual flow vol will be calculated in second solve + density = 1000 * pyunits.kg / pyunits.m**3 + feed_pressure_in = -12 * pyunits.psi # 101325 * pyunits.Pa + + m.fs.unit.control_volume.properties_in[0].pressure.setlb(None) + m.fs.unit.control_volume.properties_in[0].pressure.domain = Reals + + feed_pressure_out = 46 * pyunits.psi + feed_mass_frac_TDS = 0.0005 + feed_temperature = 298.15 + + geometric_head = pyunits.convert( + 12 * pyunits.psi / (density * 9.81 * pyunits.m / pyunits.s**2), + to_units=pyunits.m, + ) + + # Fix pump characteristics + m.fs.unit.system_curve_geometric_head.fix(geometric_head) + m.fs.unit.ref_speed_fraction.fix(1.0) + m.fs.unit.inlet.pressure[0].fix(feed_pressure_in) + m.fs.unit.inlet.temperature[0].fix(feed_temperature) + m.fs.unit.outlet.pressure[0].fix(feed_pressure_out) + + # Calculated feed conditions + feed_flow_mass = feed_flow_vol * density + feed_mass_frac_H2O = 1 - feed_mass_frac_TDS + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix( + feed_flow_mass * feed_mass_frac_TDS + ) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + calculate_scaling_factors(m) + m.fs.unit.initialize() + assert degrees_of_freedom(m) == 0 + + solver = get_solver() + test_pump_speed = 0.71 + m.fs.unit.design_speed_fraction.fix(test_pump_speed) + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].unfix() + m.fs.unit.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].unfix() + m.fs.unit.control_volume.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].fix() + calculate_scaling_factors(m) + + results = solver.solve(m) + assert_optimal_termination(results) + expected_power = 80 # kW + + assert value( + pyunits.convert(m.fs.unit.work_mechanical[0], to_units=pyunits.kW) + ) == pytest.approx(expected_power, rel=0.15) diff --git a/src/wrd/utilities.py b/src/wrd/utilities.py index 17c7126d..b4f6592e 100644 --- a/src/wrd/utilities.py +++ b/src/wrd/utilities.py @@ -333,3 +333,28 @@ def touch_flow_and_conc(b): sb = b.find_component(f"{x}_state") sb[0].flow_vol_phase sb[0].conc_mass_phase_comp + + +def ft_head_to_psi(head): + """ + Convert head (in feet) to pressure (in psi) + Default fluid density is for water at room temperature (998.2 kg/m^3) + Default gravity is 9.81 m/s^2 + """ + # head should have units of ft already + fluid_density = 62.4 * pyunits.lb / pyunits.ft**3 + gravity = 32.174 * pyunits.ft / pyunits.s**2 + pressure = fluid_density * gravity * head + return pyunits.convert(pressure, to_units=pyunits.psi) + + +def psi_to_ft_head(pressure): + """ + Convert pressure (in psi) to head (in feet) + Default fluid density is for water at room temperature (998.2 kg/m^3) + Default gravity is 9.81 m/s^2 + """ + fluid_density = 62.4 * pyunits.lb / pyunits.ft**3 + gravity = 32.174 * pyunits.ft / pyunits.s**2 + head = pyunits.convert(pressure / (fluid_density * gravity), to_units=pyunits.ft) + return head