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..fb3252e8a --- /dev/null +++ b/h2integrate/converters/combustion_machines/NGCT_thermo_model.py @@ -0,0 +1,332 @@ +import numpy as np +import pyfluids + +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: + """ + 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 + ) = None, # m**3/s, design volumetric flowrate; unit-mass analysis 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_isentropic_compression_outlet_state( + fluid_ambient, + 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( + 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_isentropic_expansion_outlet_state( + fluid_combusted, + ratio_P=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 = ThermodynamicCycleResult() + + # 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 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] + ) + + # 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..50b322b99 --- /dev/null +++ b/h2integrate/converters/combustion_machines/__init__.py @@ -0,0 +1,3 @@ +from h2integrate.converters.combustion_machines.turbine_simple_cycle import ( + SimpleCycleTurbinePerformanceModel, +) 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..092e49c4b --- /dev/null +++ b/h2integrate/converters/combustion_machines/test_NGCT_thermo_model.py @@ -0,0 +1,266 @@ +import numpy as np +import pytest +import pyfluids + +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_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, +) + + +@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_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) + 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_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) + assert exhaust.temperature < combusted.temperature + + +@pytest.mark.unit +def test_helper_sign_conventions(iso_ambient_state): + 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_isentropic_expansion_outlet_state(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 = ThermodynamicCycleResult(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/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 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..2de6a4bd0 --- /dev/null +++ b/h2integrate/converters/combustion_machines/turbine_simple_cycle.py @@ -0,0 +1,399 @@ +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 + + +@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( + "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( + "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 + [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[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 + + 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, + }, + "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 + }, + }, + }, + }, + } + 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 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() + + # 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", diff --git a/pyproject.toml b/pyproject.toml index 3346b8887..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", @@ -115,6 +116,7 @@ gis = [ "reverse_geocoder" ] ard = ["ard-nrel"] +thermo = ["pyfluids"] extras = ["h2integrate[gis,ard]"] all = ["h2integrate[develop,extras,examples]"]