diff --git a/src/watertap_contrib/reflo/analysis/multiperiod/thermal/thermal_process_steady_state_flowsheet.py b/src/watertap_contrib/reflo/analysis/multiperiod/thermal/thermal_process_steady_state_flowsheet.py new file mode 100644 index 00000000..7dbba5d0 --- /dev/null +++ b/src/watertap_contrib/reflo/analysis/multiperiod/thermal/thermal_process_steady_state_flowsheet.py @@ -0,0 +1,488 @@ +from pyomo.environ import ( + Var, + ConcreteModel, + Objective, + Param, + Expression, + Constraint, + Block, + log10, + TransformationFactory, + assert_optimal_termination, + value, + units as pyunits, + NonNegativeReals, +) + +from idaes.core.util.initialization import propagate_state as _prop_state +from idaes.core.util.model_statistics import * +from pyomo.network import Arc +from idaes.core import FlowsheetBlock +from idaes.core.solvers.get_solver import get_solver +from idaes.models.unit_models import Product, Feed +from idaes.core.util.model_statistics import * +from idaes.core.util.scaling import ( + set_scaling_factor, + calculate_scaling_factors, + constraint_scaling_transform, +) +from idaes.core import UnitModelCostingBlock +from idaes.core.util.initialization import propagate_state +from idaes.core.solvers.get_solver import get_solver +import idaes.logger as idaeslog + +from watertap_contrib.reflo.unit_models.thermal_energy_storage import ThermalEnergyStorage +from watertap_contrib.reflo.solar_models.zero_order.flat_plate_physical import FlatPlatePhysical +from watertap_contrib.reflo.core import SolarModelType + +from idaes.models.unit_models.heat_exchanger import( + HeatExchanger, + delta_temperature_lmtd_callback, + delta_temperature_amtd_callback, + HeatExchangerFlowPattern, + +) +# Using a slightly modified version tp use the same property package across all the other models +from watertap_contrib.reflo.unit_models.heat_exchanger_ntu import ( + HeatExchangerNTU, + HXNTUInitializer, +) +from watertap.core.util.model_diagnostics.infeasible import * + +from idaes.models.unit_models import Heater + +from watertap.property_models.water_prop_pack import WaterParameterBlock + +from idaes.core.util.scaling import ( + calculate_scaling_factors, + constraint_scaling_transform, + unscaled_variables_generator, + unscaled_constraints_generator, + badly_scaled_var_generator, + list_badly_scaled_variables, + constraint_autoscale_large_jac +) + + +def propagate_state(arc): + _prop_state(arc) + + +def build_thermal_flowsheet(m = None): + + ''' + This function adds the relevant components to the flowsheet + ''' + + m.fs = FlowsheetBlock(dynamic = False) + m.fs.properties = WaterParameterBlock() + + # Include FPC model + m.fs.fpc = FlatPlatePhysical( + property_package = m.fs.properties, solar_model_type=SolarModelType.physical + ) + + # Include TES model + m.fs.tes = ThermalEnergyStorage( + property_package = m.fs.properties + ) + + # HX between FPC and TES + m.fs.hx_solar = HeatExchangerNTU( + hot_side_name = 'fpc', + cold_side_name = 'tes', + fpc = {"property_package": m.fs.properties}, + tes = {"property_package": m.fs.properties} + ) + + # Grid heater to meet process inlet temperature set point + m.fs.grid_heater = Heater( + property_package = m.fs.properties) + + return m + + +def create_coupling_variables(blk): + # Create coupling variables + + blk.fs.previous_hx_solar_hot_outlet_temperature = Var( + domain = NonNegativeReals, + initialize = 35+273.15, + bounds = (20+273.15, 99+273.15), + units = pyunits.K, + doc = 'Outlet temperature from HX solar hot side from the previous time step' + ) + + blk.fs.previous_fpc_outlet_temperature = Var( + domain = NonNegativeReals, + initialize = 40+273.15, + bounds = (20+273.15, 99+273.15), + units = pyunits.K, + doc = 'Outlet temperature from FPC from the previous time step' + ) + + blk.fs.previous_tes_tank_temp = Var( + domain = NonNegativeReals, + initialize = 60+273.15, + bounds = (20+273.15, 99+273.15), + units = pyunits.K, + doc= 'Temperature of the thermal storage tank from the previous time step' + ) + + blk.fs.previous_hx_solar_cold_outlet_temperature = Var( + domain = NonNegativeReals, + initialize = 40+273.15, + bounds = (20+273.15, 99+273.15), + units = pyunits.K, + doc = 'Outlet temperature from HX solar cold side from the previous time step' + ) + + + # blk.fs.previous_process_outlet_temperature = Var( + # domain = NonNegativeReals, + # initialize = 35+273.15, + # bounds = (20+273.15, 99+273.15), + # units = pyunits.K, + # doc = 'Outlet temperature from process from the previous time step' + # ) + + blk.fs.previous_grid_duty = Var( + initialize = 0, + units = pyunits.W, + doc = 'Grid heat duty from previous step' + ) + + blk.fs.acc_grid_duty = Var( + initialize = 0, + units = pyunits.W, + doc = 'Accumulated heat duty at each step' + ) + + blk.fs.previous_acc_grid_duty = Var( + initialize = 0, + units = pyunits.W, + doc = 'Accumulated heat duty from previous step' + ) + + return blk + + +def create_feed_streams(m, + mass_fr_fpc = 0.05, + mass_fr_tes_hx_solar=0.1, + mass_fr_tes_process=0.05, + ): + + # Creating state blocks for the outlet stream from TES with the TES tank temperature from the previous step + + # Outlet stream from the FPC going back to the solar HX hot inlet + m.fs.fpc_outlet = Feed(property_package = m.fs.properties) + m.fs.fpc_outlet.properties[0].temperature.fix(m.fs.previous_fpc_outlet_temperature()) + m.fs.fpc_outlet.properties[0].flow_mass_phase_comp['Liq','H2O'].fix(mass_fr_fpc) + m.fs.fpc_outlet.properties[0].flow_mass_phase_comp['Vap','H2O'].fix(0) + m.fs.fpc_outlet.properties[0].pressure.fix(101325) + + + # Outlet stream from the TES going back to the solar HX cold inlet + m.fs.tes_hx_outlet = Feed(property_package = m.fs.properties) + m.fs.tes_hx_outlet.properties[0].temperature.fix(m.fs.previous_tes_tank_temp) + m.fs.tes_hx_outlet.properties[0].flow_mass_phase_comp['Liq','H2O'].fix(mass_fr_tes_hx_solar) + m.fs.tes_hx_outlet.properties[0].flow_mass_phase_comp['Vap','H2O'].fix(0) + m.fs.tes_hx_outlet.properties[0].pressure.fix(101325) + + + # Outlet stream from the solar HX hot outlet to the FPC inlet + m.fs.hx_solar_hot_outlet = Feed(property_package = m.fs.properties) + m.fs.hx_solar_hot_outlet.properties[0].temperature.fix(m.fs.previous_hx_solar_hot_outlet_temperature()) + m.fs.hx_solar_hot_outlet.properties[0].flow_mass_phase_comp['Liq','H2O'].fix(mass_fr_fpc) + m.fs.hx_solar_hot_outlet.properties[0].flow_mass_phase_comp['Vap','H2O'].fix(0) + m.fs.hx_solar_hot_outlet.properties[0].pressure.fix(101325) + + + # Outlet stream from the solar HX cold outlet to the TES + m.fs.hx_solar_cold_outlet = Feed(property_package = m.fs.properties) + m.fs.hx_solar_cold_outlet.properties[0].temperature.fix(m.fs.previous_hx_solar_cold_outlet_temperature()) + m.fs.hx_solar_cold_outlet.properties[0].flow_mass_phase_comp['Liq','H2O'].fix(mass_fr_tes_hx_solar) + m.fs.hx_solar_cold_outlet.properties[0].flow_mass_phase_comp['Vap','H2O'].fix(0) + m.fs.hx_solar_cold_outlet.properties[0].pressure.fix(101325) + + + # Outlet stream from process outlet to the TES + # m.fs.process_outlet = Feed(property_package = m.fs.properties) + # m.fs.process_outlet.properties[0].temperature.fix(m.fs.previous_process_outlet_temperature()) + # m.fs.process_outlet.properties[0].flow_mass_phase_comp['Liq','H2O'].fix(mass_fr_tes_hx_solar) + # m.fs.process_outlet.properties[0].flow_mass_phase_comp['Vap','H2O'].fix(0) + # m.fs.process_outlet.properties[0].pressure.fix(101325) + + + # Outlet stream from the TES going back to the grid heater inlet + m.fs.tes_process_outlet = Feed(property_package = m.fs.properties) + m.fs.tes_process_outlet.properties[0].temperature.fix(m.fs.previous_tes_tank_temp()) + m.fs.tes_process_outlet.properties[0].flow_mass_phase_comp['Liq','H2O'].fix(mass_fr_tes_process) + m.fs.tes_process_outlet.properties[0].flow_mass_phase_comp['Vap','H2O'].fix(0) + m.fs.tes_process_outlet.properties[0].pressure.fix(101325) + + return m + + +def create_arcs(m): + + m.fs.fpc_hx_solar = Arc(source = m.fs.fpc_outlet.outlet, + destination = m.fs.hx_solar.hot_side_inlet, + doc = 'Connect FPC outlet to the solar HX hotside inlet') + + m.fs.hx_solar_fpc = Arc(source = m.fs.hx_solar_hot_outlet.outlet, + destination = m.fs.fpc.inlet, + doc = 'Connect solar hx hot side outlet back to fpc') + + m.fs.tes_hx_solar = Arc(source = m.fs.tes_hx_outlet.outlet, + destination = m.fs.hx_solar.cold_side_inlet, + doc = 'Connect TES from previous time step to the solar HX cold side inlet') + + m.fs.hx_solar_tes = Arc(source = m.fs.hx_solar_cold_outlet.outlet, + destination = m.fs.tes.tes_hx_inlet, + doc = 'Connect solar HX cold side outlet back to TES inlet') + + m.fs.tes_gridHtr = Arc(source = m.fs.tes_process_outlet.outlet, + destination = m.fs.grid_heater.inlet, + doc = 'Connect TES to the grid heater inlet') + + + # m.fs.process_tes = Arc(source = m.fs.process_outlet.outlet, + # destination = m.fs.tes.tes_process_inlet, + # doc = 'Connect process outlet to the TES') + + TransformationFactory("network.expand_arcs").apply_to(m) + + +def fix_dof_and_initialize( + blk, + mass_fr_tes_hx_solar, + mass_fr_tes_process, + tank_vol, + fpc_collector_area, + process_inlet_temp, + process_outlet_temp, + outlvl=idaeslog.WARNING, +): + """Fix degrees of freedom and initialize the flowsheet + + This function fixes the degrees of freedom of each unit and initializes the entire flowsheet. + + Args: + m: Pyomo `Block` or `ConcreteModel` containing the flowsheet + outlvl: Logger (default: idaeslog.WARNING) + """ + + solver = get_solver() + optarg = solver.options + + blk.fs.previous_grid_duty.fix(0) + blk.fs.previous_acc_grid_duty.fix(0) + + # if process_outlet_temp!= None: + # blk.fs.previous_process_outlet_temperature.fix(process_outlet_temp + 273.15) + + blk.fs.fpc_outlet.properties[0].temperature.fix(blk.fs.previous_fpc_outlet_temperature()) + blk.fs.tes_hx_outlet.properties[0].temperature.fix(blk.fs.previous_tes_tank_temp()) + blk.fs.hx_solar_hot_outlet.properties[0].temperature.fix(blk.fs.previous_hx_solar_hot_outlet_temperature()) + blk.fs.hx_solar_cold_outlet.properties[0].temperature.fix(blk.fs.previous_hx_solar_cold_outlet_temperature()) + # blk.fs.process_outlet.properties[0].temperature.fix(blk.fs.previous_process_outlet_temperature()) + blk.fs.tes_process_outlet.properties[0].temperature.fix(blk.fs.previous_tes_tank_temp()) + + # Initializing and propagate the FPC + # FPC + # blk.fs.previous_hx_solar_hot_outlet_temperature.fix() + blk.fs.hx_solar_hot_outlet.initialize(optarg=optarg) + + propagate_state(blk.fs.hx_solar_fpc) + # m.fs.fpc.inlet.fix() + + # blk.fs.fpc.total_irradiance.fix(GHI) + blk.fs.fpc.collector_area.fix(fpc_collector_area) + blk.fs.fpc.outlet.pressure.fix(101325) + blk.fs.fpc.initialize() + + # Initialize and propragate for solar HX + # blk.fs.previous_fpc_outlet_temperature.fix() + blk.fs.fpc_outlet.initialize() + propagate_state(blk.fs.fpc_hx_solar) + + # blk.fs.previous_tes_tank_temp.fix() + blk.fs.tes_hx_outlet.initialize() + propagate_state(blk.fs.tes_hx_solar) + + blk.fs.hx_solar.effectiveness.fix(0.7) + blk.fs.hx_solar.area.fix(1) + blk.fs.hx_solar.hot_side_outlet.flow_mass_phase_comp[0,'Vap','H2O'].fix(0) + blk.fs.hx_solar.cold_side_outlet.flow_mass_phase_comp[0,'Vap','H2O'].fix(0) + + blk.fs.hx_solar.initialize_build() + + # Initialize and propagate from TES + # blk.fs.previous_hx_solar_cold_outlet_temperature.fix() + blk.fs.hx_solar_cold_outlet.initialize() + propagate_state(blk.fs.hx_solar_tes) + + # blk.fs.previous_process_outlet_temperature.fix() + # blk.fs.process_outlet.initialize() + # propagate_state(blk.fs.process_tes) + + # blk.fs.tes.dt.fix(dt) + # m.fs.tes.tes_hx_inlet.fix() + blk.fs.tes.tes_initial_temperature.fix(blk.fs.previous_tes_tank_temp()) + blk.fs.tes.tes_volume.fix(tank_vol) + blk.fs.tes.hours_storage.fix(6) + # m.fs.tes.heat_load.fix(0.5) + + blk.fs.tes.tes_process_inlet.temperature.fix(process_outlet_temp + 273.15) + blk.fs.tes.tes_process_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(mass_fr_tes_process) + blk.fs.tes.tes_process_inlet.flow_mass_phase_comp[0, "Vap", "H2O"].fix(0) + blk.fs.tes.tes_process_inlet.pressure.fix(101325) + + # Fix outlet vapor flow to be 0 + blk.fs.tes.tes_hx_outlet.flow_mass_phase_comp[0, "Vap", "H2O"].fix(0) + blk.fs.tes.tes_hx_outlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(mass_fr_tes_hx_solar) + blk.fs.tes.tes_hx_outlet.pressure.fix() + # m.fs.tes.tes_hx_outlet.temperature.fix(m.fs.previous_tes_tank_temp()) + + blk.fs.tes.tes_process_outlet.flow_mass_phase_comp[0, "Vap", "H2O"].fix(0) + blk.fs.tes.tes_process_outlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(mass_fr_tes_process) + blk.fs.tes.tes_process_outlet.pressure.fix() + # m.fs.tes.tes_process_outlet.temperature.fix(m.fs.previous_tes_tank_temp()) + try: + blk.fs.tes.initialize() + except ValueError: + print_infeasible_constraints(blk.fs.tes) + print(blk.fs.tes.display()) + + # Grid Heater + # m.fs.grid_heater.outlet.temperature[0].setub(99+273.15) + propagate_state(blk.fs.tes_gridHtr) + blk.fs.tes_process_outlet.initialize() + blk.fs.grid_heater.outlet.flow_mass_phase_comp[0,'Vap','H2O'].fix(0) + # m.fs.grid_heater.outlet.pressure.fix() + blk.fs.grid_heater.outlet.temperature[0].fix(process_inlet_temp + 273.15) + + blk.fs.grid_heater.initialize() + + calculate_scaling_factors(blk) + + +def create_mp_steady_state( + m= None, + dt=3600, + GHI = 900, + mass_fr_fpc = 0.05, + mass_fr_tes_hx_solar = 0.1, + mass_fr_tes_process = 0.05, +): + + if m is None: + m = ConcreteModel() + + m = build_thermal_flowsheet(m) + print('\nDegrees of freedom after build thermal flowsheet:',degrees_of_freedom(m)) + + m = create_coupling_variables(m) + print('\nDegrees of freedom after creating coupling variables:',degrees_of_freedom(m)) + create_feed_streams(m, mass_fr_fpc, mass_fr_tes_hx_solar, mass_fr_tes_process) + print('Degrees of freedom after creating feed streams:',degrees_of_freedom(m)) + create_arcs(m) + + m.fs.fpc.total_irradiance.fix(GHI) + m.fs.tes.dt.fix(dt) + + print('Degrees of freedom after creating multiperiod flowsheet:',degrees_of_freedom(m)) + + m.fs.fpc.total_irradiance.fix(GHI) + # Constraints on the TES output so that the TES tank temperature is the same as the temperature to process and solar hx + + @m.Constraint(doc='Temperature to the process and solar HX are the same') + def eq_hx_process(b,t): + return b.fs.tes.tes_process_outlet.temperature[0] == b.fs.tes.tes_hx_outlet.temperature[0] + + @m.Constraint(doc='Calculated TES temperature is the same as the outlet to solar HX') + def eq_temp_process(b,t): + return b.fs.tes.tes_temperature[0] == b.fs.tes.tes_hx_outlet.temperature[0] + + @m.Constraint(doc='Calculating the accumulated grid heat duty') + def eq_acc_grid_duty(b): + return (b.fs.acc_grid_duty == ((b.fs.previous_grid_duty+b.fs.previous_grid_duty)/2 + b.fs.previous_acc_grid_duty)) + + return m + + +def print_results(m): + + print('\nTank Variables') + print('Initial tank temperature:',m.fs.tes.tes_initial_temperature.value -273.15 ) + print('Tank temperature:', value(m.fs.tes.tes_temperature[0])-273.15) + print('Temperature of stream from process to TES:', m.fs.tes.tes_process_inlet.temperature[0]()-273.15) + print('Temperature of stream exiting to grid heater:',m.fs.tes.tes_process_outlet.temperature[0]()-273.15) + print('Temperature of stream exiting to solar hx:',m.fs.tes.tes_hx_outlet.temperature[0]()-273.15) + print('TES process outlet composition:',m.fs.tes.tes_process_outlet.display() ) + print('TES solar hx outlet composition:',m.fs.tes.tes_hx_outlet.display() ) + + # FPC + print('\nFlat plate collector') + print('Total irradiance:',m.fs.fpc.total_irradiance.value ) + print('FPC inlet temperature:',m.fs.fpc.inlet.temperature[0].value-273.15) + print('FPC outlet temperature:',m.fs.fpc.outlet.temperature[0]()-273.15) + print('GHI: ', m.fs.fpc.total_irradiance()) + print('FPC net heat gain:',m.fs.fpc.net_heat_gain[0].value) + print('FPC outlet composition:',m.fs.fpc.outlet.display() ) + + # Solar HX + print('\nSolar HX') + print('Solar HX hot side inlet:',m.fs.hx_solar.hot_side_inlet.temperature[0].value - 273.15) + print('Solar HX hot side outlet:',m.fs.hx_solar.hot_side_outlet.temperature[0].value - 273.15) + print('Solar HX cold side inlet:',m.fs.hx_solar.cold_side_inlet.temperature[0].value - 273.15) + print('Solar HX cold side outlet:',m.fs.hx_solar.cold_side_outlet.temperature[0].value - 273.15) + print('Solar HX hot side outlet composition:',m.fs.hx_solar.hot_side_outlet.display() ) + print('Solar HX cold side outlet composition:',m.fs.hx_solar.cold_side_outlet.display() ) + + # Grid inlet temperature + print('\nGrid powered heater') + print('Grid heater inlet temperature:',m.fs.grid_heater.inlet.temperature[0]() - 273.15) + print('Grid heater outlet temperature:',m.fs.grid_heater.outlet.temperature[0]() - 273.15) + print('Grid outlet composition:',m.fs.grid_heater.outlet.display() ) + + # Grid heat duty + print('Grid heater heat duty:', m.fs.grid_heater.heat_duty[0]()) + print('Accumulated grid heat duty:', m.fs.acc_grid_duty()) + print('Previous accumulated grid heat duty:', m.fs.previous_acc_grid_duty()) + print('Previous grid duty:', m.fs.previous_grid_duty()) + print('Degrees of freedom:', degrees_of_freedom(m)) + +def main(): + m = create_mp_steady_state( + dt=3600, + GHI = 0, + mass_fr_fpc = 0.05, + mass_fr_tes_hx_solar = 0.1, + mass_fr_tes_process = 0.05, + ) + + fix_dof_and_initialize(m, + mass_fr_tes_hx_solar = 0.1, + mass_fr_tes_process = 0.05, + tank_vol = 2, + fpc_collector_area=3, + process_inlet_temp=50, + process_outlet_temp=40) + + # print_results(m) + print(degrees_of_freedom(m)) + + solver = get_solver() + results = solver.solve(m) + + print(degrees_of_freedom(m)) + + return m + + +if __name__ == "__main__": + m = main() + print_results(m) diff --git a/src/watertap_contrib/reflo/analysis/multiperiod/thermal/thermal_system_mp.py b/src/watertap_contrib/reflo/analysis/multiperiod/thermal/thermal_system_mp.py new file mode 100644 index 00000000..f7a69871 --- /dev/null +++ b/src/watertap_contrib/reflo/analysis/multiperiod/thermal/thermal_system_mp.py @@ -0,0 +1,399 @@ +from pyomo.environ import ( + ConcreteModel, + Var, + value, + Objective, + assert_optimal_termination, + units as pyunits +) +from idaes.core.util.model_statistics import * +from idaes.core import FlowsheetBlock +import logging +import pandas as pd +import datetime +from idaes.apps.grid_integration.multiperiod.multiperiod import MultiPeriodModel +from watertap_contrib.reflo.analysis.multiperiod.thermal.thermal_process_steady_state_flowsheet import ( + create_mp_steady_state, + fix_dof_and_initialize, + print_results +) + +from idaes.core.solvers.get_solver import get_solver +import matplotlib.pyplot as plt +import numpy as np +from watertap.core.solvers import get_solver + + +# Can optimize storage tank volume +def unfix_dof(blk, process_outlet_temperature = None): + blk.fs.previous_hx_solar_hot_outlet_temperature.unfix() + blk.fs.previous_fpc_outlet_temperature.unfix() + blk.fs.previous_tes_tank_temp.unfix() + blk.fs.previous_hx_solar_cold_outlet_temperature.unfix() + blk.fs.previous_acc_grid_duty.unfix() + blk.fs.previous_grid_duty.unfix() + + # Checks if the process outlet temperature is fixed and provided. This is the steady-state operation of the water treatment process + if process_outlet_temperature!= None: + blk.fs.previous_process_outlet_temperature.unfix() + + return None + +def get_variable_pairs_steady_treatment(t1,t2): + """ + This function returns pairs of variables that need to be connected across two time periods + This function assumes that water treatment process is operating at steady state and the stream return to the TES has a constant temperature + + Args: + t1: current time block + t2: next time block + + Returns: + List of tuples with variable pairing + """ + return [ + + (t1.fs.hx_solar.hot_side_outlet.temperature[0], t2.fs.previous_hx_solar_hot_outlet_temperature), + (t1.fs.fpc.outlet.temperature[0], t2.fs.previous_fpc_outlet_temperature), + (t1.fs.tes.tes_temperature[0], t2.fs.previous_tes_tank_temp), + (t1.fs.hx_solar.cold_side_outlet.temperature[0], t2.fs.previous_hx_solar_cold_outlet_temperature), + + # In steady-state treatment case- the process outlet temperature is fixed + # (t1.fs.tes.tes_process_outlet.temperature[0], t2.fs.previous_process_outlet_temperature), + + (t1.fs.acc_grid_duty, t2.fs.previous_acc_grid_duty), + (t1.fs.grid_heater.heat_duty[0], t2.fs.previous_grid_duty), + + ] + +def get_variable_pairs(t1,t2): + """ + This function returns pairs of variables that need to be connected across two time periods + + Args: + t1: current time block + t2: next time block + + Returns: + List of tuples with variable pairing + """ + + return [ + + (t1.fs.hx_solar.hot_side_outlet.temperature[0], t2.fs.previous_hx_solar_hot_outlet_temperature), + (t1.fs.fpc.outlet.temperature[0], t2.fs.previous_fpc_outlet_temperature), + (t1.fs.tes.tes_temperature[0], t2.fs.previous_tes_tank_temp), + (t1.fs.hx_solar.cold_side_outlet.temperature[0], t2.fs.previous_hx_solar_cold_outlet_temperature), + (t1.fs.tes.tes_process_outlet.temperature[0], t2.fs.previous_process_outlet_temperature), + (t1.fs.acc_grid_duty, t2.fs.previous_acc_grid_duty), + (t1.fs.grid_heater.heat_duty[0], t2.fs.previous_grid_duty), + + ] + +def create_multiperiod_thermal_model( + n_time_points = 10, + # 24-hr GHI in Phoenix, AZ on June 18th (W/m2) + GHI = [10, 20, 30, 0, 0, 23, 170, 386, 596, 784, 939, 1031, + 1062, 1031, 938, 790, 599, 383, 166, 31, 0, 0, 0, 0], + initial_temperature = { + 'solar_hx_hot_outlet': 49, + 'fpc_outlet': 49, + 'tes_tank': 49, + 'solar_hx_cold_outlet':49, + 'process_outlet': 49 + }, + process_inlet_temperature = 70, + process_outlet_temperature = None, + + tank_vol = 2, + fpc_collector_area = 3 +): + + """ + This function creates a multi-period pv battery flowsheet object. This object contains + a pyomo model with a block for each time instance. + + Args: + n_time_points: Number of time blocks to create + + Returns: + Object containing multi-period vagmd batch flowsheet model + """ + + if process_outlet_temperature!= None: + + mp = MultiPeriodModel( + n_time_points = n_time_points, + process_model_func = create_mp_steady_state, + linking_variable_func = get_variable_pairs_steady_treatment, + initialization_func = fix_dof_and_initialize, + unfix_dof_func = unfix_dof, + outlvl = logging.WARNING, + ) + + else: + mp = MultiPeriodModel( + n_time_points = n_time_points, + process_model_func = create_mp_steady_state, + linking_variable_func = get_variable_pairs, + initialization_func = fix_dof_and_initialize, + unfix_dof_func = unfix_dof, + outlvl = logging.WARNING, + ) + + + + model_data = { + t: { + "dt" : 3600, + "GHI" : GHI[t], + "mass_fr_fpc" : 0.05, + "mass_fr_tes_hx_solar" : 0.1, + "mass_fr_tes_process" : 0.05, + } + for t in range(n_time_points) + } + + flowsheet_options = { + "mass_fr_fpc" : 0.05, + "mass_fr_tes_hx_solar" : 0.1, + "mass_fr_tes_process" : 0.05, + "tank_vol" : tank_vol, + "fpc_collector_area": fpc_collector_area, + "process_inlet_temp": process_inlet_temperature, # This the set point the grid heater needs to achieve + "process_outlet_temp": process_outlet_temperature # This is the fixed temperature exiting the process + } + + + # create the multiperiod object + # model_data_kwargs-Passes the inputs required for the process_model_func + # flowsheet_options-Passes the variables required for the initializtion function + + mp.build_multi_period_model( + model_data_kwargs = model_data, + flowsheet_options= flowsheet_options, + initialization_options= None, + # { + # "GHI" : 0, + # "mass_fr_fpc" : 0.05, + # "mass_fr_tes_hx_solar" : 0.1, + # "mass_fr_tes_process" : 0.05, + # }, + unfix_dof_options=None, + ) + + active_blocks = mp.get_active_process_blocks() + + # Initialize the first step + active_blocks[0].fs.previous_hx_solar_hot_outlet_temperature.fix(initial_temperature['solar_hx_hot_outlet'] + 273.15) + active_blocks[0].fs.previous_fpc_outlet_temperature.fix(initial_temperature['fpc_outlet'] + 273.15) + active_blocks[0].fs.previous_tes_tank_temp.fix(initial_temperature['tes_tank'] + 273.15) + active_blocks[0].fs.previous_hx_solar_cold_outlet_temperature.fix(initial_temperature['solar_hx_cold_outlet'] + 273.15) + + if process_outlet_temperature == None: + active_blocks[0].fs.previous_process_outlet_temperature.fix(initial_temperature['process_outlet'] + 273.15) + + active_blocks[0].fs.previous_grid_duty.fix(0) + active_blocks[0].fs.previous_acc_grid_duty.fix(0) + + # Initialize and unfix dof for each period + solver = get_solver() + for blk in active_blocks: + print('\nDegrees of freedom before initialization function \n', degrees_of_freedom(blk)) + fix_dof_and_initialize( + blk=blk, + mass_fr_tes_hx_solar=0.1, + mass_fr_tes_process= 0.05, + tank_vol = tank_vol, + fpc_collector_area = fpc_collector_area, + process_inlet_temp=process_inlet_temperature, + process_outlet_temp=process_outlet_temperature + ) + print('dof before init: ', degrees_of_freedom(blk)) + result = solver.solve(blk) + print('\nDegrees of freedom after initialization function \n', degrees_of_freedom(blk)) + unfix_dof(blk) + print('\nDegrees of freedom after unfix_dof function \n', degrees_of_freedom(blk)) + + # Initialize the first step + active_blocks[0].fs.previous_hx_solar_hot_outlet_temperature.fix(initial_temperature['solar_hx_hot_outlet'] + 273.15) + active_blocks[0].fs.previous_fpc_outlet_temperature.fix(initial_temperature['fpc_outlet'] + 273.15) + active_blocks[0].fs.previous_tes_tank_temp.fix(initial_temperature['tes_tank'] + 273.15) + active_blocks[0].fs.previous_hx_solar_cold_outlet_temperature.fix(initial_temperature['solar_hx_cold_outlet'] + 273.15) + + if process_outlet_temperature == None: + active_blocks[0].fs.previous_process_outlet_temperature.fix(initial_temperature['process_outlet'] + 273.15) + + active_blocks[0].fs.previous_grid_duty.fix(0) + active_blocks[0].fs.previous_acc_grid_duty.fix(0) + + return mp + + +if __name__ == "__main__": + + n_time_points = 72 + GHI = [0, 0, 0, 0, 0, 23, 170, 386, 596, 784, 939, 1031, + 1062, 1031, 938, 790, 599, 383, 166, 31, 0, 0, 0, 0,] + + GHI = np.tile(GHI,3) + + process_inlet_temperature = 50 + process_outlet_temperature = 45 + + fpc_collector_area = 3 #m2 + tank_vol = 3 # m3 + + # Initial temperatures in C + initial_temperature = { + 'solar_hx_hot_outlet': 49, + 'fpc_outlet': 49, + 'tes_tank': 49, + 'solar_hx_cold_outlet':49, + 'process_outlet': 49 + } + + + mp = create_multiperiod_thermal_model( + n_time_points = n_time_points, + # 24-hr GHI in Phoenix, AZ on June 18th (W/m2) + GHI = GHI, + initial_temperature = initial_temperature, + + process_inlet_temperature = process_inlet_temperature, + process_outlet_temperature = process_outlet_temperature, + + tank_vol = tank_vol, + fpc_collector_area =fpc_collector_area + + ) + + # print('\nDegrees of freedom after initialization', degrees_of_freedom(mp)) + + solver = get_solver() + results = solver.solve(mp) + print('dof after solving: ', degrees_of_freedom(mp)) + + # assert_optimal_termination(results) + + print('\nStep 1') + print_results(mp.blocks[0].process) + + print('\nStep 2') + print_results(mp.blocks[1].process) + + print('\nStep 3') + print_results(mp.blocks[2].process) + + + tes_temp = [] + grid_duty = [] + + fpc_inlet = [] + fpc_outlet = [] + + hx_hot_inlet = [] + hx_hot_outlet = [] + + hx_cold_inlet = [] + hx_cold_outlet = [] + + process_outlet = [] + grid_heater_outlet = [] + + # fig, ((ax00,ax0,ax1), (ax2, ax3)) = plt.subplots(3,3) + fig = plt.figure() + gs = fig.add_gridspec(3,2) + ax00 = fig.add_subplot(gs[0, :]) + ax0 = fig.add_subplot(gs[1, 0]) + ax1 = fig.add_subplot(gs[1, 1]) + ax2 = fig.add_subplot(gs[2, 0]) + ax3 = fig.add_subplot(gs[2, 1]) + + for n in range(0,n_time_points): + tes_temp.append(mp.blocks[n].process.fs.tes.tes_temperature[0].value-273.15) + process_outlet.append(mp.blocks[n].process.fs.tes.tes_process_inlet.temperature[0]()-273.15) + grid_heater_outlet.append(mp.blocks[n].process.fs.grid_heater.outlet.temperature[0]() - 273.15) + + fpc_inlet.append(mp.blocks[n].process.fs.fpc.inlet.temperature[0].value-273.15) + fpc_outlet.append(mp.blocks[n].process.fs.fpc.outlet.temperature[0]()-273.15) + + hx_hot_inlet.append(mp.blocks[n].process.fs.hx_solar.hot_side_inlet.temperature[0].value - 273.15) + hx_hot_outlet.append(mp.blocks[n].process.fs.hx_solar.hot_side_outlet.temperature[0].value - 273.15) + + hx_cold_inlet.append(mp.blocks[n].process.fs.hx_solar.cold_side_inlet.temperature[0].value - 273.15) + hx_cold_outlet.append(mp.blocks[n].process.fs.hx_solar.cold_side_outlet.temperature[0].value - 273.15) + + grid_duty.append(mp.blocks[n].process.fs.grid_heater.heat_duty[0].value) + + ax00.plot(range(0,n_time_points),GHI[0:n_time_points], label = 'GHI', marker = 'o') + + ax0.plot(range(0,n_time_points),fpc_inlet, label = 'FPC inlet', marker = 'o') + ax0.plot(range(0,n_time_points),fpc_outlet, label = 'FPC outlet', marker = 'o') + + ax1.plot(range(0,n_time_points),hx_hot_inlet, label = 'Solar HX hot inlet', marker = 'o') + ax1.plot(range(0,n_time_points),hx_hot_outlet, label = 'Solar HX hot outlet', marker = 'o') + + ax1.plot(range(0,n_time_points),hx_cold_inlet, label = 'Solar HX cold inlet', marker = 'o') + ax1.plot(range(0,n_time_points),hx_cold_outlet, label = 'Solar HX cold outlet', marker = 'o') + + ax2.plot(range(0,n_time_points),tes_temp, label = 'TES temp', marker = 'o') + ax2.plot(range(0,n_time_points),process_outlet, label = 'Process to TES', marker = 'o') + ax2.plot(range(0,n_time_points),grid_heater_outlet, label = 'Grid heater outlet', marker = 'o') + + ax3.plot(range(0,n_time_points),grid_duty, marker = 'o') + + ax00.set_title('GHI',fontsize = 20) + ax0.set_title('FPC',fontsize = 20) + ax1.set_title('Solar HX',fontsize = 20) + ax2.set_title('TES',fontsize = 20) + ax3.set_title('Additional grid duty',fontsize = 20) + + + ax00.set_ylabel(r'W/m${^2}$',fontsize = 20) + ax0.set_ylabel(r'Temperature (${^o}$C)',fontsize = 20) + ax1.set_ylabel(r'Temperature (${^o}$C)',fontsize = 20) + ax2.set_ylabel(r'Temperature (${^o}$C)',fontsize = 20) + ax3.set_ylabel('Heat Duty (W)',fontsize = 20) + + + ax0.set_ylim([40,70]) + ax1.set_ylim([40,70]) + ax3.set_ylim([-200,800]) + + # Function to find xmin and xmax for the axvspan function + # Get first non-zero point and first zero point + xmin_list = [] + xmax_list = [] + for n in range(0,len(GHI)-1): + if GHI[n]<=0 and GHI[n+1]>0: + xmin_list.append(n) + if GHI[n]>0 and GHI[n+1]<=0: + xmax_list.append(n) + + print(xmin_list) + print(xmax_list) + + for ax in (ax00,ax0,ax1, ax2, ax3): + ax.legend() + ax.set_xlabel('Time (h)',fontsize = 20) + ax.tick_params(axis='both', which='major', labelsize=20) + + for xmin, xmax in zip(xmin_list,xmax_list): + ax.axvspan(xmin,xmax, facecolor ='gold',alpha=0.3) + + # gs.tight_layout(fig, pad = 0.1) + plt.subplots_adjust(left=0.1, + bottom=0.05, + right=0.9, + top=0.9, + wspace=0.2, + hspace=0.5) + plt.show() + + # Calculating the average grid heat duty + grid_duty = np.array(grid_duty) + avg_grid_duty = grid_duty[grid_duty>0].mean() + + print('Average grid:',avg_grid_duty) + \ No newline at end of file diff --git a/src/watertap_contrib/reflo/solar_models/zero_order/flat_plate_physical.py b/src/watertap_contrib/reflo/solar_models/zero_order/flat_plate_physical.py index b083fc52..641dbc0f 100644 --- a/src/watertap_contrib/reflo/solar_models/zero_order/flat_plate_physical.py +++ b/src/watertap_contrib/reflo/solar_models/zero_order/flat_plate_physical.py @@ -262,9 +262,9 @@ def build(self): initialize=1, units=pyunits.kW, doc="Rated plant heat capacity in MW" ) - self.heat_annual = Var( - initialize=1, units=pyunits.MWh, doc="Annual heat generated by flat plate" - ) + # self.heat_annual = Var( + # initialize=1, units=pyunits.MWh, doc="Annual heat generated by flat plate" + # ) self.heat_loss_coeff = Expression(expr=self.FR_UL / self.FR_ta) @@ -355,15 +355,15 @@ def eq_outlet_temp(b, t): def eq_P_pump(b): return b.electricity == (b.pump_power / b.pump_eff) - @self.Constraint(self.flowsheet().config.time, doc="Useful heat generated") - def eq_heat(b, t): - return b.heat == b.net_heat_gain[t] + # @self.Constraint(self.flowsheet().config.time, doc="Useful heat generated") + # def eq_heat(b, t): + # return b.heat == b.net_heat_gain[t] - @self.Constraint(doc="Heat generated annually") - def eq_heat_annual(b): - return b.heat_annual == b.heat * pyunits.convert( - 1 * pyunits.year, to_units=pyunits.hour - ) + # @self.Constraint(doc="Heat generated annually") + # def eq_heat_annual(b): + # return b.heat_annual == b.heat * pyunits.convert( + # 1 * pyunits.year, to_units=pyunits.hour + # ) # Need to check this @self.Constraint(doc="Heat load of system") @@ -388,11 +388,12 @@ def initialize( # Create solver opt = get_solver(solver=solver, options=optarg) - self.inlet_block.initialize( + flags = self.inlet_block.initialize( outlvl=outlvl, optarg=optarg, solver=solver, state_args=state_args, + hold_state = True ) if state_args is None: @@ -421,7 +422,8 @@ def initialize( # solve unit with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(self, tee=slc.tee) - + + self.inlet_block.release_state(flags, outlvl=outlvl) init_log.info("FPC initialization status {}.".format(idaeslog.condition(res))) def calculate_scaling_factors(self): @@ -455,13 +457,13 @@ def calculate_scaling_factors(self): sf = iscale.get_scaling_factor(self.heat_load, default=1e-3) iscale.set_scaling_factor(self.heat_load, sf) - if iscale.get_scaling_factor(self.heat_annual) is None: - sf = iscale.get_scaling_factor(self.heat_annual, default=1e-6) - iscale.set_scaling_factor(self.heat_annual, sf) + # if iscale.get_scaling_factor(self.heat_annual) is None: + # sf = iscale.get_scaling_factor(self.heat_annual, default=1e-6) + # iscale.set_scaling_factor(self.heat_annual, sf) - if iscale.get_scaling_factor(self.heat) is None: - sf = iscale.get_scaling_factor(self.heat, default=1e-3) - iscale.set_scaling_factor(self.heat, sf) + # if iscale.get_scaling_factor(self.heat) is None: + # sf = iscale.get_scaling_factor(self.heat, default=1e-3) + # iscale.set_scaling_factor(self.heat, sf) if iscale.get_scaling_factor(self.electricity) is None: sf = iscale.get_scaling_factor(self.electricity, default=0.1) diff --git a/src/watertap_contrib/reflo/unit_models/heat_exchanger_ntu.py b/src/watertap_contrib/reflo/unit_models/heat_exchanger_ntu.py new file mode 100644 index 00000000..8daf5342 --- /dev/null +++ b/src/watertap_contrib/reflo/unit_models/heat_exchanger_ntu.py @@ -0,0 +1,611 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +IDAES heat exchanger model using effectiveness-NTU method + +Assumptions: + * No phase equilibrium or reactions occur within unit + +""" + +## Modified model to work with mass flow instead of mole + +# Import Pyomo libraries +from pyomo.environ import ( + Block, + check_optimal_termination, + Constraint, + Expression, + Param, + PositiveReals, + Reference, + units as pyunits, + Var, +) +from pyomo.common.config import Bool, ConfigBlock, ConfigValue, In + +# Import IDAES cores +from idaes.core import ( + ControlVolume0DBlock, + declare_process_block_class, + MaterialBalanceType, + EnergyBalanceType, + MomentumBalanceType, + UnitModelBlockData, + useDefault, +) +from idaes.models.unit_models.heat_exchanger import hx_process_config, add_hx_references +from idaes.core.util.config import is_physical_parameter_block +from idaes.core.util.tables import create_stream_table_dataframe +from idaes.core.util.math import smooth_min, smooth_max +from idaes.core.solvers import get_solver +from idaes.core.util.exceptions import InitializationError +import idaes.logger as idaeslog +from idaes.core.initialization import SingleControlVolumeUnitInitializer + +__author__ = "Paul Akula, Andrew Lee" +# Updates by Mukta Hardikar to match REFLO property packages + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +class HXNTUInitializer(SingleControlVolumeUnitInitializer): + """ + Initializer for NTU Heat Exchanger units. + + """ + + def initialization_routine( + self, + model: Block, + plugin_initializer_args: dict = None, + copy_inlet_state: bool = False, + duty=1000 * pyunits.W, + ): + """ + Common initialization routine for NTU Heat Exchangers. + + This routine starts by initializing the hot and cold side properties. Next, the heat + transfer between the two sides is fixed to an initial guess for the heat duty (provided by the duty + argument), the associated constraint deactivated, and the model is then solved. Finally, the heat + duty is unfixed and the heat transfer constraint reactivated followed by a final solve of the model. + + Args: + model: Pyomo Block to be initialized + plugin_initializer_args: dict-of-dicts containing arguments to be passed to plug-in Initializers. + Keys should be submodel components. + copy_inlet_state: bool (default=False). Whether to copy inlet state to other states or not + (0-D control volumes only). Copying will generally be faster, but inlet states may not contain + all properties required elsewhere. + duty: initial guess for heat duty to assist with initialization. Can be a Pyomo expression with units. + + Returns: + Pyomo solver results object + """ + return super(SingleControlVolumeUnitInitializer, self).initialization_routine( + model=model, + plugin_initializer_args=plugin_initializer_args, + copy_inlet_state=copy_inlet_state, + duty=duty, + ) + + def initialize_main_model( + self, + model: Block, + copy_inlet_state: bool = False, + duty=1000 * pyunits.W, + ): + """ + Initialization routine for main NTU HX models. + + Args: + model: Pyomo Block to be initialized. + copy_inlet_state: bool (default=False). Whether to copy inlet state to other states or not + (0-D control volumes only). Copying will generally be faster, but inlet states may not contain + all properties required elsewhere. + duty: initial guess for heat duty to assist with initialization. Can be a Pyomo expression with units. + + Returns: + Pyomo solver results object. + + """ + # TODO: Aside from one differences in constraint names, this is + # identical to the Initializer for the 0D HX unit. + # Set solver options + init_log = idaeslog.getInitLogger( + model.name, self.get_output_level(), tag="unit" + ) + solve_log = idaeslog.getSolveLogger( + model.name, self.get_output_level(), tag="unit" + ) + + # Create solver + solver = get_solver(self.config.solver, self.config.solver_options) + + self.initialize_control_volume(model.hot_side, copy_inlet_state) + init_log.info_high("Initialization Step 1a (hot side) Complete.") + + self.initialize_control_volume(model.cold_side, copy_inlet_state) + init_log.info_high("Initialization Step 1b (cold side) Complete.") + + # --------------------------------------------------------------------- + # Solve unit without heat transfer equation + model.energy_balance_constraint.deactivate() + + model.cold_side.heat.fix(duty) + for i in model.hot_side.heat: + model.hot_side.heat[i].set_value(-duty) + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver.solve(model, tee=slc.tee) + + init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res))) + + model.cold_side.heat.unfix() + model.energy_balance_constraint.activate() + # --------------------------------------------------------------------- + # Solve unit + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver.solve(model, tee=slc.tee) + init_log.info("Initialization Completed, {}".format(idaeslog.condition(res))) + + return res + + +@declare_process_block_class("HeatExchangerNTU") +class HeatExchangerNTUData(UnitModelBlockData): + """Heat Exchanger Unit Model using NTU method.""" + + default_initializer = HXNTUInitializer + + CONFIG = UnitModelBlockData.CONFIG(implicit=True) + + # Configuration template for fluid specific arguments + _SideCONFIG = ConfigBlock() + + _SideCONFIG.declare( + "material_balance_type", + ConfigValue( + default=MaterialBalanceType.useDefault, + domain=In(MaterialBalanceType), + description="Material balance construction flag", + doc="""Indicates what type of mass balance should be constructed, +**default** - MaterialBalanceType.useDefault. +**Valid values:** { +**MaterialBalanceType.useDefault - refer to property package for default +balance type +**MaterialBalanceType.none** - exclude material balances, +**MaterialBalanceType.componentPhase** - use phase component balances, +**MaterialBalanceType.componentTotal** - use total component balances, +**MaterialBalanceType.elementTotal** - use total element balances, +**MaterialBalanceType.total** - use total material balance.}""", + ), + ) + _SideCONFIG.declare( + "energy_balance_type", + ConfigValue( + default=EnergyBalanceType.useDefault, + domain=In(EnergyBalanceType), + description="Energy balance construction flag", + doc="""Indicates what type of energy balance should be constructed, +**default** - EnergyBalanceType.useDefault. +**Valid values:** { +**EnergyBalanceType.useDefault - refer to property package for default +balance type +**EnergyBalanceType.none** - exclude energy balances, +**EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material, +**EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase, +**EnergyBalanceType.energyTotal** - single energy balance for material, +**EnergyBalanceType.energyPhase** - energy balances for each phase.}""", + ), + ) + _SideCONFIG.declare( + "momentum_balance_type", + ConfigValue( + default=MomentumBalanceType.pressureTotal, + domain=In(MomentumBalanceType), + description="Momentum balance construction flag", + doc="""Indicates what type of momentum balance should be constructed, +**default** - MomentumBalanceType.pressureTotal. +**Valid values:** { +**MomentumBalanceType.none** - exclude momentum balances, +**MomentumBalanceType.pressureTotal** - single pressure balance for material, +**MomentumBalanceType.pressurePhase** - pressure balances for each phase, +**MomentumBalanceType.momentumTotal** - single momentum balance for material, +**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""", + ), + ) + _SideCONFIG.declare( + "has_pressure_change", + ConfigValue( + default=False, + domain=Bool, + description="Pressure change term construction flag", + doc="""Indicates whether terms for pressure change should be +constructed, +**default** - False. +**Valid values:** { +**True** - include pressure change terms, +**False** - exclude pressure change terms.}""", + ), + ) + _SideCONFIG.declare( + "property_package", + ConfigValue( + default=useDefault, + domain=is_physical_parameter_block, + description="Property package to use ", + doc="""Property parameter object used to define property calculations + **default** - useDefault. + **Valid values:** { + **useDefault** - use default package from parent model or flowsheet, + **PhysicalParameterObject** - a PhysicalParameterBlock object.}""", + ), + ) + _SideCONFIG.declare( + "property_package_args", + ConfigBlock( + implicit=True, + description="Arguments to use for constructing property package", + doc="""A ConfigBlock with arguments to be passed to + property block(s) and used when constructing these, + **default** - None. + **Valid values:** { + see property package for documentation.}""", + ), + ) + + # Create individual config blocks for hot and cold sides + CONFIG.declare("hot_side", _SideCONFIG(doc="Hot fluid config arguments")) + CONFIG.declare("cold_side", _SideCONFIG(doc="Cold fluid config arguments")) + CONFIG.declare( + "hot_side_name", + ConfigValue( + default=None, + domain=str, + doc="Hot side name, sets control volume and inlet and outlet names", + ), + ) + CONFIG.declare( + "cold_side_name", + ConfigValue( + default=None, + domain=str, + doc="Cold side name, sets control volume and inlet and outlet names", + ), + ) + + def build(self): + # Call UnitModel.build to setup model + super().build() + hx_process_config(self) + + # --------------------------------------------------------------------- + # Build hot-side control volume + self.hot_side = ControlVolume0DBlock( + dynamic=self.config.dynamic, + has_holdup=self.config.has_holdup, + property_package=self.config.hot_side.property_package, + property_package_args=self.config.hot_side.property_package_args, + ) + + # TODO : Add support for phase equilibrium? + self.hot_side.add_state_blocks(has_phase_equilibrium=False) + + self.hot_side.add_material_balances( + balance_type=self.config.hot_side.material_balance_type, + has_phase_equilibrium=False, + ) + + self.hot_side.add_energy_balances( + balance_type=self.config.hot_side.energy_balance_type, + has_heat_transfer=True, + ) + + self.hot_side.add_momentum_balances( + balance_type=self.config.hot_side.momentum_balance_type, + has_pressure_change=self.config.hot_side.has_pressure_change, + ) + + # --------------------------------------------------------------------- + # Build cold-side control volume + self.cold_side = ControlVolume0DBlock( + dynamic=self.config.dynamic, + has_holdup=self.config.has_holdup, + property_package=self.config.cold_side.property_package, + property_package_args=self.config.cold_side.property_package_args, + ) + + self.cold_side.add_state_blocks(has_phase_equilibrium=False) + + self.cold_side.add_material_balances( + balance_type=self.config.cold_side.material_balance_type, + has_phase_equilibrium=False, + ) + + self.cold_side.add_energy_balances( + balance_type=self.config.cold_side.energy_balance_type, + has_heat_transfer=True, + ) + + self.cold_side.add_momentum_balances( + balance_type=self.config.cold_side.momentum_balance_type, + has_pressure_change=self.config.cold_side.has_pressure_change, + ) + + # --------------------------------------------------------------------- + # Add Ports to control volumes + self.add_inlet_port( + name="hot_side_inlet", block=self.hot_side, doc="Hot side inlet port" + ) + self.add_outlet_port( + name="hot_side_outlet", block=self.hot_side, doc="Hot side outlet port" + ) + + self.add_inlet_port( + name="cold_side_inlet", block=self.cold_side, doc="Cold side inlet port" + ) + self.add_outlet_port( + name="cold_side_outlet", block=self.cold_side, doc="Cold side outlet port" + ) + + # --------------------------------------------------------------------- + # Add unit level References + # Set references to balance terms at unit level + self.heat_duty = Reference(self.cold_side.heat[:]) + + # Add references to the user provided aliases (if applicable). + add_hx_references(self) + + # --------------------------------------------------------------------- + # Add performance equations + # All units of measurement will be based on hot side + hunits = self.config.hot_side.property_package.get_metadata().get_derived_units + + # Common heat exchanger variables + self.area = Var( + initialize=1, + units=hunits("area"), + domain=PositiveReals, + doc="Heat transfer area", + ) + + self.heat_transfer_coefficient = Var( + self.flowsheet().time, + initialize=1, + units=hunits("heat_transfer_coefficient"), + domain=PositiveReals, + doc="Overall heat transfer coefficient", + ) + + # Overall energy balance + def rule_energy_balance(blk, t): + return blk.hot_side.heat[t] == -pyunits.convert( + blk.cold_side.heat[t], to_units=hunits("power") + ) + + self.energy_balance_constraint = Constraint( + self.flowsheet().time, rule=rule_energy_balance + ) + + # Add e-NTU variables + self.effectiveness = Var( + self.flowsheet().time, + initialize=1, + units=pyunits.dimensionless, + domain=PositiveReals, + doc="Effectiveness factor for NTU method", + ) + + # Minimum heat capacitance ratio for e-NTU method + self.eps_cmin = Param( + initialize=1e-3, + mutable=True, + units=hunits("power") / hunits("temperature"), + doc="Epsilon parameter for smooth Cmin and Cmax", + ) + + # TODO : Support both mass and mole based flows + def rule_Cmin(blk, t): + caph = ( + blk.hot_side.properties_in[t].flow_vol*blk.hot_side.properties_in[t].dens_mass_phase['Liq'] + * blk.hot_side.properties_in[t].cp_mass_phase['Liq'] + ) + capc = pyunits.convert( + blk.cold_side.properties_in[t].flow_vol*blk.cold_side.properties_in[t].dens_mass_phase['Liq'] + * blk.cold_side.properties_in[t].cp_mass_phase['Liq'], + to_units=hunits("power") / hunits("temperature"), + ) + return smooth_min(caph, capc, eps=blk.eps_cmin) + + self.Cmin = Expression( + self.flowsheet().time, rule=rule_Cmin, doc="Minimum heat capacitance rate" + ) + + def rule_Cmax(blk, t): + caph = ( + blk.hot_side.properties_in[t].flow_vol*blk.hot_side.properties_in[t].dens_mass_phase['Liq'] + * blk.hot_side.properties_in[t].cp_mass_phase['Liq'] + ) + capc = pyunits.convert( + blk.cold_side.properties_in[t].flow_vol*blk.cold_side.properties_in[t].dens_mass_phase['Liq'] + * blk.cold_side.properties_in[t].cp_mass_phase['Liq'], + to_units=hunits("power") / hunits("temperature"), + ) + return smooth_max(caph, capc, eps=blk.eps_cmin) + + self.Cmax = Expression( + self.flowsheet().time, rule=rule_Cmax, doc="Maximum heat capacitance rate" + ) + + # Heat capacitance ratio + def rule_Cratio(blk, t): + return blk.Cmin[t] / blk.Cmax[t] + + self.Cratio = Expression( + self.flowsheet().time, rule=rule_Cratio, doc="Heat capacitance ratio" + ) + + def rule_NTU(blk, t): + return blk.heat_transfer_coefficient[t] * blk.area / blk.Cmin[t] + + self.NTU = Expression( + self.flowsheet().time, rule=rule_NTU, doc="Number of heat transfer units" + ) + + # Heat transfer by e-NTU method + def rule_entu(blk, t): + return blk.hot_side.heat[t] == -( + blk.effectiveness[t] + * blk.Cmin[t] + * ( + blk.hot_side.properties_in[t].temperature + - pyunits.convert( + blk.cold_side.properties_in[t].temperature, + to_units=hunits("temperature"), + ) + ) + ) + + self.heat_duty_constraint = Constraint(self.flowsheet().time, rule=rule_entu) + + # TODO : Add scaling methods + + def initialize_build( + self, + hot_side_state_args=None, + cold_side_state_args=None, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + duty=None, + ): + """ + Heat exchanger initialization method. + + Args: + hot_side_state_args : a dict of arguments to be passed to the + property initialization for the hot side (see documentation of + the specific property package) (default = None). + cold_side_state_args : a dict of arguments to be passed to the + property initialization for the cold side (see documentation of + the specific property package) (default = None). + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None, use + default solver options) + solver : str indicating which solver to use during + initialization (default = None, use default solver) + duty : an initial guess for the amount of heat transferred. This + should be a tuple in the form (value, units), (default + = (1000 J/s)) + + Returns: + None + + """ + # Set solver options + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") + + hot_side = self.hot_side + cold_side = self.cold_side + + # Create solver + opt = get_solver(solver, optarg) + + flags1 = hot_side.initialize( + outlvl=outlvl, optarg=optarg, solver=solver, state_args=hot_side_state_args + ) + + init_log.info_high("Initialization Step 1a (hot side) Complete.") + + flags2 = cold_side.initialize( + outlvl=outlvl, optarg=optarg, solver=solver, state_args=cold_side_state_args + ) + + init_log.info_high("Initialization Step 1b (cold side) Complete.") + + # --------------------------------------------------------------------- + # Solve unit without heat transfer equation + self.energy_balance_constraint.deactivate() + + # Get side 1 and side 2 heat units, and convert duty as needed + s1_units = hot_side.heat.get_units() + s2_units = cold_side.heat.get_units() + + if duty is None: + # Assume 1000 J/s and check for unitless properties + if s1_units is None and s2_units is None: + # Backwards compatibility for unitless properties + s1_duty = -1000 + s2_duty = 1000 + else: + s1_duty = pyunits.convert_value( + -1000, from_units=pyunits.W, to_units=s1_units + ) + s2_duty = pyunits.convert_value( + 1000, from_units=pyunits.W, to_units=s2_units + ) + else: + # Duty provided with explicit units + s1_duty = -pyunits.convert_value( + duty[0], from_units=duty[1], to_units=s1_units + ) + s2_duty = pyunits.convert_value( + duty[0], from_units=duty[1], to_units=s2_units + ) + + cold_side.heat.fix(s2_duty) + for i in hot_side.heat: + hot_side.heat[i].value = s1_duty + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + + init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res))) + + cold_side.heat.unfix() + self.energy_balance_constraint.activate() + + # --------------------------------------------------------------------- + # Solve unit + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res))) + + # --------------------------------------------------------------------- + # Release Inlet state + hot_side.release_state(flags1, outlvl=outlvl) + cold_side.release_state(flags2, outlvl=outlvl) + + init_log.info("Initialization Completed, {}".format(idaeslog.condition(res))) + + if not check_optimal_termination(res): + raise InitializationError( + f"{self.name} failed to initialize successfully. Please check " + f"the output logs for more information." + ) + + return res + + def _get_stream_table_contents(self, time_point=0): + return create_stream_table_dataframe( + { + "Hot Inlet": self.hot_side_inlet, + "Hot Outlet": self.hot_side_outlet, + "Cold Inlet": self.cold_side_inlet, + "Cold Outlet": self.cold_side_outlet, + }, + time_point=time_point, + )