From b0db08f144827ff59459def3dd709d8d70e5edaa Mon Sep 17 00:00:00 2001 From: Cory Frontin Date: Fri, 17 Jul 2026 23:01:55 -0600 Subject: [PATCH 1/6] progress toward NGCT thermodynamic model in H2I --- .../combustion_machines/NGCT_thermo_model.py | 620 ++++++++++++++++++ .../combustion_machines/__init__.py | 3 + .../test_NGCT_thermo_model.py | 260 ++++++++ .../turbine_simple_cycle.py | 357 ++++++++++ h2integrate/core/supported_models.py | 1 + 5 files changed, 1241 insertions(+) create mode 100644 h2integrate/converters/combustion_machines/NGCT_thermo_model.py create mode 100644 h2integrate/converters/combustion_machines/__init__.py create mode 100644 h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py create mode 100644 h2integrate/converters/combustion_machines/turbine_simple_cycle.py diff --git a/h2integrate/converters/combustion_machines/NGCT_thermo_model.py b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py new file mode 100644 index 000000000..b7fc3ee0d --- /dev/null +++ b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py @@ -0,0 +1,620 @@ +from pprint import pprint + +from matplotlib.pylab import power +import numpy as np + +import pyfluids + +DRY_AIR_MASS_FRACTIONS = { + pyfluids.FluidsList.Nitrogen: 75.51760366879898, # generated by copilot, matches engineeringtoolbox.com data + pyfluids.FluidsList.Oxygen: 23.139561050152018, # generated by copilot, matches engineeringtoolbox.com data + pyfluids.FluidsList.Argon: 1.2881375564363484, # generated by copilot, matches engineeringtoolbox.com data + pyfluids.FluidsList.CarbonDioxide: 0.05469772461264626, # generated by copilot, matches engineeringtoolbox.com data +} + + +def humidity_ratio_to_water_mass_fraction(humidity_ratio: float) -> float: + return humidity_ratio / (1.0 + humidity_ratio) + + +def enforce_gas_phase(fluid): + phases = getattr(pyfluids, "Phases", None) + phase_gas = getattr(phases, "Gas", None) if phases is not None else None + if phase_gas is None: + return fluid + + try: + fluid.specify_phase(phase_gas) + except Exception: + # Keep behavior backward compatible if phase forcing is not supported. + pass + return fluid + + +def make_humid_air_mixture( + pressure: float, + temp_rel: float, + rel_humidity: float, +) -> pyfluids.Mixture: + + # get absolute mass humidity ratio using conditions w/ pyfluid humidair + humidity_ratio = ( + pyfluids.HumidAir() + .with_state( + pyfluids.InputHumidAir.pressure(pressure), + pyfluids.InputHumidAir.temperature(temp_rel), + pyfluids.InputHumidAir.relative_humidity(rel_humidity), + ) + .humidity + ) + + # use a mix of water and dry air + water_mass_fraction = 100.0 * humidity_ratio_to_water_mass_fraction(humidity_ratio) + dry_air_mass_fraction = 100.0 - water_mass_fraction + + # create a mixture using our dry air mix plus water vapor + fluids = list(DRY_AIR_MASS_FRACTIONS.keys()) + [pyfluids.FluidsList.Water] + fractions = [ + dry_air_mass_fraction * mass_fraction / 100.0 + for mass_fraction in DRY_AIR_MASS_FRACTIONS.values() + ] + [water_mass_fraction] + humid_air_mixture = pyfluids.Mixture(fluids, fractions) + enforce_gas_phase(humid_air_mixture) + + return humid_air_mixture + + +def compute_compressor_outlet_state_isentropic( + state_inlet: pyfluids.fluids.abstract_fluid.AbstractFluid, + ratio_P: float, +) -> pyfluids.fluids.abstract_fluid.AbstractFluid: + + # compute the outlet state + state_outlet = state_inlet.with_state( + pyfluids.Input.pressure( + ratio_P * state_inlet.pressure + ), # set by pressure ratio of the compressor + pyfluids.Input.entropy( + state_inlet.entropy + ), # by isentropic compression assumption + ) + + return state_outlet + + +def compute_work_rate_isentropic( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, # by default unity => returns unit-mass kJ/kg +): + # work done by a system should be positive + return mass_flow_fluid * ( + state_inlet.enthalpy / 1.0e3 - state_outlet.enthalpy / 1.0e3 + ) # kW = kJ/s + + +def compute_turbine_work_rate( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, + isentropic_efficiency=1.0, +): + work_rate_isentropic = compute_work_rate_isentropic( + state_inlet, + state_outlet, + mass_flow_fluid=mass_flow_fluid, + ) + return isentropic_efficiency * work_rate_isentropic # kW = kJ/s + + +def compute_compressor_work_rate( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg + isentropic_efficiency=1.0, +): + work_rate_isentropic = compute_work_rate_isentropic( + state_inlet, + state_outlet, + mass_flow_fluid=mass_flow_fluid, + ) + return 1.0 / isentropic_efficiency * work_rate_isentropic # kW + + +def compute_heat_transfer( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg +): + # heat transfer to a system should be positive + return mass_flow_fluid * ( + state_outlet.enthalpy / 1.0e3 - state_inlet.enthalpy / 1.0e3 + ) # kJ/s + + +def compute_turbine_outlet_state_isentropic(state_inlet, ratio_P): + + state_outlet = state_inlet.with_state( + pyfluids.Input.pressure( + state_inlet.pressure / ratio_P + ), # set by pressure ratio of the turbine + pyfluids.Input.entropy( + state_inlet.entropy + ), # by isentropic expansion assumption + ) + + return state_outlet + + +class NGCTResult: + + closed: bool = True + desc: str = "" + mass_flowrate: np.ndarray # kg/s + states: dict[int, pyfluids.fluids.abstract_fluid.AbstractFluid] + process_work_unit: dict[tuple[int, int], np.ndarray] + process_heat_unit: dict[tuple[int, int], np.ndarray] + state_names: dict[int, str] + process_names: dict[tuple[int, int], str] + + def __init__( + self, + desc: str = None, + closed: bool = None, + ): + if closed: + self.closed = closed + if desc: + self.desc = desc + + self.mass_flowrate = None # kg/s + self.states = {} + self.process_work_unit = {} + self.process_heat_unit = {} + self.state_names = {} + self.process_names = {} + + def add_state( + self, + index: int, + state: pyfluids.fluids.abstract_fluid.AbstractFluid, + name: str, + ): + # validations + if index in self.states: + raise ValueError(f"index {index} already in states") + if index in self.state_names: + raise ValueError(f"index {index} already in states") + + # add the state + self.states[index] = state + self.state_names[index] = name + + def add_process( + self, + index_in: int, + index_out: int, + work_unit: np.ndarray, # kJ/kg + heat_unit: np.ndarray, # kJ/kg + name: str, + ): + process_index = (index_in, index_out) + # validations + if (index_in not in self.states) or (index_in not in self.state_names): + raise ValueError( + f"index {index_in} must be in states before process " + f"({index_in}, {index_out}) can be added." + ) + if (index_out not in self.states) or (index_out not in self.state_names): + raise ValueError( + f"index {index_out} must be in states before process " + f"({index_in}, {index_out}) can be added." + ) + if process_index in self.states: + raise ValueError(f"process {process_index} already in states") + if process_index in self.state_names: + raise ValueError(f"process {process_index} already in states") + + # add the process + self.process_work_unit[process_index] = work_unit + self.process_heat_unit[process_index] = heat_unit + self.process_names[process_index] = name + + def get_net_work(self): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + return mass_flowrate * sum(self.process_work_unit.values()) # kJ/s + + def get_net_heat_input(self): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + return mass_flowrate * sum( + np.maximum(0.0, list(self.process_heat_unit.values())) + ) # kJ/s + + def get_net_heat_rejection(self): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + return mass_flowrate * sum( + np.minimum(0.0, list(self.process_heat_unit.values())) + ) # kJ/s + + def get_back_work_ratio(self): + back_work = -sum( + np.minimum(0.0, list(self.process_work_unit.values())) + ) # kJ/s + forward_work = sum( + np.maximum(0.0, list(self.process_work_unit.values())) + ) # kJ/s + return back_work / forward_work # - + + def get_efficiency(self): + net_work = self.get_net_work() # kJ/kg + net_heat_input = self.get_net_heat_input() # kJ/kg + return net_work / net_heat_input # - + + def get_fuel_mass_flow(self, LHV_fuel, heating_process=(2, 3)): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + combustor_heat_input_rate = ( + mass_flowrate * self.process_heat_unit[heating_process] + ) # [kg/s]*[kJ/kg] = kJ/s + mass_flowrate_fuel = combustor_heat_input_rate / ( + 1000.0 * LHV_fuel + ) # [J/s]/[J/kg] = kg/s + return mass_flowrate_fuel + + def get_fuel_air_ratio(self, LHV_fuel, heating_process=(2, 3)): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + mass_flowrate_fuel = self.get_fuel_mass_flow( + LHV_fuel, heating_process=heating_process + ) + return mass_flowrate_fuel / mass_flowrate + + def print_states( + self, + ): + for idx_state, state in self.states.items(): + print(f"{(idx_state):>3d}:", end="") + print(f" T={state.temperature:.2f}°C;", end="") + print(f" h={state.enthalpy/1.0e3:.2f} MJ/kg;", end="") + print(f" P={state.pressure/1.0e3:.2f} kPa;", end="") + print() + + def __str__(self): + output = f"\n{self.desc}\n" if self.desc else "\n" + output += "states:\n" + for idx_state, state in self.states.items(): + state_name = self.state_names.get(idx_state, "") + output += f"\t{idx_state} ({state_name}): T={state.temperature:.2f} C, h={state.enthalpy/1.0e3:.2f} MJ/kg, P={state.pressure/1.0e3:.2f} kPa\n" + + output += "\nprocesses:\n" + for (idx_in, idx_out), work_unit in self.process_work_unit.items(): + heat_unit = self.process_heat_unit.get((idx_in, idx_out), 0.0) + process_name = self.process_names.get((idx_in, idx_out), "") + output += f"\t{idx_in} -> {idx_out} ({process_name}): W={work_unit:.2f} kJ/kg, Q={heat_unit:.2f} kJ/kg\n" + + return output + + +class NGCT: + """ + Thermodynamic model for the performance of a natural gas combustion turbine. + + Thermodynamic model for the performance of a specific natural gas combustion turbine. + We assume a Brayton cycle with isentropic efficiency corrections for the compressor and turbine work. + The turbine is assumed to ingest the intake fluid at maximum a specific volumetric flowrate. + We also assume that the combustion chamber and turbine set the maximum designed firing temperature. + """ + + ratio_P: float # -, pressure ratio + Trel_firing: float # deg C, relative temperature + isentropic_efficiency_compressor: float # -, compressor isentropic efficiency + isentropic_efficiency_turbine: float # -, turbine isentropic efficiency + Q_fluid_max: float # m**3/s, design volumetric flowrate + + design_conditions: dict # design conditions dictionary + + def __init__( + self, + ratio_P: float, # -, pressure ratio + Trel_firing: float, # deg C, relative temperature + isentropic_efficiency_compressor: float, # -, compressor isentropic efficiency + isentropic_efficiency_turbine: float, # -, turbine isentropic efficiency + Q_fluid_max: float = None, # m**3/s, design volumetric flowrate; allows unit-mass analysis (units) if None + design_conditions=None, # dict, design conditions dictionary + ): + # set in variables + self.ratio_P = ratio_P + self.Trel_firing = Trel_firing + self.isentropic_efficiency_compressor = isentropic_efficiency_compressor + self.isentropic_efficiency_turbine = isentropic_efficiency_turbine + self.Q_fluid_max = Q_fluid_max + + if design_conditions is not None: + self.design_conditions = design_conditions + + def run_turbine_model( + self, + fluid_ambient_list: list[pyfluids.fluids.abstract_fluid.AbstractFluid] | pyfluids.fluids.abstract_fluid.AbstractFluid, + power_rated: float = np.inf, # kW, maximum rated power of the system + heatrate_fuel_capacity: float = np.inf, # kJ/s, maximum allowable heat rate + ): + + singleton = False # make a flag if there's just one ambient state to run + if isinstance( + fluid_ambient_list, + pyfluids.fluids.abstract_fluid.AbstractFluid, + ): + fluid_ambient_list = [ + fluid_ambient_list, + ] + singleton = True + + results = [] # vector to store results for each ambient fluid state + for fluid_ambient in fluid_ambient_list: # loop over ambient fluid states + + ### CALCULATE THE CYCLE STATES + # note: gas phase enforcement speeds up computation significantly + + enforce_gas_phase( + fluid_ambient + ) # ensure the ambient fluid is in the gas phase + fluid_compressed = compute_compressor_outlet_state_isentropic( + fluid_ambient, + self.ratio_P, + ) # after the compressor + enforce_gas_phase( + fluid_compressed + ) # ensure the compressed fluid is in the gas phase + fluid_combusted = fluid_compressed.with_state( + pyfluids.Input.pressure(fluid_compressed.pressure), + pyfluids.Input.temperature(self.Trel_firing), + ) # after the combustion chamber + enforce_gas_phase( + fluid_combusted + ) # ensure the combusted fluid is in the gas phase + fluid_exhaust = compute_turbine_outlet_state_isentropic( + fluid_combusted, + self.ratio_P, + ) # exhaust from the turbine + enforce_gas_phase( + fluid_exhaust + ) # ensure the exhaust fluid is in the gas phase + + ### CALCULATE THE WORK INPUTS/OUTPUTS + wdot_compressor = compute_compressor_work_rate( + fluid_ambient, + fluid_compressed, + isentropic_efficiency=self.isentropic_efficiency_compressor, + ) # now that the gas states are determined, compute the work + wdot_turbine = compute_turbine_work_rate( + fluid_combusted, + fluid_exhaust, + isentropic_efficiency=self.isentropic_efficiency_turbine, + ) # now that the gas states are determined, compute the work + + ### CALCULATE THE HEAT INPUTS/OUTPUTS + qdot_HX_combustor = compute_heat_transfer( + fluid_compressed, + fluid_combusted, + ) # now that the gas states are determined, compute the heat transfer + qdot_HX_exhaust_rejection = compute_heat_transfer( + fluid_exhaust, + fluid_ambient, + ) # now that the gas states are determined, compute the heat transfer + + ### CALCULATE THE MASS FLOWRATE W/ CONTROL IF APPLICABLE + mass_flowrate = NGCT.get_mass_flowrate( + fluid_ambient, self.Q_fluid_max, + wdot_turbine, wdot_compressor, qdot_HX_combustor, + power_rated=power_rated, heatrate_fuel_capacity=heatrate_fuel_capacity, + ) if self.Q_fluid_max is not None else None # use static method + + ### PACK THE RESULT + result = NGCTResult() + + # add mass flowrate + result.mass_flowrate = mass_flowrate + + # add states + result.add_state(1, fluid_ambient, "ambient intake air") + result.add_state(2, fluid_compressed, "compressed air") + result.add_state(3, fluid_combusted, "products of combustion") + result.add_state(4, fluid_exhaust, "exhaust air") + + # add processes + result.add_process( + 1, + 2, + wdot_compressor, + 0.0, + "isentropic compressor w/ work efficiency correction", + ) + result.add_process( + 2, + 3, + 0.0, + qdot_HX_combustor, + "constant-pressure heating", + ) + result.add_process( + 3, + 4, + wdot_turbine, + 0.0, + "isentropic turbine w/ work efficiency correction", + ) + result.add_process( + 4, + 1, + 0.0, + qdot_HX_exhaust_rejection, + "constant-pressure (atmospheric) cooling", + ) + + # return logic + if singleton: + return result + results.append(result) + + return results + + @staticmethod + def celsius_to_kelvin(degC): + return degC + 273.15 # K + + @staticmethod + def kelvin_to_celsius(degK): + return degK - 273.15 # C + + @staticmethod + def get_mass_flowrate( + fluid_ambient: pyfluids.fluids.abstract_fluid.AbstractFluid, + Q_fluid_max: float, + wdot_turbine: float, + wdot_compressor: float, + qdot_HX_combustor: float, + power_rated: float = np.inf, + heatrate_fuel_capacity: float = np.inf, + ) -> float: + """ + Calculate the mass flowrate constrained by design and operational limits. + + Determines the actual mass flowrate as the minimum of multiple limiting cases: + volumetric flow capacity, power rating, and heat rate capacity. + + Args: + fluid_ambient: Ambient fluid mixture with density property (pyfluids.fluids.abstract_fluid.AbstractFluid) + Q_fluid_max: Maximum volumetric flowrate design limit (m**3/s) + wdot_turbine: Turbine work output (kW/kg) + wdot_compressor: Compressor work input (kW/kg) + qdot_HX_combustor: Heat released by combustor (kW/kg) + power_rated: Generator power rating constraint (kW), optional + heatrate_fuel_capacity: Fuel heat capacity flow limit (kJ/s), optional + + Returns: + Mass flowrate (kg/s) limited by the most restrictive constraint + """ + # maximum mass flowrate allowed by the design + mass_flowrate_max = fluid_ambient.density * Q_fluid_max # kg/s + # make a list of potential limiting cases + mass_flowrate_candidates = [mass_flowrate_max] + + # throttle air/fuel flow to stay under given power rating (i.e. generator) + mass_flowrate_candidates.append( + power_rated / (wdot_turbine + wdot_compressor) # kg/s <= [kW = kJ/s]/[kJ/kg] + ) + + # if flowrate is limited by a MMBtu/s (or kJ/s) flow limit + mass_flowrate_candidates.append( + heatrate_fuel_capacity / qdot_HX_combustor # kg/s <= [kJ/s] / [kJ/kg] + ) + + # print(f"DEBUG!!!!! {'flowrate limited' if np.argmin(mass_flowrate_candidates) == 0 else 'power limited' if np.argmin(mass_flowrate_candidates) == 1 else 'fuel limited'}") + # if np.argmin(mass_flowrate_candidates) == 1: + # print(f"\tDEBUG!!!!! power_rated: {power_rated}") + # print(f"\tDEBUG!!!!! wdot_turbine: {wdot_turbine} wdot_compressor: {wdot_compressor}") + # print(f"\tDEBUG!!!!! power-set mass flowrate: {power_rated/(wdot_turbine + wdot_compressor)}") + + # the actual mass flowrate is the minimum of the limiting cases + return float(np.min(mass_flowrate_candidates)) + + +class GE_7FA05_NGCT(NGCT): + """ + Thermodynamic model for the GE 7FA.05 NGCT in simple-cycle performance. + + A thermodynamic model for the performance of the GE 7FA.05 natural gas combustion turbine in simple-cycle operation. + We use a mix of fact- and data-sheet performance descriptions, textbook assumptions, and reverse-engineering to characterize the system. + Reverse engineering targets matching the ISO performance at the design conditions. + """ + + def __init__(self): + ratio_P = 18.712850988834695 # -, by reverse-engineering from datasheet + Trel_firing = 1300.0 # °C, by textbook assumption + isentropic_efficiency_compressor = 0.85 # -, by textbook assumption + isentropic_efficiency_turbine = 0.90 # -, by textbook assumption + Q_fluid = ( + 477.78734437019483 # m**3/s, by reverse-engineering from datasheet datasheet + ) + + design_conditions = { + "P_ISO": 101325.0, # Pa, ISO conditions + "T_ISO": 288.15, # K, ISO conditions + "rel_humidity_ISO": 60.0, # %, ISO conditions + "eta_th_ISO": 0.385, # -, from GEfactsheet for 7F.05 + "W_net_ISO": 239.0e3, # kW, from GEfactsheet for 7F.05 + } + + # specialize the parent class + super().__init__( + ratio_P, + Trel_firing, + isentropic_efficiency_compressor, + isentropic_efficiency_turbine, + Q_fluid, + design_conditions, + ) + + +class GE_7EA_NGCT_2014(NGCT): + """ + Thermodynamic model for the GE 7EA NGCT in simple-cycle performance. + + A thermodynamic model for the performance of the GE 7EA natural gas combustion turbine in simple-cycle operation. + We use a mix of fact- and data-sheet performance descriptions, textbook assumptions, and reverse-engineering to characterize the system. + Reverse engineering targets matching the ISO performance at the design conditions. + + Matching is middling and the information on the turbine is very old. + """ + + def __init__(self): + ratio_P = 12.6 # -, from NYISO GE factsheet + Trel_firing = 1300.0 # °C, by textbook assumption + isentropic_efficiency_compressor = 0.85 # -, by textbook assumption + isentropic_efficiency_turbine = 0.90 # -, by textbook assumption + Q_fluid = 238.3673469388 # m**3/s, from NYISO GE factsheet + + design_conditions = { + "P_ISO": 101325.0, # Pa, ISO conditions + "T_ISO": 288.15, # K, ISO conditions + "rel_humidity_ISO": 60.0, # %, ISO conditions + "eta_th_ISO": 0.3275, # -, from GEfactsheet for 7F.05 + "W_net_ISO": 85.4e3, # kW, from GEfactsheet for 7F.05 + } + + # specialize the parent class + super().__init__( + ratio_P, + Trel_firing, + isentropic_efficiency_compressor, + isentropic_efficiency_turbine, + Q_fluid, + design_conditions, + ) + + +if __name__ == "__main__": + P_ambient = 101325.0 # Pa + Trel_ambient = 15.0 # Pa + rel_humidity_ambient = 60 + + working_fluid = make_humid_air_mixture( + P_ambient, Trel_ambient, rel_humidity_ambient + ) + fluid_ambient = working_fluid.with_state( + pyfluids.Input.pressure(P_ambient), + pyfluids.Input.temperature(Trel_ambient), + ) + + ngct = GE_7EA_NGCT_2014() + # ngct = GE_7FA05_NGCT() + + result = ngct.run_turbine_model(fluid_ambient) + + net_work = result.get_net_work() + net_heat_input = result.get_net_heat_input() + efficiency = net_work / net_heat_input + exhaust_temperature = result.states[4].temperature + + print(f"net work: {net_work}") + print(f"net heat input: {net_heat_input}") + print(f"efficiency: {efficiency}") + print(f"exhaust temperature: {exhaust_temperature}") diff --git a/h2integrate/converters/combustion_machines/__init__.py b/h2integrate/converters/combustion_machines/__init__.py new file mode 100644 index 000000000..22c71f460 --- /dev/null +++ b/h2integrate/converters/combustion_machines/__init__.py @@ -0,0 +1,3 @@ +from h2integrate.converters.combustion_machines.turbine_simple_cycle import ( + SimpleCycleTurbinePerformanceModel +) \ No newline at end of file diff --git a/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py new file mode 100644 index 000000000..aa973e649 --- /dev/null +++ b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py @@ -0,0 +1,260 @@ +import numpy as np +import pyfluids +import pytest + +from h2integrate.converters.combustion_machines.NGCT_thermo_model import GE_7EA_NGCT_2014 +from h2integrate.converters.combustion_machines.NGCT_thermo_model import GE_7FA05_NGCT +from h2integrate.converters.combustion_machines.NGCT_thermo_model import NGCT +from h2integrate.converters.combustion_machines.NGCT_thermo_model import NGCTResult +from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_compressor_outlet_state_isentropic +from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_compressor_work_rate +from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_heat_transfer +from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_turbine_outlet_state_isentropic +from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_turbine_work_rate +from h2integrate.converters.combustion_machines.NGCT_thermo_model import humidity_ratio_to_water_mass_fraction +from h2integrate.converters.combustion_machines.NGCT_thermo_model import make_humid_air_mixture + + +@pytest.fixture(scope="module") +def iso_conditions(): + return { + "pressure": 101325.0, + "temperature": 15.0, + "relative_humidity": 60.0, + } + + +@pytest.fixture(scope="module") +def iso_ambient_state(iso_conditions): + working_fluid = make_humid_air_mixture( + iso_conditions["pressure"], + iso_conditions["temperature"], + iso_conditions["relative_humidity"], + ) + return working_fluid.with_state( + pyfluids.Input.pressure(iso_conditions["pressure"]), + pyfluids.Input.temperature(iso_conditions["temperature"]), + ) + + +@pytest.fixture(scope="module") +def basic_ngct(): + return NGCT( + ratio_P=15.0, + Trel_firing=1300.0, + isentropic_efficiency_compressor=0.85, + isentropic_efficiency_turbine=0.90, + Q_fluid_max=1.0, + ) + + +@pytest.fixture(scope="module") +def ge_7fa05(): + return GE_7FA05_NGCT() + + +@pytest.fixture(scope="module") +def ge_7ea_2014(): + return GE_7EA_NGCT_2014() + + +@pytest.fixture( + params=[ + {"pressure": 101325.0, "temperature": 0.0, "relative_humidity": 40.0}, + {"pressure": 101325.0, "temperature": 35.0, "relative_humidity": 80.0}, + {"pressure": 80000.0, "temperature": 15.0, "relative_humidity": 60.0}, + ] +) +def varied_ambient_state(request): + conditions = request.param + working_fluid = make_humid_air_mixture( + conditions["pressure"], + conditions["temperature"], + conditions["relative_humidity"], + ) + return working_fluid.with_state( + pyfluids.Input.pressure(conditions["pressure"]), + pyfluids.Input.temperature(conditions["temperature"]), + ) + + +@pytest.mark.unit +def test_humidity_ratio_to_water_mass_fraction_roundtrip(): + humidity_ratio = 0.015 + water_mass_fraction = humidity_ratio_to_water_mass_fraction(humidity_ratio) + recovered_humidity_ratio = water_mass_fraction / (1.0 - water_mass_fraction) + + assert water_mass_fraction == pytest.approx(0.014778325123152709) + assert recovered_humidity_ratio == pytest.approx(humidity_ratio) + + +@pytest.mark.unit +def test_make_humid_air_mixture_builds_valid_state(iso_conditions): + working_fluid = make_humid_air_mixture( + iso_conditions["pressure"], + iso_conditions["temperature"], + iso_conditions["relative_humidity"], + ) + fluid_ambient = working_fluid.with_state( + pyfluids.Input.pressure(iso_conditions["pressure"]), + pyfluids.Input.temperature(iso_conditions["temperature"]), + ) + + assert np.isfinite(fluid_ambient.density) + assert np.isfinite(fluid_ambient.enthalpy) + assert np.isfinite(fluid_ambient.entropy) + assert fluid_ambient.density > 0.0 + assert fluid_ambient.pressure == pytest.approx(iso_conditions["pressure"]) + assert fluid_ambient.temperature == pytest.approx(iso_conditions["temperature"]) + + +@pytest.mark.unit +def test_compressor_outlet_state_isentropic_preserves_entropy(iso_ambient_state): + ratio_p = 12.6 + compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, ratio_p) + + assert compressed.pressure == pytest.approx(iso_ambient_state.pressure * ratio_p) + assert compressed.entropy == pytest.approx(iso_ambient_state.entropy, rel=1e-6, abs=1e-3) + assert compressed.temperature > iso_ambient_state.temperature + + +@pytest.mark.unit +def test_turbine_outlet_state_isentropic_preserves_entropy(iso_ambient_state): + combusted = iso_ambient_state.with_state( + pyfluids.Input.pressure(iso_ambient_state.pressure * 12.6), + pyfluids.Input.temperature(1300.0), + ) + exhaust = compute_turbine_outlet_state_isentropic(combusted, 12.6) + + assert exhaust.pressure == pytest.approx(combusted.pressure / 12.6) + assert exhaust.entropy == pytest.approx(combusted.entropy, rel=1e-6, abs=1e-3) + assert exhaust.temperature < combusted.temperature + + +@pytest.mark.unit +def test_helper_sign_conventions(iso_ambient_state): + compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, 10.0) + combusted = compressed.with_state( + pyfluids.Input.pressure(compressed.pressure), + pyfluids.Input.temperature(1300.0), + ) + exhaust = compute_turbine_outlet_state_isentropic(combusted, 10.0) + + wdot_compressor = compute_compressor_work_rate(iso_ambient_state, compressed, isentropic_efficiency=0.85) + wdot_turbine = compute_turbine_work_rate(combusted, exhaust, isentropic_efficiency=0.90) + qdot_combustor = compute_heat_transfer(compressed, combusted) + qdot_exhaust = compute_heat_transfer(exhaust, iso_ambient_state) + + assert wdot_compressor < 0.0 + assert wdot_turbine > 0.0 + assert qdot_combustor > 0.0 + assert qdot_exhaust < 0.0 + + +@pytest.mark.unit +def test_ngct_result_aggregates_cycle_quantities(iso_ambient_state): + result = NGCTResult(desc="synthetic cycle") + for index in range(1, 5): + result.add_state(index, iso_ambient_state, f"state {index}") + + result.mass_flowrate = 5.0 + result.add_process(1, 2, -10.0, 0.0, "compressor") + result.add_process(2, 3, 0.0, 30.0, "combustor") + result.add_process(3, 4, 20.0, 0.0, "turbine") + result.add_process(4, 1, 0.0, -5.0, "cooler") + + assert result.get_net_work() == pytest.approx(50.0) + assert result.get_net_heat_input() == pytest.approx(150.0) + assert result.get_net_heat_rejection() == pytest.approx(-25.0) + assert result.get_back_work_ratio() == pytest.approx(0.5) + assert result.get_efficiency() == pytest.approx(1.0 / 3.0) + + +@pytest.mark.integration +def test_run_turbine_model_singleton_matches_batch_result(ge_7ea_2014, iso_ambient_state): + result_single = ge_7ea_2014.run_turbine_model(iso_ambient_state) + result_batch = ge_7ea_2014.run_turbine_model([iso_ambient_state])[0] + + assert result_single.mass_flowrate == pytest.approx(result_batch.mass_flowrate) + assert result_single.get_net_work() == pytest.approx(result_batch.get_net_work()) + assert result_single.get_efficiency() == pytest.approx(result_batch.get_efficiency()) + assert result_single.states[4].temperature == pytest.approx(result_batch.states[4].temperature) + + +@pytest.mark.regression +def test_run_turbine_model_applies_minimum_mass_flow_constraint(basic_ngct, iso_ambient_state): + unconstrained = basic_ngct.run_turbine_model(iso_ambient_state) + specific_net_work = unconstrained.get_net_work() / unconstrained.mass_flowrate + specific_heat_input = unconstrained.process_heat_unit[(2, 3)] + + power_limit = 0.75 * unconstrained.get_net_work() + heat_limit = 0.60 * unconstrained.get_net_heat_input() + + constrained = basic_ngct.run_turbine_model( + iso_ambient_state, + power_rated=power_limit, + heatrate_fuel_capacity=heat_limit, + ) + + expected_mass_flowrate = min( + unconstrained.mass_flowrate, + power_limit / specific_net_work, + heat_limit / specific_heat_input, + ) + + assert constrained.mass_flowrate == pytest.approx(expected_mass_flowrate) + + +@pytest.mark.regression +def test_ge_7fa05_iso_design_point_regression(ge_7fa05, iso_ambient_state): + result = ge_7fa05.run_turbine_model(iso_ambient_state) + + assert result.get_efficiency() == pytest.approx( + ge_7fa05.design_conditions["eta_th_ISO"], + rel=0.05, + ) + assert result.get_net_work() == pytest.approx( + ge_7fa05.design_conditions["W_net_ISO"], + rel=0.08, + ) + assert result.states[4].temperature > iso_ambient_state.temperature + + +# @pytest.mark.regression +# def test_ge_7ea_2014_iso_design_point_regression(ge_7ea_2014, iso_ambient_state): +# result = ge_7ea_2014.run_turbine_model(iso_ambient_state) +# +# assert result.get_efficiency() == pytest.approx( +# ge_7ea_2014.design_conditions["eta_th_ISO"], +# rel=0.05, +# ) +# assert result.get_net_work() == pytest.approx( +# ge_7ea_2014.design_conditions["W_net_ISO"], +# rel=0.08, +# ) +# assert result.states[4].temperature > iso_ambient_state.temperature + + +@pytest.mark.integration +def test_ideal_brayton_cycle_energy_balance(iso_ambient_state): + ideal_ngct = NGCT( + ratio_P=15.0, + Trel_firing=1300.0, + isentropic_efficiency_compressor=1.0, + isentropic_efficiency_turbine=1.0, + Q_fluid_max=1.0, + ) + result = ideal_ngct.run_turbine_model(iso_ambient_state) + net_heat = result.get_net_heat_input() + result.get_net_heat_rejection() + + assert result.get_net_work() == pytest.approx(net_heat, rel=1e-9) + + +@pytest.mark.integration +def test_varied_ambient_conditions_produce_finite_outputs(ge_7ea_2014, varied_ambient_state): + result = ge_7ea_2014.run_turbine_model(varied_ambient_state) + + assert np.isfinite(result.mass_flowrate) + assert np.isfinite(result.get_net_work()) + assert np.isfinite(result.get_efficiency()) + assert 0.0 < result.get_efficiency() < 1.0 diff --git a/h2integrate/converters/combustion_machines/turbine_simple_cycle.py b/h2integrate/converters/combustion_machines/turbine_simple_cycle.py new file mode 100644 index 000000000..6e939f701 --- /dev/null +++ b/h2integrate/converters/combustion_machines/turbine_simple_cycle.py @@ -0,0 +1,357 @@ +from attrs import field, define + +import numpy as np + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.validators import gt_zero, gte_zero +from h2integrate.core.model_baseclasses import PerformanceModelBaseClass + +import h2integrate.converters.combustion_machines.NGCT_thermo_model as NGCT + + +@define(kw_only=True) +class SimpleCycleTurbinePerformanceConfig(BaseConfig): + """ + Configuration class for a simple Brayton-cycle turbine performance model. + + This configuration class handles the parameters for generic combustion + turbines assuming simple Brayton-cycle thermodynamics models with isentropic + efficiency corrections for compressor and turbine work, and compensation + for ambient conditions. + + Attributes: + ----------- + flowrate_max_fluid_cubic_m_per_s (float): maximum volumetric flowrate of working fluid through the turbine in m^3/s. + + system_capacity_mw (float): rated system capacity in MW. + + generator_efficiency (float): efficiency of the generator connected to the turbine. + + pressure_ratio (float): pressure ratio generated by the compressor. + + firing_temp_C (float): firing temperature of the turbine in degrees Celsius. + + isentropic_efficiency_compressor (float): isentropic efficiency of the compressor. + + isentropic_efficiency_turbine (float): isentropic efficiency of the turbine. + """ + + flowrate_max_fluid_cubic_m_per_s: float = field(validator=gt_zero) + system_capacity_mw: float = field(validator=gte_zero) + generator_efficiency: float = field(validator=gt_zero) + pressure_ratio: float = field(validator=gt_zero) + firing_temp_C: float = field(validator=gt_zero) + isentropic_efficiency_compressor: float = field(validator=gt_zero) + isentropic_efficiency_turbine: float = field(validator=gt_zero) + + +class SimpleCycleTurbinePerformanceModel(PerformanceModelBaseClass): + """ + Performance model for simple Brayton-cycle turbines. + + This model will compute the electricity output from a simple Brayton-cycle + combustion turbine directly connected at the shaft to a generator, using + real-time or provided ambient conditions and a working fluid and fuel of + choice. + + Inputs: + ??? + + Outputs: + ??? + """ + + _time_step_bounds = ( + 300, # 5 min. run + 3600, # 1 hr. run + ) # (min, max) time step lengths (in seconds) compatible with this model + _control_classifier = "dispatchable" + + def initialize(self, fuel_source="natural_gas"): + super().initialize() + self.commodity = "electricity" + self.commodity_rate_units = "MW" + self.commodity_amount_units = "MW*h" + + self.fuel_source = fuel_source + + def setup(self): + super().setup() + + self.config = SimpleCycleTurbinePerformanceConfig.from_dict( + merge_shared_inputs( + self.options["tech_config"]["model_inputs"], "performance" + ), + additional_cls_name=self.__class__.__name__, + ) + + # add solar resource data to access ambient conditions consistently + self.add_discrete_input( + "solar_resource_data", + val={}, + desc="Solar resource data dictionary", + ) + + # add natural gas consumed output + self.add_output( + f"{self.fuel_source}_consumed", + val=0.0, + shape=self.n_timesteps, + units="kJ/s", + desc="natural gas consumed by the plant", + ) + + # add system max rated capacity as an input w/ config value as default + self.add_input( + f"system_capacity", + val=self.config.system_capacity_mw, + units="MW", + desc="rated generator capacity in MW", + ) + + # default the electricity command value input as the rated capacity + self.add_input( + f"{self.commodity}_command_value", + val=self.config.system_capacity_mw, + shape=self.n_timesteps, + units=self.commodity_rate_units, + desc="electricity command value for natural gas plant", + ) + + # add max volumetric flowrate as an input w/ config value as default + self.add_input( + f"flowrate_max_fluid", + val=self.config.flowrate_max_fluid_cubic_m_per_s, + units="m**3/s", + desc="maximum volumetric flowrate of working fluid through the turbine", + ) + + # add fuel input, default to 0 --> set using feedstock component + self.add_input( + f"{self.fuel_source}_in", + val=0.0, + shape=self.n_timesteps, + units="kJ/s", + desc="fuel input energy to the turbine/plant in MMBtu/h", + ) + + self.add_output( + f"unmet_{self.commodity}_demand", + val=0.0, + shape=self.n_timesteps, + units=self.commodity_rate_units, + desc="unmet demand for electricity from the turbine/plant in MW", + ) + + self.ngct = NGCT.NGCT( + ratio_P=self.config.pressure_ratio, + Trel_firing=self.config.firing_temp_C, + isentropic_efficiency_compressor=self.config.isentropic_efficiency_compressor, + isentropic_efficiency_turbine=self.config.isentropic_efficiency_turbine, + Q_fluid_max=1.0, # (explicitly) unit-mass analysis + ) + + # try to run NGCT cycle here expensively, and all compute variations + # should be derivable w/o major new computation after that + P_ambient = 101325.0 # Pa + T_ambient = 15.0 # °C + self.working_fluid = NGCT.pyfluids.Fluid(NGCT.pyfluids.FluidsList.Air) + wx_data = [(T_ambient, P_ambient)] + self.wx_checksums = [hash(wxd) for wxd in wx_data] + self.ambient_fluid_list = [self.working_fluid.with_state( + NGCT.pyfluids.Input.pressure(P), + NGCT.pyfluids.Input.temperature(T), + ) for T, P in wx_data] + self.fundamental_cycle = self.ngct.run_turbine_model(self.ambient_fluid_list) + + def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): + + # working variables for inputs + system_capacity_mw = float(inputs["system_capacity"][0]) + flowrate_max_fluid_m3_per_s = float(inputs["flowrate_max_fluid"][0]) + fuel_source_heating_in_kJ_per_s = inputs[f"{self.fuel_source}_in"] + + # extract from the solar resource data + temperature_degC = discrete_inputs["solar_resource_data"]["temperature"] + pressure_Pa = discrete_inputs["solar_resource_data"]["pressure"]*100.0 # convert hPa to Pa + wx_data = list(zip(temperature_degC.tolist(), pressure_Pa.tolist())) + wx_checksums = [hash(wxd) for wxd in wx_data] + + # figure out if we can re-use the cache + if not all([a == b for a, b in zip(wx_checksums, self.wx_checksums)]): + self.wx_checksums = None # break everything that will be recomputed + self.ambient_fluid_list = None # break everything that will be recomputed + self.fundamental_cycle = None # break everything that will be recomputed + + self.wx_checksums = wx_checksums + self.ambient_fluid_list = [self.working_fluid.with_state( + NGCT.pyfluids.Input.pressure(P), + NGCT.pyfluids.Input.temperature(T), + ) for T, P in wx_data] + self.fundamental_cycle = self.ngct.run_turbine_model(self.ambient_fluid_list) + # print("DEBUG!!!!! RE-COMPUTING UNIT-MASS CYCLE B/C AMBIENT FLUID STATES CHANGED") + else: + pass + # print("DEBUG!!!!! RE-USING ORIGINAL CYCLE (NO AMBIENT FLUID CHANGE)") + + # extract the net unit mass net work + unit_mass_net_work_vec = [(result.process_work_unit[(3,4)] + result.process_work_unit[(1,2)]) for result in self.fundamental_cycle] # kJ/kg + unit_mass_net_heat_input_vec = [result.process_heat_unit[(2,3)] for result in self.fundamental_cycle] # kJ/kg + eta_th_vec = [result.get_efficiency() for result in self.fundamental_cycle] # - + + # compute the rating/flowrate/input heat limited flowrate + mass_flowrates = [NGCT.NGCT.get_mass_flowrate( + fluid_ambient=fluid_ambient, + Q_fluid_max=flowrate_max_fluid_m3_per_s, + wdot_turbine=result.process_work_unit[(3,4)], + wdot_compressor=result.process_work_unit[(1,2)], + qdot_HX_combustor=result.process_heat_unit[(2,3)], + power_rated=system_capacity_mw*1000.0, # expects kW, convert from MW + heatrate_fuel_capacity=fuel_source_in, + ) for fluid_ambient, result, fuel_source_in in zip( + self.ambient_fluid_list, self.fundamental_cycle, fuel_source_heating_in_kJ_per_s.squeeze().tolist(), + )] + # volume_flowrates = [m/fluid_ambient.density for m, fluid_ambient in zip(mass_flowrates, self.ambient_fluid_list)] # m^3/s + + net_work_vec = [m*w/1000.0 for m, w in zip(mass_flowrates, unit_mass_net_work_vec)] # MW + net_heat_input_vec = [m*q for m, q in zip(mass_flowrates, unit_mass_net_heat_input_vec)] # kJ/s + + # import matplotlib.pyplot as plt + # fig, axes = plt.subplots(5, 1, sharex=True) + # axes[0].plot(temperature_degC) + # axes[1].plot(pressure_Pa) + # axes[2].plot(net_work_vec) + # axes[2].axhline(system_capacity_mw, c='r', linestyle="--") + # axes[3].plot(volume_flowrates) + # axes[3].axhline(flowrate_max_fluid_m3_per_s, c='r', linestyle="--") + # axes[4].plot(mass_flowrates) + # plt.show() + + generator_efficiency = self.config.generator_efficiency + electricity_out = generator_efficiency*np.array(net_work_vec) # MW + + outputs["electricity_out"] = electricity_out + outputs[f"{self.fuel_source}_consumed"] = net_heat_input_vec + + outputs["rated_electricity_production"] = system_capacity_mw # QUESTION!!!!! WHAT IS RATED??? + + max_production = inputs["system_capacity"] * len(electricity_out) * self.dt / 3600.0 + + outputs["total_electricity_produced"] = np.sum(electricity_out) * self.dt / 3600.0 + outputs["capacity_factor"] = outputs["total_electricity_produced"].sum() / max_production + outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( + 1 / self.fraction_of_year_simulated + ) + outputs[f"unmet_{self.commodity}_demand"] = inputs["electricity_command_value"] - electricity_out + + +if __name__ == "__main__": + from h2integrate import H2IntegrateModel + + # Create an H2I model + driver_config = { + "name": "driver_config", + "description": "this analysis runs a brayton-cycle NG power plant", + "general": {"folder_output": "outputs"}, + } + technology_config = { + "name": "detailed NG plant", + "description": "a detailed NG plant for brayton-cycle analysis", + "technologies": { + "ng_feedstock": { + "performance_model": { + "model": "FeedstockPerformanceModel", + }, + "cost_model": { + "model": "FeedstockCostModel", + }, + "model_inputs": { + "shared_parameters": { + "commodity": "natural_gas", + "commodity_rate_units": "MMBtu/h", + }, + "performance_parameters": { + "rated_capacity": 3000.0, + }, + "cost_parameters": { + "cost_year": 2023, + "commodity_amount_units": "MMBtu", + "price": 4.2, + "annual_cost": 0.0, + "start_up_cost": 100000.0, + }, + }, + }, + "ng": { + "performance_model": { + "model": "SimpleCycleTurbinePerformanceModel", + }, + "cost_model": { + "model": "NaturalGasCostModel", + }, + "model_inputs": { + "performance_parameters": { + "system_capacity_mw": 239.0, + "firing_temp_C": 1300.0, + "pressure_ratio": 18.712850988834695, + "flowrate_max_fluid_cubic_m_per_s": 477.78734437019483, + "isentropic_efficiency_compressor": 0.85, + "isentropic_efficiency_turbine": 0.90, + "generator_efficiency": 1.0, + }, + }, + }, + }, + } + plant_config = { + "name": "brayton plant", + "description": "a brayton-cycle plant located somewhere in Texas", + "sites": { + "site": { + "latitude": 34.382308, + "longitude": -101.816607, + "resources": { + "solar_resource": { + "resource_model": "GOESAggregatedSolarAPI", + "resource_parameters": { + "resource_year": 2024, + "resource_dir": "../../../examples/11_hybrid_energy_plant/tech_inputs/weather/solar", + "resource_filename": "30.6617_-101.7096_psmv3_60_2013.csv", + }, + }, + }, + }, + }, + "plant": { + "plant_life": 30, + "simulation": { + "n_timesteps": 8760, + "dt": 3600, + }, + }, + "technology_interconnections": [ + ["ng_feedstock", "ng", "natural_gas", "pipe"], + ], + "resource_to_tech_connections": [ + ["site.solar_resource", "ng", "solar_resource_data"], + ] + } + + h2i = H2IntegrateModel({ + "name": "brayton", + "system_summary": "a brayon-cycle NG plant sim using real time ambient condition data", + "driver_config": driver_config, + "technology_config": technology_config, + "plant_config": plant_config, + }) + + # Run the model + h2i.run() + + # show N2 diagram + if False: + import openmdao.api as om + om.n2(h2i.model) + + # Post-process the results + h2i.post_process() diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index f9e43ce89..d61721cb7 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -129,6 +129,7 @@ def copy(self): "QuinnNuclearPerformanceModel": "converters.nuclear:QuinnNuclearPerformanceModel", "QuinnNuclearCostModel": "converters.nuclear:QuinnNuclearCostModel", "NaturalGasCostModel": "converters.natural_gas:NaturalGasCostModel", + "SimpleCycleTurbinePerformanceModel": "converters.combustion_machines:SimpleCycleTurbinePerformanceModel", # Transport "cable": "transporters:CablePerformanceModel", "pipe": "transporters:PipePerformanceModel", From e4da05526f8e22375e8ec6265bb7432a4f9b2209 Mon Sep 17 00:00:00 2001 From: Cory Frontin Date: Mon, 20 Jul 2026 11:27:02 -0600 Subject: [PATCH 2/6] updated ngct for cleanliness --- .../combustion_machines/NGCT_thermo_model.py | 165 ++++++++---------- .../turbine_simple_cycle.py | 150 ++++++++++------ 2 files changed, 166 insertions(+), 149 deletions(-) diff --git a/h2integrate/converters/combustion_machines/NGCT_thermo_model.py b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py index b7fc3ee0d..6eadbc3cc 100644 --- a/h2integrate/converters/combustion_machines/NGCT_thermo_model.py +++ b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py @@ -1,15 +1,13 @@ -from pprint import pprint - -from matplotlib.pylab import power import numpy as np - import pyfluids + DRY_AIR_MASS_FRACTIONS = { - pyfluids.FluidsList.Nitrogen: 75.51760366879898, # generated by copilot, matches engineeringtoolbox.com data - pyfluids.FluidsList.Oxygen: 23.139561050152018, # generated by copilot, matches engineeringtoolbox.com data - pyfluids.FluidsList.Argon: 1.2881375564363484, # generated by copilot, matches engineeringtoolbox.com data - pyfluids.FluidsList.CarbonDioxide: 0.05469772461264626, # generated by copilot, matches engineeringtoolbox.com data + # generated by copilot, matches engineeringtoolbox.com data + pyfluids.FluidsList.Nitrogen: 75.51760366879898, + pyfluids.FluidsList.Oxygen: 23.139561050152018, + pyfluids.FluidsList.Argon: 1.2881375564363484, + pyfluids.FluidsList.CarbonDioxide: 0.05469772461264626, } @@ -23,11 +21,7 @@ def enforce_gas_phase(fluid): if phase_gas is None: return fluid - try: - fluid.specify_phase(phase_gas) - except Exception: - # Keep behavior backward compatible if phase forcing is not supported. - pass + fluid.specify_phase(phase_gas) return fluid @@ -36,7 +30,6 @@ def make_humid_air_mixture( temp_rel: float, rel_humidity: float, ) -> pyfluids.Mixture: - # get absolute mass humidity ratio using conditions w/ pyfluid humidair humidity_ratio = ( pyfluids.HumidAir() @@ -53,7 +46,7 @@ def make_humid_air_mixture( dry_air_mass_fraction = 100.0 - water_mass_fraction # create a mixture using our dry air mix plus water vapor - fluids = list(DRY_AIR_MASS_FRACTIONS.keys()) + [pyfluids.FluidsList.Water] + fluids = [*list(DRY_AIR_MASS_FRACTIONS.keys()), pyfluids.FluidsList.Water] fractions = [ dry_air_mass_fraction * mass_fraction / 100.0 for mass_fraction in DRY_AIR_MASS_FRACTIONS.values() @@ -68,15 +61,12 @@ def compute_compressor_outlet_state_isentropic( state_inlet: pyfluids.fluids.abstract_fluid.AbstractFluid, ratio_P: float, ) -> pyfluids.fluids.abstract_fluid.AbstractFluid: - # compute the outlet state state_outlet = state_inlet.with_state( pyfluids.Input.pressure( ratio_P * state_inlet.pressure ), # set by pressure ratio of the compressor - pyfluids.Input.entropy( - state_inlet.entropy - ), # by isentropic compression assumption + pyfluids.Input.entropy(state_inlet.entropy), # by isentropic compression assumption ) return state_outlet @@ -127,27 +117,21 @@ def compute_heat_transfer( mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg ): # heat transfer to a system should be positive - return mass_flow_fluid * ( - state_outlet.enthalpy / 1.0e3 - state_inlet.enthalpy / 1.0e3 - ) # kJ/s + return mass_flow_fluid * (state_outlet.enthalpy / 1.0e3 - state_inlet.enthalpy / 1.0e3) # kJ/s def compute_turbine_outlet_state_isentropic(state_inlet, ratio_P): - state_outlet = state_inlet.with_state( pyfluids.Input.pressure( state_inlet.pressure / ratio_P ), # set by pressure ratio of the turbine - pyfluids.Input.entropy( - state_inlet.entropy - ), # by isentropic expansion assumption + pyfluids.Input.entropy(state_inlet.entropy), # by isentropic expansion assumption ) return state_outlet class NGCTResult: - closed: bool = True desc: str = "" mass_flowrate: np.ndarray # kg/s @@ -159,8 +143,8 @@ class NGCTResult: def __init__( self, - desc: str = None, - closed: bool = None, + desc: str | None = None, + closed: bool | None = None, ): if closed: self.closed = closed @@ -226,23 +210,15 @@ def get_net_work(self): def get_net_heat_input(self): mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - return mass_flowrate * sum( - np.maximum(0.0, list(self.process_heat_unit.values())) - ) # kJ/s + return mass_flowrate * sum(np.maximum(0.0, list(self.process_heat_unit.values()))) # kJ/s def get_net_heat_rejection(self): mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - return mass_flowrate * sum( - np.minimum(0.0, list(self.process_heat_unit.values())) - ) # kJ/s + return mass_flowrate * sum(np.minimum(0.0, list(self.process_heat_unit.values()))) # kJ/s def get_back_work_ratio(self): - back_work = -sum( - np.minimum(0.0, list(self.process_work_unit.values())) - ) # kJ/s - forward_work = sum( - np.maximum(0.0, list(self.process_work_unit.values())) - ) # kJ/s + back_work = -sum(np.minimum(0.0, list(self.process_work_unit.values()))) # kJ/s + forward_work = sum(np.maximum(0.0, list(self.process_work_unit.values()))) # kJ/s return back_work / forward_work # - def get_efficiency(self): @@ -255,16 +231,12 @@ def get_fuel_mass_flow(self, LHV_fuel, heating_process=(2, 3)): combustor_heat_input_rate = ( mass_flowrate * self.process_heat_unit[heating_process] ) # [kg/s]*[kJ/kg] = kJ/s - mass_flowrate_fuel = combustor_heat_input_rate / ( - 1000.0 * LHV_fuel - ) # [J/s]/[J/kg] = kg/s + mass_flowrate_fuel = combustor_heat_input_rate / (1000.0 * LHV_fuel) # [J/s]/[J/kg] = kg/s return mass_flowrate_fuel def get_fuel_air_ratio(self, LHV_fuel, heating_process=(2, 3)): mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - mass_flowrate_fuel = self.get_fuel_mass_flow( - LHV_fuel, heating_process=heating_process - ) + mass_flowrate_fuel = self.get_fuel_mass_flow(LHV_fuel, heating_process=heating_process) return mass_flowrate_fuel / mass_flowrate def print_states( @@ -282,13 +254,19 @@ def __str__(self): output += "states:\n" for idx_state, state in self.states.items(): state_name = self.state_names.get(idx_state, "") - output += f"\t{idx_state} ({state_name}): T={state.temperature:.2f} C, h={state.enthalpy/1.0e3:.2f} MJ/kg, P={state.pressure/1.0e3:.2f} kPa\n" + output += ( + f"\t{idx_state} ({state_name}): T={state.temperature:.2f} C, " + f"h={state.enthalpy/1.0e3:.2f} MJ/kg, P={state.pressure/1.0e3:.2f} kPa\n" + ) output += "\nprocesses:\n" for (idx_in, idx_out), work_unit in self.process_work_unit.items(): heat_unit = self.process_heat_unit.get((idx_in, idx_out), 0.0) process_name = self.process_names.get((idx_in, idx_out), "") - output += f"\t{idx_in} -> {idx_out} ({process_name}): W={work_unit:.2f} kJ/kg, Q={heat_unit:.2f} kJ/kg\n" + output += ( + f"\t{idx_in} -> {idx_out} ({process_name}): " + f"W={work_unit:.2f} kJ/kg, Q={heat_unit:.2f} kJ/kg\n" + ) return output @@ -297,10 +275,12 @@ class NGCT: """ Thermodynamic model for the performance of a natural gas combustion turbine. - Thermodynamic model for the performance of a specific natural gas combustion turbine. - We assume a Brayton cycle with isentropic efficiency corrections for the compressor and turbine work. - The turbine is assumed to ingest the intake fluid at maximum a specific volumetric flowrate. - We also assume that the combustion chamber and turbine set the maximum designed firing temperature. + Thermodynamic model for the performance of a specific natural gas combustion + turbine. We assume a Brayton cycle with isentropic efficiency corrections + for the compressor and turbine work. The turbine is assumed to ingest the + intake fluid at maximum a specific volumetric flowrate. We also assume that + the combustion chamber and turbine set the maximum designed firing + temperature. """ ratio_P: float # -, pressure ratio @@ -317,7 +297,8 @@ def __init__( Trel_firing: float, # deg C, relative temperature isentropic_efficiency_compressor: float, # -, compressor isentropic efficiency isentropic_efficiency_turbine: float, # -, turbine isentropic efficiency - Q_fluid_max: float = None, # m**3/s, design volumetric flowrate; allows unit-mass analysis (units) if None + Q_fluid_max: float + | None = None, # m**3/s, design volumetric flowrate; unit-mass analysis if None design_conditions=None, # dict, design conditions dictionary ): # set in variables @@ -332,11 +313,11 @@ def __init__( def run_turbine_model( self, - fluid_ambient_list: list[pyfluids.fluids.abstract_fluid.AbstractFluid] | pyfluids.fluids.abstract_fluid.AbstractFluid, + fluid_ambient_list: list[pyfluids.fluids.abstract_fluid.AbstractFluid] + | pyfluids.fluids.abstract_fluid.AbstractFluid, power_rated: float = np.inf, # kW, maximum rated power of the system heatrate_fuel_capacity: float = np.inf, # kJ/s, maximum allowable heat rate ): - singleton = False # make a flag if there's just one ambient state to run if isinstance( fluid_ambient_list, @@ -349,34 +330,25 @@ def run_turbine_model( results = [] # vector to store results for each ambient fluid state for fluid_ambient in fluid_ambient_list: # loop over ambient fluid states - ### CALCULATE THE CYCLE STATES # note: gas phase enforcement speeds up computation significantly - enforce_gas_phase( - fluid_ambient - ) # ensure the ambient fluid is in the gas phase + enforce_gas_phase(fluid_ambient) # ensure the ambient fluid is in the gas phase fluid_compressed = compute_compressor_outlet_state_isentropic( fluid_ambient, self.ratio_P, ) # after the compressor - enforce_gas_phase( - fluid_compressed - ) # ensure the compressed fluid is in the gas phase + enforce_gas_phase(fluid_compressed) # ensure the compressed fluid is in the gas phase fluid_combusted = fluid_compressed.with_state( pyfluids.Input.pressure(fluid_compressed.pressure), pyfluids.Input.temperature(self.Trel_firing), ) # after the combustion chamber - enforce_gas_phase( - fluid_combusted - ) # ensure the combusted fluid is in the gas phase + enforce_gas_phase(fluid_combusted) # ensure the combusted fluid is in the gas phase fluid_exhaust = compute_turbine_outlet_state_isentropic( fluid_combusted, self.ratio_P, ) # exhaust from the turbine - enforce_gas_phase( - fluid_exhaust - ) # ensure the exhaust fluid is in the gas phase + enforce_gas_phase(fluid_exhaust) # ensure the exhaust fluid is in the gas phase ### CALCULATE THE WORK INPUTS/OUTPUTS wdot_compressor = compute_compressor_work_rate( @@ -401,11 +373,19 @@ def run_turbine_model( ) # now that the gas states are determined, compute the heat transfer ### CALCULATE THE MASS FLOWRATE W/ CONTROL IF APPLICABLE - mass_flowrate = NGCT.get_mass_flowrate( - fluid_ambient, self.Q_fluid_max, - wdot_turbine, wdot_compressor, qdot_HX_combustor, - power_rated=power_rated, heatrate_fuel_capacity=heatrate_fuel_capacity, - ) if self.Q_fluid_max is not None else None # use static method + mass_flowrate = ( + NGCT.get_mass_flowrate( + fluid_ambient, + self.Q_fluid_max, + wdot_turbine, + wdot_compressor, + qdot_HX_combustor, + power_rated=power_rated, + heatrate_fuel_capacity=heatrate_fuel_capacity, + ) + if self.Q_fluid_max is not None + else None + ) # use static method ### PACK THE RESULT result = NGCTResult() @@ -477,11 +457,12 @@ def get_mass_flowrate( """ Calculate the mass flowrate constrained by design and operational limits. - Determines the actual mass flowrate as the minimum of multiple limiting cases: - volumetric flow capacity, power rating, and heat rate capacity. + Determines the actual mass flowrate as the minimum of multiple limiting + cases: volumetric flow capacity, power rating, and heat rate capacity. Args: - fluid_ambient: Ambient fluid mixture with density property (pyfluids.fluids.abstract_fluid.AbstractFluid) + fluid_ambient: Ambient fluid mixture with density property + (pyfluids.fluids.abstract_fluid.AbstractFluid) Q_fluid_max: Maximum volumetric flowrate design limit (m**3/s) wdot_turbine: Turbine work output (kW/kg) wdot_compressor: Compressor work input (kW/kg) @@ -507,12 +488,6 @@ def get_mass_flowrate( heatrate_fuel_capacity / qdot_HX_combustor # kg/s <= [kJ/s] / [kJ/kg] ) - # print(f"DEBUG!!!!! {'flowrate limited' if np.argmin(mass_flowrate_candidates) == 0 else 'power limited' if np.argmin(mass_flowrate_candidates) == 1 else 'fuel limited'}") - # if np.argmin(mass_flowrate_candidates) == 1: - # print(f"\tDEBUG!!!!! power_rated: {power_rated}") - # print(f"\tDEBUG!!!!! wdot_turbine: {wdot_turbine} wdot_compressor: {wdot_compressor}") - # print(f"\tDEBUG!!!!! power-set mass flowrate: {power_rated/(wdot_turbine + wdot_compressor)}") - # the actual mass flowrate is the minimum of the limiting cases return float(np.min(mass_flowrate_candidates)) @@ -521,9 +496,11 @@ class GE_7FA05_NGCT(NGCT): """ Thermodynamic model for the GE 7FA.05 NGCT in simple-cycle performance. - A thermodynamic model for the performance of the GE 7FA.05 natural gas combustion turbine in simple-cycle operation. - We use a mix of fact- and data-sheet performance descriptions, textbook assumptions, and reverse-engineering to characterize the system. - Reverse engineering targets matching the ISO performance at the design conditions. + A thermodynamic model for the performance of the GE 7FA.05 natural gas + combustion turbine in simple-cycle operation. We use a mix of fact- and + data-sheet performance descriptions, textbook assumptions, and + reverse-engineering to characterize the system. Reverse engineering targets + matching the ISO performance at the design conditions. """ def __init__(self): @@ -531,9 +508,7 @@ def __init__(self): Trel_firing = 1300.0 # °C, by textbook assumption isentropic_efficiency_compressor = 0.85 # -, by textbook assumption isentropic_efficiency_turbine = 0.90 # -, by textbook assumption - Q_fluid = ( - 477.78734437019483 # m**3/s, by reverse-engineering from datasheet datasheet - ) + Q_fluid = 477.78734437019483 # m**3/s, by reverse-engineering from datasheet datasheet design_conditions = { "P_ISO": 101325.0, # Pa, ISO conditions @@ -558,9 +533,11 @@ class GE_7EA_NGCT_2014(NGCT): """ Thermodynamic model for the GE 7EA NGCT in simple-cycle performance. - A thermodynamic model for the performance of the GE 7EA natural gas combustion turbine in simple-cycle operation. - We use a mix of fact- and data-sheet performance descriptions, textbook assumptions, and reverse-engineering to characterize the system. - Reverse engineering targets matching the ISO performance at the design conditions. + A thermodynamic model for the performance of the GE 7EA natural gas + combustion turbine in simple-cycle operation. We use a mix of fact- and + data-sheet performance descriptions, textbook assumptions, and + reverse-engineering to characterize the system. Reverse engineering targets + matching the ISO performance at the design conditions. Matching is middling and the information on the turbine is very old. """ @@ -596,9 +573,7 @@ def __init__(self): Trel_ambient = 15.0 # Pa rel_humidity_ambient = 60 - working_fluid = make_humid_air_mixture( - P_ambient, Trel_ambient, rel_humidity_ambient - ) + working_fluid = make_humid_air_mixture(P_ambient, Trel_ambient, rel_humidity_ambient) fluid_ambient = working_fluid.with_state( pyfluids.Input.pressure(P_ambient), pyfluids.Input.temperature(Trel_ambient), diff --git a/h2integrate/converters/combustion_machines/turbine_simple_cycle.py b/h2integrate/converters/combustion_machines/turbine_simple_cycle.py index 6e939f701..2de6a4bd0 100644 --- a/h2integrate/converters/combustion_machines/turbine_simple_cycle.py +++ b/h2integrate/converters/combustion_machines/turbine_simple_cycle.py @@ -1,13 +1,11 @@ -from attrs import field, define - import numpy as np +from attrs import field, define +import h2integrate.converters.combustion_machines.NGCT_thermo_model as NGCT from h2integrate.core.utilities import BaseConfig, merge_shared_inputs from h2integrate.core.validators import gt_zero, gte_zero from h2integrate.core.model_baseclasses import PerformanceModelBaseClass -import h2integrate.converters.combustion_machines.NGCT_thermo_model as NGCT - @define(kw_only=True) class SimpleCycleTurbinePerformanceConfig(BaseConfig): @@ -21,7 +19,8 @@ class SimpleCycleTurbinePerformanceConfig(BaseConfig): Attributes: ----------- - flowrate_max_fluid_cubic_m_per_s (float): maximum volumetric flowrate of working fluid through the turbine in m^3/s. + flowrate_max_fluid_cubic_m_per_s (float): maximum volumetric flowrate of + working fluid through the turbine in m^3/s. system_capacity_mw (float): rated system capacity in MW. @@ -79,9 +78,7 @@ def setup(self): super().setup() self.config = SimpleCycleTurbinePerformanceConfig.from_dict( - merge_shared_inputs( - self.options["tech_config"]["model_inputs"], "performance" - ), + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), additional_cls_name=self.__class__.__name__, ) @@ -103,7 +100,7 @@ def setup(self): # add system max rated capacity as an input w/ config value as default self.add_input( - f"system_capacity", + "system_capacity", val=self.config.system_capacity_mw, units="MW", desc="rated generator capacity in MW", @@ -120,7 +117,7 @@ def setup(self): # add max volumetric flowrate as an input w/ config value as default self.add_input( - f"flowrate_max_fluid", + "flowrate_max_fluid", val=self.config.flowrate_max_fluid_cubic_m_per_s, units="m**3/s", desc="maximum volumetric flowrate of working fluid through the turbine", @@ -158,14 +155,16 @@ def setup(self): self.working_fluid = NGCT.pyfluids.Fluid(NGCT.pyfluids.FluidsList.Air) wx_data = [(T_ambient, P_ambient)] self.wx_checksums = [hash(wxd) for wxd in wx_data] - self.ambient_fluid_list = [self.working_fluid.with_state( - NGCT.pyfluids.Input.pressure(P), - NGCT.pyfluids.Input.temperature(T), - ) for T, P in wx_data] + self.ambient_fluid_list = [ + self.working_fluid.with_state( + NGCT.pyfluids.Input.pressure(P), + NGCT.pyfluids.Input.temperature(T), + ) + for T, P in wx_data + ] self.fundamental_cycle = self.ngct.run_turbine_model(self.ambient_fluid_list) def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): - # working variables for inputs system_capacity_mw = float(inputs["system_capacity"][0]) flowrate_max_fluid_m3_per_s = float(inputs["flowrate_max_fluid"][0]) @@ -173,21 +172,26 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): # extract from the solar resource data temperature_degC = discrete_inputs["solar_resource_data"]["temperature"] - pressure_Pa = discrete_inputs["solar_resource_data"]["pressure"]*100.0 # convert hPa to Pa + pressure_Pa = ( + discrete_inputs["solar_resource_data"]["pressure"] * 100.0 + ) # convert hPa to Pa wx_data = list(zip(temperature_degC.tolist(), pressure_Pa.tolist())) wx_checksums = [hash(wxd) for wxd in wx_data] # figure out if we can re-use the cache - if not all([a == b for a, b in zip(wx_checksums, self.wx_checksums)]): + if not all(a == b for a, b in zip(wx_checksums, self.wx_checksums)): self.wx_checksums = None # break everything that will be recomputed self.ambient_fluid_list = None # break everything that will be recomputed self.fundamental_cycle = None # break everything that will be recomputed self.wx_checksums = wx_checksums - self.ambient_fluid_list = [self.working_fluid.with_state( - NGCT.pyfluids.Input.pressure(P), - NGCT.pyfluids.Input.temperature(T), - ) for T, P in wx_data] + self.ambient_fluid_list = [ + self.working_fluid.with_state( + NGCT.pyfluids.Input.pressure(P), + NGCT.pyfluids.Input.temperature(T), + ) + for T, P in wx_data + ] self.fundamental_cycle = self.ngct.run_turbine_model(self.ambient_fluid_list) # print("DEBUG!!!!! RE-COMPUTING UNIT-MASS CYCLE B/C AMBIENT FLUID STATES CHANGED") else: @@ -195,45 +199,71 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): # print("DEBUG!!!!! RE-USING ORIGINAL CYCLE (NO AMBIENT FLUID CHANGE)") # extract the net unit mass net work - unit_mass_net_work_vec = [(result.process_work_unit[(3,4)] + result.process_work_unit[(1,2)]) for result in self.fundamental_cycle] # kJ/kg - unit_mass_net_heat_input_vec = [result.process_heat_unit[(2,3)] for result in self.fundamental_cycle] # kJ/kg - eta_th_vec = [result.get_efficiency() for result in self.fundamental_cycle] # - + unit_mass_net_work_vec = [ + (result.process_work_unit[(3, 4)] + result.process_work_unit[(1, 2)]) + for result in self.fundamental_cycle + ] # kJ/kg + unit_mass_net_heat_input_vec = [ + result.process_heat_unit[(2, 3)] for result in self.fundamental_cycle + ] # kJ/kg + [result.get_efficiency() for result in self.fundamental_cycle] # - # compute the rating/flowrate/input heat limited flowrate - mass_flowrates = [NGCT.NGCT.get_mass_flowrate( - fluid_ambient=fluid_ambient, - Q_fluid_max=flowrate_max_fluid_m3_per_s, - wdot_turbine=result.process_work_unit[(3,4)], - wdot_compressor=result.process_work_unit[(1,2)], - qdot_HX_combustor=result.process_heat_unit[(2,3)], - power_rated=system_capacity_mw*1000.0, # expects kW, convert from MW - heatrate_fuel_capacity=fuel_source_in, - ) for fluid_ambient, result, fuel_source_in in zip( - self.ambient_fluid_list, self.fundamental_cycle, fuel_source_heating_in_kJ_per_s.squeeze().tolist(), - )] - # volume_flowrates = [m/fluid_ambient.density for m, fluid_ambient in zip(mass_flowrates, self.ambient_fluid_list)] # m^3/s - - net_work_vec = [m*w/1000.0 for m, w in zip(mass_flowrates, unit_mass_net_work_vec)] # MW - net_heat_input_vec = [m*q for m, q in zip(mass_flowrates, unit_mass_net_heat_input_vec)] # kJ/s + mass_flowrates = [ + NGCT.NGCT.get_mass_flowrate( + fluid_ambient=fluid_ambient, + Q_fluid_max=flowrate_max_fluid_m3_per_s, + wdot_turbine=result.process_work_unit[(3, 4)], + wdot_compressor=result.process_work_unit[(1, 2)], + qdot_HX_combustor=result.process_heat_unit[(2, 3)], + power_rated=system_capacity_mw * 1000.0, # expects kW, convert from MW + heatrate_fuel_capacity=fuel_source_in, + ) + for fluid_ambient, result, fuel_source_in in zip( + self.ambient_fluid_list, + self.fundamental_cycle, + fuel_source_heating_in_kJ_per_s.squeeze().tolist(), + ) + ] + # volume_flowrates = [ + # m/fluid_ambient.density + # for m, fluid_ambient in zip(mass_flowrates, self.ambient_fluid_list) + # ] # m^3/s + + net_work_vec = [ + m * w / 1000.0 for m, w in zip(mass_flowrates, unit_mass_net_work_vec) + ] # MW + net_heat_input_vec = [ + m * q for m, q in zip(mass_flowrates, unit_mass_net_heat_input_vec) + ] # kJ/s # import matplotlib.pyplot as plt # fig, axes = plt.subplots(5, 1, sharex=True) # axes[0].plot(temperature_degC) + # axes[0].set_ylabel("Temperature (°C)") # axes[1].plot(pressure_Pa) + # axes[1].set_ylabel("Pressure (Pa)") # axes[2].plot(net_work_vec) # axes[2].axhline(system_capacity_mw, c='r', linestyle="--") + # axes[2].set_ylabel("Net Work (MW)") # axes[3].plot(volume_flowrates) # axes[3].axhline(flowrate_max_fluid_m3_per_s, c='r', linestyle="--") + # axes[3].set_ylabel("Max Flowrate (m³/s)") # axes[4].plot(mass_flowrates) + # axes[4].set_ylabel("Mass Flowrate (kg/s)") + # axes[4].set_xlabel("Time Step") + # # plt.show() generator_efficiency = self.config.generator_efficiency - electricity_out = generator_efficiency*np.array(net_work_vec) # MW + electricity_out = generator_efficiency * np.array(net_work_vec) # MW outputs["electricity_out"] = electricity_out outputs[f"{self.fuel_source}_consumed"] = net_heat_input_vec - outputs["rated_electricity_production"] = system_capacity_mw # QUESTION!!!!! WHAT IS RATED??? + outputs["rated_electricity_production"] = ( + system_capacity_mw # QUESTION!!!!! WHAT IS RATED??? + ) max_production = inputs["system_capacity"] * len(electricity_out) * self.dt / 3600.0 @@ -242,7 +272,9 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): outputs["annual_electricity_produced"] = outputs["total_electricity_produced"] * ( 1 / self.fraction_of_year_simulated ) - outputs[f"unmet_{self.commodity}_demand"] = inputs["electricity_command_value"] - electricity_out + outputs[f"unmet_{self.commodity}_demand"] = ( + inputs["electricity_command_value"] - electricity_out + ) if __name__ == "__main__": @@ -286,9 +318,9 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): "performance_model": { "model": "SimpleCycleTurbinePerformanceModel", }, - "cost_model": { - "model": "NaturalGasCostModel", - }, + # "cost_model": { + # "model": "NaturalGasCostModel", + # }, "model_inputs": { "performance_parameters": { "system_capacity_mw": 239.0, @@ -299,6 +331,12 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): "isentropic_efficiency_turbine": 0.90, "generator_efficiency": 1.0, }, + "cost_parameters": { + "capex_per_kw": 1000, # $/kW - typical for NGCC; stolen from ex. 16 + "fixed_opex_per_kw_per_year": 10.0, # $/kW/year; stolen from ex. 16 + "variable_opex_per_mwh": 2.5, # $/MWh; stolen from ex. 16 + "cost_year": 2023, # stolen from ex. 16 + }, }, }, }, @@ -315,7 +353,8 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): "resource_model": "GOESAggregatedSolarAPI", "resource_parameters": { "resource_year": 2024, - "resource_dir": "../../../examples/11_hybrid_energy_plant/tech_inputs/weather/solar", + "resource_dir": "../../../examples/11_hybrid_energy_plant/" + "tech_inputs/weather/solar", "resource_filename": "30.6617_-101.7096_psmv3_60_2013.csv", }, }, @@ -334,16 +373,18 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): ], "resource_to_tech_connections": [ ["site.solar_resource", "ng", "solar_resource_data"], - ] + ], } - h2i = H2IntegrateModel({ - "name": "brayton", - "system_summary": "a brayon-cycle NG plant sim using real time ambient condition data", - "driver_config": driver_config, - "technology_config": technology_config, - "plant_config": plant_config, - }) + h2i = H2IntegrateModel( + { + "name": "brayton", + "system_summary": "a Brayton-cycle NG plant sim using real time ambient condition data", + "driver_config": driver_config, + "technology_config": technology_config, + "plant_config": plant_config, + } + ) # Run the model h2i.run() @@ -351,6 +392,7 @@ def compute(self, inputs, outputs, discrete_inputs, _discrete_outputs): # show N2 diagram if False: import openmdao.api as om + om.n2(h2i.model) # Post-process the results From 9eac59e3db9c893dd19a8563081ab56f814d22de Mon Sep 17 00:00:00 2001 From: Cory Frontin Date: Mon, 20 Jul 2026 13:17:48 -0600 Subject: [PATCH 3/6] precommit fixes --- .../combustion_machines/__init__.py | 4 +- .../test_NGCT_thermo_model.py | 324 +++++++++--------- 2 files changed, 166 insertions(+), 162 deletions(-) diff --git a/h2integrate/converters/combustion_machines/__init__.py b/h2integrate/converters/combustion_machines/__init__.py index 22c71f460..50b322b99 100644 --- a/h2integrate/converters/combustion_machines/__init__.py +++ b/h2integrate/converters/combustion_machines/__init__.py @@ -1,3 +1,3 @@ from h2integrate.converters.combustion_machines.turbine_simple_cycle import ( - SimpleCycleTurbinePerformanceModel -) \ No newline at end of file + SimpleCycleTurbinePerformanceModel, +) diff --git a/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py index aa973e649..b10c914a3 100644 --- a/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py +++ b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py @@ -1,223 +1,227 @@ import numpy as np -import pyfluids import pytest +import pyfluids -from h2integrate.converters.combustion_machines.NGCT_thermo_model import GE_7EA_NGCT_2014 -from h2integrate.converters.combustion_machines.NGCT_thermo_model import GE_7FA05_NGCT -from h2integrate.converters.combustion_machines.NGCT_thermo_model import NGCT -from h2integrate.converters.combustion_machines.NGCT_thermo_model import NGCTResult -from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_compressor_outlet_state_isentropic -from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_compressor_work_rate -from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_heat_transfer -from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_turbine_outlet_state_isentropic -from h2integrate.converters.combustion_machines.NGCT_thermo_model import compute_turbine_work_rate -from h2integrate.converters.combustion_machines.NGCT_thermo_model import humidity_ratio_to_water_mass_fraction -from h2integrate.converters.combustion_machines.NGCT_thermo_model import make_humid_air_mixture +from h2integrate.converters.combustion_machines.NGCT_thermo_model import ( + NGCT, + GE_7FA05_NGCT, + GE_7EA_NGCT_2014, + NGCTResult, + compute_heat_transfer, + make_humid_air_mixture, + compute_turbine_work_rate, + compute_compressor_work_rate, + humidity_ratio_to_water_mass_fraction, + compute_turbine_outlet_state_isentropic, + compute_compressor_outlet_state_isentropic, +) @pytest.fixture(scope="module") def iso_conditions(): - return { - "pressure": 101325.0, - "temperature": 15.0, - "relative_humidity": 60.0, - } + return { + "pressure": 101325.0, + "temperature": 15.0, + "relative_humidity": 60.0, + } @pytest.fixture(scope="module") def iso_ambient_state(iso_conditions): - working_fluid = make_humid_air_mixture( - iso_conditions["pressure"], - iso_conditions["temperature"], - iso_conditions["relative_humidity"], - ) - return working_fluid.with_state( - pyfluids.Input.pressure(iso_conditions["pressure"]), - pyfluids.Input.temperature(iso_conditions["temperature"]), - ) + working_fluid = make_humid_air_mixture( + iso_conditions["pressure"], + iso_conditions["temperature"], + iso_conditions["relative_humidity"], + ) + return working_fluid.with_state( + pyfluids.Input.pressure(iso_conditions["pressure"]), + pyfluids.Input.temperature(iso_conditions["temperature"]), + ) @pytest.fixture(scope="module") def basic_ngct(): - return NGCT( - ratio_P=15.0, - Trel_firing=1300.0, - isentropic_efficiency_compressor=0.85, - isentropic_efficiency_turbine=0.90, - Q_fluid_max=1.0, - ) + return NGCT( + ratio_P=15.0, + Trel_firing=1300.0, + isentropic_efficiency_compressor=0.85, + isentropic_efficiency_turbine=0.90, + Q_fluid_max=1.0, + ) @pytest.fixture(scope="module") def ge_7fa05(): - return GE_7FA05_NGCT() + return GE_7FA05_NGCT() @pytest.fixture(scope="module") def ge_7ea_2014(): - return GE_7EA_NGCT_2014() + return GE_7EA_NGCT_2014() @pytest.fixture( - params=[ - {"pressure": 101325.0, "temperature": 0.0, "relative_humidity": 40.0}, - {"pressure": 101325.0, "temperature": 35.0, "relative_humidity": 80.0}, - {"pressure": 80000.0, "temperature": 15.0, "relative_humidity": 60.0}, - ] + params=[ + {"pressure": 101325.0, "temperature": 0.0, "relative_humidity": 40.0}, + {"pressure": 101325.0, "temperature": 35.0, "relative_humidity": 80.0}, + {"pressure": 80000.0, "temperature": 15.0, "relative_humidity": 60.0}, + ] ) def varied_ambient_state(request): - conditions = request.param - working_fluid = make_humid_air_mixture( - conditions["pressure"], - conditions["temperature"], - conditions["relative_humidity"], - ) - return working_fluid.with_state( - pyfluids.Input.pressure(conditions["pressure"]), - pyfluids.Input.temperature(conditions["temperature"]), - ) + conditions = request.param + working_fluid = make_humid_air_mixture( + conditions["pressure"], + conditions["temperature"], + conditions["relative_humidity"], + ) + return working_fluid.with_state( + pyfluids.Input.pressure(conditions["pressure"]), + pyfluids.Input.temperature(conditions["temperature"]), + ) @pytest.mark.unit def test_humidity_ratio_to_water_mass_fraction_roundtrip(): - humidity_ratio = 0.015 - water_mass_fraction = humidity_ratio_to_water_mass_fraction(humidity_ratio) - recovered_humidity_ratio = water_mass_fraction / (1.0 - water_mass_fraction) + humidity_ratio = 0.015 + water_mass_fraction = humidity_ratio_to_water_mass_fraction(humidity_ratio) + recovered_humidity_ratio = water_mass_fraction / (1.0 - water_mass_fraction) - assert water_mass_fraction == pytest.approx(0.014778325123152709) - assert recovered_humidity_ratio == pytest.approx(humidity_ratio) + assert water_mass_fraction == pytest.approx(0.014778325123152709) + assert recovered_humidity_ratio == pytest.approx(humidity_ratio) @pytest.mark.unit def test_make_humid_air_mixture_builds_valid_state(iso_conditions): - working_fluid = make_humid_air_mixture( - iso_conditions["pressure"], - iso_conditions["temperature"], - iso_conditions["relative_humidity"], - ) - fluid_ambient = working_fluid.with_state( - pyfluids.Input.pressure(iso_conditions["pressure"]), - pyfluids.Input.temperature(iso_conditions["temperature"]), - ) - - assert np.isfinite(fluid_ambient.density) - assert np.isfinite(fluid_ambient.enthalpy) - assert np.isfinite(fluid_ambient.entropy) - assert fluid_ambient.density > 0.0 - assert fluid_ambient.pressure == pytest.approx(iso_conditions["pressure"]) - assert fluid_ambient.temperature == pytest.approx(iso_conditions["temperature"]) + working_fluid = make_humid_air_mixture( + iso_conditions["pressure"], + iso_conditions["temperature"], + iso_conditions["relative_humidity"], + ) + fluid_ambient = working_fluid.with_state( + pyfluids.Input.pressure(iso_conditions["pressure"]), + pyfluids.Input.temperature(iso_conditions["temperature"]), + ) + + assert np.isfinite(fluid_ambient.density) + assert np.isfinite(fluid_ambient.enthalpy) + assert np.isfinite(fluid_ambient.entropy) + assert fluid_ambient.density > 0.0 + assert fluid_ambient.pressure == pytest.approx(iso_conditions["pressure"]) + assert fluid_ambient.temperature == pytest.approx(iso_conditions["temperature"]) @pytest.mark.unit def test_compressor_outlet_state_isentropic_preserves_entropy(iso_ambient_state): - ratio_p = 12.6 - compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, ratio_p) + ratio_p = 12.6 + compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, ratio_p) - assert compressed.pressure == pytest.approx(iso_ambient_state.pressure * ratio_p) - assert compressed.entropy == pytest.approx(iso_ambient_state.entropy, rel=1e-6, abs=1e-3) - assert compressed.temperature > iso_ambient_state.temperature + assert compressed.pressure == pytest.approx(iso_ambient_state.pressure * ratio_p) + assert compressed.entropy == pytest.approx(iso_ambient_state.entropy, rel=1e-6, abs=1e-3) + assert compressed.temperature > iso_ambient_state.temperature @pytest.mark.unit def test_turbine_outlet_state_isentropic_preserves_entropy(iso_ambient_state): - combusted = iso_ambient_state.with_state( - pyfluids.Input.pressure(iso_ambient_state.pressure * 12.6), - pyfluids.Input.temperature(1300.0), - ) - exhaust = compute_turbine_outlet_state_isentropic(combusted, 12.6) + combusted = iso_ambient_state.with_state( + pyfluids.Input.pressure(iso_ambient_state.pressure * 12.6), + pyfluids.Input.temperature(1300.0), + ) + exhaust = compute_turbine_outlet_state_isentropic(combusted, 12.6) - assert exhaust.pressure == pytest.approx(combusted.pressure / 12.6) - assert exhaust.entropy == pytest.approx(combusted.entropy, rel=1e-6, abs=1e-3) - assert exhaust.temperature < combusted.temperature + assert exhaust.pressure == pytest.approx(combusted.pressure / 12.6) + assert exhaust.entropy == pytest.approx(combusted.entropy, rel=1e-6, abs=1e-3) + assert exhaust.temperature < combusted.temperature @pytest.mark.unit def test_helper_sign_conventions(iso_ambient_state): - compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, 10.0) - combusted = compressed.with_state( - pyfluids.Input.pressure(compressed.pressure), - pyfluids.Input.temperature(1300.0), - ) - exhaust = compute_turbine_outlet_state_isentropic(combusted, 10.0) - - wdot_compressor = compute_compressor_work_rate(iso_ambient_state, compressed, isentropic_efficiency=0.85) - wdot_turbine = compute_turbine_work_rate(combusted, exhaust, isentropic_efficiency=0.90) - qdot_combustor = compute_heat_transfer(compressed, combusted) - qdot_exhaust = compute_heat_transfer(exhaust, iso_ambient_state) - - assert wdot_compressor < 0.0 - assert wdot_turbine > 0.0 - assert qdot_combustor > 0.0 - assert qdot_exhaust < 0.0 + compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, 10.0) + combusted = compressed.with_state( + pyfluids.Input.pressure(compressed.pressure), + pyfluids.Input.temperature(1300.0), + ) + exhaust = compute_turbine_outlet_state_isentropic(combusted, 10.0) + + wdot_compressor = compute_compressor_work_rate( + iso_ambient_state, compressed, isentropic_efficiency=0.85 + ) + wdot_turbine = compute_turbine_work_rate(combusted, exhaust, isentropic_efficiency=0.90) + qdot_combustor = compute_heat_transfer(compressed, combusted) + qdot_exhaust = compute_heat_transfer(exhaust, iso_ambient_state) + + assert wdot_compressor < 0.0 + assert wdot_turbine > 0.0 + assert qdot_combustor > 0.0 + assert qdot_exhaust < 0.0 @pytest.mark.unit def test_ngct_result_aggregates_cycle_quantities(iso_ambient_state): - result = NGCTResult(desc="synthetic cycle") - for index in range(1, 5): - result.add_state(index, iso_ambient_state, f"state {index}") + result = NGCTResult(desc="synthetic cycle") + for index in range(1, 5): + result.add_state(index, iso_ambient_state, f"state {index}") - result.mass_flowrate = 5.0 - result.add_process(1, 2, -10.0, 0.0, "compressor") - result.add_process(2, 3, 0.0, 30.0, "combustor") - result.add_process(3, 4, 20.0, 0.0, "turbine") - result.add_process(4, 1, 0.0, -5.0, "cooler") + result.mass_flowrate = 5.0 + result.add_process(1, 2, -10.0, 0.0, "compressor") + result.add_process(2, 3, 0.0, 30.0, "combustor") + result.add_process(3, 4, 20.0, 0.0, "turbine") + result.add_process(4, 1, 0.0, -5.0, "cooler") - assert result.get_net_work() == pytest.approx(50.0) - assert result.get_net_heat_input() == pytest.approx(150.0) - assert result.get_net_heat_rejection() == pytest.approx(-25.0) - assert result.get_back_work_ratio() == pytest.approx(0.5) - assert result.get_efficiency() == pytest.approx(1.0 / 3.0) + assert result.get_net_work() == pytest.approx(50.0) + assert result.get_net_heat_input() == pytest.approx(150.0) + assert result.get_net_heat_rejection() == pytest.approx(-25.0) + assert result.get_back_work_ratio() == pytest.approx(0.5) + assert result.get_efficiency() == pytest.approx(1.0 / 3.0) @pytest.mark.integration def test_run_turbine_model_singleton_matches_batch_result(ge_7ea_2014, iso_ambient_state): - result_single = ge_7ea_2014.run_turbine_model(iso_ambient_state) - result_batch = ge_7ea_2014.run_turbine_model([iso_ambient_state])[0] + result_single = ge_7ea_2014.run_turbine_model(iso_ambient_state) + result_batch = ge_7ea_2014.run_turbine_model([iso_ambient_state])[0] - assert result_single.mass_flowrate == pytest.approx(result_batch.mass_flowrate) - assert result_single.get_net_work() == pytest.approx(result_batch.get_net_work()) - assert result_single.get_efficiency() == pytest.approx(result_batch.get_efficiency()) - assert result_single.states[4].temperature == pytest.approx(result_batch.states[4].temperature) + assert result_single.mass_flowrate == pytest.approx(result_batch.mass_flowrate) + assert result_single.get_net_work() == pytest.approx(result_batch.get_net_work()) + assert result_single.get_efficiency() == pytest.approx(result_batch.get_efficiency()) + assert result_single.states[4].temperature == pytest.approx(result_batch.states[4].temperature) @pytest.mark.regression def test_run_turbine_model_applies_minimum_mass_flow_constraint(basic_ngct, iso_ambient_state): - unconstrained = basic_ngct.run_turbine_model(iso_ambient_state) - specific_net_work = unconstrained.get_net_work() / unconstrained.mass_flowrate - specific_heat_input = unconstrained.process_heat_unit[(2, 3)] + unconstrained = basic_ngct.run_turbine_model(iso_ambient_state) + specific_net_work = unconstrained.get_net_work() / unconstrained.mass_flowrate + specific_heat_input = unconstrained.process_heat_unit[(2, 3)] - power_limit = 0.75 * unconstrained.get_net_work() - heat_limit = 0.60 * unconstrained.get_net_heat_input() + power_limit = 0.75 * unconstrained.get_net_work() + heat_limit = 0.60 * unconstrained.get_net_heat_input() - constrained = basic_ngct.run_turbine_model( - iso_ambient_state, - power_rated=power_limit, - heatrate_fuel_capacity=heat_limit, - ) + constrained = basic_ngct.run_turbine_model( + iso_ambient_state, + power_rated=power_limit, + heatrate_fuel_capacity=heat_limit, + ) - expected_mass_flowrate = min( - unconstrained.mass_flowrate, - power_limit / specific_net_work, - heat_limit / specific_heat_input, - ) + expected_mass_flowrate = min( + unconstrained.mass_flowrate, + power_limit / specific_net_work, + heat_limit / specific_heat_input, + ) - assert constrained.mass_flowrate == pytest.approx(expected_mass_flowrate) + assert constrained.mass_flowrate == pytest.approx(expected_mass_flowrate) @pytest.mark.regression def test_ge_7fa05_iso_design_point_regression(ge_7fa05, iso_ambient_state): - result = ge_7fa05.run_turbine_model(iso_ambient_state) + result = ge_7fa05.run_turbine_model(iso_ambient_state) - assert result.get_efficiency() == pytest.approx( - ge_7fa05.design_conditions["eta_th_ISO"], - rel=0.05, - ) - assert result.get_net_work() == pytest.approx( - ge_7fa05.design_conditions["W_net_ISO"], - rel=0.08, - ) - assert result.states[4].temperature > iso_ambient_state.temperature + assert result.get_efficiency() == pytest.approx( + ge_7fa05.design_conditions["eta_th_ISO"], + rel=0.05, + ) + assert result.get_net_work() == pytest.approx( + ge_7fa05.design_conditions["W_net_ISO"], + rel=0.08, + ) + assert result.states[4].temperature > iso_ambient_state.temperature # @pytest.mark.regression @@ -237,24 +241,24 @@ def test_ge_7fa05_iso_design_point_regression(ge_7fa05, iso_ambient_state): @pytest.mark.integration def test_ideal_brayton_cycle_energy_balance(iso_ambient_state): - ideal_ngct = NGCT( - ratio_P=15.0, - Trel_firing=1300.0, - isentropic_efficiency_compressor=1.0, - isentropic_efficiency_turbine=1.0, - Q_fluid_max=1.0, - ) - result = ideal_ngct.run_turbine_model(iso_ambient_state) - net_heat = result.get_net_heat_input() + result.get_net_heat_rejection() + ideal_ngct = NGCT( + ratio_P=15.0, + Trel_firing=1300.0, + isentropic_efficiency_compressor=1.0, + isentropic_efficiency_turbine=1.0, + Q_fluid_max=1.0, + ) + result = ideal_ngct.run_turbine_model(iso_ambient_state) + net_heat = result.get_net_heat_input() + result.get_net_heat_rejection() - assert result.get_net_work() == pytest.approx(net_heat, rel=1e-9) + assert result.get_net_work() == pytest.approx(net_heat, rel=1e-9) @pytest.mark.integration def test_varied_ambient_conditions_produce_finite_outputs(ge_7ea_2014, varied_ambient_state): - result = ge_7ea_2014.run_turbine_model(varied_ambient_state) + result = ge_7ea_2014.run_turbine_model(varied_ambient_state) - assert np.isfinite(result.mass_flowrate) - assert np.isfinite(result.get_net_work()) - assert np.isfinite(result.get_efficiency()) - assert 0.0 < result.get_efficiency() < 1.0 + assert np.isfinite(result.mass_flowrate) + assert np.isfinite(result.get_net_work()) + assert np.isfinite(result.get_efficiency()) + assert 0.0 < result.get_efficiency() < 1.0 From abeece6277e6b93adda5380f236125649809bdfb Mon Sep 17 00:00:00 2001 From: Cory Frontin Date: Mon, 20 Jul 2026 16:48:11 -0600 Subject: [PATCH 4/6] add pyfluids dependency --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 3346b8887..558ffd353 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -115,6 +115,7 @@ gis = [ "reverse_geocoder" ] ard = ["ard-nrel"] +thermo = ["pyfluids"] extras = ["h2integrate[gis,ard]"] all = ["h2integrate[develop,extras,examples]"] From 06d8f1752f80f4c650322a7cd6abbdf3a329742f Mon Sep 17 00:00:00 2001 From: Cory Frontin Date: Mon, 20 Jul 2026 16:53:00 -0600 Subject: [PATCH 5/6] pyfluids to main deps but also with optional --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 558ffd353..b54024232 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,6 +48,7 @@ dependencies = [ "openmeteo-requests", "pandas>=2.0.3", "ProFAST", + "pyfluids>=4.0.0", "Pyomo>=6.1.2", "rainflow", "requests", From 54f26b52bb4cabc8b98478b1fd494ec38c0fc594 Mon Sep 17 00:00:00 2001 From: Cory Frontin Date: Fri, 24 Jul 2026 12:02:28 -0600 Subject: [PATCH 6/6] prep for generalization of code for recip. engines, apply black reformatting --- .../combustion_machines/NGCT_thermo_model.py | 307 ++---------------- .../test_NGCT_thermo_model.py | 26 +- .../combustion_machines/thermo_tools.py | 304 +++++++++++++++++ 3 files changed, 340 insertions(+), 297 deletions(-) create mode 100644 h2integrate/converters/combustion_machines/thermo_tools.py diff --git a/h2integrate/converters/combustion_machines/NGCT_thermo_model.py b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py index 6eadbc3cc..fb3252e8a 100644 --- a/h2integrate/converters/combustion_machines/NGCT_thermo_model.py +++ b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py @@ -1,274 +1,16 @@ import numpy as np import pyfluids - -DRY_AIR_MASS_FRACTIONS = { - # generated by copilot, matches engineeringtoolbox.com data - pyfluids.FluidsList.Nitrogen: 75.51760366879898, - pyfluids.FluidsList.Oxygen: 23.139561050152018, - pyfluids.FluidsList.Argon: 1.2881375564363484, - pyfluids.FluidsList.CarbonDioxide: 0.05469772461264626, -} - - -def humidity_ratio_to_water_mass_fraction(humidity_ratio: float) -> float: - return humidity_ratio / (1.0 + humidity_ratio) - - -def enforce_gas_phase(fluid): - phases = getattr(pyfluids, "Phases", None) - phase_gas = getattr(phases, "Gas", None) if phases is not None else None - if phase_gas is None: - return fluid - - fluid.specify_phase(phase_gas) - return fluid - - -def make_humid_air_mixture( - pressure: float, - temp_rel: float, - rel_humidity: float, -) -> pyfluids.Mixture: - # get absolute mass humidity ratio using conditions w/ pyfluid humidair - humidity_ratio = ( - pyfluids.HumidAir() - .with_state( - pyfluids.InputHumidAir.pressure(pressure), - pyfluids.InputHumidAir.temperature(temp_rel), - pyfluids.InputHumidAir.relative_humidity(rel_humidity), - ) - .humidity - ) - - # use a mix of water and dry air - water_mass_fraction = 100.0 * humidity_ratio_to_water_mass_fraction(humidity_ratio) - dry_air_mass_fraction = 100.0 - water_mass_fraction - - # create a mixture using our dry air mix plus water vapor - fluids = [*list(DRY_AIR_MASS_FRACTIONS.keys()), pyfluids.FluidsList.Water] - fractions = [ - dry_air_mass_fraction * mass_fraction / 100.0 - for mass_fraction in DRY_AIR_MASS_FRACTIONS.values() - ] + [water_mass_fraction] - humid_air_mixture = pyfluids.Mixture(fluids, fractions) - enforce_gas_phase(humid_air_mixture) - - return humid_air_mixture - - -def compute_compressor_outlet_state_isentropic( - state_inlet: pyfluids.fluids.abstract_fluid.AbstractFluid, - ratio_P: float, -) -> pyfluids.fluids.abstract_fluid.AbstractFluid: - # compute the outlet state - state_outlet = state_inlet.with_state( - pyfluids.Input.pressure( - ratio_P * state_inlet.pressure - ), # set by pressure ratio of the compressor - pyfluids.Input.entropy(state_inlet.entropy), # by isentropic compression assumption - ) - - return state_outlet - - -def compute_work_rate_isentropic( - state_inlet, - state_outlet, - mass_flow_fluid=1.0, # by default unity => returns unit-mass kJ/kg -): - # work done by a system should be positive - return mass_flow_fluid * ( - state_inlet.enthalpy / 1.0e3 - state_outlet.enthalpy / 1.0e3 - ) # kW = kJ/s - - -def compute_turbine_work_rate( - state_inlet, - state_outlet, - mass_flow_fluid=1.0, - isentropic_efficiency=1.0, -): - work_rate_isentropic = compute_work_rate_isentropic( - state_inlet, - state_outlet, - mass_flow_fluid=mass_flow_fluid, - ) - return isentropic_efficiency * work_rate_isentropic # kW = kJ/s - - -def compute_compressor_work_rate( - state_inlet, - state_outlet, - mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg - isentropic_efficiency=1.0, -): - work_rate_isentropic = compute_work_rate_isentropic( - state_inlet, - state_outlet, - mass_flow_fluid=mass_flow_fluid, - ) - return 1.0 / isentropic_efficiency * work_rate_isentropic # kW - - -def compute_heat_transfer( - state_inlet, - state_outlet, - mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg -): - # heat transfer to a system should be positive - return mass_flow_fluid * (state_outlet.enthalpy / 1.0e3 - state_inlet.enthalpy / 1.0e3) # kJ/s - - -def compute_turbine_outlet_state_isentropic(state_inlet, ratio_P): - state_outlet = state_inlet.with_state( - pyfluids.Input.pressure( - state_inlet.pressure / ratio_P - ), # set by pressure ratio of the turbine - pyfluids.Input.entropy(state_inlet.entropy), # by isentropic expansion assumption - ) - - return state_outlet - - -class NGCTResult: - closed: bool = True - desc: str = "" - mass_flowrate: np.ndarray # kg/s - states: dict[int, pyfluids.fluids.abstract_fluid.AbstractFluid] - process_work_unit: dict[tuple[int, int], np.ndarray] - process_heat_unit: dict[tuple[int, int], np.ndarray] - state_names: dict[int, str] - process_names: dict[tuple[int, int], str] - - def __init__( - self, - desc: str | None = None, - closed: bool | None = None, - ): - if closed: - self.closed = closed - if desc: - self.desc = desc - - self.mass_flowrate = None # kg/s - self.states = {} - self.process_work_unit = {} - self.process_heat_unit = {} - self.state_names = {} - self.process_names = {} - - def add_state( - self, - index: int, - state: pyfluids.fluids.abstract_fluid.AbstractFluid, - name: str, - ): - # validations - if index in self.states: - raise ValueError(f"index {index} already in states") - if index in self.state_names: - raise ValueError(f"index {index} already in states") - - # add the state - self.states[index] = state - self.state_names[index] = name - - def add_process( - self, - index_in: int, - index_out: int, - work_unit: np.ndarray, # kJ/kg - heat_unit: np.ndarray, # kJ/kg - name: str, - ): - process_index = (index_in, index_out) - # validations - if (index_in not in self.states) or (index_in not in self.state_names): - raise ValueError( - f"index {index_in} must be in states before process " - f"({index_in}, {index_out}) can be added." - ) - if (index_out not in self.states) or (index_out not in self.state_names): - raise ValueError( - f"index {index_out} must be in states before process " - f"({index_in}, {index_out}) can be added." - ) - if process_index in self.states: - raise ValueError(f"process {process_index} already in states") - if process_index in self.state_names: - raise ValueError(f"process {process_index} already in states") - - # add the process - self.process_work_unit[process_index] = work_unit - self.process_heat_unit[process_index] = heat_unit - self.process_names[process_index] = name - - def get_net_work(self): - mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - return mass_flowrate * sum(self.process_work_unit.values()) # kJ/s - - def get_net_heat_input(self): - mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - return mass_flowrate * sum(np.maximum(0.0, list(self.process_heat_unit.values()))) # kJ/s - - def get_net_heat_rejection(self): - mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - return mass_flowrate * sum(np.minimum(0.0, list(self.process_heat_unit.values()))) # kJ/s - - def get_back_work_ratio(self): - back_work = -sum(np.minimum(0.0, list(self.process_work_unit.values()))) # kJ/s - forward_work = sum(np.maximum(0.0, list(self.process_work_unit.values()))) # kJ/s - return back_work / forward_work # - - - def get_efficiency(self): - net_work = self.get_net_work() # kJ/kg - net_heat_input = self.get_net_heat_input() # kJ/kg - return net_work / net_heat_input # - - - def get_fuel_mass_flow(self, LHV_fuel, heating_process=(2, 3)): - mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - combustor_heat_input_rate = ( - mass_flowrate * self.process_heat_unit[heating_process] - ) # [kg/s]*[kJ/kg] = kJ/s - mass_flowrate_fuel = combustor_heat_input_rate / (1000.0 * LHV_fuel) # [J/s]/[J/kg] = kg/s - return mass_flowrate_fuel - - def get_fuel_air_ratio(self, LHV_fuel, heating_process=(2, 3)): - mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 - mass_flowrate_fuel = self.get_fuel_mass_flow(LHV_fuel, heating_process=heating_process) - return mass_flowrate_fuel / mass_flowrate - - def print_states( - self, - ): - for idx_state, state in self.states.items(): - print(f"{(idx_state):>3d}:", end="") - print(f" T={state.temperature:.2f}°C;", end="") - print(f" h={state.enthalpy/1.0e3:.2f} MJ/kg;", end="") - print(f" P={state.pressure/1.0e3:.2f} kPa;", end="") - print() - - def __str__(self): - output = f"\n{self.desc}\n" if self.desc else "\n" - output += "states:\n" - for idx_state, state in self.states.items(): - state_name = self.state_names.get(idx_state, "") - output += ( - f"\t{idx_state} ({state_name}): T={state.temperature:.2f} C, " - f"h={state.enthalpy/1.0e3:.2f} MJ/kg, P={state.pressure/1.0e3:.2f} kPa\n" - ) - - output += "\nprocesses:\n" - for (idx_in, idx_out), work_unit in self.process_work_unit.items(): - heat_unit = self.process_heat_unit.get((idx_in, idx_out), 0.0) - process_name = self.process_names.get((idx_in, idx_out), "") - output += ( - f"\t{idx_in} -> {idx_out} ({process_name}): " - f"W={work_unit:.2f} kJ/kg, Q={heat_unit:.2f} kJ/kg\n" - ) - - return output +from h2integrate.converters.combustion_machines.thermo_tools import ( + ThermodynamicCycleResult, + enforce_gas_phase, + compute_heat_transfer, + make_humid_air_mixture, + compute_turbine_work_rate, + compute_compressor_work_rate, + compute_isentropic_expansion_outlet_state, + compute_isentropic_compression_outlet_state, +) class NGCT: @@ -297,8 +39,9 @@ def __init__( Trel_firing: float, # deg C, relative temperature isentropic_efficiency_compressor: float, # -, compressor isentropic efficiency isentropic_efficiency_turbine: float, # -, turbine isentropic efficiency - Q_fluid_max: float - | None = None, # m**3/s, design volumetric flowrate; unit-mass analysis if None + Q_fluid_max: ( + float | None + ) = None, # m**3/s, design volumetric flowrate; unit-mass analysis if None design_conditions=None, # dict, design conditions dictionary ): # set in variables @@ -313,8 +56,10 @@ def __init__( def run_turbine_model( self, - fluid_ambient_list: list[pyfluids.fluids.abstract_fluid.AbstractFluid] - | pyfluids.fluids.abstract_fluid.AbstractFluid, + fluid_ambient_list: ( + list[pyfluids.fluids.abstract_fluid.AbstractFluid] + | pyfluids.fluids.abstract_fluid.AbstractFluid + ), power_rated: float = np.inf, # kW, maximum rated power of the system heatrate_fuel_capacity: float = np.inf, # kJ/s, maximum allowable heat rate ): @@ -334,9 +79,9 @@ def run_turbine_model( # note: gas phase enforcement speeds up computation significantly enforce_gas_phase(fluid_ambient) # ensure the ambient fluid is in the gas phase - fluid_compressed = compute_compressor_outlet_state_isentropic( + fluid_compressed = compute_isentropic_compression_outlet_state( fluid_ambient, - self.ratio_P, + ratio_P=self.ratio_P, ) # after the compressor enforce_gas_phase(fluid_compressed) # ensure the compressed fluid is in the gas phase fluid_combusted = fluid_compressed.with_state( @@ -344,9 +89,9 @@ def run_turbine_model( pyfluids.Input.temperature(self.Trel_firing), ) # after the combustion chamber enforce_gas_phase(fluid_combusted) # ensure the combusted fluid is in the gas phase - fluid_exhaust = compute_turbine_outlet_state_isentropic( + fluid_exhaust = compute_isentropic_expansion_outlet_state( fluid_combusted, - self.ratio_P, + ratio_P=self.ratio_P, ) # exhaust from the turbine enforce_gas_phase(fluid_exhaust) # ensure the exhaust fluid is in the gas phase @@ -388,7 +133,7 @@ def run_turbine_model( ) # use static method ### PACK THE RESULT - result = NGCTResult() + result = ThermodynamicCycleResult() # add mass flowrate result.mass_flowrate = mass_flowrate @@ -436,14 +181,6 @@ def run_turbine_model( return results - @staticmethod - def celsius_to_kelvin(degC): - return degC + 273.15 # K - - @staticmethod - def kelvin_to_celsius(degK): - return degK - 273.15 # C - @staticmethod def get_mass_flowrate( fluid_ambient: pyfluids.fluids.abstract_fluid.AbstractFluid, diff --git a/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py index b10c914a3..092e49c4b 100644 --- a/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py +++ b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py @@ -2,18 +2,20 @@ import pytest import pyfluids -from h2integrate.converters.combustion_machines.NGCT_thermo_model import ( - NGCT, - GE_7FA05_NGCT, - GE_7EA_NGCT_2014, - NGCTResult, +from h2integrate.converters.combustion_machines.thermo_tools import ( + ThermodynamicCycleResult, compute_heat_transfer, make_humid_air_mixture, compute_turbine_work_rate, compute_compressor_work_rate, humidity_ratio_to_water_mass_fraction, - compute_turbine_outlet_state_isentropic, - compute_compressor_outlet_state_isentropic, + compute_isentropic_expansion_outlet_state, + compute_isentropic_compression_outlet_state, +) +from h2integrate.converters.combustion_machines.NGCT_thermo_model import ( + NGCT, + GE_7FA05_NGCT, + GE_7EA_NGCT_2014, ) @@ -113,7 +115,7 @@ def test_make_humid_air_mixture_builds_valid_state(iso_conditions): @pytest.mark.unit def test_compressor_outlet_state_isentropic_preserves_entropy(iso_ambient_state): ratio_p = 12.6 - compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, ratio_p) + compressed = compute_isentropic_compression_outlet_state(iso_ambient_state, ratio_p) assert compressed.pressure == pytest.approx(iso_ambient_state.pressure * ratio_p) assert compressed.entropy == pytest.approx(iso_ambient_state.entropy, rel=1e-6, abs=1e-3) @@ -126,7 +128,7 @@ def test_turbine_outlet_state_isentropic_preserves_entropy(iso_ambient_state): pyfluids.Input.pressure(iso_ambient_state.pressure * 12.6), pyfluids.Input.temperature(1300.0), ) - exhaust = compute_turbine_outlet_state_isentropic(combusted, 12.6) + exhaust = compute_isentropic_expansion_outlet_state(combusted, 12.6) assert exhaust.pressure == pytest.approx(combusted.pressure / 12.6) assert exhaust.entropy == pytest.approx(combusted.entropy, rel=1e-6, abs=1e-3) @@ -135,12 +137,12 @@ def test_turbine_outlet_state_isentropic_preserves_entropy(iso_ambient_state): @pytest.mark.unit def test_helper_sign_conventions(iso_ambient_state): - compressed = compute_compressor_outlet_state_isentropic(iso_ambient_state, 10.0) + compressed = compute_isentropic_compression_outlet_state(iso_ambient_state, 10.0) combusted = compressed.with_state( pyfluids.Input.pressure(compressed.pressure), pyfluids.Input.temperature(1300.0), ) - exhaust = compute_turbine_outlet_state_isentropic(combusted, 10.0) + exhaust = compute_isentropic_expansion_outlet_state(combusted, 10.0) wdot_compressor = compute_compressor_work_rate( iso_ambient_state, compressed, isentropic_efficiency=0.85 @@ -157,7 +159,7 @@ def test_helper_sign_conventions(iso_ambient_state): @pytest.mark.unit def test_ngct_result_aggregates_cycle_quantities(iso_ambient_state): - result = NGCTResult(desc="synthetic cycle") + result = ThermodynamicCycleResult(desc="synthetic cycle") for index in range(1, 5): result.add_state(index, iso_ambient_state, f"state {index}") diff --git a/h2integrate/converters/combustion_machines/thermo_tools.py b/h2integrate/converters/combustion_machines/thermo_tools.py new file mode 100644 index 000000000..7e10abb08 --- /dev/null +++ b/h2integrate/converters/combustion_machines/thermo_tools.py @@ -0,0 +1,304 @@ +import numpy as np +import pyfluids + + +DRY_AIR_MASS_FRACTIONS = { + # generated by copilot, matches engineeringtoolbox.com data + pyfluids.FluidsList.Nitrogen: 75.51760366879898, + pyfluids.FluidsList.Oxygen: 23.139561050152018, + pyfluids.FluidsList.Argon: 1.2881375564363484, + pyfluids.FluidsList.CarbonDioxide: 0.05469772461264626, +} + + +def celsius_to_kelvin(degC): + return degC + 273.15 # K + + +def kelvin_to_celsius(degK): + return degK - 273.15 # C + + +def humidity_ratio_to_water_mass_fraction(humidity_ratio: float) -> float: + return humidity_ratio / (1.0 + humidity_ratio) + + +def enforce_gas_phase(fluid): + phases = getattr(pyfluids, "Phases", None) + phase_gas = getattr(phases, "Gas", None) if phases is not None else None + if phase_gas is None: + return fluid + + fluid.specify_phase(phase_gas) + return fluid + + +def make_humid_air_mixture( + pressure: float, + temp_rel: float, + rel_humidity: float, +) -> pyfluids.Mixture: + # get absolute mass humidity ratio using conditions w/ pyfluid humidair + humidity_ratio = ( + pyfluids.HumidAir() + .with_state( + pyfluids.InputHumidAir.pressure(pressure), + pyfluids.InputHumidAir.temperature(temp_rel), + pyfluids.InputHumidAir.relative_humidity(rel_humidity), + ) + .humidity + ) + + # use a mix of water and dry air + water_mass_fraction = 100.0 * humidity_ratio_to_water_mass_fraction(humidity_ratio) + dry_air_mass_fraction = 100.0 - water_mass_fraction + + # create a mixture using our dry air mix plus water vapor + fluids = [*list(DRY_AIR_MASS_FRACTIONS.keys()), pyfluids.FluidsList.Water] + fractions = [ + dry_air_mass_fraction * mass_fraction / 100.0 + for mass_fraction in DRY_AIR_MASS_FRACTIONS.values() + ] + [water_mass_fraction] + humid_air_mixture = pyfluids.Mixture(fluids, fractions) + enforce_gas_phase(humid_air_mixture) + + return humid_air_mixture + + +def compute_isentropic_compression_outlet_state( + state_inlet: pyfluids.fluids.abstract_fluid.AbstractFluid, + ratio_P: float | None = None, + ratio_compression: float | None = None, +) -> pyfluids.fluids.abstract_fluid.AbstractFluid: + if not (ratio_P or ratio_compression) or (ratio_P and ratio_compression): + raise ValueError( + "compression process must specify exactly one of pressure or compression ratio" + ) + # compute the outlet state + state_outlet = state_inlet.with_state( + ( + pyfluids.Input.pressure(ratio_P * state_inlet.pressure) + if ratio_P + else pyfluids.Input.specific_volume( + 1.0 / ratio_compression * state_inlet.specific_volume + ) + ), # set by pressure/compression ratio of the compressor + pyfluids.Input.entropy(state_inlet.entropy), # by isentropic compression assumption + ) + + return state_outlet + + +def compute_work_rate_isentropic( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, # by default unity => returns unit-mass kJ/kg +): + # work done by a system should be positive + return mass_flow_fluid * ( + state_inlet.enthalpy / 1.0e3 - state_outlet.enthalpy / 1.0e3 + ) # kW = kJ/s + + +def compute_turbine_work_rate( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, + isentropic_efficiency=1.0, +): + work_rate_isentropic = compute_work_rate_isentropic( + state_inlet, + state_outlet, + mass_flow_fluid=mass_flow_fluid, + ) + return isentropic_efficiency * work_rate_isentropic # kW = kJ/s + + +def compute_compressor_work_rate( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg + isentropic_efficiency=1.0, +): + work_rate_isentropic = compute_work_rate_isentropic( + state_inlet, + state_outlet, + mass_flow_fluid=mass_flow_fluid, + ) + return 1.0 / isentropic_efficiency * work_rate_isentropic # kW + + +def compute_heat_transfer( + state_inlet, + state_outlet, + mass_flow_fluid=1.0, # by default unity, returns unit-mass kJ/kg +): + # heat transfer to a system should be positive + return mass_flow_fluid * (state_outlet.enthalpy / 1.0e3 - state_inlet.enthalpy / 1.0e3) # kJ/s + + +def compute_isentropic_expansion_outlet_state( + state_inlet, + ratio_P: float | None = None, + ratio_compression: float | None = None, +): + if not (ratio_P or ratio_compression) or (ratio_P and ratio_compression): + raise ValueError( + "compression process must specify exactly one of pressure or compression ratio" + ) + # compute the outlet state + state_outlet = state_inlet.with_state( + ( + pyfluids.Input.pressure(1.0 / ratio_P * state_inlet.pressure) + if ratio_P + else pyfluids.Input.specific_volume(ratio_compression * state_inlet.specific_volume) + ), # set by pressure/compression ratio of the turbine + pyfluids.Input.entropy(state_inlet.entropy), # by isentropic expansion assumption + ) + + return state_outlet + + +class ThermodynamicCycleResult: + closed: bool = True + desc: str = "" + mass_flowrate: np.ndarray # kg/s + states: dict[int, pyfluids.fluids.abstract_fluid.AbstractFluid] + process_work_unit: dict[tuple[int, int], np.ndarray] + process_heat_unit: dict[tuple[int, int], np.ndarray] + state_names: dict[int, str] + process_names: dict[tuple[int, int], str] + + def __init__( + self, + desc: str | None = None, + closed: bool | None = None, + ): + if closed: + self.closed = closed + if desc: + self.desc = desc + + self.mass_flowrate = None # kg/s + self.states = {} + self.process_work_unit = {} + self.process_heat_unit = {} + self.state_names = {} + self.process_names = {} + + def add_state( + self, + index: int, + state: pyfluids.fluids.abstract_fluid.AbstractFluid, + name: str, + ): + # validations + if index in self.states: + raise ValueError(f"index {index} already in states") + if index in self.state_names: + raise ValueError(f"index {index} already in states") + + # add the state + self.states[index] = state + self.state_names[index] = name + + def add_process( + self, + index_in: int, + index_out: int, + work_unit: np.ndarray, # kJ/kg + heat_unit: np.ndarray, # kJ/kg + name: str, + ): + process_index = (index_in, index_out) + # validations + if (index_in not in self.states) or (index_in not in self.state_names): + raise ValueError( + f"index {index_in} must be in states before process " + f"({index_in}, {index_out}) can be added." + ) + if (index_out not in self.states) or (index_out not in self.state_names): + raise ValueError( + f"index {index_out} must be in states before process " + f"({index_in}, {index_out}) can be added." + ) + if process_index in self.states: + raise ValueError(f"process {process_index} already in states") + if process_index in self.state_names: + raise ValueError(f"process {process_index} already in states") + + # add the process + self.process_work_unit[process_index] = work_unit + self.process_heat_unit[process_index] = heat_unit + self.process_names[process_index] = name + + def get_net_work(self): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + return mass_flowrate * sum(self.process_work_unit.values()) # kJ/s + + def get_net_heat_input(self): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + return mass_flowrate * sum(np.maximum(0.0, list(self.process_heat_unit.values()))) # kJ/s + + def get_net_heat_rejection(self): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + return mass_flowrate * sum(np.minimum(0.0, list(self.process_heat_unit.values()))) # kJ/s + + def get_back_work_ratio(self): + back_work = -sum(np.minimum(0.0, list(self.process_work_unit.values()))) # kJ/s + forward_work = sum(np.maximum(0.0, list(self.process_work_unit.values()))) # kJ/s + return back_work / forward_work # - + + def get_efficiency(self): + net_work = self.get_net_work() # kJ/kg + net_heat_input = self.get_net_heat_input() # kJ/kg + return net_work / net_heat_input # - + + def get_fuel_mass_flow(self, LHV_fuel, heating_process=(2, 3)): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + combustor_heat_input_rate = ( + mass_flowrate * self.process_heat_unit[heating_process] + ) # [kg/s]*[kJ/kg] = kJ/s + mass_flowrate_fuel = combustor_heat_input_rate / (1000.0 * LHV_fuel) # [J/s]/[J/kg] = kg/s + return mass_flowrate_fuel + + def get_fuel_air_ratio(self, LHV_fuel, heating_process=(2, 3)): + mass_flowrate = self.mass_flowrate if self.mass_flowrate else 1.0 + mass_flowrate_fuel = self.get_fuel_mass_flow(LHV_fuel, heating_process=heating_process) + return mass_flowrate_fuel / mass_flowrate + + def print_states( + self, + ): + for idx_state, state in self.states.items(): + print(f"{(idx_state):>3d}:", end="") + print(f" T={state.temperature:.2f}°C;", end="") + print(f" h={state.enthalpy/1.0e3:.2f} MJ/kg;", end="") + print(f" P={state.pressure/1.0e3:.2f} kPa;", end="") + print(f" rho={state.density:.2f} kg/m**3;", end="") + print() + + def __str__(self): + output = f"\n{self.desc}\n" if self.desc else "\n" + output += "states:\n" + for idx_state, state in self.states.items(): + state_name = self.state_names.get(idx_state, "") + output += ( + f"\t{idx_state} ({state_name}): " + f"T={state.temperature:.2f} C ({celsius_to_kelvin(state.temperature):.2f} K), " + f"h={state.enthalpy/1.0e3:.2f} MJ/kg, " + f"P={state.pressure/1.0e3:.2f} kPa\n" + f"rho={state.density:.2f} kg/m**3\n" + ) + + if len(self.process_names): + output += "\nprocesses:\n" + for (idx_in, idx_out), work_unit in self.process_work_unit.items(): + heat_unit = self.process_heat_unit.get((idx_in, idx_out), 0.0) + process_name = self.process_names.get((idx_in, idx_out), "") + output += ( + f"\t{idx_in} -> {idx_out} ({process_name}): " + f"W={work_unit:.2f} kJ/kg, Q={heat_unit:.2f} kJ/kg\n" + ) + + return output